Learning Statistics with Python
Chapter 5 · Learning Statistics with Python
Priors, likelihoods and posteriors for risk and reward, credible intervals, Bayesian regression with predictive checks, robust regression with fat tails, dynamic updating.
Prof. Xuhu Wan
ISOM, HKUST Business School · 2026 Edition
Why this matters
Financial data are short, noisy and fat-tailed. Bayesian methods let you say “the probability that this strategy’s Sharpe ratio exceeds 1 is 99 %” — a sentence no p-value can produce — and let you carry beliefs forward as markets change.
| Section | Concept | Tool |
|---|---|---|
| 5.1 | Prior × likelihood → posterior; Sharpe ratio of a Student-t return model; credible intervals | grid approximation, scipy.stats.t.logpdf, Metropolis in numpy |
| 5.2 | Regression as a normal distribution; prior predictive and posterior predictive checks | conjugate normal-linear posterior, rng.multivariate_normal |
| 5.3 | Fat tails break OLS; Student-t likelihood repairs it | scipy.optimize.minimize, robust beta of GOOG on the S&P 500 |
| 5.4 | Sequential updating, forgetting, the local-level model | closed-form normal update, sm.tsa.UnobservedComponents |
You will build the posterior distribution of GOOG’s Sharpe ratio three ways — on a grid, then with a sampler you write yourself — and discover why the classical estimate is too small. Predict every number before you run it.
\[\underbrace{p(\theta \mid X)}_{\text{posterior}} \;\propto\; \underbrace{p(X \mid \theta)}_{\text{likelihood}}\;\underbrace{p(\theta)}_{\text{prior}}\]
Classical statistics fixes \(\theta\) and asks how the data would vary. Bayes fixes the data and asks how \(\theta\) varies.
After seeing 1 258 daily returns, which object lets you say “P(Sharpe > 1) = 0.99”?
GOOG daily returns, October 2015 to October 2020 — the notebook’s data. The annualised Sharpe ratio is \(\sqrt{252}\,\bar r / s\).
Mean daily return ≈ 0.00080, daily sd ≈ 0.01644. What does the classical Sharpe ratio print (3 dp)?
0.771
Sharpe 0.771, excess kurtosis 6.72. A normal has kurtosis 0 — these returns have far fatter tails. The model must say so.
\[r_t \sim \text{StudentT}(\nu,\ \mu,\ \sigma)\]
Three parameters: location \(\mu\), scale \(\sigma\), and degrees of freedom \(\nu\) that set the tail weight.
Which choice of \(\nu\) produces the fattest tails?
On a log scale the normal’s tail is a parabola falling off a cliff; the t(3) tail is almost a straight line. A 4-sigma event: 6.3e-05 under the normal, 0.028 under t(3) — about 440 times more likely.
The notebook’s choices, one per parameter:
| Parameter | Prior | Why |
|---|---|---|
| \(\mu\) | \(N(\bar r,\ s)\) | centred on the sample mean, very wide (sd = one daily sd) |
| \(\sigma\) | \(\text{Uniform}(s/1000,\ 1000\,s)\) | positive, essentially flat |
| \(\nu\) | \(2 + \text{Exponential}(\lambda = 1/29)\) | mean 31, but mass near 2 keeps fat tails possible; \(\nu > 2\) guarantees a finite variance |
What does the prior \(\nu = 2 + \text{Exp}(1/29)\) force?
Check it: expon.cdf(8, scale=29) = 0.241 — the prior gives only 24 % to \(\nu < 10\). It allows fat tails without insisting on them.
In Colab
PyMC is not available in the browser. This is the notebook’s model; the next slides reproduce its posterior with numpy and scipy.
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)Four lines of priors, one likelihood, one derived quantity. Everything else is the sampler — which we now build by hand.
With one or two parameters you do not need a sampler at all: evaluate prior × likelihood on a grid and normalise. This is Bayes’ rule made literal.
Same five years, S&P 500. Treat \(\sigma\) as known (the sample sd). Prior \(\mu \sim N(0,\ 0.001^2)\). For a normal likelihood the posterior has a closed form — predict what the grid should reproduce.
Posterior precision = prior precision + \(n/\sigma^2\). With \(n = 1258\), \(\sigma = 0.01186\), prior sd 0.001, what is the posterior sd (6 dp)?
0.000317
Grid 0.000497 0.000317, closed form 0.000497 0.000317 — identical to six decimals.
log_post - log_post.max() before exp is not cosmetic: 1 258 log-densities sum to about 3 800, and exp of that overflows.0.000497 sits between the prior mean 0 and the sample mean 0.00055: a precision-weighted average. The data’s precision (\(n/\sigma^2 \approx 8.9\) million) dwarfs the prior’s (1 million), so the data win — mostly.Two parameters? A 20 × 20 grid is still only 400 points. Let’s do \(\mu\) and \(\sigma\) for GOOG under a Student-t likelihood with \(\nu = 4\).
Posterior mean σ = 0.0108 — well below the sample sd 0.0164. The t-scale is not the standard deviation; fat tails inflate the sd. And Sharpe 2.035 with P(Sharpe > 1) = 0.972. Hold that thought.
A 400-point grid per parameter and a 3-parameter model (\(\mu, \sigma, \nu\)). How many likelihood evaluations?
The curse of dimensionality is why we need Markov chain Monte Carlo: instead of visiting every point, wander through parameter space so that time spent in a region is proportional to its posterior probability.
Propose a random step, accept it with probability \(\min(1, p^\star / p_t)\), otherwise stay. That is the whole algorithm — and it never needs the normalising constant.
At \(\theta_t\) the posterior density is \(p_t\); you propose \(\theta^\star\) with density \(p^\star\).
\(p^\star = 0.3\,p_t\). What does Metropolis do?
In code: if log(u) < log_post(prop) - log_post(current). Work in logs; the difference of two logs is a ratio of densities, so the unknown normaliser \(p(X)\) cancels.
Sample \(\theta = (\mu, \log\sigma, \log(\nu - 2))\) so the chain can never leave the allowed region. The extra ls and le terms are the Jacobians of those transforms.
3438.4 at the sample sd, 3443.2 at half of it. The likelihood already prefers a scale smaller than the sample sd — the tails are being handled by \(\nu\), not by \(\sigma\).
Acceptance 0.47 — inside the healthy 0.2–0.5 band. \(\nu \approx 2.7\): GOOG’s tails are heavier than a t(3) (the notebook’s PyMC run found 2.79). \(\sigma = 0.0097\), six-tenths of the sample sd.
What does a healthy trace look like?
Posterior mean 2.38, 95 % credible interval [1.30, 3.59], P(Sharpe > 1) = 0.99. The notebook’s PyMC run gave 2.35 and [1.31, 3.42] — the same answer from a twelve-line sampler. The classical 0.77 lies outside the credible interval. “The Sharpe ratio is under-estimated by classical statistics.”
The posterior mean Sharpe is 2.38 against a classical 0.77. The main reason?
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 Sharpe on that basis gives a posterior mean of about 1.18. The “under-estimation” is a statement about ordinary-day risk. Decide which risk you are pricing before you quote a Sharpe.
Lo (2002): the classical standard error of an annualised Sharpe is \(\sqrt{252\,(1 + \tfrac12 \widehat{SR}_d^2)/n}\).
Which interval lets you say “there is a 95 % probability the true Sharpe lies in here”?
[-0.11, 1.65][1.30, 3.59]Note the classical interval includes zero: a t-test cannot reject “GOOG earned nothing”. The posterior says P(Sharpe > 0.5) = 0.999.
The notebook’s practice. The starter runs the sampler on AAPL (2 515 returns). Switch r_p to NVDA from nv (2023–2024, 500 returns) and compare posterior mean and classical Sharpe for each.
Next: §5.2 — regression as a normal distribution, and how to check a model before fitting it.
A regression is a statement that \(y\) is normal with a mean that moves with \(x\). You will simulate from the prior before touching the data, discover which prior is absurd, then simulate from the posterior to see whether the model can reproduce what it was fitted to.
\[y_i \sim N(\mu_i,\ \sigma^2), \qquad \mu_i = \beta_0 + \beta_1 x_i\]
The notebook’s artificial data: \(\beta_0 = 1\), \(\beta_1 = 2\), \(\sigma = 0.5\).
In the Bayesian regression, which quantities get a prior?
y mean 2.041, sd 0.759. The sd of \(y\) (0.76) exceeds \(\sigma\) (0.5) because part of the spread is the line itself. Keep that distinction: spread of y ≠ noise.
In Colab
The notebook fits this with NUTS. We will get the same posterior in closed form on the next slides.
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)The posterior draws come back as a table: one row per draw, one column per parameter. Every plot in this section is built from that table.
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 found out for free.
\(x \in [-10, 10]\), \(y = 0.5 + 3x + N(0, 6^2)\), 200 points.
\(x\) has sd 5.8 and the slope is 3, noise sd 6. Roughly what is the sd of \(y\)? (\(\sqrt{3^2 \cdot 5.8^2 + 6^2}\))
18.4
Draw 100 lines with \(\beta_0, \beta_1 \sim N(0, 10)\). At \(x = 10\), what range covers 95 % of them?
Prior I: [-173.8, 209.7]. Prior II: [-20.8, 25.2]. The data at \(x = 10\) sit near 30. Prior II covers plausible slopes without covering nonsense — but it puts the true slope 3 at 2.5 prior sd. With 200 points the data will drag it there anyway; with 10 points they might not.
Normal likelihood, normal prior on \(\beta\): the posterior is normal and you can write it down. The formula is the same precision-weighting you saw on the grid, now in matrix form.
\[\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)\]
OLS on (xt, yt) gives slope 3.011. Prior II is \(N(0, 1.2^2)\) on the slope, and the data precision on the slope is about \(\sum x^2/\sigma^2 \approx 6700/32.5 \approx 206\). Will the posterior mean slope be nearer 3.011 or 0?
≈ 3.00 — data precision 206 vs prior precision 0.69
OLS [0.385, 3.011], posterior mean [0.331, 3.000], posterior sd [0.374, 0.069]. The intercept is shrunk a little toward 0 (its prior sd 1 competes with a data sd 0.4); the slope barely moves. With \(n = 200\) the prior is a rounding error.
The 50 red lines are 50 rows of post_df. Slope credible interval [2.861, 3.139] — the same information a classical CI gives, but readable as a probability.
Now go one step further: for each posterior draw, simulate a whole replicated data set. If the replicas do not look like the data, the model is missing something.
For posterior draw \(s\): \(\mu_i^{(s)} = b_0^{(s)} + b_1^{(s)}x_i\), then \(y_i^{\text{rep},(s)} \sim N(\mu_i^{(s)}, \sigma)\) for every \(i\).
yrep has shape (1000, 200). What is yrep[7]?
Thirty cyan replicas and the black data share one shape. If the data had a bump the replicas never produce, the normal likelihood would be wrong — this is how the §5.3 outliers will announce themselves.
A Bayesian p-value: the fraction of replicas whose statistic exceeds the observed statistic. Near 0.5 is good; near 0 or 1 means the model cannot produce what you saw.
p_mean is done. Fix p_max so it compares the maximum of each replica row with the observed maximum yt.max(). Expect about 0.13 — the normal model slightly under-produces the largest value.
p_max = 0.13 is the first hint that normal tails are thin.Next: §5.3 — put three outliers in and watch the normal likelihood surrender.
Three planted outliers will drag an OLS slope from 5 to 4.4. You will predict the direction, then swap the likelihood for a Student-t and watch the slope come back — first on simulated data, then on GOOG’s market beta.
\(y = 0.5 + 5x + N(0, 0.5^2)\) on 100 points, plus three points at \((0.1, 8), (0.15, 6), (0.2, 9)\) — far above the line on the left.
How will OLS respond to three high points at small \(x\)?
OLS [1.033, 4.403] — intercept doubled, slope off by 0.6. Under a normal likelihood a residual of 7 costs \(7^2/2 = 24.5\) log-units; the fit bends to reduce that at everyone else’s expense.
In Colab
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"])
# robust version: replace the last line by
likelihood = pm.StudentT("y", mu=mu, sigma=std, nu=3, observed=data2["y"].values)Same priors, one word changed: Normal → StudentT. Let’s see what that word does.
No sampler needed to see the effect: minimise the negative t log-likelihood with scipy.optimize.minimize, \(\nu = 3\).
Robust [0.492, 5.167], scale 0.457 — the truth is 0.5, 5, 0.5. The green dashed line goes through the cloud; the red one is tugged up-left by three points.
The t log-density’s derivative gives each point an implicit weight \(w_i = \dfrac{\nu + 1}{\nu + r_i^2}\), where \(r_i\) is the standardised residual.
The outlier at \((0.2, 9)\) has standardised residual ≈ 16.5 under the robust fit. With \(\nu = 3\), what is its weight (3 dp)?
0.015
Outliers get weights 0.017, 0.036, 0.015; a typical point gets 1.13. The t-likelihood is a self-weighting least squares — it down-weights exactly the points that do not fit. Nobody chose which points to drop.
Same sampler as §5.1, new log-posterior: Prior II from §5.2, \(\sigma \sim \text{Exponential}(1)\), t(3) likelihood.
Slope posterior mean 5.06, credible interval [4.64, 5.41]. OLS’s 4.40 lies outside it. The model is telling you the outliers are not from the same process.
Daily returns 2015-10 to 2020-10, 1 258 days including the March 2020 crash. Two likelihoods, one question: what is \(\beta\)?
OLS gives \(\beta \approx 1.06\). What will a Student-t likelihood give?
OLS β 1.064, robust β 1.094, \(\nu = 3.35\) — close to the 2.7 the return model found in §5.1. The point estimates differ by 0.03; the difference in inference is larger: a normal model with kurtosis-11 residuals reports standard errors that are fiction.
Twelve days with \(|r_{SP}| > 5\%\), nine of them in March 2020. Those twelve points — 1 % of the sample — carry 37 % of \(\sum x^2\), the leverage that drives OLS. The t model weights each by its surprise, not by its size.
250-day windows, stepped every 25 days. The OLS beta is done; replace the placeholder so b_rob is the t(3) slope from minimize (hint: mimic nll_t with xx, yy). Then compare the two paths.
Normal → StudentT) is the whole robust model; minimize or Metropolis fits it in under a second.Next: §5.4 — stop refitting from scratch; let yesterday’s posterior become today’s prior.
Every posterior you computed so far used all the data at once. Markets arrive one day at a time. You will update a belief day by day, discover that a static parameter cannot adapt, add one line to fix it — and find you have written the Kalman filter.
Normal mean, known \(\sigma\). One observation \(r_t\) updates \((m, v)\) — the posterior mean and variance of \(\mu\):
\[\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)\]
After 253 daily updates with \(\sigma = 0.02\) and prior sd 0.002, what is the posterior sd of \(\mu\)?
253 days, mean 0.00277, sd 0.02939. We fix \(\sigma = 0.02\) (known variance is the price of a closed form) and start from a prior \(\mu \sim N(0, 0.002^2)\).
Day 1 sd 0.00199, day 253 mean 0.00199, sd 0.00106 — and the batch formula gives exactly 0.00199, 0.00106. Sequential = batch. The posterior is a sufficient summary of the past.
During the March 2020 crash the posterior sd of μ…
The fix is one line: before each update, admit that \(\mu\) may have moved overnight.
\[\mu_t = \mu_{t-1} + \eta_t, \quad \eta_t \sim N(0, q) \qquad\Longrightarrow\qquad v_{t-1} \leftarrow v_{t-1} + q\]
With \(q > 0\) added every day, what happens to the posterior variance \(v_t\) in the long run?
Observation equation \(r_t = \mu_t + \varepsilon_t\), state equation \(\mu_t = \mu_{t-1} + \eta_t\): this is the local-level model, and the recursion you just wrote is the Kalman filter for it.
UnobservedComponents estimates both variances by maximum likelihood. First on the returns themselves.
Level variance 8.8e-13 — effectively zero. The data say the mean of AAPL’s daily return does not wander: a constant is as good as it gets. That is a finding, not a failure. Risk is a different story.
\(\sqrt{\pi/2}\,|r_t|\) is an unbiased daily-volatility proxy under normality. Give it a wandering level.
Filtered vol 1.2 % on 2 Jan 2020, 8.8 % on 16 Mar, back to 1.7 % by year-end. Signal-to-noise q = 0.025: the level moves, but each day’s |r| is mostly noise.
In steady state the filter’s update is \(m_t = m_{t-1} + K\,(y_t - m_{t-1})\) with gain \(K = \tfrac{-q + \sqrt{q^2 + 4q}}{2}\).
With \(q = 0.0251\), what is the steady-state gain \(K\) (3 dp)?
0.146
Your five-line recursion reproduces filtered_state to 5e-07. \(K = 0.146\) is an EWMA with a 12.7-day span — RiskMetrics’ \(\lambda = 0.94\) (span 32) is the same filter with a smaller \(q\). Chapter 6 adds trends, seasons and regressors to this state.
The notebook’s practice: 20 days set the prior, then update every 10 days with the posterior of the previous block as the new prior. The starter resets the prior every block. Move the reset out of the loop so beliefs carry forward. sd_final should then equal \(\sigma/\sqrt{253}\) = 0.00126: a prior built from 20 days carries exactly 20 observations of information, so block updates end where a one-shot update would.
0.00199 ± 0.00106 exactly — the posterior is a sufficient summary of the past.UnobservedComponents to better than \(10^{-6}\); the gain \(K = 0.146\) is an EWMA with a 12.7-day span.Next: Chapter 6 — the same state-space idea with trends, ARIMA and the full Kalman filter.
Warning
On the eve of the 8 November 2016 US presidential election the Princeton Election Consortium gave Clinton a 99 % chance of winning and the HuffPost model 98 %; FiveThirtyEight said 71 %. National polls had her about 3 points ahead and she won the popular vote by 2.1 — an ordinary polling miss. But Wisconsin, Michigan and Pennsylvania, polling at Clinton +2 to +6, went to Trump by under a point, and the Electoral College 306–232.
The 99 % models treated the fifty state polling errors as independent, so a 2-point lead became near-certainty as the errors averaged away. FiveThirtyEight’s likelihood gave the errors a common component — one demographic miss shared across the Rust Belt — and its posterior stayed wide.
Lesson: a posterior’s width is only as honest as the dependence structure in the likelihood. Errors modelled as independent when they are really one error shrink the credible interval by \(\sqrt{n}\). As with Student-t tails, run the predictive check before you print 99 %: could this model have generated the miss you are about to see?
The posterior of the Sharpe ratio is not the deliverable. A recommendation with a probability attached is.
To: CIO, Multi-Strategy Fund From: <Your name>, quantitative analyst Subject: Allocation to the GOOG single-name sleeve Date: 2026-09-15
Recommendation: Allocate 5 % of risk budget; review after 250 trading days.
Evidence (Student-t model, 1 258 daily returns 2015–2020, Metropolis posterior): - Posterior mean annualised Sharpe 2.38; 95 % credible interval [1.30, 3.59]; classical CI [−0.11, 1.65] includes zero. - P(Sharpe > 1) = 0.99; P(Sharpe > 0.5) = 1.00. - Tail index ν ≈ 2.7: a 4σ day is 440 × more likely than a normal model says.
Caveats: - Sharpe uses the t-scale σ; on the standard-deviation basis the posterior mean is ≈ 1.18. - The local-level filter (§5.4, AAPL) finds a constant mean but a volatility that moved from 1.2 % to 8.8 % in Q1 2020 — size positions on filtered vol, not the sample sd. - Single name; no cross-sectional shrinkage applied (Chapter 5 of the book).
Next step: re-run the posterior every 10 days with the previous posterior as prior; trigger a review if P(Sharpe > 0.5) falls below 0.8.
Three prompts that work with this chapter’s tools — and the pitfalls each one hides.
“Write a Metropolis sampler for this log-posterior and report the acceptance rate.” Good. Then check the trace yourself: an LLM will happily call a 0.98 acceptance rate “excellent” when it means the step size is far too small and the chain has explored nothing.
“Which prior should I use for ν?” It will suggest Exponential(1/29) + 2 or Gamma(2, 0.1) — both fine — but ask it to draw prior predictive samples and show you P(ν < 5). The prior you accept should be one you have seen, not one you were told.
“Is the Bayesian Sharpe of 2.38 better than the classical 0.77?” The copilot will say yes. It is answering a different question. Ask instead: “which σ does each number divide by, and which one matches my risk limit?” — the answer in this chapter was “they measure different risks”, and no model chooses your risk definition for you.
| Concept | Tool |
|---|---|
| Posterior ∝ likelihood × prior | np.exp(log_post - log_post.max()) on a grid |
| Student-t return model, Sharpe posterior | scipy.stats.t.logpdf, Metropolis loop with rng.uniform() |
| Credible interval | np.percentile(draws, [2.5, 97.5]) |
| Prior / posterior predictive check | rng.normal(...) on parameter draws; Bayesian p-value |
| Conjugate normal-linear posterior | Sn = inv(inv(S0) + X.T @ X / σ²), rng.multivariate_normal |
| Robust regression | minimize(nll_t) with student_t.logpdf; weights \((\nu+1)/(\nu+r^2)\) |
| Sequential update / local level | precision recursion; sm.tsa.UnobservedComponents(y, level="local level") |
In Colab the same models are four lines of PyMC each — now you know what pm.sample does.
Next: Chapter 6 — Time Series Models for the Mean.
Your fund reports Sharpe ratios on a standard-deviation basis to investors but sizes positions on the t-scale σ. Which posterior should go in the marketing deck, and what has to be disclosed alongside it?
The prior predictive check rejected Prior I as “absurd”, yet with 200 points Prior I and Prior II give nearly identical posteriors. When does the choice of prior actually change a decision — and how would you demonstrate that to a sceptical PM?
The Student-t fit gave the March 2020 days weights near 0.02. A risk manager objects: “those are exactly the days I care about”. Is robust regression the wrong tool for risk, the right tool for alpha, or both?
The local-level MLE said AAPL’s mean return does not wander but its volatility does. What would you expect for a fund’s strategy returns, and how would you set \(q\) if you had only 60 days of live track record?
Prof. Xuhu Wan · HKUST ISOM · Learning Statistics with Python