Chapter 3: Reshaping Statistics

Chapter Introduction

Chapter 2 left you with clean frames: date-typed, audited for missing and duplicated values, reshaped by groupby, merge and pivot_table. This chapter turns those frames into statistics — and the first thing it reshapes is not a table but a habit of mind. Most introductory courses train you to compute a mean and a standard deviation, run a $t$-test and report a number. That training is adequate for a tame world. It is dangerous in a world where the extremes do most of the work, where one observation in a thousand drives the year’s profit and loss, and where the next observation may not resemble any of the previous ones. In August 2007 the chief financial officer of Goldman Sachs explained a week of losses by saying that his funds had seen “25-standard-deviation moves, several days in a row”. Under the Normal model his risk systems assumed, one 25σ event has a probability of order $10^{-137}$. The models were not unlucky; they were wrong about the shape of the distribution, and everything downstream — the intervals, the tests, the correlations, the loss limits — inherited the error.

“Reshaping statistics” therefore means moving from averages to distributions. Every number you compute is a draw from a distribution; every comparison is a comparison of distributions; every “relationship” between two variables is a property of their joint distribution; and every honest statement carries uncertainty. A population parameter is a fact. A sample statistic is an estimate with a distribution of its own, and the five sections are five ways of taking that seriously. §3.1 builds the empirical distribution from five numbers and then from 4 309 days of Apple returns, smooths it into a density with a bandwidth and with a formula, learns to read a Student-$t$ off a Q-Q plot, and then lets the Shapiro–Wilk and Kolmogorov–Smirnov tests put numbers on what the eye saw. §3.2 asks how sure you can be of any of those numbers: the bootstrap gives a standard error and a confidence interval for any statistic — even the kurtosis, which has no textbook formula worth trusting — and the classical tests of a mean — one-sample, permutation, two-sample — are recast as decisions, one of which becomes a rolling regime detector. §3.3 reverses the direction: you design the experiment — how many users, when to stop, and when a bandit that earns while it learns beats a test altogether. §3.4 takes two variables at once and ranks Pearson, Spearman, Kendall and distance correlation by the shape of dependence each can see, then asks whether the association between two stocks survives on the worst days. §3.5 goes to the tail itself: extreme value theory on the Dow Jones around Black Monday, a fitted GEV, a peaks-over-threshold GPD, the Hill estimator, return levels, and the extreme VaR — the number that names your highest pain threshold.

The datasets are the ones the slides use, and the printed numbers are the same. Apple’s daily closes from 2005 to 2022 (apple_2005_2022.csv) carry §3.1 and §3.2; a panel of world index moves (indicepanel.csv) and ten years of S&P 500 and Tesla returns (returns.csv) carry §3.4; the Dow from 1985 to 1990 (dji.csv) carries §3.5, chosen so that 19 October 1987 is in the sample. §3.3 is simulation by necessity — you cannot rerun a market or an A/B test on demand — with fixed seeds. Wherever a cell prints something, the prose says what the number is and what it means; wherever a number could have been predicted before the cell ran, the text asks you to predict it. That habit — decide what the distribution should look like, compute it, and treat any surprise as a claim to be checked — is the chapter’s real subject.

By the end you will be able to describe a variable by its empirical CDF, estimate its density two ways, bootstrap any statistic, run and read the standard tests, size an experiment, compute four measures of association, and fit the limit laws of the extremes into a loss limit a risk committee can argue with. Chapter 4 puts these tools inside predictive models, Chapter 5 replaces the frequentist interval with a posterior, Chapter 6 adds time, and Chapter 7 returns to the multiplicity problem of §3.3 when it puts forty-two trading rules on trial.

Population vs Sample and Distributions in Practice
Bootstrap Confidence Intervals and Hypothesis Testing
Experimental Design: A/B Tests and Bandits
Association: Linear, Monotonic, and Nonlinear
Extreme Value Theory: GEV, Block Maxima, POT, Extreme VaR

Population vs Sample and Distributions in Practice

A distribution is the most honest summary of a variable you can produce. A mean throws away information; a standard deviation throws away more; a histogram throws away very little; the empirical CDF throws away nothing except the order in which the observations arrived. The first habit of a serious analyst is to look at the shape of every variable before doing anything else with it, because every method that follows — the $t$-test, the regression of Chapter 4, the Gaussian noise inside a Kalman filter in Chapter 6 — is a statement about a shape, and a wrong shape corrupts every number downstream. This section builds the three ways of looking: the empirical distribution (histogram and ECDF), the smoothed empirical distribution (a kernel density estimate), and the theoretical distribution (a parametric family fitted by maximum likelihood). It ends with the diagnostics that tell them apart.

Population vs sample: why pandas divides by n − 1

The lecture’s sixty-three numbers are a whole population; nothing is hidden. Draw ten of them without replacement and you have a sample. The population mean and variance are descriptive facts, computed with ddof=0 — divide by $N$. A sample variance that divides by $n$ underestimates the population variance on average, because the sample mean sits closer to the sample points than the true mean does; dividing by $n - 1$ (ddof=1, pandas’ default for std and var) removes that bias. NumPy’s default is ddof=0; the difference between the two libraries is the source of a thousand off-by-a-little discrepancies between a notebook and a textbook.

The population of 63 has mean 47.746 and standard deviation 28.505. The sample of 10 drawn with seed 1 has mean 49.7 and standard deviation 23.372 with ddof=1; with ddof=0 the same ten numbers give 22.172, five per cent smaller, and the smaller number is the biased one. Change the seed and both sample numbers move — that movement is the sampling distribution, and every confidence interval and hypothesis test from here on is a statement about how far it can move. One more method matters here: sample(frac=0.8, replace=True) draws with replacement, so a value can appear twice and another not at all. That is the primitive of the bootstrap: resample your own data as if it were the population, recompute the statistic, and read its spread off the replicates. When no closed-form standard error exists — for a Sharpe ratio, a quantile, a maximum drawdown — resampling is the universal solvent, and §3.2 builds it from this one call.

Univariate exploration: Apple, 2005–2022

The variable for the rest of the section is Apple’s daily log return, $r_t = \ln S_t - \ln S_{t-1}$, over 4 309 trading days from January 2005 to February 2022. Log returns add across time and are symmetric under sign reversal, which makes them the natural object for distribution fitting; they differ from simple returns by a term of order $r^2$, negligible on ordinary days and visible only on the crashes that §3.5 is about. agg with a list of four names gives the shape in one table. Before running, decide: is the excess kurtosis of daily returns positive or negative?

The mean daily return is 0.12 % with a standard deviation of 2.08 %. The skew is −0.28 (a slightly longer left tail) and the excess kurtosis is 5.99. Pandas reports Fisher kurtosis, so a Normal scores 0; a value near 6 means far more mass in the tails than a Normal with the same variance would place there. The price column tells a different story — skew 1.98 and kurtosis 3.36 — but the shape of a price series is driven by growth and means little; it is the return whose distribution matters. The two red lines mark the quartiles: half of all days lie in a band only 2.03 % wide, yet the axis runs from −19.8 % to +13.0 %. That gap between the crowded middle and the sparse extremes is what the rest of the section, and all of §3.5, quantifies. Notice what the table does not tell you: skew and kurtosis are ratios of sample moments, and the fourth moment of a fat-tailed sample is dominated by a handful of days — §3.2 will bootstrap the kurtosis of one year and find a standard error a third of its own size. A shape summary is a starting point, not a conclusion.

From theoretical to empirical: build an ECDF from five numbers

Textbooks hand you a density. Real work hands you a sample. The bridge between them is one line of arithmetic, worth building by hand before any library does it for you. Take the sample [3, 1, 4, 1, 5] and define $\hat F(x)$ as the fraction of the sample that is at most $x$. What fraction is $\le 3$? Three of five, 0.6. What fraction is $\le 1$? Two of five, 0.4. You have just defined the empirical cumulative distribution function. Given an i.i.d. sample $X_1, \dots, X_n$,

\[\hat F_n(x) = \frac{1}{n}\sum_{i=1}^{n}\mathbf 1\{X_i \le x\},\]

a step function that rises by $1/n$ at every observation. It assumes nothing about shape, it is the non-parametric maximum-likelihood estimate of the true CDF, and a sorted sample’s $i$-th value sits at height $i/n$. The second half of the cell applies that to Apple: studentise the 4 309 returns — subtract the mean, divide by the standard deviation — sort them, and read the staircase at three heights. Predict first: the ECDF height at the 2 155th of 4 309 sorted points is $2155/4309 = 0.500$, and its $x$-coordinate is the median.

The five-number ECDF prints 0.4, 0.6 and 1.0 at $x = 1, 3, 5$. On Apple, the median studentised return is −0.008, essentially zero. The 10 % and 90 % points sit at −1.110 and +1.086, inside the Normal’s ±1.282: the body is narrower than a Normal with the same variance. And yet 76.9 % of days lie within one standard deviation, against the Normal’s 68.3 %, while only 98.51 % lie within three, against 99.73 %. Both departures point the same way. A taller, narrower peak and fatter tails is what a variance-matched fat-tailed distribution must look like — the variance is fixed, so the extra mass in the tails is paid for by mass pulled in from the shoulders. That is excess kurtosis 5.99, read off a staircase rather than a fourth moment, and the staircase is the more robust reading: it cannot be dragged around by three days in 2008.

Why doesn’t a histogram do the job? The ECDF of a bimodal mixture

A histogram’s shape depends on the bin edges and the bin width; a wrong width can merge two clusters into one hump or split one hump into a comb. The ECDF has no tuning knob, so it cannot lie about shape. The lecture’s demonstration is a two-component mixture (seed 7): with probability one half a draw comes from a wide Normal centred at +5, otherwise from a narrow Normal centred at −10. Predict what its ECDF looks like — a single smooth S, or something else?

The curve climbs steeply through −10, flattens between about −5 and 0, then climbs again around +5. The quartiles say the same thing: the first quartile is −10.00, the median −5.18, the third quartile 5.03 — a median that sits in the empty valley between the two clusters, the surest sign that “the middle” of this variable is not a place where any observation lives. The flat stretch is the bimodality. A histogram with bins 10 wide would have hidden it; the ECDF shows it with no choice made at all.

How fast does it converge? Glivenko–Cantelli and the DKW bound

The ECDF is built from $n$ random points, so it is itself random. Does it land on the truth as $n$ grows? The Glivenko–Cantelli theorem (1933) says yes, and uniformly: $\sup_x |\hat F_n(x) - F(x)| \to 0$ almost surely — the worst gap anywhere on the real line goes to zero. The Dvoretzky–Kiefer–Wolfowitz inequality is the finite-$n$ version,

