%%{init: {"flowchart": {"rankSpacing": 20}, "themeVariables": {"fontSize": "20px"}}}%%
flowchart LR
A["Chapter 2<br/>aligned features"] --> B["Training-only<br/>Bayesian fit"]
B --> C["Model weights<br/>predictive mixture"]
C --> D["Frozen rule<br/>weights and cash"]
D --> E["Next outcomes<br/>costs and NAV"]
E --> F["Chapter 3<br/>evidence review"]
Chapter 5: Bayesian Methods for Investment Research
Chapter Introduction
Bayesian methods can change the signals a researcher estimates, the information shared across assets and the portfolio decisions made from those signals. The chapter follows that investment research process while retaining the statistical foundations needed to explain it. Bayes’ theorem combines a prior and a likelihood into a posterior; a predictive distribution integrates over uncertain parameters; a decision rule connects that distribution to an action and its costs.
The original return examples develop priors, likelihoods, normalisation, grid calculations, Metropolis sampling, credible intervals and Student-t robustness. Keep the distinction between a t-distribution’s scale and its standard deviation visible: a reward-to-scale ratio of 2.38 and a standard-deviation Sharpe posterior mean of 1.18 are different quantities, not proof that one statistical school discovers three times as much profit. Posterior conclusions remain conditional on the data definition, likelihood, prior and computational checks.
Sections 5.1–5.4 retain those foundations, regression predictive checks and sequential learning. The dynamic risk example now separates a retrospective fit from a forward filter with parameters estimated on earlier data. Sequential updates must add new evidence; refitting an expanding history from the original prior is a different valid calculation. Carrying an old posterior and applying the same historical likelihood again counts observations twice.
Section 5.5 connects the chapter to the nine-stock feature panel from Chapter 2. A hierarchical normal regression shares information across stock-specific effects through a proper prior. You will calculate shrinkage, interpret feature coefficients in their training-standardised units, compare prior strengths and distinguish uncertainty about a conditional mean from uncertainty about next month’s realised return. Fixed hyperparameters make the posterior transparent; their uncertainty is an explicit extension rather than something silently claimed to have been estimated.
Section 5.6 combines three declared prior models using integrated training evidence, builds posterior predictive mixtures and applies a fixed portfolio rule. Three annual refits use labels available before each freeze. Forecast RMSE, a simple baseline, exact-mixture density scores, exposure, costs and reconciled net accounts are separate pieces of evidence. The public survivor sample has already been explored in earlier chapters; this is a reproducible historical teaching experiment, not newly untouched evidence of a live trading edge.
For self-study, read each definition and assumption, predict an output, run an example, change one declared input in that cell and explain the result. The downloadable investment notebook mirrors §§5.5–5.6 and includes worked repairs. The three-hour class selects experiments and discussion from the full chapter; it does not define a limit on the statistical material available for independent study.
How to use the numbered learning path
Follow the numbered parts in order within each section. Each part states a question and a result. Read the concept, predict the example, run it, and explain the output before attempting the check.
The 3.5-hour class selects examples, group discussion and notebook work. The complete Slides and Book support self-study.
Worked examples run automatically when the page opens or refreshes. Each cell includes required earlier imports, data and definitions under Previously introduced setup. Run Code includes that setup; Copy complete code copies the whole cell. Practice and repair cells run when you choose. Make input changes in the current cell before the operation you are testing.
Table of Contents
Bayesian Foundations and Risk Estimates
| Part | Question |
|---|---|
| 5.1.1 Prior, likelihood and posterior | How do assumptions and data determine a belief? |
| 5.1.2 Posterior calculation on a grid | How is a posterior normalised numerically? |
| 5.1.3 Sampling and uncertainty in a risk statistic | Can posterior draws support a defensible risk summary? |
Prior, likelihood and posterior
Question: How do assumptions and data determine a belief?
Result: Write a return model with an explicit likelihood and prior.
Begin with the definitions and assumptions; then trace the worked example.
Which object does a Bayesian actually compute?
Bayes’ theorem, in its statistical form, reads
\[ \underbrace{p(\theta \mid X)}_{\text{posterior}} \;=\; \frac{\overbrace{p(X \mid \theta)}^{\text{likelihood}}\;\overbrace{p(\theta)}^{\text{prior}}}{\underbrace{p(X)}_{\text{evidence}}} \;\propto\; p(X \mid \theta)\, p(\theta). \]
Three objects appear on the right. The prior \(p(\theta)\) is what you believed about the parameter before this data set arrived — from theory, from other data, or from the simple knowledge that a daily equity return does not have a mean of 50 %. The likelihood \(p(X \mid \theta)\) is exactly the function that classical statistics maximises: the probability of the observed data as a function of the parameter. Their product, normalised so that it integrates to one, is the posterior \(p(\theta \mid X)\) — the distribution of the parameter given the data you actually saw. The denominator \(p(X)\), the marginal probability of the data, is the integral of the numerator over all \(\theta\); it is a constant with respect to \(\theta\), which is why the proportionality sign carries the whole argument. Almost every computational trick in this chapter is a way of never having to compute it.
It helps to state the reversal of roles bluntly. Classical statistics fixes \(\theta\) at some hypothesised value and asks how the data would vary across imaginary repeated samples; the p-value and the confidence interval are statements about that imaginary repetition. Bayesian statistics fixes the data — there is only one October 2015 to October 2020 — and asks how \(\theta\) varies given it. Neither the likelihood nor the prior alone is the answer to a question about \(\theta\). After you have seen 1 258 daily returns, the object that lets you say “the probability that reward-to-t-scale exceeds 1 is 0.99” is the posterior, and only the posterior. A likelihood is not a probability distribution over \(\theta\) at all (it does not integrate to one in \(\theta\)), and a p-value is a probability about data, not about parameters.
Once you have the posterior, three outputs follow immediately and each has a plain meaning. The posterior mean, median or mode are different point estimates; choose one to match the decision or loss function, and inspect skew rather than assuming they coincide. A credible interval at level \(1-\alpha\) is any region holding \(1-\alpha\) of the posterior probability; unlike a confidence interval it does mean “the parameter lies in here with probability \(1-\alpha\)”, because the parameter is the random object. The posterior predictive distribution of a new observation, \(p(\tilde x \mid X) = \int p(\tilde x \mid \theta)\,p(\theta \mid X)\,d\theta\), averages the likelihood of the new point over the posterior and so automatically carries parameter uncertainty into predictions — something a plug-in classical prediction ignores. Section 5.2 makes the predictive the centre of model checking.
A common objection is that the prior is subjective and therefore suspect. Two replies. First, with a data set as large as 1 258 observations and a prior as wide as the ones used below, the prior’s influence on the posterior is numerically negligible — you will see this measured in precision units on the grid. Second, a flat prior is also a prior, and an unreasonable one for many parameters: it says a daily mean return of 0.5 is as plausible as 0.0005. Encoding what you actually believe produces estimates that are more stable in small samples and more honest in large ones. The cost is that you must write the prior down; the benefit is that you cannot pretend you did not have one.
Start with the classical number
The data are GOOG’s daily closes from 23 October 2015 to 22 October 2020 — the notebook’s file, 1 259 prices and hence 1 258 returns. Here risk-free return is set to zero and square-root annualisation assumes no serial dependence. The classical annualised Sharpe ratio is \(\sqrt{252}\,\bar r / s\) with \(\bar r\) the mean daily return and \(s\) the sample standard deviation. The mean is about 0.00080 and the standard deviation about 0.01644; commit to a value for the Sharpe ratio to three decimals before you run the cell.
Mean daily return ≈ 0.00080, daily sd ≈ 0.01644, 252 trading days a year. What does the classical Sharpe ratio print (approximately), and is the excess kurtosis above or below zero?
The Sharpe ratio prints 0.771: \(\sqrt{252} \times 0.00080 / 0.01644 = 15.87 \times 0.0485\). The second line is the one to dwell on. Excess kurtosis is 6.72. A normal distribution has excess kurtosis exactly 0, so these returns have far heavier tails than a normal — a few days contribute a disproportionate share of the fourth moment, and, as it will turn out, of the second. Any model that describes these returns as normal is wrong in a way that matters for exactly the quantity we are after, because the Sharpe ratio’s denominator is a measure of dispersion and the tails are where dispersion is concentrated. The model must say so.
Why Student-t, not normal?
The notebook’s model treats each daily return as a draw from a Student-t distribution with three parameters:
\[ r_t \sim \text{StudentT}(\nu,\ \mu,\ \sigma). \]
Here \(\mu\) is the location (the centre of the distribution), \(\sigma\) is the scale — the width of the central body — and \(\nu\), the degrees of freedom, controls the weight of the tails. As \(\nu \to \infty\) the Student-t converges to the normal with mean \(\mu\) and standard deviation \(\sigma\); at \(\nu = 30\) the two are practically indistinguishable; at \(\nu = 3\) the tails are heavy enough that a five-scale-unit day is rare but perfectly possible; at \(\nu = 1\) the distribution is the Cauchy, which has no mean at all. The degrees of freedom is therefore the parameter that answers the kurtosis of 6.72, and estimating it from the data — rather than fixing it at infinity by choosing a normal — is the single most consequential decision in the model.
One fact about the Student-t will be needed repeatedly, so it is worth stating now. The standard deviation of a \(\text{StudentT}(\nu, \mu, \sigma)\) is not \(\sigma\); it is \(\sigma\sqrt{\nu/(\nu-2)}\), which is finite only for \(\nu > 2\) and which for \(\nu = 3\) is \(\sigma\sqrt{3} \approx 1.73\,\sigma\). The scale describes ordinary days; the standard deviation includes the tail days. Keep the two apart — the whole “Bayesian Sharpe versus classical Sharpe” argument at the end of this section turns on it.
The plot is on a log scale on purpose. The normal density is a parabola in log space — it falls off a cliff, and by four scale units it is at \(10^{-4}\) and still accelerating downward. The t(3) density is almost a straight line in the tails: it decays as a power of \(|x|\) rather than as an exponential of \(x^2\). The two probabilities make this concrete. A four-scale-unit event has probability 6.3e-05 under the normal and 0.028 under t(3) — a ratio of roughly 440. On a normal model a 4σ day is a once-in-sixty-years event; on a t(3) model it happens about seven times a year. GOOG’s kurtosis of 6.72 says the second description is much closer to the truth.
The last line previews the prior on \(\nu\) that the next subsection introduces; hold it in mind.
How to choose priors
The notebook chooses one prior per parameter, and each choice is worth a sentence of justification because each is a decision you will have to make yourself on other data.
| Parameter | Prior | Why |
|---|---|---|
| \(\mu\) | \(N(\bar r,\ s)\) | centred on the sample mean, with standard deviation equal to one daily sd — enormously wide relative to the plausible range of a daily mean |
| \(\sigma\) | \(\text{Uniform}(s/1000,\ 1000\,s)\) | positive, and essentially flat over six orders of magnitude |
| \(\nu\) | \(2 + \text{Exponential}(\lambda = 1/29)\) | mean \(2 + 29 = 31\), but with density highest near 2 so that fat tails remain possible; the shift by 2 guarantees a finite variance |
The prior on \(\mu\) is calibrated from the same sample: this is an empirical choice, not independent prior information. The uniform scale prior is proper only because its bounds are finite, and those bounds are assumptions. The exponential prior on \(\nu-2\) gives probabilities 0.241 for \(\nu<10\) and 0.098 for \(\nu<5\). It restricts \(\nu>2\) so variance exists, while allowing heavy tails. Sensitivity must be assessed rather than declaring a prior irrelevant. In forward research, calibrate any data-dependent prior using eligible training data only and preserve that choice during evaluation.
In Colab — the notebook’s model in PyMC
PyMC is not available in the browser. This is the notebook’s model, verbatim; the rest of this section reproduces its posterior with numpy and scipy.
rmean, rstd = google["return"].mean(), google["return"].std()
with pm.Model() as sr_model:
mean = pm.Normal("mean", mu=rmean, sigma=rstd)
std = pm.Uniform("std", lower=rstd/1000, upper=rstd*1000)
df = pm.Exponential("df", 1/29, initval=5) + 2.0
returns = pm.StudentT("returns", nu=df, mu=mean, sigma=std,
observed=google["return"])
pm.Deterministic("sharpe", np.sqrt(252) * mean / std)
trace = pm.sample(tune=1500, draws=1000, chains=2, random_seed=8888)
az.summary(trace)The notebook’s az.summary reports, for its quantity named sharpe (actually reward-to-t-scale), a posterior mean of 2.349 with sd 0.560 and a 94 % HDI of [1.310, 3.417]; for df (the exponential part, before the shift by 2) a mean of 0.787, so \(\nu \approx 2.79\); for std a mean of 0.010. Effective sample sizes are 1 069 to 1 291 from 2 000 draws and every \(\hat R\) is 1.00. Four lines of priors, one likelihood, one derived quantity — everything else is the sampler, which we now build by hand.
The posterior \(p(\theta \mid X)\). It is a probability distribution over the parameter given the data actually observed. The likelihood is the same function classical statistics maximises but is not a distribution over \(\theta\); the prior is what you believed before; a p-value is a probability about data under a hypothesis, not about the parameter.
Posterior calculation on a grid
Question: How is a posterior normalised numerically?
Result: Compute one- and two-parameter grids and explain their scaling limit.
Builds on Prior, likelihood and posterior.
Warm-up: a posterior on a grid
With one or two parameters you do not need a sampler at all. Evaluate the prior times the likelihood at every point of a fine grid, divide by the sum, and you have a posterior approximation whose accuracy depends on grid resolution and bounds. This is Bayes’ rule made literal, and it is the right first computation because it removes every source of mystery: the posterior is two arrays multiplied and normalised.
Start with one parameter, on the S&P 500 over the same five years. Treat \(\sigma\) as known (set it to the sample standard deviation, 0.01186) and put a prior \(\mu \sim N(0,\ 0.001^2)\) on the daily mean — centred on zero, with a standard deviation of a tenth of a percent, which is a genuinely informative prior for a daily mean. Because a normal likelihood with known variance and a normal prior on the mean are a conjugate pair, the posterior is normal and has a closed form:
\[ \mu \mid y \sim N\!\left(\frac{\mu_0/\tau_0^2 + n\bar y/\sigma^2}{1/\tau_0^2 + n/\sigma^2},\; \Big(\frac{1}{\tau_0^2} + \frac{n}{\sigma^2}\Big)^{-1}\right). \]
Read this in precision units — precision being one over variance. The posterior precision is the prior precision plus the data precision \(n/\sigma^2\); the posterior mean is the average of the prior mean and the sample mean, each weighted by its precision. Conjugacy is what makes the arithmetic possible: when the prior and likelihood belong to matching families the posterior is in the prior’s family and the update is a formula rather than an integral. The Beta–Binomial pair (a Beta prior on a probability updated by counts of successes and failures, so that the prior’s two parameters behave like pseudo-counts) is the other workhorse; Section 5.4 will use the normal–normal update 253 times in a row. The grid is the check that the formula is right.
Posterior precision = prior precision + \(n/\sigma^2\). With \(n = 1258\), \(\sigma = 0.01186\) and prior sd 0.001, what is the posterior sd to six decimals — and will the grid match it?
Grid 0.000497 0.000317, closed form 0.000497 0.000317 — identical to six decimals. Two details of the code carry the lesson.
First, log_post - log_post.max() before the exp is not cosmetic. Each of the 1 258 log-densities is a number around 3, so their sum is about 3 800, and \(e^{3800}\) overflows a double-precision float (the limit is about \(e^{709}\)). Subtracting the maximum makes the largest value exactly \(e^0 = 1\) and every other value a ratio to it; since the posterior is normalised afterwards, the constant subtracted does not matter. Every MCMC sampler in this chapter works in log densities for the same reason.
Second, look at where the posterior mean sits. The sample mean of the S&P 500 over these years is 0.00055; the prior mean is 0; the posterior mean is 0.000497 — between them, and much closer to the data. The prior precision is \(1/0.001^2 = 10^6\); the data precision is \(1258/0.01186^2 \approx 8.9 \times 10^6\). The data have about nine times the weight of the prior, so the posterior mean is roughly 90 % of the way from the prior to the sample mean: \(0.9 \times 0.00055 = 0.000497\). That is the precision-weighted average, and it is the template for every conjugate update in this chapter. A prior only matters when its precision is comparable to the data’s; here even a deliberately informative prior with sd 0.1 % per day loses nine to one.
Two parameters is still easy: a 20 × 20 grid is 400 points. Return to GOOG and put \(\mu\) and \(\sigma\) on a grid under a Student-t likelihood with \(\nu\) fixed at 4, with a flat prior on the grid. The interesting output is not \(\mu\) or \(\sigma\) alone but the Sharpe ratio \(\sqrt{252}\,\mu/\sigma\) evaluated at every grid point and averaged under the posterior weights.
The posterior mean of \(\sigma\) is 0.01082 — well below the sample standard deviation of 0.01644. This is the first appearance of the fact flagged above: the t-scale is not the standard deviation. Under a t(4) the standard deviation is \(\sigma\sqrt{4/2} = 1.41\,\sigma \approx 0.0153\), close to the sample sd; the extra dispersion lives in the tails, and \(\nu\) accounts for it, so \(\sigma\) can be smaller. The posterior mean of \(\mu\) is 0.00139, larger than the sample mean of 0.0008, because the t-likelihood is less impressed by a few extreme negative days than the sample mean is. Put together, the posterior mean reward-to-t-scale is 2.035 and P(reward-to-t-scale > 1) = 0.972. Hold that thought — the full model with \(\nu\) estimated will move the numbers further in the same direction, and the last subsection will come back to ask what they mean.
Why not a grid for everything?
A grid with \(g\) points per dimension and \(k\) parameters needs \(g^k\) likelihood evaluations. Four hundred points per axis and the three-parameter model \((\mu, \sigma, \nu)\) is \(400^3 = 64\) million sums over 1 258 observations — feasible but slow; ten parameters is hopeless. The cost grows exponentially in the number of parameters, and the fraction of the grid that carries any posterior mass shrinks just as fast, so most of the evaluations are wasted on regions of negligible probability. This is the curse of dimensionality, and it is why grids were abandoned in the early 1990s.
The alternative is to stop visiting every point and instead wander through parameter space in a way that spends time in each region in proportion to its posterior probability. If that can be arranged, then the collection of visited points is a sample from the posterior, and every posterior quantity — mean, credible interval, \(P(\text{Sharpe} > 1)\) — is a sample average. Arranging it is what a Markov chain Monte Carlo algorithm does, and the simplest one, Metropolis, is short enough to write and understand in full.
Sampling and uncertainty in a risk statistic
Question: Can posterior draws support a defensible risk summary?
Result: Inspect sampling diagnostics and distinguish t-scale from standard-deviation Sharpe.
Builds on Posterior calculation on a grid.
Metropolis: the sampler behind PyMC
The Metropolis algorithm, published by Metropolis, Rosenbluth, Rosenbluth, Teller and Teller in 1953, is a rule for taking a random walk whose long-run occupancy matches a target density. From the current point \(\theta_t\), with posterior density \(p_t\), propose a candidate \(\theta^\star = \theta_t + \epsilon\) where \(\epsilon\) is a small random step — a normal with some proposal standard deviation. Compute the density \(p^\star\) at the candidate. Then:
- If \(p^\star \ge p_t\) — the proposal is uphill — accept it: \(\theta_{t+1} = \theta^\star\).
- If \(p^\star < p_t\) — downhill — accept it with probability \(p^\star/p_t\); otherwise stay: \(\theta_{t+1} = \theta_t\).
That is the whole algorithm: accept with probability \(\min(1,\ p^\star/p_t)\). Two properties make it work. Uphill moves are always taken, so the chain is drawn toward high-density regions; downhill moves are sometimes taken, in proportion to how much lower the density is, so the chain does not collapse onto the mode but explores the whole distribution. Always accepting would ignore the posterior and give a pure random walk; never going downhill would give an optimiser that finds the mode and stops. The compromise is exactly calibrated — one can prove that the posterior is the stationary distribution of the resulting Markov chain — so that time-in-region is proportional to posterior mass. Suppose \(p^\star = 0.3\,p_t\): Metropolis accepts with probability 0.3, and that is the correct answer regardless of what the prior says about either point, because the prior is already inside \(p\).
The detail that makes the algorithm usable is that only the ratio \(p^\star/p_t\) appears. The normalising constant \(p(X)\) — the integral nobody can compute — is the same in numerator and denominator and cancels. In code you work in logs: if log(u) < log_post(prop) - log_post(current) for a uniform \(u\), because the difference of two log-densities is the log of the ratio, and the log of a uniform is the log of the acceptance threshold. Everything in modern Bayesian computation — Stan, PyMC, NumPyro — is an elaboration of this loop. The proposal may be cleverer, the chains may run in parallel, the diagnostics may be richer, but the accept–reject core is what follows.
One practical matter before the code. The three parameters live on different domains: \(\mu\) on the whole real line, \(\sigma > 0\) and \(\nu > 2\). A random-walk proposal on \(\sigma\) itself would sometimes propose negative scales, which the likelihood cannot evaluate. The fix is to sample the transformed vector \(\theta = (\mu,\ \log\sigma,\ \log(\nu - 2))\), which lives on \(\mathbb{R}^3\), so the chain can never leave the allowed region. The price is a Jacobian: a density on \(\sigma\) becomes, on \(\log\sigma\), the same density times \(\sigma\), which adds \(\log\sigma\) — the ls term — to the log-prior; and likewise \(\log(\nu - 2)\) — the le term — for the exponential prior on \(\nu - 2\). PyMC does this transformation and its Jacobian automatically; when you write the sampler yourself you must not forget it.
The two diagnostic evaluations come first. At the sample standard deviation and \(\nu = 3\) the log-posterior is 3438.4; at half the sample standard deviation it is 3443.2 — higher. Before a single step of the sampler has been taken, the likelihood is already telling you that it prefers a scale smaller than the sample sd, because the tails are being handled by \(\nu\) rather than by \(\sigma\). This is the grid result again, from a different angle.
Then twelve lines of Metropolis: 3 000 iterations, each proposing a three-dimensional normal step with per-parameter standard deviations of 0.0004, 0.03 and 0.15, and accepting on the log rule. The first 500 draws are discarded as burn-in — the chain starts at an arbitrary point (here the sample mean, the sample sd and \(\nu = 3\)) and needs time to find the high-density region; draws made before it arrives are not from the posterior. The acceptance rate is 0.473 for this proposal. Very small steps can produce high acceptance yet poor exploration; very large steps can produce long rejection stretches. Acceptance alone is not a convergence or precision test. Check several dispersed chains, rank-normalised R-hat, effective sample size and Monte Carlo error before relying on a reported tail probability.
The posterior means are \(\mu\) 0.00145, \(\sigma\) 0.00974 and \(\nu\) 2.72. Take them in reverse order. A tail index of 2.7 says GOOG’s daily returns have heavier tails than a t(3); the notebook’s PyMC run found 2.79, and the posterior 95 % credible interval for \(\nu\) runs from about 2.2 to 3.4, with 84 % of the mass below 3. The scale \(\sigma\) is six-tenths of the sample standard deviation — the t model’s implied standard deviation \(\sigma\sqrt{\nu/(\nu-2)}\) averages about 0.020, above the sample sd of 0.0164, so the model has not lost any dispersion; it has moved it from the body to the tails. And \(\mu\) is nearly double the sample mean, for the reason already given: a Student-t likelihood discounts extreme negative days when locating the centre.
Read the trace before you trust it
A sampler always produces numbers; the question is whether they are draws from the posterior. The first and most important diagnostic is the trace plot — the sequence of draws against iteration number.
A converged chain looks like a fuzzy caterpillar: a flat band with no trend, no long horizontal stretches, and the same width from the start of the trace to the end. Each pathology has a signature. A smooth drift means the burn-in was too short and the chain is still travelling toward the posterior — discard more. Long flat segments with occasional jumps mean the step is too large and almost every proposal is rejected. A trace that moves but very slowly, with the caterpillar’s fuzz stretched into long waves, means the step is too small and successive draws are highly correlated; the chain will get there, but 2 500 draws may carry the information of 50 independent ones. Both traces here are flat bands from the first retained draw to the last, and the \(\nu\) trace makes visible the asymmetry of its posterior: a hard floor near 2.2 and occasional excursions above 3.5.
PyMC runs two or more chains from different starting points and reports \(\hat R\), the ratio of between-chain to within-chain variance, which should be 1.00; and the effective sample size, the number of independent draws the autocorrelated chain is worth. With one hand-rolled chain you have only the eye, which is why the trace comes before the histogram.
What NUTS does that Metropolis does not
The notebook’s pm.sample does not run the loop you just wrote. It runs the No-U-Turn Sampler (NUTS), an adaptive form of Hamiltonian Monte Carlo, and the difference is worth understanding because it explains the numbers in the notebook’s output.
Metropolis proposes blindly. The step is a random direction of a fixed size, and the algorithm learns nothing about the shape of the posterior from one proposal to the next. In three dimensions this is fine — acceptance 0.47, a clean trace. In thirty dimensions it is a disaster: a random step of any useful size is almost certain to leave the thin shell where the posterior mass lives, so it must be made tiny to be accepted, and then the chain crawls. The effective sample size per iteration collapses as dimension grows — roughly as \(1/d\) — which is why the random-walk sampler is not what any serious library uses by default.
Hamiltonian Monte Carlo replaces the blind step with a physical simulation. Treat the negative log-posterior as a potential energy surface, give the current point a random momentum, and let it roll along the surface for a while under Hamilton’s equations, using the gradient of the log-posterior to steer. The trajectory follows the contours of the posterior instead of cutting across them, so a proposal can travel a long way and still land in a high-density region; acceptance rates of 0.8 to 0.95 with distant proposals are routine, and the notebook’s acceptance_rate column shows values like 0.90, 0.87, 0.98 for exactly this reason. The cost is the gradient, which PyMC obtains by automatic differentiation of the model you wrote — one reason the model has to be expressed in PyMC’s own tensor language rather than in plain numpy.
NUTS adds the two pieces of tuning that HMC otherwise leaves to the user. It chooses the length of each trajectory automatically by extending it until the path starts to double back on itself (the “U-turn”), and during the tune=1500 warm-up iterations it adapts the step size to hit a target acceptance rate and estimates a mass matrix that rescales the parameters so that \(\mu\) (on a scale of \(10^{-3}\)) and \(\nu\) (on a scale of 1) are explored equally well. This is why the notebook needed no step vector: the 1 500 tuning draws are discarded like burn-in, but they also did the hand-craft this section did by trial. The result is visible in the notebook’s summary — effective sample sizes of 1 069 to 1 291 from 2 000 retained draws, so that each NUTS draw is worth more than half an independent draw, and \(\hat R = 1.00\) on every parameter. The 2 500 retained Metropolis draws are correlated. Their count alone cannot establish precision comparable to NUTS; calculate effective sample size and Monte Carlo error.
None of this changes what is being computed. NUTS and Metropolis target the same posterior when implemented correctly. Their reported means are numerically similar; separate convergence and Monte Carlo error checks are needed to explain their difference. What NUTS buys is efficiency and, through its divergence diagnostics (a trajectory whose simulated energy is not conserved signals a region the sampler cannot handle — the divergences column in the notebook’s trace), a warning system Metropolis lacks. When you read pm.sample, read it as “run NUTS, an adaptive Hamiltonian Monte Carlo method with gradient-guided proposals”.
The posterior of the reward-to-t-scale ratio
The 2 500 retained draws of \((\mu, \sigma)\) give 2 500 draws of the Sharpe ratio \(\sqrt{252}\,\mu/\sigma\) — the array SR already plotted as a trace; this is reward-to-t-scale, not return-sd Sharpe. This is the object the section was built to produce. A derived quantity’s posterior costs nothing extra: any function of the parameter draws is a draw from that function’s posterior, which is why the notebook could define sharpe as a pm.Deterministic and get its distribution for free.
Posterior mean 2.376, 95 % credible interval [1.304, 3.594], and P(Sharpe > 1) = 0.988. The notebook’s PyMC run gave 2.35 and [1.31, 3.42] — the same answer from a twelve-line sampler as from NUTS. The red line is the conventional Sharpe estimate of 0.77. It cannot be tested against a credible interval for a differently defined ratio. These are differently defined ratios. Their difference does not prove downward bias in classical Sharpe estimation. Use the return-standard-deviation basis for a like-defined conventional Sharpe comparison.
Why is the Bayesian number three times the classical one? It is not the prior on \(\mu\): with a prior standard deviation of one daily sd and 1 258 observations the prior has essentially no pull, and the posterior mean of \(\mu\) (0.00145) is driven by the t-likelihood, not by where the prior is centred. It is not a bias in MCMC — the grid, with no sampler at all, gave 2.04 under \(\nu = 4\). It is not a forgotten \(\sqrt{252}\). It is the denominator. The classical Sharpe ratio divides by the sample standard deviation, which a handful of crash days inflate. The t model splits dispersion into a scale \(\sigma\) for ordinary days and a tail parameter \(\nu\) for extreme days; \(\sigma\) is 0.0097 against a sample sd of 0.0164, so \(\mu/\sigma\) is larger even before \(\mu\) itself has moved.
Be honest about what σ means
The t-distribution’s standard deviation is \(\sigma\sqrt{\nu/(\nu-2)}\), which for \(\nu \approx 2.7\) is almost twice \(\sigma\). Re-computing the Sharpe ratio on that basis gives a posterior mean of about 1.18 — the cell below does it. The t-scale ratio is a distinct estimand; it should be labelled reward-to-scale rather than conventional Sharpe. A fund that reports Sharpe ratios on a standard-deviation basis to investors but sizes positions on the t-scale is quoting two different risks; decide which one you are pricing before you quote a Sharpe ratio, and say which it is.
Report a standard-deviation Sharpe on a common basis
For Student-t returns with \(\nu>2\), the return sd is \(\sigma\sqrt{\nu/(\nu-2)}\). The ratio \(\sqrt{252}\mu/\sigma\) is a reward-to-scale ratio; it is not the conventional Sharpe ratio.
The mean is 1.178, credible interval [0.575, 1.852], and probability above one about 0.70. The difference from the sample ratio does not prove that classical statistics systematically understates performance. Risk-free return is set to zero; annualisation assumes independent, stable daily returns.
A trace is a diagnostic, not proof of convergence
A visually stable trace and an acceptance rate of 0.47 do not prove that all relevant posterior regions were explored. Run dispersed chains, inspect rank-normalised R-hat, bulk/tail effective sample size, Monte Carlo standard errors and sampler warnings. Check the derived investment quantity as well as individual parameters.
The single-chain browser sampler teaches Metropolis mechanics. It is not a validated computational pipeline for allocating capital. Normal conjugate calculations in §§5.5–5.6 avoid Monte Carlo error under their stated assumptions.
Credible interval versus confidence interval
The classical procedure has a standard error for the annualised Sharpe ratio too. Under i.i.d. normal returns, the usual approximation is that the standard error of the daily Sharpe ratio \(\widehat{SR}_d\) is \(\sqrt{(1 + \tfrac12 \widehat{SR}_d^2)/n}\), so the annualised standard error is \(\sqrt{252\,(1 + \tfrac12 \widehat{SR}_d^2)/n}\). The cell computes the 95 % confidence interval, sets it beside the credible interval, and adds the standard-deviation-basis posterior promised in the box above.
The classical 95 % confidence interval is [-0.107, 1.648]. It includes zero: the stated approximation does not reject a zero population mean. This does not imply the realised investment earned nothing. The reward-to-t-scale credible interval is [1.304, 3.594], and the posterior probability that reward-to-t-scale exceeds 0.5 is 0.999. The two intervals are not just numerically different; they are different kinds of statement. A confidence interval describes a procedure: in repeated samples from the same process, 95 % of the intervals so constructed would cover the fixed true value — this particular interval either does or does not, and the 95 % is not a probability about it. A credible interval is a probability statement about the parameter given this data set. Only the latter licenses the sentence “there is a 95 % probability the true Sharpe ratio lies in here”, under the specified likelihood and prior. The conventional Sharpe comparison must use the sd-basis posterior below, because the t-scale interval concerns a different estimand.
The last line is the honest number. On the standard-deviation basis the posterior mean Sharpe ratio is 1.178, with a 95 % credible interval of about [0.58, 1.85] and P(Sharpe > 1) = 0.70. That is still well above the classical 0.77 — the Student-t likelihood’s different estimate of \(\mu\) contributes to the remaining difference — but it is a different order of claim from 2.38. Both numbers are correct answers to different questions. The reward-to-t-scale ratio divides by a distribution scale; it is not a complete measure of ordinary-day or tail risk; the sd-basis Sharpe ratio measures it against all days including the crashes. A risk manager who cares about drawdowns wants the second; position sizing also requires the predictive distribution and tail scenarios; a marketing deck that quotes the first without saying so is misleading.
Your turn: AAPL and NVDA Sharpe posteriors
The notebook’s practice asks for the same posterior on other stocks. Wrap the sampler in a function and run it on AAPL (ten years, 2 515 returns) and on NVDA over 2023–2024 (500 returns). The step sizes are scaled to the data — the proposal sd for \(\mu\) is \(1.5\,s/\sqrt{n}\), roughly the posterior sd — so the same code works on both.
The following posterior numbers use reward-to-t-scale. AAPL: conventional Sharpe 0.922, reward-to-scale posterior mean 1.531, credible interval [0.877, 2.209], P(Sharpe > 1) = 0.922. NVDA: classical 2.509, posterior mean 2.609, credible interval [1.105, 4.139], P(Sharpe > 1) = 0.985. Two things to notice. The gap between classical and Bayesian is large for AAPL (whose posterior \(\nu\) is about 3.4 — fat tails, so the t-scale is well below the sd) and small for NVDA (posterior \(\nu\) about 5.6 over a period without a crash, so the t-scale and the sd nearly coincide). And NVDA’s credible interval is three times wider than AAPL’s despite the higher mean, because 500 observations pin a Sharpe ratio down far less tightly than 2 515: a two-year track record, however spectacular, cannot rule out a Sharpe ratio near 1. That last sentence is the kind of statement the posterior exists to make.
It accepts with probability 0.3 (uphill moves are always accepted; downhill ones with probability \(p^\star/p_t\)). Only the ratio of densities appears, so \(p(X)\) cancels; in code the test is log(u) < log_post(prop) − log_post(current).
The classical Sharpe divides by the sample sd, which crash days inflate; the t model absorbs those days in \(\nu \approx 2.7\) and its scale \(\sigma\) (0.0097) is six-tenths of the sample sd (0.0164). Dividing instead by \(\sigma\sqrt{\nu/(\nu-2)}\) gives a posterior mean of about 1.18 — still above 0.77, because the t-likelihood also gives a different estimate of \(\mu\), but a different claim.
Only the credible interval. A confidence interval is a statement about a procedure across repeated samples — this one either covers the fixed truth or does not. A credible interval is a probability statement about the parameter given this data. Here the classical interval could not even exclude zero.
It uses the gradient of the log-posterior to propose along the contours (Hamiltonian dynamics), so distant proposals are accepted 80–95 % of the time even in high dimension; and it tunes step size, trajectory length and a mass matrix automatically during warm-up, then reports \(\hat R\), effective sample size and divergences. Both target the same posterior.
Bayesian Regression and Predictive Checks
| Part | Question |
|---|---|
| 5.2.1 Regression as a generative model | What observations could the assumed regression generate? |
| 5.2.2 Prior predictive checks | Do the priors imply plausible observations? |
| 5.2.3 Posterior coefficients | How do prior and data precision combine? |
| 5.2.4 Posterior predictive checks | Does the fitted model reproduce relevant data features? |
Regression as a generative model
Question: What observations could the assumed regression generate?
Result: Separate coefficients, residual noise and generated observations.
Begin with the definitions and assumptions; then trace the worked example.
Regression from the normal-distribution perspective
Chapter 4 wrote a regression as \(y_i = \beta_0 + \beta_1 x_i + \varepsilon_i\) with \(\varepsilon_i\) a zero-mean error. A Bayesian writes the same model as a statement about the distribution of \(y\):
\[ y_i \sim N(\mu_i,\ \sigma^2), \qquad \mu_i = \beta_0 + \beta_1 x_i . \]
Nothing has changed numerically, but the emphasis has. The response is a normal random variable whose mean moves with \(x\) and whose spread \(\sigma\) does not. There are three unknowns — \(\beta_0\), \(\beta_1\) and \(\sigma\) — and a Bayesian puts a prior on each of them. The \(x\) values are treated as fixed design points; they get no prior. The residuals \(\varepsilon_i\) are random but are not parameters — they are integrated out by the likelihood, which is the product over \(i\) of the normal density of \(y_i\) at \(\mu_i\). Once the model is phrased this way, everything from Section 5.1 applies without modification: the posterior of \((\beta_0, \beta_1, \sigma)\) is prior times likelihood, it can be sampled by Metropolis or NUTS, and any function of the draws — the fitted line at a new \(x\), the probability that the slope is positive — has a posterior for free.
The notebook uses artificial data with known parameters, and so does this section, because the point of the exercise is to watch a model recover truths you planted. The first data set has \(\beta_0 = 1\), \(\beta_1 = 2\) and \(\sigma = 0.5\) on 200 points with \(x\) evenly spaced on \([0, 1]\).
The mean of \(y\) is 2.041 — the true line’s average height, \(1 + 2 \times 0.5 = 2\), plus noise — and the standard deviation of \(y\) is 0.759. That second number is larger than \(\sigma = 0.5\), and the reason is a distinction that will matter twice in this section: part of the spread of \(y\) is the line itself. The variance of \(y\) is the variance of \(\mu(x)\) across the design points, \(2^2 \times \text{Var}(x) = 4 \times 1/12 = 0.333\), plus the noise variance \(0.25\), giving \(\sqrt{0.583} = 0.76\). The spread of \(y\) is not the noise. A model that confuses the two — a prior on \(\sigma\) centred on the sd of \(y\), say — will be wrong by a predictable amount.
In Colab — the notebook’s regression in PyMC
The notebook fits this with NUTS; the posterior comes back as a table with one row per draw and one column per parameter, and every plot in this section is built from that table.
with pm.Model() as linear_regression1:
std = pm.HalfCauchy("sigma", beta=10, initval=1)
beta0 = pm.Normal("intercept", 0, sigma=20)
beta1 = pm.Normal("slope", 0, sigma=20)
likelihood = pm.Normal("y", mu=beta0 + beta1 * data["x"], sigma=std,
observed=data["y"])
trace1 = pm.sample(tune=6000, draws=2500, chains=2)
trace1.to_dataframe().head()The notebook’s first five posterior rows have intercepts 0.985, 0.972, 0.975, 0.961, 0.985, slopes 2.07, 2.03, 2.05, 2.03, 2.02 and sigmas 0.508, 0.503, 0.509, 0.513, 0.515 — the planted 1, 2 and 0.5, each with its uncertainty. The half-Cauchy prior on \(\sigma\) is a weakly informative choice for a scale: positive, heavy-tailed so that large values are not ruled out, but with most of its mass at moderate values.
Prior predictive checks
Question: Do the priors imply plausible observations?
Result: Compare simulated lines with the problem’s measurement scale.
Builds on Regression as a generative model.
Prior predictive check
Before fitting, draw parameters from the prior and plot the regression lines they imply. If the lines are absurd, the prior is absurd — and you have found out for free, before the data have had any say. The coefficient fan is a prior check on regression means. A full prior predictive check additionally specifies a likelihood, draws its noise scale and simulates observations; it needs no fitted posterior. Do not confuse the two.
The notebook’s second data set makes the point vividly: \(x \in [-10, 10]\), \(y = 0.5 + 3x + N(0, 6^2)\), 200 points. The sd of \(x\) is 5.8, so the line alone contributes \(3 \times 5.8 = 17.4\) to the sd of \(y\) and the noise contributes 6; the sd of \(y\) should be about \(\sqrt{17.4^2 + 36} = 18.4\). The data at \(x = 10\) sit near \(0.5 + 30 = 30.5\).
Prior I puts \(N(0, 10^2)\) on both coefficients — a “flat” prior in the sense that it is meant to be uninformative. Draw 100 pairs \((\beta_0, \beta_1)\) from it and plot the lines. Before running, think about the line at \(x = 10\): its height is \(\beta_0 + 10\,\beta_1\), a normal with standard deviation \(\sqrt{10^2 + 100 \times 10^2} \approx 100\). What range will cover 95 % of the lines there?
Slope sd 10, and \(x\) reaches 10. At \(x = 10\), roughly what interval covers 95 % of the prior lines — about ±3, ±10, ±30 or ±200?
The data have sd 18.37, as predicted. The prior lines at \(x = 10\) span [-173.8, 209.7] — roughly ±200, where the data live inside ±35. Prior I is not wrong; the true line is inside its support and with 200 points the posterior will find it. But it spends the overwhelming majority of its probability on lines no analyst would believe for a moment: slopes of −15, intercepts of +20. A prior is supposed to encode what you know, and you know more than this. The cyan fan in the figure, most of which is nowhere near the black points, is the prior telling you it has not been thought about.
Prior II tightens both: \(N(0, 1)\) on the intercept and \(N(0, 1.2^2)\) on the slope.
Now the 95 % range at \(x = 10\) is [-20.8, 25.2]. The fan covers slopes of plausible size without covering nonsense. This is a regularising prior, and it is worth being candid about its one weakness: the true slope of 3 sits at \(3/1.2 = 2.5\) prior standard deviations, near the edge of what the prior considers likely. With 200 points the data will drag the posterior there regardless — you will see the slope come out at 3.000 in a moment. With 10 points they might not, and the posterior slope would be shrunk toward zero by an amount you could compute from the precisions. Neither prior is “the right one”; the check has told you what each of them assumes, which is the purpose. A prior you have seen is a prior you can defend to a sceptical portfolio manager.
Posterior coefficients
Question: How do prior and data precision combine?
Result: A coefficient posterior with its covariance retained.
Builds on Prior predictive checks.
Posterior in closed form
With a normal likelihood and a normal prior on the coefficient vector \(\beta = (\beta_0, \beta_1)\), the posterior is normal and can be written down — the multivariate version of the conjugate update from Section 5.1. Take \(\sigma\) as known for the moment (plug in the OLS residual sd; the notebook’s NUTS run estimates it jointly), let \(X\) be the design matrix with a column of ones and a column of \(x\), and let the prior be \(\beta \sim N(\mu_0, \Sigma_0)\). Then
\[ \Sigma_n = \left(\Sigma_0^{-1} + \tfrac{1}{\sigma^2} X^\top X\right)^{-1}, \qquad \mu_n = \Sigma_n\left(\Sigma_0^{-1}\mu_0 + \tfrac{1}{\sigma^2} X^\top y\right). \]
The structure is exactly the scalar formula in matrix dress. \(\Sigma_0^{-1}\) is the prior precision matrix and \(X^\top X/\sigma^2\) is the data precision matrix — the inverse of the classical covariance of \(\hat\beta_{\text{OLS}}\) — and they add. The posterior mean is the precision-weighted average of the prior mean \(\mu_0\) and the OLS estimate (note that \(X^\top y/\sigma^2 = (X^\top X/\sigma^2)\,\hat\beta_{\text{OLS}}\)). When the prior precision is negligible next to the data’s, \(\Sigma_n \to \sigma^2 (X^\top X)^{-1}\) and \(\mu_n \to \hat\beta_{\text{OLS}}\): Bayes with a flat prior is OLS, with the classical covariance matrix reinterpreted as a posterior covariance. When the prior precision is comparable, the posterior is pulled toward \(\mu_0\) — this is ridge regression, which Chapter 4 introduced as a penalty and which is, exactly, the posterior mode under a normal prior centred at zero.
Before computing, predict. OLS on (xt, yt) gives a slope of 3.011. Prior II’s precision on the slope is \(1/1.2^2 = 0.69\). The data’s precision on the slope is \(\sum x_i^2/\sigma^2\), and \(\sum x_i^2 = 6\,734\) while \(\sigma^2 \approx 32.5\), so about 207. Two hundred and seven against 0.69: the posterior mean slope will be within a rounding error of 3.011.
OLS gives [0.385, 3.011] with \(\hat\sigma\) = 5.705 (the truth is 0.5, 3, 6); the posterior mean is [0.331, 3.000] with posterior sd [0.374, 0.069]. The slope barely moved — 3.011 to 3.000, a shift of one sixth of its posterior sd, which is the 0.69-versus-207 arithmetic made visible. The intercept moved more, from 0.385 to 0.331, because its prior is relatively stronger: the data precision on the intercept is \(n/\sigma^2 = 200/32.5 = 6.1\) against a prior precision of 1, so the prior holds about 14 % of the weight and pulls the intercept 14 % of the way toward 0. With \(n = 200\) the prior is a rounding error on the slope and a small correction on the intercept; with \(n = 20\) it would be a tenth of the story on both.
Sampling the posterior: one row is one line
The closed form gives a mean vector and a covariance matrix. To do anything beyond quoting those two — to plot lines, to compute the probability that the slope exceeds 2.9, to build the predictive check that follows — you want draws. Because the posterior is a bivariate normal, rng.multivariate_normal(mn, Sn, size=1000) produces them directly; each row of the resulting array is one pair \((b_0, b_1)\), which is to say one regression line. This is the same object the notebook’s trace.to_dataframe() returns, obtained without a sampler because the model is conjugate.
The first five rows of post_df are five plausible lines: slopes of 2.972, 2.917, 3.178, 3.008, 3.021 with intercepts scattered around 0.3. The 50 red lines on the figure are 50 such rows, and they form a narrow bow-tie pinched at the centre of the data — the posterior is most certain about the line’s height where the design points are densest and least certain at the extremes, where the slope uncertainty has the most leverage. The 95 % credible interval for the slope is [2.861, 3.139]: the same information a classical confidence interval \(3.011 \pm 1.96 \times 0.069\) would give, but readable as a probability — there is a 95 % chance the slope lies between 2.86 and 3.14, and the probability it is below 2.9 is about 5 %.
Posterior predictive checks
Question: Does the fitted model reproduce relevant data features?
Result: A predictive check including both parameter and observation uncertainty.
Builds on Posterior coefficients.
Posterior predictive check
Now go one step further. A posterior draw is a set of parameters; a set of parameters plus the likelihood is a recipe for generating data. For each posterior draw \(s\), compute \(\mu_i^{(s)} = b_0^{(s)} + b_1^{(s)} x_i\) at every design point and then draw a whole replicated data set \(y_i^{\text{rep},(s)} \sim N(\mu_i^{(s)}, \sigma)\) for every \(i\). Do this for all 1 000 draws and you have 1 000 fake data sets, each of which the fitted model considers as plausible as the real one. If the replicas do not look like the data — in their spread, their extremes, their skew, some feature you care about — the model is missing something, and the posterior predictive check has found it.
The shape of the resulting array is the thing to get straight. yrep will be (1000, 200): rows index posterior draws, columns index observations. yrep[7] is not 1 000 predictions for the seventh observation; it is one complete simulated data set of 200 points, generated from posterior draw number 7 — one cyan curve in the notebook’s az.plot_ppc figure. yrep[:, 7] would be the 1 000 replicated values of the seventh observation. Getting the axes wrong is the most common bug in a hand-rolled PPC, and the debug exercise below is built around it.
Thirty cyan replica histograms and the black data histogram share one shape — wide, roughly flat-topped (because \(y\) is dominated by the uniformly spaced line rather than by the noise), tailing off at about ±35. The mean of replica 0 is 0.82 against an observed mean of 0.39; the difference is noise of the size the model expects, since the sd of a 200-point mean of \(y\) is \(18.4/\sqrt{200} = 1.3\). If the data had a bump the replicas never produce, or a tail they never reach, the normal likelihood would be wrong — and this is exactly how the outliers of Section 5.3 will announce themselves.
Turn the picture into a number
A picture with thirty curves is persuasive but not quantitative. The standard summary is the Bayesian p-value: choose a test statistic \(T\) — the mean, the maximum, the standard deviation, anything — compute it for every replica and for the observed data, and report the fraction of replicas whose \(T\) exceeds the observed value. A value near 0.5 says the model produces data like yours on that dimension about as often as not; a value near 0 or 1 says the observed statistic sits in a tail of what the model can generate, and the model cannot reproduce that feature of the data.
The cell below computes the p-value for the mean correctly and the p-value for the maximum incorrectly. It runs without error and prints a number. Read it before you run it, decide what is wrong, and predict what the wrong line will print.
The second line compares the mean of each replica row with the maximum of the observed data. A replica’s mean is around 0.4; the observed maximum is 42.2; no replica mean will ever exceed it, so p_max prints 0.0 — a value that, taken at face value, would say the model is catastrophically unable to reproduce the data’s maximum. The bug is a copy-and-paste of .mean(axis=1) where .max(axis=1) was needed: the statistic must be computed on each replica in the same way it is computed on the observed data. The fix is one word, and the corrected version — together with two further statistics — follows.
The p-value for the mean is 0.454 and for the standard deviation 0.457 — the model reproduces the centre and the spread of the data about as often as not, which is what a correctly specified normal model should do. The p-value for the minimum is 0.583, fine. The p-value for the maximum is 0.135: only 13.5 % of replicas produce a maximum as large as the observed 42.21, whose typical replica counterpart is 38.49. That is not a rejection — a p-value of 0.13 is well inside the range chance can produce — but it is the first hint of a pattern. A normal likelihood generates extremes reluctantly, and the observed data, even these normally generated data, happen to have a maximum that the fitted model finds slightly surprising. On real financial data the same statistic would come out near 0.00, and that is the signal that the tails need a heavier distribution. Section 5.3 supplies one.
\(\beta_0\), \(\beta_1\) and \(\sigma\) receive priors. The \(x_i\) are fixed design points and the residuals \(\varepsilon_i\) are not parameters — they are integrated out by the likelihood.
yrep[7] is one complete simulated data set of 200 points generated from posterior draw number 7 — one replica curve. yrep[:, 7] is the 1 000 replicated values of the seventh observation. A Bayesian p-value compares a statistic of each row with the same statistic of the observed data.
Robust Regression with Fat Tails
| Part | Question |
|---|---|
| 5.3.1 Outlier sensitivity and likelihood choice | Why can a few observations dominate a Gaussian fit? |
| 5.3.2 Robust posterior uncertainty | How uncertain is the robust slope? |
| 5.3.3 Market beta and rolling robustness | When do ordinary and robust market associations differ? |
Outlier sensitivity and likelihood choice
Question: Why can a few observations dominate a Gaussian fit?
Result: Compare squared-error sensitivity with a Student-t likelihood.
Begin with the definitions and assumptions; then trace the worked example.
Plant the outliers
OLS assumes normal errors, and a normal error distribution has a specific opinion about large residuals: they are almost impossible. That opinion is what breaks the estimator when the data disagree. Under a normal likelihood the log-density of a residual of size \(r\) standardised units is \(-r^2/2\) plus a constant, so a residual of 7 costs \(7^2/2 = 24.5\) log-units — as much as sixty ordinary residuals of size 0.9. The fit will bend a long way to reduce that one term, at everyone else’s expense. When the errors really are normal this never arises, because residuals of 7 do not occur. When the errors are fat-tailed — as they are for every asset return series in this book — a handful of extreme observations can dominate the estimate, pull \(\hat\beta\) around and, worse, produce standard errors that describe a world that does not exist.
The notebook makes this visible by planting outliers where you can see them. The base data are \(y = 0.5 + 5x + N(0, 0.5^2)\) on 100 points with \(x\) on \([0, 1]\). Three more points are added at \((0.1, 8)\), \((0.15, 6)\) and \((0.2, 9)\) — far above the line, at small \(x\). Before running, decide how OLS will respond. The outliers are high and to the left; to reduce their residuals the fitted line must be higher there, and the cheapest way to do that with a straight line is to raise the intercept and tilt the slope downward. So: intercept up, slope down.
OLS gives [1.033, 4.403]. The intercept has doubled, from 0.5 to 1.03, and the slope is off by 0.6, from 5 to 4.40. Three points out of 103 — under 3 % of the sample — have moved the slope by 12 %, and nothing in the OLS output will tell you so: the \(R^2\) is still high, the t-statistics are still enormous, and the residual plot will show three large positive residuals that a careless reader could mistake for noise.
In Colab — normal and Student-t likelihoods in PyMC
with pm.Model() as model:
beta0 = pm.Normal("b0", 0.0, 1.0); beta1 = pm.Normal("b1", 0, 1.2)
std = pm.Exponential("sd", 1.0)
output = pm.Normal("y", mu=beta0 + beta1 * data2["x"], sigma=std, observed=data2["y"])
with pm.Model() as robust_model:
beta0 = pm.Normal("b0", 0.0, 1.0); beta1 = pm.Normal("b1", 0, 1.2)
std = pm.Exponential("sd", 1.0)
mu = beta0 + beta1 * data2["x"].values
likelihood = pm.StudentT("y", mu=mu, sigma=std, nu=3, observed=data2["y"].values)
robust_trace = pm.sample(tune=2000)Same priors (Prior II from Section 5.2, and an exponential on the scale), one word changed: Normal → StudentT with nu=3. The notebook’s first posterior rows under the normal likelihood have slopes of 4.02, 4.24, 3.47, 3.40, 3.35; under the Student-t they are 5.06, 5.06, 4.40, 4.39, 4.56. The rest of this section shows what that one word does.
Student-t likelihood by maximum likelihood
You do not need a sampler to see the effect of the likelihood; you need only the likelihood itself. Replace the normal density with a Student-t with \(\nu = 3\), write the negative log-likelihood as a function of \((\beta_0, \beta_1, \log\sigma)\), and minimise it with scipy.optimize.minimize. This is maximum likelihood under a t error model — the posterior mode under a flat prior — and it takes a few milliseconds.
The robust fit gives intercept 0.492, slope 5.167 and scale 0.457; the truth is 0.5, 5 and 0.5. The green dashed line goes through the cloud, essentially on top of the dotted truth; the red OLS line is tugged up and to the left by three points. The slope is now off by 0.17 rather than 0.6, and in the other direction — with only 100 clean points and noise sd 0.5 the sampling error in the slope is about \(0.5/\sqrt{100 \times 1/12} = 0.17\), so the robust estimate is within one standard error of the truth, as it should be. Nobody deleted a row, nobody decided which points were outliers; the likelihood did it.
Why does the t-likelihood ignore them?
The mechanism is worth deriving because it turns “robust” from a slogan into a formula. The log-density of a Student-t residual \(r\) (in standardised units) is \(-\frac{\nu+1}{2}\log(1 + r^2/\nu)\) plus a constant. Differentiate with respect to the parameters and the score contribution of observation \(i\) is \(w_i\,r_i\,\partial r_i/\partial\beta\) with
\[ w_i = \frac{\nu + 1}{\nu + r_i^2}. \]
Compare the normal: its log-density is \(-r^2/2\), its score is \(r_i\,\partial r_i/\partial\beta\), and so its weight is identically 1. The t-likelihood is therefore a self-weighting least squares: it solves the OLS normal equations with each observation weighted by \(w_i\), where \(w_i\) is near \((\nu+1)/\nu\) — slightly above 1 — for a residual of ordinary size and falls toward zero like \((\nu+1)/r_i^2\) as the residual grows. The weights depend on the residuals, which depend on the fit, so the solution is an iteration (this is what minimize is doing, and what the classical “iteratively reweighted least squares” M-estimators of Huber do explicitly), but the interpretation is immediate: the t-likelihood down-weights exactly the points that do not fit, by an amount that grows with their surprise. A smaller \(\nu\) makes the weights fall faster; \(\nu \to \infty\) recovers OLS.
The outlier at \((0.2, 9)\) has a standardised residual of about 16.5 under the robust fit. With \(\nu = 3\), what weight does it receive (approximately) — and roughly what weight does a typical point receive?
The outliers receive weights 0.017, 0.036 and 0.015 — for the point at \((0.2, 9)\), \(4/(3 + 16.5^2) = 0.015\) — while a typical point receives 1.128. Each planted point counts for about a sixtieth of an ordinary one. That is why the robust slope is 5.17 instead of 4.40: the three points are still in the data, still in the likelihood, but their vote has been scaled by the model’s assessment of how implausible they are under the fitted line. A risk manager will object at this point that “those are exactly the days I care about”, and the objection is right for one purpose and wrong for another; the GOOG example below returns to it.
Robust posterior uncertainty
Question: How uncertain is the robust slope?
Result: Posterior slope draws with the likelihood and prior stated.
Builds on Outlier sensitivity and likelihood choice.
Posterior of the robust slope by Metropolis
Maximum likelihood gives a point; the posterior gives the point and its uncertainty. The same twelve-line sampler from Section 5.1 works with a new log-posterior: Prior II on the coefficients, \(\sigma \sim \text{Exponential}(1)\) (whose log-density on \(\log\sigma\) is \(-\sigma + \log\sigma\), the Jacobian again), and the t(3) likelihood. Start the chain at the maximum-likelihood point, which removes most of the burn-in problem.
Acceptance 0.337, posterior mean slope 5.059, 95 % credible interval [4.644, 5.409]. The posterior mean is a little below the maximum-likelihood 5.167 because Prior II pulls the slope toward zero — the true slope of 5 sits at 4.2 prior standard deviations, and with 103 points the prior still has a small say. More important is what the interval excludes: the OLS slope of 4.40 lies outside it. Under the t model, given these data, the probability that the slope is as low as OLS says is below 2.5 %. The model is not merely giving a different point estimate; it is telling you that the OLS answer is inconsistent with the bulk of the data, and that the t likelihood reduces the influence of their large residuals. A posterior does not establish whether real observations are errors or belong to a different process.
Market beta and rolling robustness
Question: When do ordinary and robust market associations differ?
Result: Compare fitted betas and a rolling extension on matched returns.
Builds on Robust posterior uncertainty.
Real data: GOOG’s beta on the S&P 500
Simulated outliers are a demonstration; the March 2020 crash is the real thing. Regress GOOG’s daily return on the S&P 500’s over the same 1 258 days used in Section 5.1 — a market-model regression whose slope is the stock’s beta — first by OLS, then with a Student-t likelihood with \(\nu\) estimated from the data rather than fixed at 3.
Predict the direction first. OLS will give \(\beta \approx 1.06\). Will the robust beta be higher or lower? The residuals here will turn out to have excess kurtosis above 11 — extremely fat tails — so the two fits will differ. But whether the robust beta is higher or lower depends on where the extreme days sit relative to the line, which you cannot know in advance: if the crash days lie above the OLS line the robust fit will steepen, if below it will flatten. And with 1 258 days the shift will be modest either way. Heavy residual tails motivate checking inference assumptions; they do not by themselves prove the OLS standard error is wrong. Classical coefficient standard errors require an appropriate variance model, and normality matters for exact small-sample t inference. Use Chapter 4’s heteroskedasticity/dependence diagnostics and covariance estimators as appropriate.
OLS: alpha 0.0002, beta 1.0636, \(R^2\) 0.589, and the excess kurtosis of the residuals is 11.4. The market explains 59 % of GOOG’s daily variance, which is typical of a mega-cap; the residual kurtosis is not typical of anything normal. A normal distribution has excess kurtosis 0; a t(4) has infinite kurtosis; 11.4 in a finite sample says the residuals are as fat-tailed as the raw returns were in Section 5.1, and motivates checking homoskedasticity, dependence and finite-sample inference before interpreting the OLS standard error on beta (0.025, giving a t-statistic of 42).
The robust fit gives beta 1.0937 with \(\nu\) = 3.35, against OLS 1.0636. The estimated tail index is close to the 2.7 that the return model found in Section 5.1 — the market factor explains some of GOOG’s variance but none of its kurtosis, and the residuals are as heavy-tailed as the returns. The two betas differ by 0.03, which is about one OLS standard error: modest, as predicted. The direction — robust beta higher — comes from the geometry of the crash: on the biggest down days GOOG fell somewhat less than the market, so those points sit above the line on the far left, and OLS, which must accommodate them, tilts the slope down; the t model gives them small weights and lets the slope be set by the other 1 246 days.
The table shows the twelve days with \(|r_{SP}| > 5\,\%\): one in December 2018, nine in March 2020, one in April and one in June 2020. They are 0.95 % of the sample, and they carry 0.365 — 36.5 % — of \(\sum x_i^2\), which is the quantity that determines leverage in OLS. One per cent of the days supplies more than a third of the information OLS uses to set the slope. The t model weights each of those days by its surprise relative to the fitted line rather than by its size: 12 March 2020, when GOOG fell 8.3 % against the market’s 9.6 %, gets a weight of 0.32; 16 March, when both fell 11 %, sits close to the line and gets 0.89, nearly full weight. The robust regression does not throw the crash away; it stops the crash from being the only thing that matters.
This is where the risk manager’s objection has to be answered honestly. For estimating the systematic exposure — the beta you would hedge with — the robust estimate is better: it is the slope of the relationship on the days that constitute 99 % of the sample, and its standard error is not fiction. For estimating what happens in a crash, the robust fit is the wrong tool by construction: it has been told to discount precisely those days. The two questions want two models, and the mistake is to answer the second with the first. A t-likelihood is the right tool for alpha and hedge ratios; for tail risk you go back to Chapter 2’s extreme value theory, or to the conditional volatility models of Chapter 6.
Rolling robust regression
The notebook’s practice asks for a rolling comparison. Take 250-day windows stepped every 25 days — 41 windows across the sample — fit the OLS slope and the t(3) slope in each, and compare the paths. The cell wraps nll_t so that it takes the window’s data as arguments.
In the 2020 windows the two paths are close: on 3 February 2020, before the crash entered the window, OLS says 1.232 and the robust fit 1.163; by 15 April, with the crash inside the window, OLS has fallen to 0.963 and the robust fit sits at 1.004; by October both are near 0.96. The largest gap in the whole sample is 0.112, on 24 May 2018 — a window containing the February 2018 volatility spike and a few GOOG-specific earnings days. Two lessons. Beta is not a constant: it moves from 1.2 to 0.95 within a year on either estimate, which is a preview of Section 5.4’s argument that a parameter should be allowed to wander. And the gap between the estimators is largest in exactly the windows a hedger cares about — the ones with a shock in them — so the choice of likelihood is not an academic nicety but a difference of 0.05 to 0.1 in a hedge ratio.
\(w_i = (\nu+1)/(\nu + r_i^2)\); for \(\nu = 3\) and \(r_i = 16.5\) it is \(4/275 = 0.015\). A typical point (\(r_i \approx 0.9\)) gets about 1.1. The normal likelihood gives every point weight 1, which is why OLS chases outliers.
Both, for different questions. For the systematic exposure (hedge ratio, alpha) the robust beta is better: its slope and standard error describe the 99 % of days that are not a crash. For what happens in a crash the robust fit is the wrong tool by construction; tail risk needs EVT (Chapter 2) or conditional volatility (Chapter 6).
Sequential Learning and Adaptive Risk
| Part | Question |
|---|---|
| 5.4.1 Updating a fixed unknown quantity | How does genuinely new evidence update a belief? |
| 5.4.2 Allowing a latent quantity to change | What if the underlying quantity changes over time? |
| 5.4.3 Chronological adaptive-risk evaluation | Was the risk estimate available before the decision? |
Updating a fixed unknown quantity
Question: How does genuinely new evidence update a belief?
Result: A sequential posterior that counts each observation once.
Begin with the definitions and assumptions; then trace the worked example.
Yesterday’s posterior is today’s prior
Every posterior so far used all the data at once. Markets do not arrive that way; they arrive one day at a time, and refitting from scratch and sequential updating are alternative calculations with different computational and data-bookkeeping requirements. Bayes’ theorem has a property that solves this directly. Because the likelihood of two independent observations is the product of their likelihoods, the posterior after seeing \(r_1\) and \(r_2\) is proportional to \(p(r_2 \mid \theta)\,p(r_1 \mid \theta)\,p(\theta)\) — which is the likelihood of \(r_2\) times the posterior after \(r_1\). Yesterday’s posterior is today’s prior. The order of the observations does not matter, and updating one at a time gives exactly the same answer as updating all at once.
For a normal mean with known \(\sigma\) the update is the conjugate formula of Section 5.1 applied to a single observation. Let \((m_{t-1}, v_{t-1})\) be the posterior mean and variance of \(\mu\) after \(t-1\) days; the new return \(r_t\) updates them to
\[ \frac{1}{v_t} = \frac{1}{v_{t-1}} + \frac{1}{\sigma^2}, \qquad m_t = v_t\left(\frac{m_{t-1}}{v_{t-1}} + \frac{r_t}{\sigma^2}\right). \]
Precision adds \(1/\sigma^2\) per observation; the mean is the precision-weighted average of the old mean and the new return. After \(n\) observations the precision is the prior precision plus \(n/\sigma^2\) regardless of the order in which the observations arrived, so after 253 daily updates with \(\sigma = 0.02\) and a prior sd of 0.002 the posterior sd will be \(\sqrt{1/(1/0.002^2 + 253/0.02^2)} = 0.00106\). Not \(0.002/253\), and not \(0.02/\sqrt{253} = 0.00126\) either — the latter is what a flat prior would give, and the informative prior tightens it a little.
The data are AAPL’s daily returns in 2020, a year that contains the March crash and the recovery: 253 trading days with mean 0.00277 and sd 0.02939. Fix \(\sigma = 0.02\) (a known variance is the price of a closed form; the true sd is higher because of the crash, and the consequence will be visible) and start from \(\mu \sim N(0, 0.002^2)\).
The loop updates 253 times; the batch formula updates once with the sum. Before running, decide whether the two final posteriors will agree exactly, approximately, or not at all.
After day 1 the posterior sd is 0.00199 — barely below the prior’s 0.002, because one observation with sd 0.02 carries a hundredth of the prior’s precision. By day 50 it is 0.00163; by day 253 the posterior is mean 0.00199, sd 0.00106. And the batch formula, which never saw the sequence, gives exactly 0.00199, 0.00106. Sequential equals batch, to every printed digit. The practical content of this equality is that the posterior is a sufficient summary of the past: two numbers, \((m, v)\), carry everything 253 days of returns had to say about \(\mu\), and tomorrow’s update needs only those two numbers and tomorrow’s return. You never need to store or revisit the history.
Notice also that the posterior mean, 0.00199, is well below the sample mean of 0.00277. The prior with precision \(1/0.002^2 = 250\,000\) and the data with precision \(253/0.02^2 = 632\,500\) share the weight 28 : 72, so the posterior mean is 72 % of the way from 0 to 0.00277. Unlike the S&P example in Section 5.1, here the prior is a real participant, because the daily sd of 0.02 makes each observation so uninformative about the mean that a year of them is worth only two and a half priors. This is the honest arithmetic of estimating expected returns: even a full year of daily data leaves the mean uncertain to ±0.1 % per day, which annualises to ±27 %.
Watch the belief tighten
Plot the posterior mean with a ±2 sd band through the year, and look at what happens in March.
During the March 2020 crash, did the posterior sd of \(\mu\) widen, stay constant, or keep shrinking?
The posterior mean dips negative in March — on 13 March it is −0.00014 — and then climbs steadily through the recovery to 0.002 by December. The band, meanwhile, does nothing but narrow: the posterior sd on 13 March is 0.00163, on its way monotonically from 0.002 to 0.00106. It kept shrinking through the crash. Look at the variance recursion and you will see why: \(1/v_t = 1/v_{t-1} + 1/\sigma^2\) never sees the data. It counts observations and nothing else. Every day reduces \(v\) by the same amount regardless of whether that day’s return was 0.1 % or −12.9 %, and regardless of whether the world has changed. A month in which AAPL’s daily sd was three times its usual value, and in which any reasonable observer would have said “I am less sure than I was what the mean return is”, produced a model that became more confident every single day. This is exactly the wrong response to a regime change, and it is not a bug in the code — it is the model. A static parameter cannot adapt.
Allowing a latent quantity to change
Question: What if the underlying quantity changes over time?
Result: Explain process noise, observation noise and a local-level update.
Builds on Updating a fixed unknown quantity.
A static parameter cannot adapt — so let it move
The fix is one line. Before each update, admit that \(\mu\) may have moved overnight:
\[ \mu_t = \mu_{t-1} + \eta_t, \qquad \eta_t \sim N(0, q) \qquad\Longrightarrow\qquad v_{t-1} \leftarrow v_{t-1} + q . \]
The parameter is now a random walk with step variance \(q\), and the belief about it must be inflated by \(q\) before the day’s observation tightens it again. With \(q > 0\) added every day the variance no longer shrinks to zero; it settles at a positive steady state where the daily inflation by \(q\) exactly balances the daily precision gain from one observation. At the steady state the update \(m_t = m_{t-1} + K\,(r_t - m_{t-1})\) has a constant gain \(K\), which means the estimate has a constant effective memory: it is an exponentially weighted moving average of past returns, with the decay rate set by \(q\). A large \(q\) says the parameter moves fast, so forget quickly; \(q = 0\) says it never moves, so remember everything — which is the static model you just saw fail.
Write the two equations together and you have a named model. The observation equation \(r_t = \mu_t + \varepsilon_t\) says each return is the current level plus noise; the state equation \(\mu_t = \mu_{t-1} + \eta_t\) says the level wanders. This is the local-level model, the simplest state-space model, and the recursion you have already written — inflate by \(q\), then apply the conjugate update — is the Kalman filter for it. The Kalman filter is nothing more than sequential Bayesian updating of a normal belief about a state that follows a linear dynamic; Chapter 6 adds trends, seasonal components and regressors to the state, but the two-step rhythm of predict-then-update is already all here.
Let the data choose q: the local-level model
The value of \(q\) is a modelling choice, and the honest way to make it is to estimate both variances — the observation noise \(\sigma^2_\varepsilon\) and the state noise \(\sigma^2_\eta\) — by maximum likelihood. statsmodels’ UnobservedComponents does this: it runs the Kalman filter for a candidate pair of variances, computes the likelihood of the data from the filter’s one-step-ahead prediction errors, and maximises. Apply it first to the returns themselves, over 2019–2020, to see whether the mean of AAPL’s return actually wanders.
The irregular variance is 0.000566 (a daily sd of 2.4 %); the level variance is 8.8e-13 — for all practical purposes zero. The data say the mean of AAPL’s daily return does not wander: a constant is as good as it gets, and any \(q > 0\) you imposed would be your opinion rather than the data’s. This is a finding, not a failure. Expected returns are notoriously hard to time, and the local-level MLE has just quantified that: over 505 days there is no detectable drift in the drift. Risk is a different story.
Adaptive risk: a local level for |r|
The mean does not move, but anyone who lived through March 2020 knows that the volatility did. Give it a wandering level. Under normality, \(\sqrt{\pi/2}\,|r_t|\) is an unbiased proxy for the daily standard deviation (because \(E|Z| = \sqrt{2/\pi}\) for a standard normal), so a local-level model on that series estimates a time-varying volatility, with the filtered state as the estimate.
Now the level variance is 8.9e-06, decidedly not zero, against an irregular variance of 0.0003566; their ratio, the signal-to-noise ratio \(q = \sigma^2_\eta/\sigma^2_\varepsilon\), is 0.0251. The level moves, but each day’s \(|r|\) is mostly noise — a single day’s absolute return is a very noisy reading of that day’s volatility, which is why the filtered line in the figure is smooth while the raw series is a thicket. The filtered volatility is 1.21 % per day on 2 January 2020, 8.84 % on 16 March at the height of the crash, 4.93 % at the end of March, and back to 1.74 % by 31 December. That is a seven-fold rise and a five-fold fall within a year — the adaptive risk estimate the static model could not produce, obtained from the same conjugate update plus one line.
The bridge to the Kalman filter
In steady state the filter’s gain solves a quadratic in \(q\): \(K = \tfrac{1}{2}\left(-q + \sqrt{q^2 + 4q}\right)\). With \(q = 0.0251\) that is \(K = 0.146\). Rather than take the formula on trust, run the hand-rolled recursion — inflate \(v\) by \(\sigma^2_\eta\), then the conjugate update with \(\sigma^2_\varepsilon\) — on the \(|r|\) series and compare it with statsmodels’ filtered_state.
The five-line recursion reproduces filtered_state to 5.2e-07 — the residual difference is the diffuse initialisation, which statsmodels handles exactly and the loop approximates with \(v_0 = 1\), and which has washed out completely after twenty days. The steady-state gain is 0.146: each day the volatility estimate moves 14.6 % of the way from where it was toward that day’s absolute return. A gain of \(K\) corresponds to an exponentially weighted moving average with span \(2/K - 1 = 12.7\) days, or a decay factor \(\lambda = 1 - K = 0.854\). RiskMetrics’ celebrated \(\lambda = 0.94\) (a span of 32 days) is the same filter with a smaller \(q\) — about 0.004 — chosen in 1994 by fitting across many assets rather than to this one series. The EWMA volatility that every risk desk computes is, it turns out, the Kalman filter for a local-level model, with the smoothing parameter standing in for a signal-to-noise ratio that could have been estimated.
Chronological adaptive-risk evaluation
Question: Was the risk estimate available before the decision?
Result: A training-boundary check and a corrected block-update experiment.
Builds on Allowing a latent quantity to change.
Estimate risk parameters before the period being forecast
The earlier 2019–2020 local-level fit is retrospective: even its filtered states use variances estimated from that whole interval. For a forward exercise, estimate parameters on 2019 and carry the filter into 2020 without refitting.
There are 252 training and 253 forward days. The levels are 1.15%, 5.91% and 1.78% at the displayed dates. A close-based state is available only after that close; a trading rule uses an appropriate execution lag. Filtered states and smoothed states have different information sets.
Sequential updates require genuinely new evidence
Using a previous posterior as the next prior is valid when the next likelihood contributes new data, with the conditional dependence model specified.
- Expanding-window refit: original prior + all eligible observations.
- Sequential update: previous posterior + newly available observations.
- Rolling window: remove expired information explicitly; carrying the entire old posterior does not implement a fixed window.
Do not multiply yesterday’s posterior by a likelihood containing yesterday’s observations again. A rolling calendar alone does not prevent double-counting or future-labelled training rows.
Rolling block updates: a debug exercise
The notebook’s practice for this section asks for block updating: use the first 20 days to build a prior, then update every 10 days, with the posterior of each block becoming the prior of the next. The starter below is meant to do that but contains a bug that makes the loop forget everything between blocks. It runs without error and prints plausible-looking numbers. Find the bug before running; then predict what the correct final sd should be, given that a prior built from 20 days carries exactly 20 observations of information.
The line m, v = m0, v0 sits inside the loop, so every block starts again from the 20-day prior and the posterior after the last block reflects only 20 + 3 = 23 days of information: sd 0.00417, mean 0.00307 (the 20-day prior mean updated by the final three-day block). Move the assignment above the loop and beliefs carry forward. The corrected final posterior has mean 0.00277 — the full-year sample mean, as it must, since a prior built from the first 20 days with variance \(\sigma^2/20\) is exactly what a flat prior would produce after 20 observations — and sd 0.00126, which is \(\sigma/\sqrt{253} = 0.02/\sqrt{253}\). Block updates end where a one-shot update would; the block size is irrelevant to the answer. That is the sufficiency property from the start of the section, applied to a prior that was itself a posterior.
In Colab — the notebook’s dynamic Student-t model
The notebook goes one step beyond the local level: a Student-t observation model whose location, log-scale and log-tail-index all follow random walks, so that mean, volatility and tail thickness vary through time. This is the model of Section 5.1 with the state equation of Section 5.4 attached to each parameter, fitted by NUTS.
with pm.Model() as dyn_t_model:
sigma_mu_rw = pm.Exponential("sigma_mu_rw", 1.0)
sigma_ls_rw = pm.Exponential("sigma_ls_rw", 1.0)
sigma_eta_rw = pm.Exponential("sigma_eta_rw", 1.0)
mu = pm.GaussianRandomWalk("mu", sigma=sigma_mu_rw, shape=T)
log_sigma = pm.GaussianRandomWalk("log_sigma", sigma=sigma_ls_rw, shape=T)
eta = pm.GaussianRandomWalk("eta", sigma=sigma_eta_rw, shape=T)
y_obs = pm.StudentT("y_obs", nu=2 + pm.math.exp(eta), mu=mu,
sigma=pm.math.exp(log_sigma), observed=returns)
trace = pm.sample(tune=1000, draws=1000, target_accept=0.9)With \(T\) days of data this model has \(3T + 3\) parameters — the curse of dimensionality that killed the grid in Section 5.1 — and it is precisely the case where NUTS’ gradient-guided proposals are indispensable: a random-walk Metropolis sampler in 1 500 dimensions would not converge in any practical time. The GaussianRandomWalk on mu is the local-level state equation; the exponential hyperpriors on the walk scales are the Bayesian version of estimating \(q\) by maximum likelihood.
\(\sqrt{1/(1/0.002^2 + 253/0.02^2)} = 0.00106\). Order is irrelevant: precision adds \(1/\sigma^2\) per observation, so sequential and batch updates coincide exactly — the posterior is a sufficient summary of the past.
Add process noise before each update: \(v_{t-1} \leftarrow v_{t-1} + q\), i.e. let \(\mu_t = \mu_{t-1} + \eta_t\). The variance then settles at a positive steady state with constant gain \(K\). Observation \(r_t = \mu_t + \varepsilon_t\) plus this state equation is the local-level model, and the recursion is its Kalman filter.
\(K = \tfrac12(-q + \sqrt{q^2 + 4q}) = 0.146\), an EWMA with span \(2/K - 1 = 12.7\) days (\(\lambda = 0.854\)). RiskMetrics’ \(\lambda = 0.94\) (span 32) is the same filter with \(q \approx 0.004\). The local-level MLE on the returns themselves gave \(q \approx 0\): the mean does not wander, the volatility does.
Bayesian Signals and Partial Pooling
| Part | Question |
|---|---|
| 5.5.1 Targets and shrinkage across stocks | How much should a noisy stock estimate move towards its group? |
| 5.5.2 Feature effects and posterior covariance | How do shared features and stock effects combine? |
| 5.5.3 Predictive uncertainty and prior sensitivity | Is a positive expected return the same as a positive realised return? |
Use the aligned features built in Chapter 2 to predict next-month returns. Learn what a hierarchical prior shares across stocks, how it changes noisy effects, and which probability is relevant to a decision.
Targets and shrinkage across stocks
Question: How much should a noisy stock estimate move towards its group?
Result: Compare no pooling, complete pooling and partial pooling.
Begin with the definitions and assumptions; then trace the worked example.
An investment question needs a target and a clock
The question is: which candidates have favourable expected next-month returns given their available characteristics?
Use three features: momentum, realised volatility and sector-relative momentum. Preserve stock/date columns. Keep future returns outside the feature matrix and retain pending candidates.
The nine-stock classroom panel has survivor membership, static sectors and assumed monthly availability. Its previous predictions are ignored. Its realised results have already been explored in earlier chapters.
A target is an outcome over a declared horizon. Here it is next-month simple return, rather than price level, a cumulative multi-month return or a causal treatment effect. A stock can be eligible for prediction while its outcome remains unavailable. That distinction carries through fitting, scoring and accounting.
The public panel has 639 candidate rows, 630 valid mature labels and nine pending labels. Some raw pending return values happen to be present in the file; the declared validity flags still exclude them from fitting and evaluation. Treating those raw numbers as realised outcomes would violate the exercise contract. The earlier stored prediction column is also excluded: this chapter constructs new predictions from its own training samples.
Partial pooling sits between two strong assumptions
No pooling estimates each stock separately. Complete pooling assumes a common effect. Partial pooling lets stocks differ while sharing a distribution:
\[a_i\sim N(0,\tau^2),\qquad \bar e_i\mid a_i\sim N(a_i,\sigma^2/n_i).\]
\[E[a_i\mid\bar e_i]=w_i\bar e_i,\qquad w_i=\frac{\tau^2}{\tau^2+\sigma^2/n_i}.\]
Here \(\bar e_i\) is a residual mean after a declared baseline. A large noisy residual mean need not imply a large persistent stock effect.
A stock-specific effect represents the part of the conditional mean not described by the supplied characteristics and sector effects. It is not automatically a risk-adjusted alpha: that interpretation requires an appropriate factor benchmark and return definition. The exercise calls it a stock deviation to avoid borrowing an economic interpretation that the model has not earned.
The prior standard deviation \(\tau\) describes heterogeneity across stocks. A small value expresses a belief that persistent differences are small; a large value permits more dispersion. The sampling variance \(\sigma^2/n_i\) describes how noisy each individual estimate is. Their ratio determines the weight on that stock’s own evidence. This is the same precision-weighting used in the conjugate normal update, now applied across assets.
Predict which apparent winner shrinks most
Self-check: Stock A has only three monthly residuals averaging +8%. Stock C has sixty averaging −1%. Under the same prior and residual variance, which receives less weight on its own sample mean?
- A, because three observations carry less precision
- C, because negative estimates must be removed
- Both receive the same weight
- Neither: Bayesian estimates must equal sample means
Feedback: The weight depends on sampling precision and prior variance. Fewer observations imply greater pooling; the sign alone does not determine the weight.
Compute the shrinkage instead of declaring a winner
Stock A shrinks from 8% to 0.86%, with own-data weight 0.1071. Stock B becomes 1.14%; C becomes −0.71%. The largest raw sample mean is no longer the largest posterior mean.
The fixture is simulated to vary sample sizes transparently. It demonstrates a mechanism, rather than proving that the pooled ranking will earn money.
Change the prior and explain what moves
A simulated stock has three monthly residuals averaging 8%. Residual sd is fixed at 10%; the stock effect has prior mean zero.
Try: compare τ = 0.5%, 2% and 10%. Explain both the posterior mean and interval. Prior sensitivity is part of the result.
Feature effects and posterior covariance
Question: How do shared features and stock effects combine?
Result: Fit a scaled feature model while retaining coefficient covariance.
Builds on Targets and shrinkage across stocks.
A feature model and stock effects can be estimated together
\[r_{i,t+1}=b_0+z_{it}^{\top}\beta+s_i^{\top}\gamma+a_i+\varepsilon_{it},\qquad \varepsilon_{it}\sim N(0,\sigma^2).\]
\(z\) contains three train-standardised features; \(s\) contains two sector indicators with Energy as reference. Nine stock effects share \(a_i\sim N(0,\tau^2)\).
Fixed exercise scales: intercept/features 5%, sector effects 3%, stock effects τ = 2%, residual sd σ = 10%. These are assumptions for sensitivity analysis, rather than universally appropriate market values.
This is a hierarchical normal model conditional on fixed hyperparameters. There are two levels: the observed returns depend on stock effects, and those effects share a prior distribution. A fully Bayesian extension would assign hyperpriors to \(\tau\) and \(\sigma\), estimate them and integrate over their uncertainty. The browser example fixes them so that the posterior is available in closed form and every matrix operation can be inspected.
The design contains an intercept, sector indicators and all nine stock indicators. The corresponding unpenalised design is redundant, but the proper Gaussian priors give a positive-definite posterior precision matrix. Consequently a conventional OLS decomposition cannot be copied over mechanically. Individual stock deviations are interpreted conditional on this prior and baseline. Predictions combine all terms. New stocks would need an explicit unseen-identity policy, such as integrating a new effect drawn from the shared prior.
Standardisation belongs inside the fitting boundary
The design has 15 columns. The function receives the training mean and standard deviation; it never estimates a new scale from prediction rows.
Ticker indicators represent this declared teaching universe. They do not establish historical membership. A production universe and sector history need separate point-in-time evidence.
The normal posterior makes the calculation transparent
For \(\theta\sim N(0,V_0)\) and \(y\mid\theta\sim N(X\theta,\sigma^2I)\):
\[V_n=(V_0^{-1}+X^{\top}X/\sigma^2)^{-1},\qquad m_n=V_nX^{\top}y/\sigma^2.\]
The same posterior mean minimises squared error plus a quadratic penalty. A prior changes parameter estimates and predictions, while the posterior covariance adds uncertainty around them.
The next function also calculates the normalised marginal likelihood needed in §5.6.
Prior precision and data precision add. The stock-effect block of \(V_0\) determines pooling strength; the feature and sector blocks express separate prior scales. Standardising features on training rows makes a one-unit coefficient correspond to one training standard deviation of that characteristic. It does not make economically different prior choices equivalent.
The likelihood assumes constant residual variance and independent rows. Stocks observed in the same month often share shocks; monthly returns can also exhibit temporal dependence. The posterior intervals in this exercise therefore describe a simplified model, not dependence-robust guarantees. A richer model can include common factors or structured residual covariance. Checking this assumption is as necessary as choosing a prior.
Fit parameters and preserve their covariance
The function uses linear solves rather than explicitly inverting the precision matrix. All prior variances are strictly positive. Its evidence calculation includes normalising constants; an unnormalised posterior density cannot replace it.
Read coefficients in their transformed units
There are 207 training rows over 23 outcome periods before the 1 January 2021 freeze. Posterior feature means are −0.0129, +0.0331 and +0.0050 per training standard deviation.
A coefficient’s sign is conditional on the included predictors, sector effects, stock effects and prior. It does not establish a causal trading mechanism.
The difference between 207 rows and 23 monthly snapshots matters. Nine stocks in one month do not necessarily supply nine independent market experiments. The code’s numerical solution is exact for its Gaussian assumptions; the economic evidence still depends on those assumptions and the source data.
An estimate of −0.0129 on standardised momentum means that, within this model, increasing momentum by one training standard deviation changes expected next-month simple return by about −1.29 percentage points, holding the other columns fixed. It is neither a raw-unit slope nor a statement that increasing historical momentum causes returns to fall.
Predictive uncertainty and prior sensitivity
Question: Is a positive expected return the same as a positive realised return?
Result: Separate parameter from predictive probability and test prior sensitivity.
Builds on Feature effects and posterior covariance.
Parameter uncertainty and future-return uncertainty differ
For a candidate row \(x_*\):
\[E[r_*\mid D]=x_*^\top m_n,\quad \operatorname{Var}(E[r_*\mid\theta]\mid D)=x_*^\top V_nx_*.\] \[\operatorname{Var}(r_*\mid D)=\sigma^2+x_*^\top V_nx_*.\]
A credible interval for the conditional mean leaves out the new observation’s residual shock. A predictive interval includes it.
For the January 2021 candidates, the median posterior standard deviation of the conditional mean is 1.83%, while the predictive standard deviation is 10.17%. Most of the next month’s uncertainty comes from the residual return, not from uncertainty about the fitted mean. A narrow interval for expected return is compatible with a wide range of realised returns.
This is why a probability that the conditional mean is positive cannot be described as the probability that next month’s trade will win. Both can be computed, but they answer different questions and use different variances. A portfolio-level prediction also needs cross-stock covariance; adding individual predictive variances while assuming zero covariance can severely understate aggregate risk.
A probability of positive mean is not a win probability
For AAPL, P(conditional mean > 0.2%) is 0.974, whereas P(next return > 0.2%) is 0.638.
The 0.2% screen is an assumed round-trip hurdle for this experiment, rather than a prediction of actual turnover costs. Cash earns zero. Both probabilities are conditional on the declared likelihood and prior.
Self-check: distinguish the two probabilities
Self-check: A posterior gives P(conditional mean > cost) = 0.97. What follows?
- The next trade has a 97% probability of making money
- The conditional mean exceeds the hurdle with 97% posterior probability under this model
- The stock has a causal alpha of 97%
- The strategy will be profitable after every rebalance
Feedback: Realised return includes a new residual shock. Its posterior predictive probability generally differs from the probability about the conditional mean.
Check whether a stronger prior changes the decision
Across τ = 0.5%, 2% and 8%, the maximum absolute stock deviation changes from 0.11% to 1.26% to 4.35%. AAPL remains the highest forecast in this snapshot.
An unchanged top name does not imply that predictions, confidence or portfolio weights are insensitive. Report what changed and what stayed stable.
Repair exercise: uncertainty needs the observation noise
Predictive variance must include both parameter uncertainty and residual variation. Fix the assignment; keep the conditional-mean variance separately. Repair task: Repair the omitted uncertainty term.
Worked repair: reported_predictive_var = parameter_var + investment_model["residual_sd"]**2. The starter reports the wrong information set; the repair preserves the separate quantities.
What you discovered
- Partial pooling changes a noisy signal rather than merely attaching an error bar.
- A hierarchical prior shares information across stocks; fixed hyperparameters keep the exercise transparent.
- Posterior covariance describes parameter uncertainty. Predictive variance also includes future residual noise.
- Prior sensitivity, dependence and timing remain part of the research result.
Next: §5.6 combines competing models and takes forecasts with stock and date columns through a declared portfolio policy.
Model Uncertainty and Portfolio Evidence
| Part | Question |
|---|---|
| 5.6.1 Model evidence and predictive mixtures | How should competing models contribute to a forecast? |
| 5.6.2 Decision objectives and walk-forward scoring | How does uncertainty enter a dated decision rule? |
| 5.6.3 Holdings, accounts and incremental evidence | Does the resulting policy improve the complete account? |
Bayesian modelling can produce trading signals and portfolio decisions. This section compares prior models, builds predictive mixtures and evaluates complete net accounts. Walk-forward supplies the time boundary; the Bayesian model supplies the forecast distribution.
Model evidence and predictive mixtures
Question: How should competing models contribute to a forecast?
Result: Evidence weights and a mixture that includes model disagreement.
Begin with the definitions and assumptions; then trace the worked example.
One research chain connects the chapters
The process must be auditable at each fitting cutoff. Posterior calculations do not repair future information in an input, an overlapping label or an incorrectly reconstructed trading account.
Model probabilities use integrated evidence
For a declared candidate set \(M_1,\ldots,M_K\):
\[p(M_k\mid D)=\frac{p(D\mid M_k)p(M_k)}{\sum_jp(D\mid M_j)p(M_j)},\qquad p(D\mid M_k)=\int p(D\mid\theta_k,M_k)p(\theta_k\mid M_k)d\theta_k.\]
The evidence averages fit over the prior; it is not the likelihood at the best-fitting parameter. Proper priors and their normalisation matter.
Our three models differ only in stock-effect prior sd τ = 0.5%, 2%, 8%, with equal prior model probabilities.
A maximised likelihood rewards fit without integrating parameter uncertainty. Marginal likelihood instead accounts for the probability of the observations across the model’s prior parameter distribution. Priors that spread mass over many implausible parameter values can have low evidence even if a small region fits well. Consequently an arbitrary change in prior width can change model weights. Bayes factors are ratios of evidences; they inherit this sensitivity.
All three models here condition on the same training observations, feature construction and fixed residual standard deviation. Their evidence weights are probabilities within this finite declared model set. They are not probabilities that any model is universally true, or probabilities of a profitable backtest. The candidate set omits heavy tails and residual correlation, so a dominant weight does not establish that those omissions are harmless.
Compute an evidence-weighted predictive mixture
This section is standalone: each cell includes the setup functions it needs. Models are recomputed from the current training sample; future portfolio performance never determines their weights.
Refit the same normal posterior inside each training window
The next cell defines prediction without relying on a previous section.
Mixture uncertainty includes model disagreement
\[p(r_*\mid D)=\sum_kw_k\,p(r_*\mid D,M_k),\qquad \bar m=\sum_kw_km_k.\] \[\operatorname{Var}(r_*\mid D)=\sum_kw_k(v_k+m_k^2)-\bar m^2.\]
A weighted average of standard deviations is not this variance. Averaging only \(v_k\) also misses disagreement between the component means.
Mixture tail probabilities average component tail probabilities; the mixture is generally not Normal.
There are now three uncertainty sources: uncertain fitted parameters within each model, new return noise and disagreement across models. Each component predictive variance already includes the first two; the law of total variance adds the third through the squared component means. This decomposition prevents the false precision that arises when the best-looking model is treated as certain.
The mean and variance summarise a mixture without determining its complete distribution. In particular, a Normal density with these matching moments is an approximation. For log predictive scoring we evaluate the actual weighted component densities. Predictive stacking, which estimates combination weights for predictive performance, is a different procedure from the evidence weights used here; random-row leave-one-out is also not automatically appropriate for a chronological investment question.
Use logs to normalise model weights safely
logsumexp prevents underflow when multiplying many small densities. The returned probabilities distinguish a conditional-mean hurdle from an actual return hurdle.
Read model weights within their candidate set
The model weights are 0.7728, 0.2264 and 0.0008. The training sample favours stronger pooling within these assumptions.
This does not mean a 77.28% probability that the tightest-prior strategy will make money.
Self-check: the highest evidence is not a profit guarantee
Self-check: A model has posterior model probability 0.90 within a three-model comparison. What does 0.90 describe?
- Its probability of positive net profit next year
- Its relative posterior probability within the declared model set and priors
- Its probability of being the true causal mechanism
- Its expected annual return
Feedback: Model probability depends on the candidate set, prior probabilities and integrated likelihood. Future trading profit is a different random quantity.
Decision objectives and walk-forward scoring
Question: How does uncertainty enter a dated decision rule?
Result: Chronological refits and predictive scores on mature outcomes.
Builds on Model evidence and predictive mixtures.
A posterior supports a decision only through an objective
A general Bayesian decision chooses an action \(w\) to maximise:
\[E[U(w^{\top}r_{t+1}-C(w,w_{t^-}))\mid D_t].\]
Utility \(U\), predictive joint returns, transaction costs \(C\), constraints and drifted holdings all matter.
This exercise uses a simpler declared rule: rank posterior mean returns; reserve three 0.33 slots; retain a slot only if P(conditional mean > 0.2%) ≥ 0.95. Otherwise keep cash. This is an illustrative policy, not a solved optimal portfolio.
The probability threshold is an assumption to investigate, not a result delivered by Bayes’ theorem. A mean-return screen ignores risk aversion unless risk enters elsewhere in the policy. A portfolio optimiser would require the joint predictive covariance and explicit exposure, turnover and concentration constraints. The small normal model here is designed to teach the chain from inference to action without claiming to solve that larger allocation problem.
The 95% threshold and 0.2% hurdle are fixed for the reported classroom comparison; neither is presented as a tuned optimum. The supplied data have already been explored in the course, so this exercise is an historical illustration rather than a newly untouched performance test. Changing the threshold after viewing its account is a sensitivity experiment and must be labelled as such. Repeated strategy searches need an additional evaluation design and a record of the attempted rules.
Each annual refit sees only labels already available
Training row counts are 315, 423 and 531. The three annual blocks produce 315 predictions: 306 evaluable, nine pending.
Features and priors are transformed inside the training window. January’s freeze uses labels strictly before 1 January. Refit from the original prior on the complete expanding history, avoiding double-counting.
Refitting and sequential updating count evidence differently
Self-check: You retain last year’s posterior and also refit on every historical observation, including last year’s data. What is the problem?
- No problem: more repetitions always improve precision
- Old observations are counted again through both the prior and likelihood
- The posterior must become uniform
- Walk-forward prohibits any refitting
Feedback: Use the original prior with all eligible history, or the previous posterior with genuinely new likelihood contributions. Combining the previous posterior with the same old data counts evidence twice.
Score predictions before reading the portfolio curve
The mixture has RMSE 0.0952, versus 0.0930 for a training-mean baseline. Its exact-mixture mean log predictive density is 0.9315 in decimal-return units; higher is better for comparable models on the same outcomes.
This forecast does not beat the simple baseline on RMSE. A favourable portfolio ranking would not reverse that finding.
RMSE evaluates the posterior predictive mean under squared-error loss. Log predictive density evaluates the complete distribution, including dispersion and tails. Density scores depend on measurement units, so the value 0.9315 is not an accuracy percentage and cannot be compared with a score computed after rescaling returns into percentages. A baseline log score under the same units and observations is needed to interpret predictive-density improvements.
An individual stock prediction is not a portfolio allocation. A model can lose on average squared error while a ranking-based rule happens to select a profitable subset, and the converse is also possible. Report the forecast metrics and the economic policy metrics separately. Excluding pending labels is necessary for either evaluation; excluding candidates before prediction because their future outcome is missing would be a different error.
Holdings, accounts and incremental evidence
Question: Does the resulting policy improve the complete account?
Result: Reconciled policy comparisons with sensitivity and leakage checks.
Builds on Decision objectives and walk-forward scoring.
Build a holdings rule without changing the candidate spine
A permitted name retains its original 0.33 allocation. Blocked slots stay in cash. The equal-weight comparator invests 0.99 across all nine names.
Exposure, costs and drawdown must accompany returns. A lower-exposure portfolio is not directly evidence of better stock selection.
Carry drifted holdings and cash through the ledger
The loop below represents changing portfolio state. Cash earns zero, trades cost 10 bps per dollar bought or sold, opening NAV is 100, and there is no terminal liquidation. The 1% target cash reserve funds the specified fees.
At every month: closing NAV − opening NAV = security P&L − fees. A positive prediction never substitutes for a realised return.
Turnover is calculated from current holdings after the previous period’s gains and losses. Subtracting the previous target weights would omit that drift. Targets in this exercise are fractions of opening NAV before the current trade fees; cash then absorbs fees and the unused budget. The account assumes the supplied simple returns can be earned on those target positions. Execution timing, spread, capacity and source accuracy need evidence beyond the accounting identity.
Only mature, valid monthly outcomes enter the realised ledger. Nine pending predictions remain in the forecast output, while their unknown holding period stays outside realised accounting. The final NAV is marked at the end of the last observed holding period without an assumed sale of all holdings. That endpoint convention matters if a different study includes terminal transaction costs.
Compare all three complete accounts
With opening NAV 100, the net ending values are 183.90 equal-weight, 204.77 ranking alone and 171.27 Bayesian screening. Maximum drawdowns are 18.15%, 34.26% and 34.26%, respectively. Exposure in this table averages the same 34 realised months as the ledger; the earlier signal summary covers all 35 prediction months.
The Bayesian threshold is evaluated as part of the whole policy. It does not inherit a performance claim merely from using posterior probabilities.
Read the incremental result, including a null result
Compare the threshold policy with ranking alone using paired monthly net returns, then apply Chapter 3’s dependence-aware uncertainty methods. Thirty-four periods and a survivor sample provide limited evidence.
Here the screen subtracts about 0.55 percentage points per month from ranking alone on paired net returns. It also changes exposure. Forecast RMSE 0.0952 is worse than the training-mean baseline 0.0930. These results support further investigation or rejection of this rule, rather than a claim that the Bayesian screen adds value. The examined historical sample is not untouched evidence.
Changing future outcomes must not change the frozen forecast
The altered future does change later training samples when those outcomes become eligible. It must not change predictions from the earlier frozen fit.
This implementation check does not verify historical source availability or validate the chosen likelihood.
Repair exercise: keep pending labels outside the score
The starter includes every prediction in an RMSE denominator. Fix the evaluation rows without removing pending prediction candidates. Repair task: Keep all forecasts; restrict evaluation to valid outcomes.
Worked repair: evaluation_rows = investment_walk.loc[investment_walk.label_valid].copy(). The starter reports the wrong information set; the repair preserves the separate quantities.
A research memo must separate signal, model and policy
Report the target and clock; likelihood and priors; pooling assumptions; model set and evidence weights; parameter versus predictive uncertainty; forecast baseline; frozen holding rule; exposure, costs and net accounts.
State which result would change your decision. Include prior sensitivity and unsuccessful comparisons. A posterior probability describes the specified model, rather than certifying economic profitability.
What you discovered
- Bayesian inference changes signals through pooling and combines uncertainty through predictive mixtures.
- Evidence weights are conditional on the declared model set and proper priors.
- Walk-forward trains and evaluates that model process at explicit time boundaries.
- Prediction quality, portfolio profit and risk are separate results to report.
Next: Chapter 6 develops time-series dependence and dynamic risk models under the same availability discipline.
Chapter Wrap-up
You can now explain and run a Bayesian investment research process, from the original single-asset statistical foundations to hierarchical signals, model uncertainty and net accounts. Examples use public returns and clearly labelled simulated fixtures.
You can turn a prior and a likelihood into a posterior — on a grid, where the posterior mean of the S&P 500’s daily return came out at 0.000497 and matched the closed form to six decimals, and with a twelve-line Metropolis sampler that gave acceptance 0.47, a tail index \(\nu \approx 2.7\) for GOOG, and a reward-to-t-scale posterior with mean 2.38 and credible interval [1.30, 3.59] against a classical 0.77 whose confidence interval could not exclude zero. You know why the gap exists — the t-scale is six-tenths of the sample sd — and that on a standard-deviation basis the posterior mean is 1.18, so you can say which risk a Sharpe ratio divides by before you quote it. You know what NUTS does that Metropolis does not: gradient-guided proposals, automatic tuning, and diagnostics.
You can check a model before and after fitting it. The prior check on regression means exposed Prior I spraying lines over ±200 where the data lived in ±35; the conjugate normal-linear posterior turned an OLS slope of 3.011 into a posterior mean of 3.000 with sd 0.069 through precision-weighting; the posterior predictive check produced Bayesian p-values of 0.45 for the mean and 0.135 for the maximum, the first hint that normal tails are thin.
You can fit a regression that reduces the influence of large residuals under a heavy-tailed likelihood. One word — Normal → StudentT — gave a slope of 5.17 where OLS gave 4.40 on data with three planted outliers, through implicit weights \((\nu+1)/(\nu + r_i^2)\) of 0.02 on the outliers; on GOOG’s beta the same change gave 1.094 against 1.064 with \(\nu = 3.35\), and you saw that twelve crash days carrying 36.5 % of \(\sum x^2\) were the reason.
You can carry beliefs forward. Sequential updating reproduced the batch posterior 0.00199 ± 0.00106 exactly; a static parameter’s variance kept shrinking through the March 2020 crash; one line of process noise turned the recursion into the Kalman filter of a local-level model, whose maximum-likelihood fit said AAPL’s mean does not wander but its volatility moved from 1.2 % to 8.8 % in ten weeks, and whose steady-state gain of 0.146 is an EWMA with a 12.7-day span.
Mistakes library: LTCM — model confidence and liquidity risk
The Federal Reserve’s May 1999 testimony identified excessive leverage and weak counterparty risk management as central lessons of LTCM’s 1998 crisis. Counterparties had not adequately understood its large positions in illiquid markets. See the Federal Reserve testimony.
Course application: a narrow posterior conditional on a return model cannot describe every funding or liquidation scenario. Stress joint losses, changing correlations, liquidity and leverage alongside posterior predictive checks. This historical case does not establish that LTCM used the Bayesian models taught here.
Partial pooling is now worked in §5.5 and evaluated through predictive mixtures and net accounts in §5.6.
The next chapter keeps the state-space idea and adds time in earnest: autocorrelation, stationarity, ARIMA, the full Kalman filter with trends and regressors, and the volatility models — ARCH and GARCH — that make the adaptive-risk filter of Section 5.4 rigorous. The slide deck for this chapter runs every cell above with the same seeds; the two are meant to be read side by side.
← Chapter 4 · Contents · Chapter 6: Time Series Models for Trading and Risk →
Download the investment research notebook for §§5.5–5.6.
Sources and Further Study
- SciPy Student-t reference: distribution parameterisation and scale.
- Stan posterior analysis: multiple chains, R-hat, effective sample size and Monte Carlo error.
- Feng and He, Factor Investing: A Bayesian Hierarchical Approach: investment research with information sharing and estimation risk. Our fixed-hyperparameter classroom model is a simpler example, not a replication.