Learning Statistics with Python
Chapter 3 · Learning Statistics with Python
From averages to distributions: empirical and theoretical distributions, bootstrap confidence intervals and hypothesis tests, experimental design with A/B tests and bandits, linear and nonlinear association, and extreme value theory for the tails.
Prof. Xuhu Wan
ISOM, HKUST Business School · 2026 Edition
Why this matters
Every quant decision is a statement about a distribution — its centre, its spread, its shape, its tail. Chapter 2 gave you clean frames; this chapter gives you the habit of predicting what the distribution should look like before you compute it, so that a 25-standard-deviation “surprise” is recognised for what it is: a wrong model, not bad luck.
| Section | Concept | Tool |
|---|---|---|
| 3.1 | Population vs sample, ECDF, KDE, parametric fit, Q-Q, normality tests | sample, np.sort, KernelDensity, stats.t.fit, probplot, shapiro, kstest, bimodality coefficient |
| 3.2 | Bootstrap SE and CIs, one-sample and two-sample tests, permutation, regime detection | rng.integers(0, n, (B, n)), percentile and BCa, ttest_1samp, rng.shuffle, ks_2samp, rolling \(\hat t\) |
| 3.3 | Power, sequential testing, multi-armed bandits | \(z\)-test, \(n = 2\sigma^2(z_\alpha+z_\beta)^2/\tau^2\), peeking simulation, ε-greedy, UCB1, Thompson |
| 3.4 | Linear, monotonic, nonlinear and tail association | pearsonr, spearmanr, kendalltau, distance correlation via pdist, tail co-exceedance |
| 3.5 | Extreme value theory: GEV, block minima, POT, extreme VaR | genextreme, genpareto, Hill estimator, return levels, EVaR |
A population parameter is a fact; a sample statistic is an estimate with a distribution of its own. You will build the empirical CDF from five numbers, read Apple’s fat tails off it, estimate a density two ways (bandwidth vs a formula), learn to see a Student-t in a Q-Q plot, and then let Shapiro–Wilk and Kolmogorov–Smirnov put numbers on what the eye saw.
The notebook’s 63 numbers are the whole population. We draw 10 without replacement (seed 1).
Why does pandas compute std() with ddof=1 by default?
ddof=0Log return \(r_t = \ln S_t - \ln S_{t-1}\). Before running: is the excess kurtosis of daily returns positive or negative?
The excess kurtosis of Apple’s daily returns will be…
Excess kurtosis 5.99 and skew −0.28. Half the days sit between the two red lines — a band about 2 % wide — yet the axis runs from roughly −20 % to +14 %. That gap between the middle and the extremes is what the next slides quantify.
Textbooks hand you a density. Real work hands you a sample. The bridge between them is one line of arithmetic — you will build it from five numbers and then read Apple’s tails off it.
You have the sample [3, 1, 4, 1, 5]. Define \(\hat F(x)\) = (fraction of the sample \(\leq x\)). Before any code, predict two values:
The object you built has a name. Given an i.i.d. sample \(X_1,\dots,X_n\), the empirical CDF is \[\hat F_n(x) = \frac{1}{n}\sum_{i=1}^{n}\mathbf 1\{X_i \le x\}.\] A step function. No assumptions about shape. It is the non-parametric maximum-likelihood estimate of the true CDF — and you computed it by hand.
A sorted sample’s \(i\)-th value sits at height \(i/n\). Studentise the 4 309 returns (\(z = (r - \bar r)/s\)) and sort them.
The ECDF height at the 2 155th of 4 309 sorted points — what is \(2155/4309\) to three decimals?
0.5
The median studentised return is −0.008, essentially zero. But 76.9 % of days lie within one standard deviation (a Normal says 68.3 %) and only 98.5 % within three (a Normal says 99.7 %): a taller peak and fatter tails. That is excess kurtosis 5.99, read off a staircase.
The notebook’s mixture sample (seed 7) is two Normals stuck together, one at −10 and one at +5. Predict what its ECDF looks like before you draw it.
The data is bimodal (two clusters). What will the ECDF curve look like?
The curve climbs steeply through −10, then flattens between about −5 and 0 (the first quartile is −10.0, the median −5.2, the third quartile 5.0), then climbs again around +5. The rule you discovered: the ECDF has no tuning knob, so it cannot lie about shape. The flat stretch is the bimodality.
The ECDF is built from \(n\) random points, so it is itself random. Will it land on the truth as \(n\) grows? Predict the worst-case gap for Apple’s \(n = 4\,309\).
The DKW bound on \(\Pr(\sup_x|\hat F_n - F| > 0.02)\) is \(2e^{-2n(0.02)^2}\). Evaluate it at \(n = 4309\), rounded to 3 decimals.
0.064
Glivenko–Cantelli (1933): \(\sup_x|\hat F_n(x) - F(x)| \xrightarrow{a.s.} 0\) — the ECDF converges uniformly to the truth. The DKW inequality is the finite-\(n\) version: with 4 309 days, the whole Apple staircase is within 0.02 of the true CDF with probability at least 93.6 %. This theorem is what the bootstrap (§3.2), the KS test and the KDE all lean on.
The ECDF is a staircase. To get a smooth density, smear each point into a small bump and add them up. The width of the smear, \(h\), is the only choice that matters — and you will discover how much it matters.
A kernel density estimate puts a small bump on every observation and adds them up. The bandwidth \(h\) is the only knob. Predict its effect on the notebook’s two-component mixture.
For a two-mode sample, what does a very small bandwidth do?
Silverman’s rule of thumb gives a starting bandwidth: \(h^\star \approx 1.06\,\hat\sigma\,n^{-1/5}\). Predict it for the studentised Apple returns (\(\hat\sigma = 1\), \(n = 4309\)).
Compute Silverman’s \(h\) for \(z\) (sd 1, \(n = 4309\)), rounded to 3 decimals.
0.199
\(h^\star \approx 0.199\) on the studentised scale — a fifth of a standard deviation — or 0.0041 on raw returns. But Silverman assumes one bump: on the mixture above it would over-smooth. Always plot a sweep; never trust one number.
Studentise the returns, then let scipy.stats.t.fit choose the degrees of freedom. Predict the meaning of the number before you see it.
If the fitted degrees of freedom come out near 3, what does that say?
The notebook fits the same \(t\) to TSLA, CVX and BAC (df between 3 and 6 for every large stock). Do it for Tesla in 2020.
Studentised 2020 Tesla returns are in z_tsla. Set nu_tsla to the degrees of freedom returned by stats.t.fit (its first element).
A Q-Q plot puts sample quantiles against a reference distribution’s quantiles. A perfect match traces the diagonal. Apple’s returns have df ≈ 3.3 tails; predict what that does to the ends against a Normal reference.
The sample is heavier-tailed than Normal. The Q-Q plot will…
A Q-Q plot compares sample quantiles with a reference distribution’s quantiles. Straight line = same shape.
Against the normal the tails peel away from the line on both sides — the classic S. Against the fitted \(t\) the points hug the line until the very last few, where the right tail is shorter than a \(t_{3.3}\) predicts.
The Q-Q plot showed where the departure is; a test puts a number on it. The notebook runs kstest(returns, stats.norm.cdf) and gets \(D \approx 0.47\). Predict why it is that large.
Why is the KS statistic 0.47 when raw returns are compared with stats.norm.cdf?
Reading the tests
Null hypothesis: the population is normal. A tiny p-value rejects it: Shapiro–Wilk \(p \approx 10^{-38}\); KS on the studentised returns \(D = 0.075\), \(p \approx 10^{-21}\). Against the fitted \(t_{3.3}\) the KS test does not reject (\(D = 0.015\), \(p = 0.27\)) — that is what “the t fits” means in numbers. Shapiro–Wilk is the more powerful normality test; KS is the more general one — any CDF, or another sample (§3.2 uses the two-sample form).
Before fitting anything, the bimodality coefficient \(b = (g_1^2 + 1)/(g_2 + 3)\) flags trouble when \(b > 5/9 \approx 0.555\). Predict the verdict on the mixture, then on Apple.
For the seed-7 mixture, \(g_1 = 0.381\) and \(g_2 = -1.214\). Is \(b\) above 0.555?
True
Mixture: \(b = 0.641\), flagged — negative kurtosis (−1.21) is the fingerprint of two separated humps. Apple: \(b = 0.120\), not flagged — one peak, fat tails. The detector tells you which problem you have before you choose a density.
ddof=1 corrects the downward bias of the sample variance.norm on raw returns tests scale, not shape — studentise first (\(D = 0.075\)); the fitted \(t_{3.3}\) is not rejected (\(D = 0.015\), \(p = 0.27\)).Next: §3.2 — how sure are you of any of those numbers? Resample, and test.
You have one sample. The bootstrap resamples it to reveal how much any statistic would wobble — no formula needed. Then the classical tests: was a year’s mean return zero, was one year different from the next, do two stocks share a distribution — and a two-sample test recomputed every day that becomes a regime detector.
Each section runs standalone. Reload Apple 2005–2022, studentise, refit the \(t\).
You want the sampling distribution of some statistic but cannot draw new samples. Classic statistics needs a formula for the standard error. The bootstrap (Efron, 1979) needs none — and you will see why by reinventing it.
You have a sample \(X_1,\dots,X_n\) and want to know how much \(\hat\theta\) would wobble across new samples — but you cannot draw new samples. Predict the bootstrap’s move.
With no access to the population, how do we mimic drawing a fresh sample?
The ECDF \(\hat F_n\) stands in for the unknown \(F\). Resampling with replacement is sampling from \(\hat F_n\). That single substitution is the whole idea — and Glivenko–Cantelli is why it works.
The recipe in code: build a \(B \times n\) matrix of resample indices, compute the statistic on each row, take the spread. A small skewed sample — 120 lognormal draws — and a statistic with no textbook SE: the median.
Bootstrap SE of the median, 120 lognormal draws (seed 11), \(B = 4000\). Order of magnitude first: 0.02, 0.2, or 2?
0.198
The trick you found: the B x n index matrix replaces every loop and every SE formula. There is no closed-form SE for a median — yet you just measured it: 0.198.
For the mean there is a formula: \(\mathrm{SE} = s/\sqrt n\). If the bootstrap is legitimate, its SE must match. Confirm it on Apple’s 251 returns of 2017.
Bootstrap the mean of r17 (\(B = 4000\), seed 5) into se_boot, and compute the textbook se_formula. They should agree within 10 %.
Both come out at about 0.0007 — so the bootstrap is not magic; it reproduces the answer you already trusted, then extends to medians, ratios and kurtosis, where no formula exists.
The simplest CI: take the 2.5 % and 97.5 % quantiles of \(\hat\theta^*\) directly. Predict whether the interval for the 2017 mean daily return contains zero.
The 95 % percentile CI for the 2017 mean daily return…
\([0.00014,\ 0.00294]\): zero is outside, barely. Clean and correct when the bootstrap distribution is symmetric and unbiased. Predict the next slide: when is that assumption dangerous?
For which statistic is the plain percentile CI least trustworthy?
BCa (bias-corrected and accelerated) fixes two blind spots: a bias correction \(z_0\) (how often \(\hat\theta^* < \hat\theta\) — should be 50 %) and an acceleration \(a\) from the jackknife (handles skew). Same cost as percentile.
Excess kurtosis is the statistic §3.1 warned you about — dominated by a handful of days. Bootstrap it (B = 2000) and compare the two intervals.
Resampling with replacement drops some of the extreme days, so most bootstrap kurtoses fall below the sample value. Is \(z_0 = \Phi^{-1}(\text{share below})\) positive or negative?
True
Kurtosis 4.58 with a bootstrap SE of 1.55 — a third of its own size. 59.4 % of resamples fall below it, so \(z_0 = +0.24\), and BCa pushes the interval from \([1.47, 7.31]\) up to \([2.23, 9.15]\). Rule discovered: for medians, quantiles, ratios and kurtosis the bias and skew are real, and BCa moves the endpoints to where they belong — for free.
A confidence interval says how much a number wobbles; a test asks whether a specific claim survives — the null / statistic / p-value logic of §3.1, now applied to a mean. Three claims: 2017’s mean return was zero, 2017 and 2018 share a mean, MSFT and AAPL share a distribution — then one test recomputed every day.
\(H_0: \mu = 0\). Predict the statistic \(\hat t = \bar x / (s/\sqrt n)\) for 2017.
With \(\bar x \approx 0.00151\), \(s \approx 0.0111\) and \(n = 251\), what is \(\hat t\) to two decimals?
2.16
The cell printed a one-sided \(p\) (Normal approximation) and ttest_1samp’s two-sided \(p\). The business question — “did Apple earn a positive return?” — is directional.
For a positive \(\hat t\), the one-sided \(p\) compared with the two-sided \(p\) is…
Both reject at 5 % — the bootstrap CI excluded zero for the same reason. The interval, the \(t\)-statistic and the \(p\)-value are three views of one fact: \(\bar x\) is 2.16 standard errors from zero.
The \(t\)-test assumed a Normal sampling distribution for \(\bar x\). For a sharp null — “the two groups have identical distributions” — you need no such assumption: build the null distribution by shuffling labels.
Under \(H_0: F_A = F_B\), why can we randomly relabel the pooled data?
No Normality, no equal-variance assumption, no large-\(n\) asymptotics. The price is \(B\) shuffles — trivial on a laptop.
Apple’s mean daily return was +0.151 % in 2017 and −0.028 % in 2018. Label the 502 returns by year, shuffle the labels 5 000 times, and ask how often a gap that large appears by chance.
The two-sided permutation \(p\) for the 2017 − 2018 gap will be roughly…
\(p = 0.18\) from shuffles, 0.18 from Welch’s \(t\) — the classical test was fine here because \(n\) is large. The permutation test’s advantage is that it would still be right for 20 days of fat-tailed returns, where the \(t\) is not. Its Monte-Carlo error shrinks like \(1/\sqrt B\), independent of \(n\).
The permutation test compared two means. ks_2samp compares two whole empirical CDFs and needs no reference distribution at all. The notebook’s MSFT-vs-AAPL check, on the catalog files’ common dates.
Null: the two samples come from the same distribution. On 778 common days \(D = 0.039\), \(p = 0.61\) — no evidence that MSFT and AAPL daily returns are distributed differently. Keep the p-value in mind — the same two-sample idea, run on a rolling window, becomes a regime detector on the next slide.
Compare the last 10 days’ returns with the 10 days ending 20 days earlier: \(\hat t = \dfrac{\bar r_{\text{now}} - \bar r_{\text{lag}}}{\sqrt{(s^2_{\text{now}} + s^2_{\text{lag}})/10}}\). Days with \(|\hat t| > 2\) are flagged.
Twenty of 249 days are flagged; the red dots cluster at the September 2012 peak and the April 2013 trough — the points where the return distribution shifted. A test statistic, recomputed every day, is a trading signal.
For each 20-day window, store the p-value of stats.ks_2samp(w_now, w_lag) in ks_p (the test compares the current 20 returns with the 20 returns lagged by 20 days). Then count windows with p < 0.2.
ks_2samp cannot tell MSFT from AAPL (\(D = 0.039\), \(p = 0.61\)). A rolling two-sample \(t\) turns “same distribution?” into a regime signal: 20 flagged days around the 2012 peak and 2013 trough.Next: §3.3 — tests you design: how many users, when to stop, and when not to test at all.
So far the data arrived and you tested it. Now you design the experiment: how many units to catch a one-point lift, what peeking at the dashboard does to your false-positive rate, and when earning while learning (a bandit) beats a fixed test. These are simulations by necessity — you cannot rerun a market.
Two arms, fixed \(N\), a hypothesis fixed before the data. The design is a century old (Fisher’s agricultural plots) and still the gold standard when you can afford it. You will discover its single most-violated requirement.
Control outcome \(\sim N(0.10, 1)\), treatment \(\sim N(0.18, 1)\), 800 each. The \(z\)-statistic is \(\hat\tau / \widehat{\mathrm{SE}}\), reject if \(|z| > 1.96\). Predict whether this experiment detects the true 0.08 lift.
There is a real 0.08 lift. Does \(|z| > 1.96\) (significant)?
A genuine effect, declared “not significant.” The discovery: the failure is not the effect — it is the sample size. Significance is a property of power, and power is something you size before you run.
For a baseline rate \(p = 0.10\) and an absolute lift of \(0.01\) (10 % → 11 %), the formula \(n \approx 2\sigma^2(z_\alpha + z_\beta)^2/\tau^2\) at 80 % power. Predict the order of magnitude.
Required \(N\) per arm for \(p = 0.10\), lift \(= 0.01\), \(z_\alpha = 1.96\), \(z_\beta = 0.84\). Integer (use math.ceil).
14112
About 14 000 per arm to catch a one-point lift. Most “failed” experiments never had the sample size to succeed. The previous slide’s \(n = 800\) was doomed before it ran.
The required \(N\) scales with \(1/\tau^2\): doubling the detectable lift should cut \(N\) by four.
Set n2 to the required \(N\) per arm (integer, math.ceil) for the same baseline \(p = 0.10\) but a lift of 0.02. Check that it is one quarter of 14 112.
Real teams peek — they check the dashboard at lunch and stop when “it looks significant.” You will discover, by simulation, exactly how much damage that does to the false-positive rate.
We simulate under a true null (no effect). At each new observation we recompute \(z\) and stop the first time \(|z| > 1.96\). Predict the resulting false-positive rate over 1 000 observations.
True null, peek-and-stop at every observation up to \(N = 1000\). False-positive rate?
Over half the null experiments produce a “significant” result. The rule you discovered: every look spends Type-I error. Naive peeking turns a 5 % test into a coin flip.
If repeated looks spend error, the cure is to budget it across looks. Predict which spending rule lets you peek hardest early.
Which sequential boundary is stringent early, relaxed later — the clinical-trial default?
Note
The cost of a properly-designed sequential test: ~5–15 % more samples than a fixed-\(N\) test of equal power. The gain: a principled option to stop early when the effect is large. Lan–DeMets \(\alpha\)-spending generalises this to any look schedule.
When the goal flips from “estimate the effect precisely” to “earn the most reward while learning,” A/B is the wrong tool. Bandits trade exploration against exploitation. You will race three policies head-to-head.
Three arms with true conversion rates [0.08, 0.10, 0.12]; arm 2 is best. Over 4 000 pulls we race \(\varepsilon\)-greedy, UCB1, and Thompson sampling. Predict which earns the most.
Which policy collects the most total reward here?
A perfect oracle expects \(0.12 \times 4000 = 480\); the gap is regret. Thompson lands closest — its pull counts pile onto arm 2 fastest. \(\varepsilon\)-greedy keeps wasting 10 % of pulls on losers forever. That is the exploration–exploitation trade-off, made concrete.
The design is decided before the data, and it depends on what the business actually pays for. Predict the right tool for each scenario.
You are picking ad creatives; the revenue during the test counts. Best tool?
Fixed-\(N\) A/B when you need a defensible point estimate (board, regulator) and the cost of a wrong call is one-time and large. Sequential when stakeholders will peek or units are expensive. Bandit when learning is earning.
You now know peeking inflates Type-I error. Predict another silent killer of experiments.
You measure 12 KPIs and declare victory on whichever turns significant. The bug?
Note
Other traps: SUTVA violations (treatment leaks between units → cluster-randomise / switchbacks); novelty effects (a new UI shines for a week → pre-register a 30-day hold-out); CUPED variance reduction can cut required \(N\) by 30–60 % at zero cost. Decisions made after seeing the numbers are p-hacking by another name — the same lesson as §7.6’s 42 rules on trial.
Next: §3.4 — two variables at once: which correlation coefficient sees a parabola?
Pearson’s \(r\) measures only the linear part of a relationship. You will build data that is strongly related yet fools Pearson completely, find which coefficient survives a monotone bend, which survives any bend at all — and then ask whether the association between two stocks is the same on the worst days as on ordinary ones.
We make \(y = (x-5)^2 + \text{noise}\) — a clean, strong, non-linear relationship. Predict what Pearson’s \(r\) reports.
Pearson \(r\) between \(x \sim U(0,10)\) and \(y = (x-5)^2 + \varepsilon\) (seed 61, \(n = 300\)). Near 0, near 0.5, or near 1?
0.04
\(r \approx 0.04\) on a relationship you can see with your eyes. Pearson is blind to anything non-monotone. Its requirements — finite variance, linearity, outlier sensitivity — are assumptions, not guarantees.
\(y = \sin x + \varepsilon\) is perfectly determined by \(x\) up to noise. Predict Pearson’s \(r\).
For y = sin(x) + noise over \(x \in [0, 100]\), Pearson’s \(r\) is…
In the index panel spy is tomorrow’s SPY change and the other columns are today’s moves: every correlation in the spy row is near zero. Today’s Hang Seng is not a linear predictor of tomorrow’s SPY. “No association” here means no linear association.
Spearman (\(\rho_S\), Pearson on ranks) and Kendall (\(\tau\), concordant minus discordant pairs) catch any monotone relationship. Predict whether they catch a parabola.
For a symmetric parabola, Spearman and Kendall will report…
Spearman’s \(\rho\) is Pearson’s \(r\) on the ranks. Predict it for \(y = \log x\).
\(\log\) is strictly increasing, so the rank of \(y\) equals the rank of \(x\) for every point. What is Spearman’s \(\rho\)?
1.0
Spearman rescues monotone nonlinearity (log: Pearson 0.92 → Spearman 1.00) and nothing else: the parabola (0.08 / 0.06) and the cosine (0.04 / 0.03) are invisible to both coefficients.
Kendall counts pairs: \(\tau = (\text{concordant} - \text{discordant}) / \binom{n}{2}\). It is more robust to ties and outliers than Spearman, and just as blind to a bend.
Log: \(\tau = 1.000\); parabola: \(0.039\); cosine: near zero. The rule discovered: rank measures fix outliers and curvature-but-monotone, not non-monotonicity. For “any dependence at all” you need distance correlation (Székely, 2007), which is zero if and only if \(X \perp Y\).
Székely’s distance correlation is zero if and only if the variables are independent. The notebook uses the dcor package; we build it from pairwise distances.
Which measure will be clearly positive for the parabola?
Distance correlation: log 0.97, parabola 0.50, cosine 0.37 — the notebook’s dcor gives 0.97 / 0.49 / 0.35 on its own draws. Only this column sees the two non-monotone shapes.
feat holds next-day return (target), today’s return (lag1) and the 5-day sum (ret5) for 2020–2021. Set dc_lag1 to the distance correlation between lag1 and target, and compare with Pearson.
returns.csv holds daily S&P 500 and Tesla returns, 2015–2024. One Pearson number describes ten years. Predict what happens conditionally on the S&P’s worst 5 % of days.
On the S&P’s worst 5 % of days, what share are also among Tesla’s worst 5 %?
Pearson 0.465 over 2 515 days; on the 126 worst S&P days Tesla was in its own worst 5 % on 32.5 % of them — 6.5 times the independence rate. (Pearson within those days, 0.355, is lower — conditioning on a slice of \(x\) truncates its variance and mechanically shrinks \(r\).) The tail is where diversification is tested, and §3.5 is about the tail.
pdist.Next: §3.5 — the tails alone: what is the worst 20-day loss you should plan for?
Averages have the central limit theorem; maxima have Fisher–Tippett–Gnedenko. You will first see that matching mean and variance does not match the tail, then collect the worst day of every 20-trading-day block on the Dow (1985–1990, so Black Monday is in the sample), fit a GEV by maximum likelihood, check it with a Q-Q plot, score each block’s rarity, and read off the extreme VaR — the highest pain threshold.
Draw 4 309 Normal returns with Apple’s mean and standard deviation (seed 71). They match on the first two moments. Predict the 0.1 % quantile comparison.
With mean and sd matched, Apple’s 0.1 % quantile versus the Normal’s is…
Same mean, same sd — and a 0.1 % quantile of −10.5 % against −5.8 %, a 1 % CVaR of −7.8 % against −5.1 %, and 25 days beyond four standard deviations where the Normal has none. Variance does not measure tail risk. Two portfolios with identical Sharpe ratios can have wildly different blow-up probabilities.
If the standardised maximum \(M_n\) of \(n\) i.i.d. draws converges at all, it converges to the GEV family:
\[F(s) = \begin{cases} \exp\!\big(-e^{-s}\big) & \xi = 0 \\[4pt] \exp\!\big(-(1 + \xi s)^{-1/\xi}\big) & \xi \ne 0 \end{cases} \qquad Q(p) = \begin{cases} -\ln(-\ln p) & \xi = 0 \\[4pt] \dfrac{1}{\xi}\big((-\ln p)^{-\xi} - 1\big) & \xi \ne 0 \end{cases}\]
with \(s = (x - \mu)/\sigma\). Three sub-families, indexed by the tail index \(\xi\).
Which sign of \(\xi\) means a heavy, unbounded tail?
scipy’s sign convention
scipy.stats.genextreme uses shape c with c = −ξ. Every fit below reports ξ = −c.
With \(\mu = 0, \sigma = 1\) the support ends at \(\mu - \sigma/\xi\): \(= 3\) for \(\xi = -1/3\) (right end), \(= -3\) for \(\xi = 1/3\) (left end). The Fréchet curve is the one that keeps going.
Which date is the worst single day in 1985–1990, and roughly what log return? (Think October 1987.)
1987-10-19, −0.2563
rolling(20).min() gives the running 20-day minimum; taking every 20th row (iloc[::20]) makes the blocks non-overlapping. Predict the number of blocks.
Roughly 1 500 trading days in 1985–1990. How many 20-day block minima?
One red dot per block, whatever happened in it: a calm month contributes a minimum of −0.5 %, October 1987 contributes −25.6 %. Block maxima treat both as one observation each — the method is rate-limited by design.
POT throws away the calendar and keeps every exceedance. Predict the order of magnitude: how many of the 1 496 days fell below −1.5 %?
How many days in 1985–1990 have a log return below −0.015?
The dots now cluster in October 1987 and October 1989 — POT captures clusters of stress that block minima flatten to one point. Its price is choosing \(u\): too high and the sample is tiny, too low and the GPD approximation fails.
Before real data, check the theorem on a case where you know the answer (the notebook’s i.i.d. Normal experiment).
10 000 standard-Normal draws, 100 blocks of 100, take each block’s maximum, fit a GEV. Which \(\xi\) should you expect?
Minima become maxima by a sign flip: fit to \(-\text{Min20}\), the block loss. Predict the support bound.
For \(\xi > 0\) the GEV support starts at \(\mu - \sigma/\xi\). With \(\mu \approx 0.0136\), \(\sigma \approx 0.0070\), \(\xi \approx 0.354\), what is the bound?
−0.0063
\(\xi \approx 0.35 > 0\): Fréchet. The 20-day worst loss of the Dow has a polynomial tail; moments of order above \(1/\xi \approx 2.8\) do not exist.
One point will sit far above the 45° line at the top right. Which block is it?
Score \(= 1 - F_{\text{GEV}}(\text{loss}) = P(\text{a block is worse than this one})\). Small score = rare block.
Which date has the smallest score, and how many blocks score below 0.10?
1987-10-19, 6
Black Monday’s score is 0.0007: the fitted GEV says a worse 20-day block arrives once in about 1 400 blocks — 110 years. Two others (January 1988, October 1989) are 1-in-50 events.
Exceedances over a high threshold follow a GPD (Pickands–Balkema–de Haan): \(\;P(X - u > y \mid X > u) \approx \big(1 + \xi y/\beta\big)^{-1/\xi}\), with the same \(\xi\) as the GEV. The POT quantile is
\[\text{VaR}_q = u + \frac{\beta}{\xi}\Big[\Big(\tfrac{n}{N_u}(1-q)\Big)^{-\xi} - 1\Big].\]
91 exceedances give \(\xi = 0.44\) against the GEV’s 0.35 — the same heavy-tail verdict from a different sample; agreement between the two estimators is your sanity check. The GPD’s 99 % one-day VaR (3.3 %) sits just above the empirical 1 % quantile (3.1 %): with 1 496 days the empirical quantile is fine at 1 %. The GPD earns its keep at 0.1 %, where there is nothing left to count.
A Fréchet tail is a power law, \(P(L > x) \sim x^{-\alpha}\) with \(\alpha = 1/\xi\). Hill (1975) reads \(\alpha\) off the top \(k\) order statistics: \(\hat\alpha_k = \big[\tfrac{1}{k}\sum_{i\le k}\ln(L_{(i)}/L_{(k+1)})\big]^{-1}\). Predict which \(k\) agrees with the GPD’s \(\xi = 0.44\).
Where will the Hill estimate of \(\xi\) form a plateau near 0.44?
\(\hat\xi\) = 0.47, 0.42, 0.44 for \(k\) = 25, 50, 100 — the plateau — then 0.63 at \(k = 200\), where the threshold (0.8 %) is an ordinary day. Picking \(k\) is the craft, exactly as picking \(u\) was for POT: too few points is noise, too many is the body of the distribution.
The \(T\)-day return level solves \(P(L > x_T) = 1/T\). Two estimators: the empirical quantile, and the fitted GPD extrapolated: \(x_T = u + \frac{\beta}{\xi}\big[(\tfrac{k}{n}T)^{\xi} - 1\big]\). Predict the empirical method’s fatal limit.
Why can’t an empirical quantile give a 1-in-2 500-day loss from 1 496 days?
At \(T = 250\) (one year) the two agree: 4.9 % against 4.7 %. At \(T = 1000\) the GPD says 8.9 % while the empirical 7.8 % rests on the two worst days in the sample; at \(T = 2500\) (ten years) only the GPD answers — 13.2 %, and Black Monday’s 25.6 % is the reminder that even that is one draw from the tail. Extrapolation with a fitted shape is what turns six years of data into a ten-year estimate.
Invert the fitted GEV at probability \(\alpha\): the loss that only a fraction \(\alpha\) of 20-day blocks will exceed.
\[\text{EVaR}_\alpha = -\Big(\mu + \frac{\sigma}{\xi}\big[(-\ln(1-\alpha))^{-\xi} - 1\big]\Big)\]
With \(\xi = 0.354\), \(\mu = 0.0136\), \(\sigma = 0.0070\), what is EVaR at \(\alpha = 5\,\%\) (as a return)?
−0.0506
A one-day 5 % VaR says “−1.6 %”. The extreme VaR says: in one 20-day block out of twenty, expect a day worse than −5.1 % — and 4.1 % of the 74 blocks did. That is the number a risk limit should be sized to.
Compute the 1 % extreme VaR from the fitted GEV into evar01 (use EVaR), then count how many of the 74 block minima fell below it.
c = −ξ.Three prompts that make an LLM useful for this chapter — and the pitfall each one guards against.
n, the standard deviation, and the excess kurtosis; if kurtosis is above 3, tell me which of the tests you propose assume Normality and give me a bootstrap or permutation alternative.” A copilot will run ttest_1samp on 20 fat-tailed returns and report a p-value to four decimals. The number is real; the assumption behind it is not.kstest against norm, confirm the data are standardised first (or pass args=(mean, sd)), and when you fit genextreme, state the sign convention and report ξ = −c.” Both are places where the code runs, returns a number, and the number is wrong — D = 0.47 on raw returns is a scale mismatch, not a shape test.Warning
In the second week of August 2007, quantitative equity funds run by Goldman Sachs, Renaissance, AQR and others lost between 10 % and 30 % in a few days as crowded long-short factor positions unwound together. Goldman’s Global Equity Opportunities fund fell about 30 % in a week and received a US$3 billion injection. Explaining it, CFO David Viniar told the Financial Times (13 August 2007): “We were seeing things that were 25-standard-deviation moves, several days in a row.”
Under a Normal model a 25σ event has probability of order \(10^{-137}\) — it should not happen once in the life of the universe, let alone on consecutive days. The models were not unlucky; they were mis-specified. A Fréchet tail with \(\xi \approx 0.35\), the number you fitted to the Dow, assigns such moves probabilities measured in years, not eons.
Lesson for this chapter: counting sigmas presumes a Gaussian. The tail index, block maxima and the extreme VaR exist precisely because “how many σ” is the wrong question about the tail.
To: Chief Risk Officer From: <Your name>, risk analytics Subject: Replace the daily-VaR-based stop with an extreme-VaR limit Date: 2026-09-15
Recommendation: Size the index book’s hard loss limit to the 5 % extreme VaR of the 20-day block minimum: −5.1 % of notional per block, reviewed quarterly.
Evidence: - GEV fitted by MLE to 74 block minima (DJI 1985–1990): ξ = 0.35, μ = 0.0136, σ = 0.0070; Q-Q diagnostic straight apart from 19 Oct 1987. - EVaR(5 %) = −5.1 % vs the daily 5 % quantile of −1.6 %: the current stop is sized to an ordinary bad day, not to an ordinary bad month. - Six blocks in six years scored below 0.10; the limit would have been breached 3 times (4.1 % of blocks), all in identified stress episodes.
Caveats: - 74 observations: the ξ confidence interval is wide (bootstrap it, §3.2, before sign-off). - Historical window ends 1990; re-estimate on 2008 and 2020 to test stability of ξ. - GPD on 91 exceedances gives ξ = 0.44 and the Hill plateau 0.42–0.47 vs the GEV’s 0.35; agreement in sign and size is the sanity check, not proof.
Next step: Re-fit on rolling 6-year windows to 2024; report EVaR(5 %) and EVaR(1 %) with bootstrap bands.
| Concept | Tool |
|---|---|
| Population vs sample | sample, std(ddof=1) |
| Empirical distribution | ECDF via np.sort and arange(1, n+1)/n, DKW bound, KernelDensity, Silverman’s \(h\) |
| Parametric fit and diagnostics | stats.t.fit, probplot, bimodality coefficient |
| Bootstrap | rng.integers(0, n, (B, n)), percentile CI, BCa (\(z_0\), \(a\) from the jackknife) |
| Tests | shapiro, kstest (studentise first), ks_2samp, ttest_1samp(alternative=), permutation via rng.shuffle, rolling \(\hat t\) |
| Experimental design | \(n = 2\sigma^2(z_\alpha+z_\beta)^2/\tau^2\), peeking simulation, \(\alpha\)-spending, ε-greedy / UCB1 / Thompson |
| Association | pearsonr, spearmanr, kendalltau, distance correlation via pdist, tail co-exceedance |
| Extreme values | genextreme.fit (ξ = −c), block minima, genpareto.fit, Hill estimator, return levels, extreme score, EVaR |
Next: Chapter 4 — Statistical Predictive Models.
Prof. Xuhu Wan · HKUST ISOM · Learning Statistics with Python