\[P\Big(\sup_x |\hat F_n(x) - F(x)| > \varepsilon\Big) \le 2e^{-2n\varepsilon^2},\]

and it is remarkable for what it does not depend on: the true $F$. Whatever Apple’s return distribution is, the whole staircase is within $\varepsilon$ of it with a probability you can compute from $n$ and $\varepsilon$ alone. Predict the bound at $\varepsilon = 0.02$ for $n = 4309$: $2e^{-2 \cdot 4309 \cdot 0.0004} = 2e^{-3.45}$.

Evaluate $2e^{-2n\varepsilon^2}$ at $n = 4309$, $\varepsilon = 0.02$ to three decimals — and decide whether the bound is even informative at $n = 100$.

At $n = 4309$ the bound is 0.064: with probability at least 93.6 %, no point of the Apple staircase is more than 0.02 from the true CDF. Inverting the inequality gives a uniform confidence band with half-width 0.0207 at 95 %: every value you read off the staircase is within about two percentage points of the truth, simultaneously for all $x$. At $n = 100$ the bound is 1.846 — larger than one, hence vacuous — and at $n = 1000$ it is 0.899, barely informative; the useful regime starts in the low thousands. This theorem is the foundation of the chapter. The bootstrap of §3.2 works because $\hat F_n$ stands in for $F$; the Kolmogorov–Smirnov test measures exactly the $\sup$ gap that DKW bounds; the kernel density estimate below is a smoothed derivative of $\hat F_n$. When any of them misbehaves, the first question is whether $n$ is in the regime where the staircase can be trusted.

Estimating a density without a formula: kernel density estimation

A staircase has no density — its derivative is a row of spikes at the observations. To get a smooth density, smear each spike into a small bump of width $h$ and add the bumps up:

\[\hat f_h(x) = \frac{1}{nh}\sum_{i=1}^{n} K\!\left(\frac{x - x_i}{h}\right),\]

with $K$ a Gaussian kernel. The bandwidth $h$ is the only knob, and it is the bias–variance trade-off in its purest form. A small $h$ makes narrow bumps; the estimate follows every sampling wiggle and looks spiky (high variance). A large $h$ merges the bumps; real features such as two modes are smeared into one (high bias). Silverman’s rule of thumb gives a starting value, $h^\star \approx 1.06\,\hat\sigma\,n^{-1/5}$, derived by minimising the integrated squared error when the truth is a single Normal. Predict it for the studentised Apple returns ($\hat\sigma = 1$, $n = 4309$): $1.06 \times 4309^{-0.2} \approx 0.199$.

With $h = 0.2$ the curve is a comb of spikes chasing individual observations. With $h = 3$ the two modes are still visible but the narrow left component has been widened to roughly the bandwidth — the estimate says “a bump near −10” but no longer says how sharp. $h = 1$ is a reasonable compromise. Silverman’s rule gives $h^\star = 0.199$ on the studentised Apple scale — a fifth of a standard deviation — or 0.00413 on raw returns; for the mixture it gives 1.97, close to the curve that over-smooths, because the rule assumes one hump and the mixture’s overall standard deviation is dominated by the distance between the humps. That failure is the lesson: Silverman is a starting point for a sweep, never a substitute for one. Two mechanical notes: scikit-learn’s KernelDensity expects a two-dimensional array, hence mix[:, None]; and score_samples returns the log-density, so you exponentiate before plotting.

Estimating a density with a formula: fit a Student-t

The alternative to a bandwidth is a family. Studentise the returns and let scipy.stats.t.fit choose the degrees of freedom $\nu$, location and scale by maximum likelihood. The Student-$t$ is the right family to try first because its single shape parameter controls the tails: $\nu \to \infty$ recovers the Normal, and small $\nu$ gives tails that decay like a power law, $x^{-(\nu+1)}$. Moments exist only up to order $\nu$: the variance is infinite for $\nu \le 2$ and the kurtosis for $\nu \le 4$. Predict what a fitted $\nu$ near 3 would mean — and whether one year of a single volatile stock is enough to see it.

The fitted degrees of freedom are 3.29, with location 0.008 and scale 0.663. Read the scale first: a $t_{3.3}$ with unit variance needs a scale well below one, because so much of its variance lives in the tails. Then read $\nu$: at 3.3 the fourth moment barely exists, which is why the sample kurtosis of returns is so unstable — add one crash day and it jumps. The dashed Normal is too low at the centre and too low in the tails, compensating by being too high at ±1.5; the $t$ hugs the histogram everywhere. The log-likelihood comparison is the rigorous version of “the $t$ looks better”: −5 731.9 against −6 113.7, a gap of 382 units for one extra parameter, when a gap of 2 or 3 would already justify it under AIC. Tesla in 2020 — 253 trading days that included a five-for-one split and a 700 % rise — fits with $\nu = 3.31$, almost identical to Apple’s seventeen-year value, and Shapiro–Wilk rejects normality with $p = 1.8 \times 10^{-6}$ on 253 observations. One year is enough to see the tails. The lecture repeats the fit for thirteen large stocks and finds $\nu$ between 3 and 6 for every one: heavy tails are the normal condition of daily equity returns.

In Colab

import yfinance as yf
from scipy import stats
for s in ["AAPL", "TSLA", "NVDA", "AMD", "QCOM", "AMZN", "GOOG", "BABA", "JD", "CVX", "AMN", "MRO", "BAC"]:
    px = yf.download(s, start="2020-11-18", end="2024-11-18", progress=False)["Close"].pct_change().dropna()
    zz = (px - px.mean()) / px.std(ddof=1)
    print(s, "df =", round(stats.t.fit(zz)[0], 2), " Shapiro p =", stats.shapiro(zz).pvalue)

Four years of thirteen tickers: df from 3.13 (BABA) to 5.60 (NVDA), every Shapiro–Wilk $p$ below $10^{-8}$; none is anywhere near the Normal’s infinity.

Testing for normality: look first with a Q-Q plot

A Q-Q plot puts the sample’s quantiles against a reference distribution’s quantiles. If the two have the same shape the points lie on a straight line; a heavier-tailed sample has more extreme values than the reference, so the largest quantiles overshoot the line at the top right and the smallest undershoot at the bottom left — the S-shape every finance student learns to recognise. stats.probplot takes the distribution name and, for the $t$, its shape parameter via sparams. Predict the two panels: against a Normal, and against the fitted $t_{3.3}$.

Against the Normal, the tails peel away from the line on both sides — the classic S. Against the fitted $t_{3.3}$ the points hug the line until the very last few, where the right tail is shorter than the $t$ predicts — the lecture’s marginal note, “right-side short tail compared with the $t$ benchmark”. A symmetric $t$ cannot capture the asymmetry that the skew of −0.28 already hinted at; a skewed-$t$ would, at the cost of one more parameter. The Q-Q plot is the eye’s kurtosis test, and it has one advantage over any formal test: it shows where the departure is. The formal tests come next, and they show that 4 309 observations reject the Normal by a margin of $10^{-38}$ and cannot reject the $t$ at all.

Testing for normality: Shapiro–Wilk and Kolmogorov–Smirnov

The Q-Q plot showed where the departure is; a hypothesis test puts a number on whether it could be chance. The logic is always three steps. State a null hypothesis $H_0$ that represents “nothing interesting” — the population is Normal, two groups have the same distribution, the mean is zero. Compute a test statistic whose distribution under $H_0$ is known. Report the p-value, the probability under $H_0$ of a statistic at least as extreme as the one observed. Then say all three before drawing any conclusion. For normality there are two standard tests. The Shapiro–Wilk statistic $W$ measures how well an ordered sample matches the expected order statistics of a Normal; it is the most powerful normality test for moderate samples. The Kolmogorov–Smirnov statistic $D$ is the largest vertical gap between the empirical CDF and a reference CDF — exactly the $\sup$ that DKW bounds — and its reference can be any distribution, or a second sample. KS is the more general tool, and it has a trap the lecture notebook fell into. kstest(returns, "norm") compares the data with the standard Normal. Raw daily returns have a standard deviation of about 0.02, so nearly all of them sit within ±0.05 — a huge distance from $N(0,1)$’s CDF — and the statistic comes out at $D \approx 0.47$, which says nothing about shape. The cell below runs, returns a number, and the number is wrong.

This cell tests whether Apple’s returns are Normal. It prints a KS statistic of 0.47 and $p = 0$ — decisive, apparently. What is wrong with the conclusion?

The bug is a scale mismatch, not a shape finding. "norm" with no arguments is $N(0, 1)$; the returns are roughly $N(0.001, 0.02)$. The 0.47 is the gap between the CDF of a variable that is almost always in $[-0.05, 0.05]$ and the CDF of one that is almost never there — it would be just as large for perfectly Normal returns with a standard deviation of 0.02. The fix is to compare shape: studentise first and test z against "norm", or pass the fitted mean and standard deviation via args=(mean, sd).

Shapiro–Wilk gives $W = 0.9415$ and $p \approx 3 \times 10^{-38}$: the null of a Normal population is rejected as decisively as a test can reject anything. KS on the studentised returns gives $D = 0.075$, $p \approx 2 \times 10^{-21}$ — also a rejection, and now an honest one, with a statistic six times smaller than the scale-corrupted 0.468. Against the fitted $t_{3.3}$ the KS test does not reject: $D = 0.015$, $p = 0.274$. That is what “the $t$ fits” means in numbers — the largest gap between the empirical CDF and the fitted $t$ CDF is 1.5 % of probability, inside the 2.1 % DKW band computed above. Only when all three tests are stated does “Apple’s returns are not Normal” become “not Normal with $D = 0.075$, but indistinguishable from a $t_{3.3}$ with $D = 0.015$”.

What a p-value is and is not

The p-value is $P(\text{data at least this extreme} \mid H_0)$, not $P(H_0 \mid \text{data})$ — Chapter 5 computes the second quantity, and it requires a prior. Run twenty tests on noise and expect one to pass at 5 %, the multiplicity problem §3.3 returns to. And a non-rejection is not a proof of the null: the $t$ “fits” in the sense that 4 309 observations cannot tell it apart from the data, not in the sense that Apple’s returns are $t$-distributed.

A cheap bimodality detector

Before committing to a density, it helps to know which problem you have: one hump with fat tails, or two humps. The bimodality coefficient $b = (g_1^2 + 1)/(g_2 + 3)$, built from the sample skewness $g_1$ and excess kurtosis $g_2$, flags trouble when $b > 5/9 \approx 0.555$, the value a uniform distribution attains. The intuition is in the denominator: two well-separated humps produce negative excess kurtosis — flatter than a Normal, with the mass pushed out to two shoulders — and a small denominator makes $b$ large. Predict the verdict on the mixture from its moments, $g_1 = 0.381$ and $g_2 = -1.214$: $(0.145 + 1)/1.786 = 0.641$, above the threshold. Then predict Apple.

