Chapter 5: Rethinking Statistics with Bayesian Methods

Chapter Introduction

Everything you did in Chapters 2 and 3 treated a parameter as a fixed, unknown number and the data as the random object. A regression coefficient is something; the sample gives you an estimate of it; a confidence interval tells you how the estimate would bounce around if you could draw the sample again. That frame is so familiar that it is easy to forget it is a choice. This chapter makes the other choice. A Bayesian treats the parameter as the uncertain object and the data as the thing that is fixed — you have the 1 258 daily returns you have, and the question is what those returns tell you about the parameter. The answer is not a number but a distribution: the posterior. Once you hold the posterior in your hands, every question you could ask — a point estimate, an interval, a prediction, a probability that a strategy is worth funding — becomes a question about that distribution, and each of them has a direct answer.

The difference is not philosophical decoration. Financial data are short, noisy and fat-tailed, and those three facts punish classical methods in specific ways. A five-year daily history of GOOG gives a classical annualised Sharpe ratio of 0.77 with a confidence interval that includes zero — a t-test cannot reject the hypothesis that the stock earned nothing. The same data, run through a Student-t model with the tails estimated rather than assumed away, yields a posterior whose 95 % credible interval is [1.30, 3.59] and which puts 99 % probability on the Sharpe ratio exceeding 1. Part of that gap is real information the classical procedure throws away; part of it is a change in what “risk” means in the denominator — and this chapter insists that you know which is which before you quote either number. Bayesian methods do not make data say more than they can; they make explicit what you assumed, and let you carry beliefs forward as the world changes.

For most of the twentieth century the posterior was a beautiful object nobody could compute. Bayes’ theorem needs a normalising integral over every parameter, and outside a handful of textbook cases that integral has no closed form. Markov chain Monte Carlo — invented for the hydrogen bomb at Los Alamos in 1953, rediscovered by statisticians in 1990 — broke the barrier by sidestepping the integral altogether: instead of computing the posterior, you wander through parameter space in such a way that the time spent in each region is proportional to its posterior probability. Today every systematic fund uses this machinery somewhere. Bridgewater fits Bayesian regime-switching factor models; AQR’s shrinkage estimators of expected returns are hierarchical Bayes in everything but name; Black–Litterman, the most widely used allocation model in institutional asset management, is a conjugate normal update. Libraries such as PyMC and Stan make the sampler a one-line call. This chapter writes the sampler by hand, in twelve lines of numpy, so that you know what the one line does.

The chapter is built around four moves, one per section, and every one of them runs in your browser on real data. Section 5.1 builds the posterior of GOOG’s Sharpe ratio — first on a grid, where Bayes’ rule is literally two arrays multiplied, then with a Metropolis sampler you write yourself — and confronts the classical estimate with the Bayesian one. Section 5.2 recasts regression as a statement that \(y\) is normally distributed around a line, introduces the prior predictive check that lets you reject an absurd prior before you have touched the data, and the posterior predictive check that lets the fitted model confess what it cannot reproduce. Section 5.3 replaces the normal likelihood with a Student-t and watches three planted outliers lose their grip on the slope, then estimates GOOG’s market beta through the March 2020 crash with both likelihoods. Section 5.4 stops refitting from scratch: yesterday’s posterior becomes today’s prior, a static parameter is shown to be unable to adapt, one line of process noise fixes it, and the recursion you have written turns out to be the Kalman filter.

Two disciplines run through the whole chapter and are worth naming now. The first is predict before you run: every numbered result is preceded by a question, and you should commit to an answer before the cell prints. The second is interpret every printed number: a posterior mean of 2.38 is not a result until you can say what it divides by, what the credible interval around it means, and why it differs from 0.77. The PyMC versions of every model appear alongside the numpy implementations in “In Colab” boxes; the slide deck for this chapter runs the same cells with the same seeds, so the numbers you see here and in the lecture will agree to the last printed digit.

A Bayesian Model of Risk and Reward: Priors, Likelihood, Posterior of the Sharpe Ratio
Bayesian Regression and Prior/Posterior Predictive Checks
Robust Regression with Fat Tails
Dynamic Bayesian Models: Rolling Updates and Adaptive Risk