Mixture: $b = 0.641$, flagged — the negative kurtosis of −1.21 is the fingerprint of two separated humps. Apple: $b = 0.120$, not flagged — one peak, fat tails, and the large positive kurtosis drives $b$ down. A flagged variable wants a mixture model or a small-bandwidth KDE; an unflagged fat-tailed one wants a Student-$t$ or, for the tail alone, the extreme value distributions of §3.5. It is not a test — no p-value, and it is fooled by skewed unimodal shapes — but as a five-line screen over a hundred columns it earns its place.

The sample variance that divides by $n$ underestimates the population variance on average, because the sample mean sits closer to the sample points than the true mean does; dividing by $n - 1$ (ddof=1) removes the bias. NumPy’s default is ddof=0. On the ten-number sample the two give 23.372 and 22.172.

The probability that the largest gap anywhere between the ECDF and the true CDF exceeds $\varepsilon$: $P(\sup_x|\hat F_n - F| > \varepsilon) \le 2e^{-2n\varepsilon^2}$, independent of $F$. At $n = 4309$, $\varepsilon = 0.02$ it is 0.064, so the whole Apple staircase is within 0.02 of the truth with probability at least 93.6 %. At $n = 100$ the bound exceeds 1 and says nothing.

Tails so heavy that moments above order 3.3 do not exist — the kurtosis is not finite, which is why the sample kurtosis is so unstable. The ECDF shows it as 76.9 % of days within ±1 sd (Normal: 68.3 %) but only 98.5 % within ±3 (Normal: 99.7 %); the Q-Q plot against a Normal shows it as an S-curve. The lecture finds df between 3 and 6 for every large stock.

No — it is a scale mismatch. 'norm' is $N(0,1)$ and raw returns have sd ≈ 0.02, so almost all of them sit near zero where $N(0,1)$ has little mass. Studentise first ($D = 0.075$, still rejected at $p \approx 10^{-21}$) or pass args=(mean, sd). Against the fitted $t_{3.3}$ the test does not reject ($D = 0.015$, $p = 0.27$).

Bootstrap Confidence Intervals and Hypothesis Testing

§3.1 produced numbers: a mean of 0.12 %, a kurtosis of 5.99, a fitted $\nu$ of 3.29. Each is a statistic computed from one sample, and each would come out differently on another. The textbook answer to “how different?” is a standard error, and for the mean there is a formula, $s/\sqrt n$. For the median the formula depends on the unknown density at the median; for the kurtosis it depends on the eighth moment, which for a $t_{3.3}$ does not exist; for a Sharpe ratio, a quantile or a maximum drawdown the formulas are unknown or wrong. The bootstrap, introduced by Bradley Efron in 1979, gives a standard error and a confidence interval for any statistic by resampling the data itself, and the first half of this section reinvents it, checks it against the one case with a known answer, and uses it where no formula exists. The second half turns to hypothesis tests — structured decisions about whether what you see could have been produced by chance, the null / statistic / p-value logic §3.1 applied to the shape of a distribution — and runs four of them on Apple’s mean: was a year’s mean return zero, was one year different from the next, do two stocks share a distribution, and, recomputed every day, has the distribution changed.

The same sample, reloaded

Every section runs standalone. Reload Apple 2005–2022, studentise, refit the $t$, and pull out 2017 — the year the one-sample test will examine.

Apple in 2017 had 251 trading days, a mean daily log return of 0.151 % and a standard deviation of 1.109 % — a calm, rising year, exactly the kind in which “was the mean really positive?” is a fair question.

The bootstrap idea

You have a sample and a statistic $\hat\theta = T(X_1, \dots, X_n)$, and you want to know how much $\hat\theta$ would wobble across new samples from the same population — but you cannot draw new samples. The bootstrap’s move is the one §3.1 prepared: the empirical CDF $\hat F_n$ is your best estimate of the unknown $F$, uniformly close to it by Glivenko–Cantelli, so drawing a fresh sample from the population is replaced by drawing a sample of size $n$ from $\hat F_n$. Sampling from $\hat F_n$ means picking each observation with probability $1/n$, independently, $n$ times — which is resampling the data with replacement. Some values appear twice, some not at all, and that is the whole mechanism. The recipe: draw $n$ observations with replacement, compute $\hat\theta^*$ on the resample, repeat $B$ times; the standard deviation of the $B$ replicates is the bootstrap standard error, and their 2.5 % and 97.5 % quantiles are the percentile confidence interval. No formula, and no assumption about the shape of $F$ beyond the one Glivenko–Cantelli licenses.

Reinvent the bootstrap SE, then check it against a case you know

The code is more compact than the recipe. Build a $B \times n$ matrix of resample indices with rng.integers(0, n, size=(B, n)), index the data with it to get $B$ resamples at once, and apply the statistic along axis 1. No loop. The first test case is a small skewed sample — 120 lognormal draws, seed 11 — and a statistic with no usable textbook standard error: the median. Predict the order of magnitude of its bootstrap SE: 0.02, 0.2 or 2? The second test case is the one statistic that does have a formula. For the mean of Apple’s 251 returns of 2017, $\mathrm{SE} = s/\sqrt n = 0.01109/\sqrt{251} \approx 0.0007$; if the bootstrap is legitimate, it must match.

The sample median is 1.659 and its bootstrap standard error is 0.198 — a number no formula gave you, measured in a millisecond. The bootstrap bias, $\bar\theta^* - \hat\theta = +0.053$, says the median of a resample tends to sit above the sample median: the lognormal is right-skewed, and resampling a skewed sample produces a skewed distribution of medians. Keep that bias in mind; it is what the BCa interval below corrects. The second block is the sanity check: the bootstrap SE of the 2017 mean is 0.0007, the formula gives 0.0007, and the ratio is 1.000 to three decimals. The bootstrap is not magic. It reproduces the answer you already trusted and then extends to medians, ratios, quantiles and kurtosis, where there is nothing to trust it against.

Percentile interval for Apple’s 2017 mean return

The simplest confidence interval is the percentile method: the 2.5 % and 97.5 % quantiles of the $B$ bootstrap means. Predict whether the interval for the 2017 mean contains zero: the mean is 0.00151 with SE 0.0007 — about 2.2 standard errors above zero — so a 95 % interval should just exclude it.

$[0.00014,\ 0.00294]$: zero is outside, barely. The textbook interval $\bar x \pm 1.96\, s/\sqrt n$ is $[0.00014, 0.00288]$ — identical at the lower end, a hair narrower at the upper end because the bootstrap distribution of a mean of slightly skewed returns is itself slightly skewed. Remember this interval: the one-sample $t$-test below reaches the same verdict with a two-sided p-value of about 0.03, because the interval and the test are the same statement — “$\bar x$ is 2.16 standard errors from zero” — read from two directions. The percentile interval is correct when the bootstrap distribution is symmetric and unbiased. For a mean of 251 observations it is. Predict when that assumption is dangerous.

When the percentile interval goes wrong: BCa for the kurtosis of 2017

For skewed statistics — the variance, the kurtosis, a ratio, an extreme quantile — or at small $n$, the bootstrap distribution is neither symmetric nor centred on $\hat\theta$, and the percentile interval is shifted in the wrong direction. Efron’s BCa interval (bias-corrected and accelerated, 1987) fixes two blind spots at the same computational cost. The bias correction $z_0 = \Phi^{-1}(\text{share of } \hat\theta^* < \hat\theta)$ measures how far the bootstrap distribution’s median sits from $\hat\theta$; a share of 50 % gives $z_0 = 0$. The acceleration $a$ measures how fast the standard error changes with the true value — the skewness of the estimator — and is estimated from the jackknife: leave one observation out at a time, recompute the statistic, and take a standardised third moment of the leave-one-out values. The two constants adjust the quantile levels at which the interval is read:

\[\alpha_{\text{lo, hi}} = \Phi\!\left(z_0 + \frac{z_0 + z_{\alpha}}{1 - a\,(z_0 + z_{\alpha})}\right), \qquad z_\alpha = \Phi^{-1}(0.025),\ \Phi^{-1}(0.975).\]

The statistic is the excess kurtosis of 2017 — the one §3.1 warned you about. Predict the sign of $z_0$ first: resampling with replacement drops some of the extreme days from most resamples, so most bootstrap kurtoses fall below the sample value, the share below exceeds 50 %, and $z_0$ is positive.

The excess kurtosis of 2017 is 4.58 with a bootstrap standard error of 1.55 — a third of its own size. Any risk model that takes “the kurtosis” as a fixed input is taking a number that could plausibly be 2 or 8. As predicted, 59.4 % of resamples fall below the sample value, so $z_0 = +0.237$; the acceleration is $a = +0.121$, large by the standards of the statistic. Together they move the quantile levels from (0.025, 0.975) to (0.117, 0.999), and the interval from $[1.47, 7.31]$ up to $[2.23, 9.15]$ — right by about 0.8 at the bottom and 1.8 at the top, in the direction the theory predicts, because the resamples systematically lost the fat-tail days. The rule is general. For a mean, the two intervals agree. For medians, quantiles, ratios and kurtosis the bias and skew are real, and BCa moves the endpoints to where they belong at the cost of $n$ jackknife evaluations.

When the bootstrap fails

When the statistic depends on the order of the data — a serial correlation, a GARCH parameter — resampling destroys the order; Chapter 6 uses a block bootstrap of contiguous stretches. For a sample maximum or minimum, $\hat F_n$ cannot extrapolate beyond the largest value seen — which is why §3.5 fits a parametric tail rather than bootstrapping the worst day. And at very small $n$, $\hat F_n$ is a poor stand-in for $F$.

One-sample t-test: was Apple’s 2017 mean return zero?

A confidence interval says how much a number wobbles; a test asks whether a specific claim survives. The one-sample $t$-test asks whether a sample mean is consistent with a hypothesised value. $H_0: \mu = 0$; the statistic is $\hat t = \bar x / (s/\sqrt n)$, the sample mean in units of its own standard error. Its null distribution is Student’s $t$ with $n - 1$ degrees of freedom if the data are Normal, and approximately standard Normal for large $n$ whatever the data are, by the central limit theorem — which is why the test works on returns that are nothing like Normal, provided $n$ is in the hundreds. Predict $\hat t$ to two decimals, and decide whether it clears the one-sided 5 % critical value of 1.645.

$\bar x = 0.00151$, $s = 0.0111$, $n = 251$. Compute $\hat t = \bar x / (s/\sqrt{n})$ in your head to two decimals — then decide how the one-sided p-value relates to the two-sided one.

$\hat t = 2.16$. The one-sided p-value against the Normal approximation is 0.0154, and the critical value is 1.645, so a one-sided test rejects $\mu = 0$: Apple’s mean daily return in 2017 was positive at a level a chance year would produce about one time in sixty-five. scipy.stats.ttest_1samp reports the same statistic with a two-sided p-value of 0.0318 on 250 degrees of freedom, and with alternative="greater" a one-sided p-value of 0.0159 — half the two-sided value, because a one-sided test counts only the tail in the hypothesised direction. Both reject at 5 %, and the bootstrap interval excluded zero for the same reason: the interval, the statistic and the p-value are three views of one fact. Sidedness is a design decision made before seeing the data. The business question “did Apple earn a positive return?” is directional, and a one-sided test is legitimate if chosen in advance; choosing it after seeing that $\hat t$ is positive is p-hacking, and two-sided is the honest default. Notice finally what the test does not say: it says nothing about 2018. A year in which the mean return was distinguishable from zero is a description, not a forecast.

Permutation test: was 2017 a different year from 2018?

The $t$-test assumed a Normal sampling distribution for $\bar x$, which the central limit theorem delivers at $n = 251$ but would not at $n = 20$. For the sharp null “the two groups have identical distributions” you need no such assumption, because under $H_0: F_A = F_B$ the group labels carry no information: every assignment of labels to the pooled values is equally likely — the observations are exchangeable. So build the null distribution by shuffling the labels. Pool the two samples, shuffle, split at the original group size, recompute the difference in means, repeat $B$ times; the permutation p-value is the share of shuffles whose difference is at least as large in absolute value as the one observed. No Normality, no equal-variance assumption, no large-$n$ asymptotics. Apple’s mean daily return was +0.151 % in 2017 and −0.028 % in 2018. Predict the two-sided p-value: with 251 days a side and a daily standard deviation of 1–2 %, the standard error of the difference is about 0.13 %, so a gap of 0.18 % is 1.3 standard errors and should be unremarkable.

The observed gap is 0.179 % a day; the standard deviation of the gap across 5 000 shuffles is 0.133 %, so the gap is 1.35 standard errors; the permutation p-value is 0.184, and Welch’s two-sample $t$-test — the version that does not assume equal variances — gives 0.183. The classical test was fine here because $n$ is large; the permutation test would still be right for twenty days of fat-tailed returns, where the $t$ is not, and its Monte-Carlo error shrinks like $1/\sqrt B$, independent of $n$. Two years whose mean returns look different — one up, one flat — cannot be told apart at any conventional level, because a year of daily returns is a noisy estimate of its own mean; this is why annual rankings of fund managers are so unstable, and why Chapter 5’s Bayesian treatment of the Sharpe ratio begins by admitting how wide the posterior is.

A rolling two-sample test as a regime detector

Here is the lecture’s most useful trick. Compare the last ten days’ returns with the ten days that ended twenty days earlier, using a two-sample $t$ statistic with the pooled standard error:

\[\hat t = \frac{\bar r_{\text{now}} - \bar r_{\text{lag}}}{\sqrt{\left(s^2_{\text{now}} + s^2_{\text{lag}}\right)/10}}.\]

Recompute it every day. Days with $|\hat t| > 2$ are days on which the recent return distribution differs from the one a month earlier — a regime switch in the mean, flagged by a test statistic. shift(20) lags the series; rolling(10).mean() and .var() do the rest without a loop.

Twenty of 249 days are flagged, and not uniformly: one in December 2012 during the slide from the September 2012 peak, two in late March 2013 on the way into the April trough, and seventeen in May, June and July 2013 as the price turned and climbed away from it — the points at which the distribution of recent returns had shifted relative to the month before. The flags sit around the turning points rather than on them, because a ten-day window needs several days of the new regime before its mean moves. The third line is the sobering one: over the full seventeen years, 5.2 % of days carry $|\hat t| > 2$ — almost exactly the 4.6 % a Normal null would produce. On the whole sample the detector fires about as often as chance says it should, and only the clustering of the flags — six in May and six in June — carries information. A test statistic recomputed every day is a trading signal, but it inherits the false-positive rate of the test (§3.3 shows what repeated looks do to that rate), and a window is a bandwidth: ten days is the $h$ of this estimator, trading responsiveness against false alarms exactly as the KDE did.

The KS version replaces the $t$ statistic with ks_2samp on twenty-day windows, so that a change in shape — volatility, tails — is detected even when the mean has not moved.

272 windows are tested and 33 have $p < 0.2$; the smallest p-value any twenty-versus-twenty comparison can reach is 0.034, which is why the lecture’s threshold is a loose 0.2 rather than 0.05. With twenty observations a side the KS test has little power — $D$ can only take values in multiples of $1/20$ — and that is the honest limit of the method: a regime detector built on short windows can see large shifts quickly or small shifts slowly, never both. Chapter 6 revisits regime detection with models that pool information across the whole history.

Because the empirical CDF $\hat F_n$ stands in for the unknown $F$, and sampling from $\hat F_n$ means picking each observation with probability $1/n$, independently, $n$ times — which is resampling with replacement. Glivenko–Cantelli (uniform convergence of $\hat F_n$ to $F$) is why the substitution is legitimate. It reproduced the mean’s known SE (0.0007) and then measured what no formula gives: 0.198 for a median, 1.55 for a kurtosis of 4.58 — and BCa moved the kurtosis interval from $[1.5, 7.3]$ to $[2.2, 9.1]$.

Under the sharp null $H_0: F_A = F_B$ the group labels carry no information — every relabelling of the pooled values is equally likely (exchangeability) — so shuffling labels and recomputing the statistic samples the null distribution directly, with no Normality, no equal-variance and no large-$n$ assumption. For Apple’s 2017 versus 2018 mean return the gap of 0.179 % a day was 1.35 shuffle-SEs: $p = 0.184$ from 5 000 shuffles, 0.183 from Welch’s $t$.

Experimental Design: A/B Tests and Bandits

So far the data arrived and you tested it. A hypothesis test interprets data you have already collected; experimental design asks the prior question — how do I plan the data collection so that the eventual test will be informative? This is where every product, marketing and clinical experiment lives, and it applies with one twist inside a quantitative fund: a backtest is an experiment whose treatment is a trading rule, and a fund that “peeks” at a strategy’s running performance and stops when it looks good is making exactly the error this section quantifies. Three parts: the classical fixed-$N$ design and why a real effect can read “not significant”; sequential testing and what checking the dashboard at lunch does to the false-positive rate; and multi-armed bandits, for when the goal is to earn the most reward while learning which arm is best. All three are simulations by necessity, with fixed seeds.

Classical A/B set-up: a real effect that reads “not significant”

Two arms, a fixed sample size $N$ per arm, and a hypothesis fixed before the data. Randomise units to control or treatment, wait for $N$ to accumulate, compute the difference in means $\hat\tau$ and its standard error, and reject $H_0: \tau = 0$ if $|z| = |\hat\tau / \widehat{\mathrm{SE}}| > 1.96$. The design is a century old — Fisher’s agricultural plots — and still the gold standard when you can afford it. Its single most-violated requirement is the one this cell demonstrates. Control outcomes are $N(0.10, 1)$, treatment outcomes $N(0.18, 1)$: a true lift of 0.08 standard deviations, 800 units per arm, seed 41. Predict whether the experiment detects the lift.

The estimate is $\hat\tau = +0.086$, close to the truth of 0.08, with a standard error of 0.051, so $z = 1.70$ — just short of the bar. A genuine effect, declared “not significant”. The second line explains why: the power of this design — the probability of rejecting $H_0$ when the lift really is 0.08 — is 0.36. Run the experiment a hundred times and sixty-four of them will fail to find the effect that is there. 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.

How big must N be?

Three quantities determine the design of any two-sample comparison: the minimum detectable effect $\tau$; the significance level $\alpha$, conventionally 0.05 two-sided, $z_\alpha = 1.96$; and the power $1 - \beta$, conventionally 0.80, $z_\beta = 0.84$. The required sample size per arm is

\[n \;\approx\; \frac{2\sigma^2\,(z_\alpha + z_\beta)^2}{\tau^2},\]

and the $\tau^2$ in the denominator is the whole story: halving the effect you want to detect quadruples the sample. For a conversion rate the outcome is Bernoulli, $\sigma^2 = p(1-p)$. Take a baseline of $p = 0.10$ and an absolute lift of one percentage point; predict the order of magnitude of $n$ per arm: $2 \times 0.09 \times 2.8^2 / 0.0001$.

About 14 000 per arm — 14 112 — to catch a one-point lift on a 10 % baseline with 80 % power. A two-point lift needs 3 528, one quarter as many; a five-point lift needs 565. (statsmodels.stats.power.NormalIndPower gives 14 745 for the same design on Cohen’s arcsine-transformed $h$; either is the right order of magnitude, which is what matters.) The previous cell’s $n = 800$ was doomed before it ran. Most “failed” experiments never had the sample size to succeed, and the failure is invisible in the output: a non-significant result from an under-powered test looks exactly like one from a well-powered test of a true null. Only the power calculation tells them apart, which is why a pre-registered experiment states $\tau$, $\alpha$, $1 - \beta$ and $N$ before the first unit is randomised.

The A/B protocol

Define the metric and pre-register the primary one; secondary metrics are exploratory.
Estimate the baseline and the minimum detectable effect; compute $N$.
Randomise units 50/50 (or another planned split).
Wait for $N$ to accumulate. Do not peek.
Test with the planned test — a $t$-test for means, a proportions $z$-test for binary outcomes.
Report the point estimate, the confidence interval and the p-value. Ship or kill.

Step 4 is the one most teams break, and the next cell measures the cost.

Sequential tests: the cost of peeking

Real teams peek. They check the dashboard at lunch and stop the experiment the first time the result “looks significant”. Simulate this under a true null. At each new observation recompute the running $z$-statistic, $z_k = \sum_{i \le k} y_i / \sqrt{k}$, and stop the first time $|z_k| > 1.96$. A fixed-$N$ test looks once and has a 5 % false-positive rate by construction. Predict the rate when the test may look after every one of 1 000 observations.