A Bayesian Model of Risk and Reward: Priors, Likelihood, Posterior of the Sharpe Ratio

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 the Sharpe ratio 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 point estimates, in decreasing order of sensitivity to skew. 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. 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 (3 dp), 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 a mild cheat — it looks at the data to choose its centre — but a prior with standard deviation 0.0164 around a mean of 0.0008 is so diffuse that its centre is irrelevant; you will measure exactly how irrelevant on the grid below. The uniform prior on \(\sigma\) is the honest expression of “I know it is positive and I know nothing else”. The prior on \(\nu\) is the one that repays thought. An exponential density is highest at zero and decays; with mean 29 it is nearly flat over the range 2 to 60 where the answer will lie, so the data decide. What the prior forces is only that \(\nu > 2\). It does not force \(\nu \approx 31\) — that is merely the prior mean — and it does not force normal tails. The cell above checked this: the prior gives probability 0.241 to \(\nu < 10\) and 0.098 to \(\nu < 5\). It allows fat tails without insisting on them, which is the right posture when you have 1 258 observations to settle the question.

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 the Sharpe ratio, 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.

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 the posterior to any precision you like. 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 Sharpe ratio is 2.035 and P(Sharpe > 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.

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, inside the healthy band of roughly 0.2 to 0.5 for a random-walk sampler. The band matters: an acceptance rate near 1 means the steps are so small that every proposal is nearly the current point and the chain explores nothing, while an acceptance rate near 0 means the steps are so large that almost everything is rejected and the chain sits still. Tuning the step vector so that acceptance lands in the band is the one piece of hand-craft the algorithm demands.

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 Metropolis chain above needs 2 500 draws to deliver comparable precision on the Sharpe ratio, and would need far more in a larger model.

None of this changes what is being computed. NUTS and Metropolis target the same posterior and converge to the same answer; the notebook’s Sharpe posterior mean of 2.35 and the hand-rolled 2.38 differ by Monte Carlo noise. 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, which is Metropolis with gradient-guided proposals and automatic tuning”.

The posterior of the Sharpe 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 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, the classical Sharpe ratio of 0.77, lies outside the credible interval: under this model, given these data, the probability that the Sharpe ratio is as low as the classical estimate is under 2.5 %. This is the notebook’s headline — “the Sharpe ratio is under-estimated by classical statistics” — and it is true in a precise sense that the next paragraphs pin down.

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 “under-estimation” is a statement about ordinary-day risk. 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.

Credible interval versus confidence interval

The classical procedure has a standard error for the annualised Sharpe ratio too. Lo (2002) shows that under i.i.d. returns 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: a t-test on the mean return cannot reject the hypothesis that GOOG earned nothing over five years. The Bayesian credible interval is [1.304, 3.594], and the posterior probability that the Sharpe ratio 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”, which is the sentence every analyst wants to say and which classical statistics never permits.

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 better estimate of \(\mu\) accounts for the remaining difference — but it is a different order of claim from 2.38. Both numbers are correct answers to different questions. The t-scale Sharpe ratio measures reward against the risk of an ordinary day; the sd-basis Sharpe ratio measures it against all days including the crashes. A risk manager who cares about drawdowns wants the second; a trader sizing a position for tomorrow may want the first; 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.

AAPL: classical 0.922, 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 locates \(\mu\) better, 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 Prior/Posterior Predictive Checks

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 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. This is the prior predictive check, and it is the cheapest piece of model criticism in statistics: it needs no likelihood, no sampler and no data beyond the range of \(x\).

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 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 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

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 (3 dp) — 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.

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 three points responsible are not from the same process as the rest.

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. What you can predict is that the OLS standard error will be wrong, because it assumes tails that the residuals do not have.

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 that the OLS standard error on beta (0.025, giving a t-statistic of 42) is computed under an assumption the data reject.

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).

Dynamic Bayesian Models: Rolling Updates and Adaptive Risk

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 a model that must be refitted from scratch each morning is both wasteful and, more subtly, unable to say what it learned overnight. 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.

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.

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.

Chapter Wrap-up

You can now do four things you could not do at the start of the chapter, and each of them was done on real data in the browser.

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 Sharpe-ratio 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 predictive check exposed Prior I spraying regression 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 make a regression ignore what does not fit. 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: Clinton 99 % — a posterior that forgot the states move together (November 2016)

On the eve of the 8 November 2016 US presidential election the Princeton Election Consortium put Clinton’s win probability at 99 % and the HuffPost forecast at 98 %; FiveThirtyEight said 71 %. The polls were not badly wrong: the national average had Clinton ahead by about 3 points and she won the popular vote by 2.1. What decided the election was that Wisconsin, Michigan and Pennsylvania, polling at Clinton +2 to +6, went to Trump by less than a point, giving him the Electoral College 306–232.

The 99 % forecasts treated the fifty state polling errors as independent draws. Under that likelihood a modest lead becomes near-certainty, because the chance that enough independent errors line up the same way is tiny — the credible interval shrinks by \(\sqrt{n}\) as states are added. FiveThirtyEight gave the state errors a shared component, so one under-sampling of non-college white voters moved every Rust-Belt state together, and its posterior stayed wide enough to contain the outcome.

The width of a posterior is only as honest as the dependence structure in the likelihood. Section 5.3 made the point for tails; this is the point for correlation. A posterior predictive check — could this model produce a three-state miss the same way? — would have flagged the 99 % before election night did.

Where this leads: hierarchical models and partial pooling

Every model in this chapter had one asset. The most important Bayesian idea for the cross-section is partial pooling. Suppose you want the expected return of each of 500 stocks. Estimating each separately (no pooling) is noisy — Section 5.4 showed that a year of daily data leaves a single mean uncertain to ±27 % annualised. Assuming all 500 are equal (complete pooling) is biased. The hierarchical model \(\mu_i \sim N(\mu_0, \tau^2)\), \(r_{it} \mid \mu_i \sim N(\mu_i, \sigma^2)\) makes each stock’s mean a draw from a common distribution whose parameters are themselves estimated, and the posterior for each \(\mu_i\) is a shrinkage estimator: the stock’s own mean pulled toward the cross-sectional mean by a factor \(\sigma^2/n\) over \((\sigma^2/n + \tau^2)\) — the same precision-weighting as every update in this chapter, with the cross-section as the prior. A biased estimator with lower mean-squared error than the unbiased one: this is Stein’s paradox, Jorion’s Bayes–Stein estimator (1986), and the mechanism inside Black–Litterman (1991). Chapter 7 applies the same logic to the covariance matrix.

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 →