0.518. Over half of the null experiments produce a “significant” result somewhere in their first thousand observations. Every look is a fresh chance to cross the 1.96 line by luck, and the running $z$-statistic is a random walk that, given enough time, crosses any fixed line with probability approaching one — the law of the iterated logarithm guarantees it. Naive peeking turns a 5 % test into a coin flip, and the p-value of 0.04 the marketing manager saw on day three is not a p-value of 0.04 once the daily looks are accounted for. This also answers the question the rolling regime detector of §3.2 raised: 249 overlapping tests at $|\hat t| > 2$ are 249 peeks, and a detector that fires on 5 % of ordinary days will find a “regime change” somewhere in almost every year of noise.

The fix: spend your alpha

If repeated looks spend Type-I error, the cure is to budget it across the looks: fix the number of interim analyses $K$ in advance and raise the boundary at each so that the total false-positive probability across all $K$ looks is 5 %. Bonferroni splits $\alpha$ equally, $\alpha/K$ per look, and is conservative because the looks are positively correlated. Pocock (1977) uses a constant boundary chosen so that the joint crossing probability is exactly $\alpha$; for $K = 5$ equally spaced looks it is 2.413. O’Brien–Fleming (1979) is stringent early and relaxes toward 1.96 at the last look, preserving most of a fixed-$N$ test’s power; it is the clinical-trial default, generalised by Lan and DeMets to an $\alpha$-spending function for any look schedule. The cell reruns the null simulation with five looks — at 200, 400, 600, 800 and 1 000 observations — under each boundary. Predict which lands closest to 0.05.

One look at 1.96 gives 0.053 — the nominal 5 % up to simulation noise. Five looks at 1.96 give 0.139: even five peeks nearly triple the error, and continuous peeking took it to 0.518. Bonferroni’s boundary of 2.576 brings the rate to 0.030, below target — safe but wasteful, because it treats five correlated looks as independent. Pocock’s 2.413 gives 0.046, on target, which is what a boundary derived from the joint distribution of the five statistics should do. The cost of a properly designed sequential test is 5–15 % more samples than a fixed-$N$ test of equal power; the gain is a principled option to stop early when the effect is large, or for futility when it is clearly absent. Modern platforms — Optimizely, Microsoft’s ExP, Netflix — use always-valid p-values (the mixture sequential probability ratio test) that remain valid when checked every minute, which is $\alpha$-spending taken to its continuous limit.

Multi-armed bandits: earning while learning

When the goal flips from “estimate the effect precisely” to “earn the most reward while learning”, A/B is the wrong tool. A fixed-$N$ test deliberately sends half its traffic to the losing arm for the entire run; if the traffic is revenue — ad creatives, headlines, a recommender — every pull of the losing arm is money lost. A multi-armed bandit shifts traffic toward the apparently best arm as evidence accumulates, trading exploration (pulling arms you are unsure about) against exploitation (pulling the arm that looks best). Three classical policies. $\varepsilon$-greedy pulls the best-looking arm with probability $1 - \varepsilon$ and a random arm with probability $\varepsilon$, here 10 %, for ever. UCB1 (Auer, Cesa-Bianchi and Fischer, 2002) pulls the arm with the highest upper confidence bound $\hat\mu_k + \sqrt{2\ln t / n_k}$, so that rarely pulled arms get an optimism bonus that shrinks as they are tried. Thompson sampling (1933) keeps a Beta posterior for each arm’s success rate, draws one sample from each posterior, and pulls the arm with the highest draw — exploring exactly as much as its uncertainty warrants. Three arms with true conversion rates 0.08, 0.10 and 0.12, 4 000 pulls, seed 51; an oracle that always pulled arm 2 would expect $0.12 \times 4000 = 480$ successes. Predict which policy lands closest.

Thompson sampling earns 469 — a regret of 11 against the oracle’s 480 — and its pull counts show why: 3 724 of 4 000 pulls went to the best arm, with only 136 and 140 spent finding out that the other two were worse. UCB1 earns 430, regret 50; it kept exploring, with 897 and 1 119 pulls on the losing arms, because its optimism bonus decays only logarithmically. $\varepsilon$-greedy earns 411, regret 69, and its pull counts expose a second failure: 3 472 pulls on arm 1, the middle arm. An early run of luck made arm 1 look best, the greedy step locked onto it, and a fixed 10 % of random exploration was too little to overturn the mistake in 4 000 pulls — $\varepsilon$-greedy explores at the same rate whether it is sure or not, too much when it is right and too little when it is wrong. That is the exploration–exploitation trade-off made concrete, and why Thompson sampling is the production default at LinkedIn, Netflix and every recommender that experiments on live traffic. It is also a first sight of Chapter 5: the Beta posterior is Bayesian updating, and the policy is nothing more than “act on a draw from the posterior”.

When to use which, and what goes wrong

The design is chosen before the data, and it depends on what the business pays for. A fixed-$N$ A/B test when you need a defensible point estimate — for a board, a regulator, a journal — and the cost of a wrong call is one-time and large. A group-sequential test when stakeholders will peek anyway, or when units are expensive and early stopping saves real money. A bandit when learning is earning and nobody needs an unbiased estimate of the losers. Beyond peeking, three traps kill experiments silently. Multiplicity across metrics: measure twelve KPIs and declare victory on whichever turns significant, and you will almost always find one — the same error as peeking, spread across metrics instead of time. SUTVA violations: when treatment leaks between units (treated sellers competing with control sellers), the arms are not independent; cluster-randomise or use switchback designs. Novelty effects: a new interface shines for a week because it is new; pre-register a 30-day hold-out. And one free lunch: CUPED, regressing the outcome on its pre-experiment value, cuts the required $N$ by 30–60 % at zero cost.

The multiplicity trap has a general cure, worth seeing once in code because it returns in Chapter 4 (which coefficients are “significant” in a wide regression) and Chapter 7 (which of forty-two trading rules is real). If you test $m$ hypotheses at $\alpha = 0.05$ and all the nulls are true, you expect $0.05\,m$ false positives — ten from two hundred. Bonferroni controls the family-wise error rate, the probability of even one false positive, by requiring $p_i < \alpha/m$; it is strict, and loses power when many nulls are false. Benjamini–Hochberg (1995) controls the false discovery rate, the expected share of false positives among the rejections: order the p-values $p_{(1)} \le \dots \le p_{(m)}$, find the largest $k$ with $p_{(k)} \le k\alpha/m$, reject the first $k$. Two hundred candidate signals, each observed for 150 periods, ten of them real with a mean of 0.4 standard deviations, the rest pure noise, seed 9:

The naive screen at $p < 0.05$ rejects 19: all ten real signals and nine false ones — just below the 9.5 expected from 190 nulls — so nearly half of the “discoveries” are noise, and nothing in the list tells you which half. Bonferroni, at $p < 0.05/200 = 0.00025$, rejects nine, all real, and misses one real signal whose p-value of 0.00067 was not extreme enough: the price of controlling the chance of any false positive. Benjamini–Hochberg rejects exactly the ten real signals and no false ones — the sweet spot here, though its guarantee is only that the expected share of false discoveries is at most 5 %. A research desk reports all three columns. Decisions made after seeing the numbers — which metric, which subset, which sidedness, which stopping time — are p-hacking by any other name, and this section is a catalogue of the ways it happens by accident.

Nothing went wrong with the effect; the design had power 0.36, so it was more likely to miss than to find. $n \approx 2\sigma^2(z_\alpha + z_\beta)^2/\tau^2$ gives 14 112 per arm for 10 % → 11 % at $\alpha = 0.05$, 80 % power; a two-point lift needs a quarter of that, 3 528.

Peeking: under a true null, continuous peek-and-stop over 1 000 observations produced a 51.8 % false-positive rate; even five looks at 1.96 gave 13.9 %. Fix by spending $\alpha$ across pre-planned looks — Pocock’s constant boundary (2.413 for five looks gave 4.6 %), O’Brien–Fleming’s early-stringent boundary, or always-valid p-values — or by not peeking. The day-three “0.04” is not a p-value of 0.04.

Association: Linear, Monotonic, and Nonlinear

Two variables can be related in many shapes, and a single correlation coefficient summarises only one of them. Research desks keep a hierarchy of association measures and apply them in sequence: Pearson for linear, Spearman and Kendall for monotone, distance correlation for arbitrary dependence — and then, because a portfolio is tested on its worst days rather than its average ones, a conditional measure of what happens in the tail. The instructor’s rule is one sentence: “no association” means independence, and the only coefficient in the list that is zero if and only if two variables are independent is the last one. Everything before it can report zero for variables that are deterministically related. This section builds data that fools Pearson completely, finds which coefficient survives a monotone bend and which survives any bend, applies the hierarchy to a real feature, and asks whether the association between the S&P 500 and Tesla is the same on the market’s worst 5 % of days as on ordinary ones.

Can Pearson see a U-shape?

Pearson’s $r = \operatorname{Cov}(X,Y)/(\sigma_X \sigma_Y)$ is the cosine of the angle between the centred $X$ and $Y$ vectors. It measures linear co-movement and nothing else; its requirements — finite variances, a linear relationship, no dominant outliers — are assumptions, not guarantees. Make $y = (x - 5)^2 + \varepsilon$ for $x$ uniform on $[0, 10]$: a clean, strong, non-linear relationship you can see with your eyes. Predict what Pearson reports — near 0, near 0.5, or near 1?

$y$ is a deterministic function of $x$ plus noise with standard deviation 2, on a parabola symmetric about $x = 5$. Is Pearson’s $r$ near 0, 0.5 or 1 — and what would it be for $|x - 5|$ against $y$?

$r = 0.040$ on a relationship that explains almost all of $y$’s variance. The parabola is symmetric about its vertex: as $x$ rises, $y$ falls and then rises by the same amount, so the best straight line through the cloud is flat and the covariance is zero. Spearman is no better at 0.025, for a reason the next cells make precise. The third number shows what Pearson can see: fold the parabola by taking $|x - 5|$ and the correlation with $y$ is 0.943, because the folded relationship is monotone and nearly linear. Pearson is blind to anything non-monotone, and a feature that enters a model through a U-shape — volatility against return, leverage against default risk, temperature against energy demand — will be reported as “unrelated” by a correlation screen and thrown away.

Pearson sees straight lines only: a sine wave and an index panel

The lecture’s second demonstration is $y = \sin x + \varepsilon$ over $x \in [0, 100]$: $y$ is a deterministic function of $x$ up to noise, yet a sine wave rises as often as it falls, so $r$ is near zero. The second half of the cell applies the same coefficient to real data: a panel of daily index moves, 2 677 days from 2008 to 2018, in which spy is tomorrow’s change in the SPY ETF and every other column — spy_lag1, the S&P 500, the Nasdaq, the Hang Seng — is today’s move. Every correlation in the spy row asks whether today’s move in some market is a linear predictor of tomorrow’s SPY.

$r = -0.033$ for the sine. In the panel, every correlation in the spy row is near zero — −0.05 with today’s SPY, −0.06 with the S&P 500, −0.03 with the Nasdaq, +0.02 with the Hang Seng — while today’s moves are strongly correlated with each other: S&P 500 with Nasdaq 0.71, with today’s SPY 0.78, Nasdaq with today’s SPY 0.88. The contemporaneous block is the ordinary fact that US indices move together. The spy row says that no index’s move today is a linear predictor of tomorrow’s SPY — what an efficient market should look like, and also what a market with a purely non-linear dependence would look like. Whether today’s Hang Seng is any other kind of predictor is exactly the question Pearson cannot answer, and the reason the hierarchy exists.

Ranks instead of values: Spearman and Kendall

Spearman’s $\rho_S$ is Pearson’s $r$ computed on the ranks. Any strictly increasing transformation leaves ranks unchanged, so Spearman detects every monotone relationship perfectly — straight, curved or stepped — and is robust to outliers, because a rank cannot be extreme: the largest value is rank $n$ whether it is 3 or 3 000. Kendall’s $\tau$ counts pairs instead: of the $\binom{n}{2}$ pairs of observations, the share that are concordant (both variables move the same way between the two points) minus the share that are discordant. It is slower, slightly more robust to ties and outliers, and has a direct probabilistic reading — $\tau = P(\text{concordant}) - P(\text{discordant})$ — that Spearman lacks. Both reward monotonicity and only monotonicity. Predict them for three shapes: $y = \log x$, strictly increasing, so the rank of $y$ equals the rank of $x$ at every point; the parabola $(x - 2.5)^2$, symmetric about the centre of the $x$ range; and $\cos 30u$ for $u \sim N(0, 0.1)$, which oscillates several times across the sample.

Rank measures rescue monotone nonlinearity and nothing else. For the log, Pearson’s 0.921 — imperfect because a log is curved — becomes Spearman 1.000 and Kendall 1.000, exactly, because every pair is concordant. The parabola scores 0.075, 0.059 and 0.039; the cosine 0.036, 0.033 and 0.014. Both shapes are deterministic functions of $x$ and both are invisible to all three coefficients, because in each case $y$ goes up and comes back down, and every rise is cancelled by a fall. The rule: rank measures fix outliers and curvature-but-monotone, which covers a great many real relationships — a diminishing return, a threshold, a saturating response — but they cannot fix non-monotonicity. For “any dependence at all” you need a measure built on a different principle.

Distance correlation sees any dependence

Distance correlation, introduced by Székely, Rizzo and Bakirov in 2007, is zero if and only if the two variables are independent — the property the instructor’s rule demands. Form the matrix of pairwise distances $a_{ij} = |x_i - x_j|$ and likewise $b_{ij}$ for $y$; double-centre each (subtract row and column means, add back the grand mean) so that $A_{ij}$ measures how unusually far apart $x_i$ and $x_j$ are; the distance covariance is the mean of $A_{ij}B_{ij}$, normalised as Pearson normalises the ordinary covariance. If $X$ and $Y$ are dependent in any way, pairs of observations that are close in $x$ are systematically close or far in $y$ — the two distance patterns co-vary — and this holds for a parabola or a cosine just as much as for a line; direction plays no role. The lecture uses the dcor package; in the browser we build it from scipy’s pdist in eight lines.

In Colab

!pip install dcor --quiet
import dcor
dcor.distance_correlation(parabo["x"], parabo["y"])      # ≈ 0.49, 12 ms on 10 000 points

The hand-written version below forms an $n \times n$ matrix and is fine for a thousand points; dcor uses an $O(n \log n)$ algorithm and took 12 ms on the notebook’s 10 000-point parabola where the naive version took 4 s.

Distance correlation scores the log at 0.966, the parabola at 0.495 and the cosine at 0.366 — the notebook’s dcor gives 0.97, 0.49 and 0.35 on its own random draws. Only this column sees the two non-monotone shapes. Two cautions. Distance correlation is not on Pearson’s scale — 0.495 for a perfectly deterministic parabola shows that its values are compressed toward the middle — and it is never negative, since direction is not part of the construction. And it is a screen, not a model: a large distance correlation with a small Pearson says “there is structure here that a line does not capture — look at the scatter”, not what the structure is. A research desk that checks only Pearson would conclude the parabola’s two variables are unrelated and throw the signal away; the permanent lesson is to compute at least one nonlinear measure before saying “no relationship”.

Does today’s return predict tomorrow’s? The hierarchy on a real feature

Apply the hierarchy to a feature a quant would actually build. feat holds Apple’s next-day log return as the target, today’s return (lag1) and the five-day sum (ret5) as candidate predictors, for 2020–2021: 505 days. The notebook’s practice asks for the distance correlation of lag1 with the target, alongside Pearson.

Pearson between today’s and tomorrow’s return is −0.214 — a short-horizon reversal, a known and tradable pattern in daily equity returns, unusually strong in these two years because 2020 contained the March crash and its violent rebounds — and the distance correlation is 0.158. The two are not on the same scale, but the comparison still answers a question: the dependence between consecutive returns is essentially the linear reversal Pearson already sees. The five-day sum tells the opposite story at smaller magnitude — Pearson −0.067, distance correlation 0.131 — a weak linear reversal and a somewhat larger nonlinear dependence, the signature to investigate with a binned plot before deciding whether it is a feature or noise. The reading rule: when distance correlation is large and Pearson is small, look at the scatter; when both are small, there is nothing there; when both are moderate and comparable, the linear model of Chapter 4 is probably enough. And every one of these numbers is a sample statistic — a Pearson of −0.214 on 505 observations has a standard error of about 0.044 — while §3.3’s multiplicity warning applies with full force to a desk that computes this table for two thousand features.

Does association survive on the worst days?

returns.csv holds daily S&P 500 and Tesla returns, 2 515 days from 2015 to 2024. One Pearson number describes the decade. But a portfolio is tested on the days the market falls hardest, and the question that matters for diversification is conditional: on the S&P’s worst 5 % of days, how often was Tesla also having one of its own worst 5 % of days? Under independence the answer is 5 %. Predict the observed share.

Pearson is 0.465 over 2 515 days and Spearman 0.459 — the two agree, so the relationship is close to linear and not driven by a few outliers. The 5 % quantiles are −1.69 % for the S&P and −5.14 % for Tesla, which is three times as volatile. 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. That number, the tail co-exceedance, is the association that matters for a portfolio, and Pearson, computed on all days, does not report it. The upper tail is weaker: on the S&P’s best 5 % of days Tesla joined its own best 5 % on 23.8 % of them — stocks fall together more reliably than they rise together, the empirical basis for every “correlations go to one in a crisis” remark. Pearson within the 126 worst days is 0.355, lower than the overall 0.465; that is not evidence that the association weakens in the tail but a mechanical effect of conditioning on a slice of $x$, which truncates its variance and shrinks $r$. The conditional exceedance probability, not the conditional correlation, is the honest tail statistic — and 126 days carry the sampling uncertainty §3.2 taught you to bootstrap. The tail is where diversification is tested, and §3.5 is about the tail alone.

Only distance correlation (0.495 for $(x - 2.5)^2$; Pearson 0.075, Spearman 0.059, Kendall 0.039). Pearson sees linear co-movement, Spearman and Kendall see monotone co-movement, and a symmetric parabola has neither — every rise is cancelled by a fall. Distance correlation is built from pairwise distances, ignores direction, and is zero if and only if the variables are independent.

On the S&P’s worst 5 % of days Tesla was in its own worst 5 % on 32.5 % of them — 6.5 times the 5 % that independence implies (23.8 % in the upper tail: co-movement is asymmetric). Pearson within those 126 days is 0.355, lower than 0.465, because conditioning on a slice of $x$ truncates its variance and mechanically shrinks $r$; the conditional exceedance probability, not the conditional correlation, is the tail statistic to report.

Extreme Value Theory: GEV, Block Maxima, POT, Extreme VaR

Averages have the central limit theorem: whatever the parent distribution, the standardised mean of $n$ independent draws converges to a Normal, which is why the $t$-test of §3.2 works on returns that are nothing like Normal. The tails have a theorem of their own, and it is the reason a risk desk can say anything at all about events it has never seen. For most decisions that matter — bank capital, the height of a sea wall, the size of a trading limit — the question is not “what is the mean?” but “how bad is the worst case?”, and the worst case is a maximum, not a sum. Extreme value theory (EVT) is the statistics of maxima. 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 from 1985 to 1990 — so that Black Monday, 19 October 1987, is in the sample — fit the limit law, check it with a Q-Q plot, and score each block’s rarity; estimate the same tail index two more ways; and read off return levels and the extreme VaR, the loss that only one block in twenty should exceed — the number that names your highest pain threshold.

Same mean, same sd — same tail?

Draw 4 309 Normal returns with Apple’s mean and standard deviation (seed 71). They match the real returns on the first two moments exactly, and a mean–variance investor would regard them as the same asset. Predict the comparison at the 0.1 % quantile — the loss exceeded on one day in a thousand — and the number of days beyond four standard deviations.

Same mean, same standard deviation — and a 0.1 % quantile of −10.5 % against −5.8 %, nearly double. The 1 % quantiles are closer, −5.8 % against −4.6 %, because one day in a hundred is not yet deep in the tail; the 1 % conditional VaR — the average loss on the worst 1 % of days — is −7.8 % against −5.1 %; and Apple has 25 days beyond four standard deviations where the Normal has none (it would expect 0.27). Variance does not measure tail risk. Two portfolios with identical Sharpe ratios can have wildly different blow-up probabilities, and the difference lives in a region the standard deviation, an average over all days, barely registers. The rest of the section is about that region alone.

The Fisher–Tippett–Gnedenko theorem

Let $M_n = \max(X_1, \dots, X_n)$ for i.i.d. draws. Fisher and Tippett (1928) and Gnedenko (1943) proved that if there are sequences $a_n > 0$ and $b_n$ such that $(M_n - b_n)/a_n$ converges in distribution to anything non-degenerate, then the limit must belong to one family, the generalised extreme value (GEV) distribution:

\[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$ for a location $\mu$ and scale $\sigma$; $Q(p)$ is the quantile function, from which the extreme VaR will be computed. The single shape parameter $\xi$ — the tail index — sorts the family into three types. $\xi > 0$ is the Fréchet case: the tail decays like a power, $1 - F(x) \sim x^{-1/\xi}$, there is no upper bound, and moments of order above $1/\xi$ do not exist; asset returns, insurance claims and earthquake magnitudes live here. $\xi = 0$ is the Gumbel case, with an exponential tail, where maxima of Normal, lognormal and exponential parents fall. $\xi < 0$ is the Weibull case, bounded above at $\mu - \sigma/\xi$. The theorem’s power is that you do not need to know the parent distribution to know that its maxima, suitably standardised, are GEV — just as you do not need to know it to know that its means are Normal. The Student-$t$ of §3.1 is a bridge: a $t_\nu$ parent has a Fréchet tail with $\xi = 1/\nu$, so Apple’s fitted $\nu = 3.3$ already predicts $\xi \approx 0.30$ for the maxima of its returns, and the Dow’s block losses below land in the same neighbourhood.

scipy’s sign convention

scipy.stats.genextreme uses a shape parameter c with c = −ξ. Every fit in this section reports ξ = −c. scipy.stats.genpareto, by contrast, uses c = +ξ. The lecture notebook flags this with a row of exclamation marks, for good reason: a fit that reports $c = -0.35$ has found a heavy tail, and reading it as light turns a Fréchet into a Weibull. The debug-yourself below is exactly that error.

With $\mu = 0$ and $\sigma = 1$, the support ends at $\mu - \sigma/\xi$: at $+3$ on the right for $\xi = -1/3$, and at $-3$ on the left for $\xi = 1/3$ (a Fréchet has a lower bound and an unbounded right tail). The Weibull curve stops dead at 3; the Gumbel decays exponentially; the Fréchet curve is the one that keeps going — visibly above the other two from about $x = 3$ onward. That gap is the entire subject of the section.

The sample: Dow Jones log returns, 1985–1990

Six years of daily log returns, 1 496 trading days. Think about which single day was the worst and roughly how bad it was — and note the .copy() on the slice, which Chapter 2 explained: the lecture notebook wrote Min20 onto a slice of dji and drew a SettingWithCopyWarning.

The worst day is 19 October 1987 with a log return of −0.2563 — a 22.6 % fall in the index, which in log terms is −25.6 %; this is the one day in the chapter on which the difference between simple and log returns is visible. The daily standard deviation over the six years is 1.28 %, so Black Monday was a 20-sigma day. Under a Normal model the probability of a 20-sigma move is of order $10^{-89}$. The plot shows why the Normal is the wrong model: one spike dwarfs everything else, and two clusters of smaller spikes — October 1987 and October 1989 — sit far outside the band of ordinary days. EVT starts from the observation that the spike is data, not an error, and asks what distribution could have produced it.

Block minima: the worst day of every 20-day block

The first way to collect extremes is block maxima: divide the sample into blocks of $n$ observations and keep the maximum of each. For losses, keep the minimum of each block and flip the sign. rolling(20).min() gives the running 20-day minimum; taking every 20th row with iloc[::20] makes the blocks non-overlapping. With 1 496 days, minus 19 lost to the rolling window, there are 74 blocks. Predict that number before you run.

Seventy-four blocks, each contributing one number. The first three minima — −1.08 %, −1.01 %, −0.60 % — come from calm months in early 1985; the median block minimum is −1.65 %, so an ordinary month’s worst day is about 1.3 daily standard deviations down; the worst block contains 19 October 1987 at −25.63 %, fifteen times the median. The red dots make the method’s character visible: one dot per block whatever happened in it. A quiet month gives a minimum of half a per cent; the crash month gives −25.6 %; both count as one observation. Block maxima are rate-limited by design — the method cannot see that October 1987 contained three terrible days rather than one, and it wastes the information in every calm block’s second-worst day. The compensating virtue is that the 74 observations are nearly independent (a month apart) and the theorem applies to them directly.

Peaks over threshold: keep every day below −1.5 %

The second way throws away the calendar and keeps every observation beyond a threshold $u$. Set $u = 1.5\,\%$ and collect every day whose log return fell below $-u$. Predict the order of magnitude: of 1 496 days, how many?

Ninety-one days — 6.1 % of the sample — closed more than 1.5 % down. The yearly counts show what block minima flattened: one exceedance in 1985, then 16, 29, 14, 10 and 21 — the exceedances cluster in 1987 and 1990, the two stress episodes. POT sees the clusters; block minima see one point per month. POT’s price is the choice of $u$. Too high and the sample is tiny; too low and the observations are no longer “extreme”, and the limit law for exceedances no longer applies. Clustered exceedances also violate the independence the theory assumes, which practitioners handle by declustering — keeping only the worst day of each run — at the cost of throwing some data back.

Warm-up: GEV on the block maxima of a Normal

Before fitting to real data, check the machinery on a case whose answer is known. The Normal is in the Gumbel domain of attraction: block maxima of Normal draws converge to $\xi = 0$. Draw 10 000 standard Normals, cut them into 100 blocks of 100, take each block’s maximum, and fit. A fitted $\xi$ near zero — slightly negative is common in finite samples — confirms the method; $\xi = 1/3$ would be a bug.

$\xi = -0.058$, effectively zero, with location 2.33 and scale 0.40. The mean of the 100 block maxima is 2.54, close to the theoretical expected maximum of 100 standard Normals of about 2.5, and the GEV’s mode $\mu$ sits a little below the mean, as it should for a right-skewed distribution. The machinery works — and it reminds you what Normal maxima look like: a 100-day worst case near 2.3 standard deviations, with a scale of 0.4. Compare that with what the Dow is about to show.

Fit the GEV to the Dow’s block losses

Minima become maxima by a sign flip: fit to $-\text{Min20}$, the block loss. gev.fit takes optional starting values for loc and scale; they only seed the optimiser. For $\xi > 0$ the support starts at $\mu - \sigma/\xi$. Predict that bound from the parameters the lecture reports — $\mu \approx 0.0136$, $\sigma \approx 0.0070$, $\xi \approx 0.354$ — before the cell reveals it.

$\xi = 0.354$, $\mu = 0.0136$, $\sigma = 0.0070$. The tail index is positive: Fréchet. The 20-day worst loss of the Dow has a polynomial tail, and moments of order above $1/\xi \approx 2.8$ do not exist — the variance of the block loss is finite, the skewness is not. The support begins at $-0.0063$, slightly below zero, which is harmless (a block “loss” of −0.6 % would be a block in which every day rose; none occurred, but the fitted family allows it). The histogram shows why a symmetric family would fail: the bulk of the 74 losses sit between 1 % and 2 %, and then a thin, long right tail runs out to 25.6 %. The red curve follows the body and keeps going into the tail, as a Fréchet should. Set the two fits side by side. Normal block maxima: $\xi \approx 0$, scale 0.40 in units of the daily standard deviation. Dow block losses: $\xi = 0.35$, and the scale of 0.70 % is more than half the daily standard deviation of 1.28 %. The Normal world’s worst month is a modest, predictable multiple of an ordinary day; the Dow’s worst month has no such ceiling.

A colleague fits the same GEV and reports that the Dow’s block losses have a bounded tail — good news for the risk limit. The cell runs without error. What went wrong?

The bug is the sign convention. gev.fit returns c = −0.354, and the cell reads c as $\xi$, concludes $\xi < 0$ — Weibull — and even computes an “upper endpoint” of $\mu - \sigma/c = 0.0335$, announcing that no 20-day block can lose more than 3.4 %. The sample contains a block that lost 25.6 %. The fix is one character: xi = -c. This is the most consequential sign error in applied risk management, because a Weibull verdict is the one a risk committee wants to hear, and a fit that appears to support it will not be questioned. Always print $\xi$ alongside the worst observed loss; if the fitted support excludes an observation you have, the model is wrong, not the observation.

Diagnostic: Q-Q plot against the fitted GEV

A fitted distribution always produces parameters; whether it fits is a separate question, and the Q-Q plot answers it. probplot accepts a frozen scipy distribution, so the comparison is against the fitted GEV itself. Decide which point will sit far above the 45° line at the top right.

The largest block loss — Black Monday’s 25.6 % — sits alone above the line at the top right: it is more extreme than even the fitted Fréchet expects the worst of 74 blocks to be. Everything else follows the line. The fit is good for the body of the tail and conservative in the wrong direction for one event: it underpredicts the crash. That is the honest state of the art. With 74 observations the tail index has a wide confidence interval, and the single most extreme point in any sample is exactly the point the model is least sure about. A risk analyst reports the fit, reports the outlier, and does not pretend the model has domesticated it.

Extreme scores: how rare was each block?

Once the GEV is fitted, every block can be scored: $\text{score} = 1 - F_{\text{GEV}}(\text{loss}) = P(\text{a block is worse than this one})$. A small score is a rare block, and the lecture’s cut-off of 0.10 defines a “rare event”. Predict which date has the smallest score and how many blocks fall below 0.10 — if the model were exactly right, 10 % of 74, about seven.

Six of 74 blocks score below 0.10 — a little fewer than the 7.4 expected if the model were exactly right, which is fine. Black Monday’s score is 0.0007: the fitted GEV says a worse 20-day block arrives once in roughly 1 400–1 500 blocks (the cell prints 1 464 from the unrounded score; the slides quote 1 400 from the rounded one), and at 12.3 blocks a year that is once in a little over a century. Two others — 8 January 1988 (score 0.021) and 13 October 1989, the “Friday the 13th mini-crash” (0.021) — are one-in-fifty-block events, about one in four years. Three of the six cluster in the twelve months after the crash, when volatility stayed high, and the price chart reads as a calendar of stress. A score is a probability with a model attached, and it converts “that looked bad” into “that was a 1-in-1 400 block under a Fréchet with $\xi = 0.35$” — a sentence a risk committee can argue with.

POT with the generalised Pareto distribution

Exceedances have their own limit law. The Pickands–Balkema–de Haan theorem (1974–75) says that for a high enough threshold $u$, the excess $Y = X - u$ given $X > u$ is approximately generalised Pareto:

\[P(X - u > y \mid X > u) \approx \Big(1 + \frac{\xi y}{\beta}\Big)^{-1/\xi},\]

with the same tail index $\xi$ as the GEV of the block maxima and a scale $\beta$ that depends on $u$. Because the GPD describes only the excess, the unconditional quantile is reassembled from the exceedance rate $N_u / n$:

\[\text{VaR}_q = u + \frac{\beta}{\xi}\Big[\Big(\frac{n}{N_u}(1 - q)\Big)^{-\xi} - 1\Big].\]

genpareto.fit with floc=0 fixes the location at the threshold, as the theory requires; note that genpareto’s shape is $+\xi$, the opposite convention from genextreme.

Ninety-one exceedances give $\xi = 0.436$ against the GEV’s 0.354. Two estimators, two samples — one number per month versus every bad day — and the same verdict: a heavy tail with an index between a third and a half. Agreement in sign and rough size is the sanity check; disagreement would send you back to the block length or the threshold. The fitted GPD’s mean excess, $\beta/(1 - \xi) = 1.14$ %, matches the observed mean excess of 1.19 % — a second check, since for a GPD the mean excess over the threshold is exactly that ratio. The GPD’s 99 % one-day VaR is 3.27 %, just above the empirical 1 % quantile of 3.05 %. With 1 496 days the empirical quantile is fine at 1 % — fifteen observations lie beyond it. The GPD earns its keep at 0.1 %, where the empirical quantile is the single worst day and there is nothing left to count.

Hill’s estimator: the tail index a third way

A Fréchet tail is a power law, $P(L > x) \sim x^{-\alpha}$ with $\alpha = 1/\xi$, and a power law is a straight line on log–log axes. Hill (1975) reads its slope off the top $k$ order statistics of the losses, $L_{(1)} \ge L_{(2)} \ge \dots$:

\[\hat\alpha_k = \Big[\frac{1}{k}\sum_{i=1}^{k}\ln\frac{L_{(i)}}{L_{(k+1)}}\Big]^{-1},\]

the reciprocal of the average log-distance of the top $k$ losses from the $(k+1)$-th. It needs no distributional fit, only the assumption that the top $k$ are in the power-law regime — and $k$ plays exactly the role $u$ played for POT. Predict which values of $k$ agree with the GPD’s $\xi = 0.44$: at $k = 25$ the threshold is a loss of about 2.5 %, deep in the tail; at $k = 200$ it is 0.8 %, an ordinary down day.

$\hat\xi$ = 0.468, 0.423 and 0.441 for $k$ = 25, 50 and 100 — the plateau, bracketing the GPD’s 0.436 — then 0.629 at $k = 200$, where the threshold of 0.82 % is inside the body of the distribution and the power-law assumption has failed: the extra points are contamination from the centre, not information about the tail. Three estimators, three samples, one verdict: GEV 0.35 on 74 block losses, GPD 0.44 on 91 exceedances, Hill 0.42–0.47 on the top 25–100 days. The spread between them is the honest uncertainty in the tail index. Picking $k$ is the craft, exactly as picking $u$ was: too few points is noise, too many is the body, and the plateau in between — $\hat\xi$ plotted against $k$, the Hill plot — is what a practitioner looks for.

Return levels: the 1-in-T-day loss

The $T$-day return level is the loss $x_T$ exceeded on average once every $T$ days: $P(L > x_T) = 1/T$. Hydrologists call it the 100-year flood; a risk desk calls it the one-in-a-thousand-day loss. Two estimators: the empirical $(1 - 1/T)$-quantile of the losses, and the fitted GPD extrapolated from the threshold, $x_T = u + \frac{\beta}{\xi}\big[(\frac{k}{n}T)^{\xi} - 1\big]$, which is the POT VaR formula with $1 - q = 1/T$. Predict the empirical method’s fatal limit: at $T = 2500$, ten years of trading days, what can 1 496 days of data say?

At $T = 250$ (one year) the two agree: 4.86 % against 4.71 %, because a one-in-250-day loss has six observations beyond it. At $T = 1000$ the GPD says 8.87 % while the empirical 7.77 % rests on an interpolation between the second and third worst days (8.38 % and 7.16 %) — two observations, and the answer would change materially if either were removed. At $T = 2500$ only the GPD answers: 13.2 %; the empirical quantile is undefined, because the $(1 - 1/2500)$-quantile of 1 496 numbers does not exist. And Black Monday’s 25.6 % is the reminder that even the GPD’s ten-year figure is one draw from a tail whose index is itself uncertain. Extrapolation with a fitted shape is what turns six years of data into a ten-year estimate; that is what EVT buys, and what the empirical distribution — which cannot go beyond its largest value, the bootstrap’s failure mode from §3.2 — cannot.

Extreme VaR: the highest pain threshold

Invert the fitted GEV at probability $\alpha$: the loss that only a fraction $\alpha$ of 20-day blocks will exceed. From the quantile function above, with the sign flipped back to a return,

\[\text{EVaR}_\alpha = -\Big(\mu + \frac{\sigma}{\xi}\big[(-\ln(1-\alpha))^{-\xi} - 1\big]\Big).\]

Predict the 5 % extreme VaR from the fitted parameters — and compare it with the 5 % quantile of daily returns, which is about −1.6 %.

$\xi = 0.354$, $\mu = 0.0136$, $\sigma = 0.0070$, $\alpha = 0.05$. Is the extreme VaR closer to −2 %, −5 % or −10 % — and what would a Gumbel ($\xi = 0$) with the same $\mu$, $\sigma$ say?

The 5 % extreme VaR is −5.06 %. A one-day 5 % VaR says “on one day in twenty, expect to lose more than 1.6 %”. The extreme VaR says something different in kind: in one 20-day block out of twenty, expect a single day worse than −5.1 %. In the six years of the sample, 4.1 % of the 74 blocks — three of them: the crash block, February 1988 and October 1989 — did breach it, close to the nominal 5 %. Tighten $\alpha$ to 1 % and the threshold is −9.5 %, breached once (Black Monday); at 0.1 % it is −22.3 %, which the crash still exceeded — the Q-Q plot’s outlier expressed in money. The final line is the cost of the tail index: a Gumbel with the same location and scale would put the 5 % limit at −3.45 %, a third less severe, and seven blocks rather than three would have breached it. The whole difference is $\xi = 0.35$ against $\xi = 0$ — one parameter, and the one that only the extremes can estimate.

The extreme VaR is the number a hard loss limit should be sized to. The daily VaR is sized to an ordinary bad day; the extreme VaR is sized to an ordinary bad month, and a month is the horizon over which a position that has gone wrong is unwound. The lecture’s decision memo makes exactly this recommendation — replace a daily-VaR stop with a 5 % EVaR limit of −5.1 % per block, reviewed quarterly — with the caveats any such memo must carry: 74 observations give $\xi$ a wide confidence band (bootstrap it before sign-off); the window ends in 1990 and should be re-estimated across 2008 and 2020; and the GPD and Hill estimates against the GEV’s 0.35 are a sanity check, not a proof.

Mistakes library: “25-standard-deviation moves” (August 2007)

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 aeons. Black Monday, a 20σ day under the same arithmetic, scored 0.0007 under the GEV: rare, but a number.

Lesson: counting sigmas presumes a Gaussian. The tail index, block maxima and the extreme VaR exist precisely because “how many σ?” is the wrong question to ask about a tail.

$\xi > 0$ — Fréchet — a polynomial tail with no upper bound; moments above order $1/\xi$ do not exist. scipy.stats.genextreme uses c = −ξ, so a fitted c = −0.35 is $\xi = +0.35$ (heavy); reading c as $\xi$ turns the Dow’s tail into a Weibull with an “upper endpoint” of 3.4 % in a sample containing a 25.6 % loss. genpareto uses c = +ξ.

In one 20-day block out of twenty, expect a single day worse than −5.1 % (4.1 % of the 74 blocks did). The daily VaR describes an ordinary bad day; the extreme VaR describes an ordinary bad month — the horizon over which a losing position is unwound, so the horizon a hard loss limit should be sized to. GEV on block losses gave $\xi = 0.35$, GPD on 91 exceedances 0.44, Hill on the top 25–100 days 0.42–0.47: one verdict, heavy tail, with the spread as the honest uncertainty.

Chapter Wrap-up

You can now describe a variable by its distribution rather than its average, and you carry a set of numbers that anchor the habit. The empirical CDF of 4 309 Apple returns showed 76.9 % of days within one standard deviation against the Normal’s 68.3 % and only 98.5 % within three against 99.7 %, and DKW said the whole staircase is within 0.02 of the truth with probability 93.6 %; a Student-$t$ fitted by maximum likelihood gave $\nu = 3.29$. The bootstrap reproduced the mean’s known standard error of 0.0007 and then measured what no formula gives — 0.198 for a median, 1.55 for a kurtosis of 4.58 — and BCa moved the kurtosis interval from $[1.5, 7.3]$ to $[2.2, 9.1]$. Shapiro–Wilk rejected normality at $10^{-38}$, a scale-corrected KS test could not reject the $t$ ($p = 0.27$), the one-sample $t$-test found 2017’s mean positive at $\hat t = 2.16$, and a permutation test found 2017 and 2018 indistinguishable at $p = 0.18$. An under-powered A/B test missed a real effect with power 0.36, a one-point lift needed 14 112 per arm, peeking turned 5 % into 51.8 %, and Thompson sampling beat $\varepsilon$-greedy by 469 to 411. Distance correlation saw the parabola (0.50) that Pearson, Spearman and Kendall could not, and Tesla joined the S&P’s worst 5 % of days a third of the time. On the Dow, three estimators agreed on a Fréchet tail — $\xi$ = 0.35, 0.44, 0.42–0.47 — from which a 5 % extreme VaR of −5.1 % per 20-day block followed, three times the daily 5 % quantile.

The through-line is one discipline: before computing a statistic, predict its distribution; after computing it, attach its uncertainty; before believing it, ask how many other statistics you looked at. Every later chapter inherits that discipline and adds a model. Chapter 4 puts the fat-tailed, associated variables of this chapter into predictive models — multiple regression, model selection, trees and boosting, logistic classification — where the residual diagnostics are the Q-Q plot of §3.1 applied to errors, and “which coefficients are significant?” is the multiplicity question of §3.3 in a wide design matrix. Chapter 5 replaces the confidence interval with a posterior and the Thompson sampler’s Beta distribution with a full Bayesian model of risk and reward. Chapter 6 adds time: the rolling regime detector becomes a Markov-switching model, the fat tails become GARCH innovations, and the bootstrap becomes a block bootstrap. Chapter 7 puts forty-two trading rules on trial with the false-discovery machinery this chapter introduced, and uses conformal prediction to give the intervals of §3.2 a finite-sample guarantee.

The slides for this chapter, at https://statpython.pages.dev/topic3.html, rehearse every section here with a prediction before every cell — same datasets, same seeds, same printed numbers; use them to test yourself, and use this chapter to understand why the predictions came out as they did.

Predict the distribution — its centre, spread, shape and tail — before you compute the statistic; attach its uncertainty after; and count how many other statistics you looked at before you believe it. A 25-standard-deviation “surprise” is a wrong model, not bad luck.

← Chapter 2: DataFrames · Contents · Chapter 4: Statistical Predictive Models →