Chapter 3: Reshaping Statistics
Chapter Introduction
Statistics supports decisions across retail, operations, customer service and finance. Sections 3.1–3.4 use customer waiting times, deliveries, branch process changes, product experiments and demand rankings to explain the general methods. Section 3.5 applies those foundations to portfolio losses and extreme tails.
Define the population and observation unit before summarising a distribution; distinguish uncertainty about an estimate from variation among observations; preserve the comparison and all tests tried; and evaluate association on the appropriate matched pairs and cohorts. Every new small business dataset is an invented teaching illustration, labelled as such. Existing financial examples remain as transfer applications with their original assumptions and numerical results.
Financial transfer case. The classroom handoff contains 639 candidates, 630 valid labels, 315 predictions with stock and date columns and 34 realised monthly portfolio returns. Nine predictions have pending labels. The fixed Ridge baseline selects the three highest scores, targets 99% equal-weight long-only investment and charges 10 basis points per dollar traded. The small current-survivor sample uses assumed data clocks and static sectors. It illustrates diagnosis and reproducibility; it does not evaluate a full MyStrategy long-short implementation.
The Apple, index, Dow, A/B and Bandit worked examples remain in their original sections. Study them alongside the strategy diagnostics; class timing sets discussion priorities.
Study route: business evidence and a financial transfer case
Study every assigned section and use the animations to predict, reveal and compare a result before running the complete Python. The primary questions concern customer and operating measurements, experimental comparisons and later prediction outcomes. The existing financial diagnostic notebook demonstrates transfer of the same methods; the class-specific teaching plan provides the actual discussion and submission instructions.
Report the question, population, observation unit, method, uncertainty and limitation. Random assignment, paired matching and resampling have different roles: matching does not itself remove confounding, and resampling does not create a randomised experiment.
How to use the numbered learning path
Numbered parts preserve the teaching order. Examples run automatically after opening or refreshing; practice and repair cells run manually. Edit inputs in the current cell; refresh clears edits.
Table of Contents
Feature Distributions and Diagnostics
| Part | Question |
|---|---|
| 3.1.1 Observations and distribution concepts | What population and distribution are we describing? |
| 3.1.2 Diagnose features across a fixed split | Does the later sample resemble development data? |
Observations and distribution concepts
Question: What population and distribution are we describing?
Result: Distinguish a sample, an ECDF, a density and a fitted distribution.
Population, sample and the quantity we estimate
A population distribution describes the random variable; a sample supplies observed units: customer visit, delivery, branch-week or trading day. Repeated observations within a branch or day can be dependent.
\[\mu=E[X],\qquad \bar x=\frac{1}{n}\sum_i x_i,\qquad s^2=\frac{1}{n-1}\sum_i(x_i-\bar x)^2.\]
The parameter \(\mu\) is unknown; \(\bar x\) is an estimator that changes across samples. Under i.i.d. sampling, \(s^2\) is unbiased for variance; \(s\) need not be unbiased for standard deviation.
Explain: Why can 639 rows describe a distribution without providing 639 independent observations?
CDF, ECDF and quantiles
\[F(x)=P(X\le x),\qquad \widehat F_n(x)=\frac{1}{n}\sum_i\mathbf1\{x_i\le x\}.\]
\[q_p=\inf\{x:F(x)\ge p\}.\]
The median is \(q_{0.5}\). Record the software quantile convention: interpolation can give an unobserved value. A 95th percentile is a location, not a confidence interval.
Density estimation: histogram, KDE and bandwidth
For a continuous variable, probability is area, not the height of its density:
\[P(a<X\le b)=\int_a^b f(x)\,dx,\qquad \widehat f_h(x)=\frac{1}{nh}\sum_i K\!\left(\frac{x-x_i}{h}\right).\]
Histogram bins and KDE bandwidth \(h\) affect apparent shape. Small \(h\) follows detail and noise; large \(h\) can merge modes.
Predict: Does a narrow density peak imply more probability than a wider, lower peak? Compare their areas.
Distribution assumptions: Normal, Student-t and diagnostics
Normal tails are thin; smaller Student-t \(\nu\) gives heavier tails. Its variance exists only for \(\nu>2\).
A Q–Q plot pairs sample and reference quantiles. Tail departures question the model; a straight centre does not validate extreme tails.
The KS distance is \(D=\sup_x|\widehat F_n(x)-F_0(x)|\). Its standard calibration assumes independent observations and a specified continuous \(F_0\). Fitting parameters on the same sample changes the null distribution; use a suitable calibrated test or simulation.
Check: A descriptive distribution fit is not proof of independent or stationary returns.
Diagnose features across a fixed split
Question: Does the later sample resemble development data?
Result: A distribution comparison with fixed training-only cutoffs.
Builds on Observations and distribution concepts.
Business application: service and delivery distributions
Use customer visits and delivery times to define the population, inspect distributions and compare a fixed later sample. The values below are invented for teaching; no company performance is claimed.
Define the population before counting the sample
A service manager wants the mean wait of all customers. Would complaints alone represent that population?
Invented teaching data.
One row is one customer visit; waits are minutes. Complaint-only selection over-represents long waits in this fixture. A larger selected sample need not remove selection bias.
Count deliveries below a promised time
Four delivery times are 10, 20, 20 and 40 minutes. What fraction arrive within 20 minutes?
Invented teaching data.
The ECDF counts observations at or below the chosen cutoff, including ties. A sample quantile reverses the question; linear interpolation may give an unobserved time.
Use histogram area to compare service times
Two calls fall in each of three unequal-width bins. Do the bins need equal density heights?
Invented teaching data.
Each bin has probability 2/6. Density divides that probability by bin width, so all bar areas equal 1/3 while their heights differ. KDE bandwidth changes smoothing, not the observed calls.
The Chapter 2 handoff
Finance transfer example: The following existing market-data case applies the same general methods. Its sample and limitations remain explicit.
Chapter 2 supplied aligned features, mature labels, walk-forward predictions and a reconciled account. Read the exported tables rather than rebuild a model in this chapter. The nine-stock teaching sample contains 639 candidates, 630 valid labels and 315 predictions, including nine predictions whose labels are still pending. A missing future label must not remove a current candidate.
Stock/date columns and training cutoffs are checked before any statistic. The manifest records the source and file hashes. This public current-survivor sample uses assumed availability clocks and static sectors; it demonstrates methods, not historical investability.
Feature distributions across a fixed split
Use 2019–2021 as the development period and 2022 onward as the held-out diagnostic period. Inspect the distribution before choosing a transform. The split is fixed for teaching; looking at the later period must not trigger repeated tuning presented as untouched OOS evidence.
ECDF and a training-only clipping rule
Derive the diagnostic clipping rule from development data and freeze it for later observations. Keep the raw feature and report the fraction affected; clipping is not a universal cleaning policy.
The ECDF is descriptive under dependence. The i.i.d. DKW band cannot automatically certify these stock-months. Use density, Student-t and Q–Q diagnostics to investigate the shape.
Fit cutoffs on development dates and apply them unchanged. Later outcomes must not determine preprocessing settings.
A large development-to-held-out quantile shift. Distribution shifts need diagnosis; the shift alone does not identify a causal or profitable effect.
Worked examples: ECDF, KDE and distribution diagnostics
Scope of these examples. The DKW bound and nominal test probabilities below rely on sampling assumptions. Financial returns need not be i.i.d.; parameters fitted on the same observations require goodness-of-fit calibration that refits them in each simulated sample. Treat the printed KS p-values against fitted distributions as uncalibrated diagnostic output, not valid nominal evidence of acceptance. A fitted Student-t tail parameter is a model estimate, not a proof of the true moment structure.
Compare three representations: the empirical distribution (histogram and ECDF), a smoothed density (KDE), and a fitted parametric family. The ECDF preserves the observed values and frequencies but drops time order; none of these summaries establishes representativeness or independence.
Population vs sample: why pandas divides by n − 1
Treat the 63 fixed values as an illustrative finite population. Its descriptive mean and variance use ddof=0. Draw ten without replacement and compare sample summaries. Under i.i.d. sampling, dividing squared deviations by \(n-1\) makes sample variance unbiased for the population variance; standard deviation need not be unbiased. For sampling without replacement from a finite population, ddof=1 targets the finite-population variance defined with denominator \(N-1\), rather than the displayed ddof=0 variance. pandas defaults to ddof=1; NumPy defaults to ddof=0. The imported stats is SciPy’s distribution and statistical-test module; its diagnostic tests appear later in §3.1.
The population SD is 28.505 with ddof=0. Seed 1 gives sample mean 49.7 and SDs 23.372 (ddof=1) and 22.172 (ddof=0), about 5% apart. Changing the seed changes the estimator. sample(frac=0.8, replace=True) draws with replacement: observations may repeat or be omitted. Recomputing a statistic on such resamples is the bootstrap primitive developed in §3.2; the resampling unit must match the dependence structure.
Univariate exploration: Apple, 2005–2022
Apple’s daily log return is \(r_t=\ln S_t-\ln S_{t-1}\), with 4,309 observed returns from January 2005 to February 2022. Log returns add over time; simple return is \(e^{r_t}-1\), differing by terms of order \(r_t^2\). The diagram traces the first prices, and agg summarises the full sample. Predict whether daily-return excess kurtosis is positive or negative.
Mean daily log return is 0.12%, sample SD 2.08%, skew −0.28 and Fisher excess kurtosis 5.99 (Normal: 0). Price skew 1.98 and excess kurtosis 3.36 describe the sampled price levels; they need not describe return risk. Q1–Q3 spans 2.03%, containing the middle half of observed returns, while the sample extremes range from −19.8% to +13.0%. Fourth-moment summaries can be dominated by a few extreme days; §3.2 examines their sampling uncertainty.
From theoretical to empirical: build an ECDF from five numbers
For the original illustrative sample [3, 1, 4, 1, 5], count observations at or below a threshold and divide by five. The two values equal to 1 produce a jump of \(2/5\). 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\},\]
Each observation adds \(1/n\); ties add together. The ECDF makes no parametric shape assumption and is the non-parametric maximum-likelihood estimate of the sampling CDF. Sort the 4,309 studentised Apple returns and read the staircase. Before running, predict its height at sorted observation 2,155: the fraction is \(2155/4309\) and that observation is the median.
The five-number ECDF gives 0.4, 0.6 and 1.0 at 1, 3 and 5. Apple’s median studentised return is −0.008; the cell’s selected 10%/90% order statistics are −1.110/+1.086, within the Normal’s ±1.282. Together with the range counts, this shows a concentrated centre and more extreme observations than a variance-matched Normal. The ECDF avoids bin tuning, but the mean and SD used to studentise it remain sensitive to extremes.
Why doesn’t a histogram do the job? The ECDF of a bimodal mixture
Histogram edges and widths can change the visible number of humps. The ECDF has no bin or bandwidth choice, but still reflects sampling variation. The original illustrative mixture (seed 7) draws with equal probability from \(N(5,5^2)\) or \(N(-10,2^2)\). Predict where its staircase rises steeply and where it rises slowly.
The quartiles are −10.00, −5.18 and 5.03. The staircase rises quickly near each component’s centre and more slowly between them. The valley is sparse rather than empty: a median between the two centres need not be typical of either component. Compare histogram widths while the sample and ECDF remain fixed.
How fast does it converge? Glivenko–Cantelli and the DKW bound
Under i.i.d. sampling, Glivenko–Cantelli gives uniform almost-sure convergence: \(\sup_x|\hat F_n(x)-F(x)|\to0\). The Dvoretzky–Kiefer–Wolfowitz inequality bounds the finite-sample worst gap,
\[P\Big(\sup_x |\hat F_n(x) - F(x)| > \varepsilon\Big) \le 2e^{-2n\varepsilon^2},\]
The bound is distribution-free under those assumptions: it does not depend on the common CDF \(F\). It is not automatically valid for serially dependent Apple returns. Predict \(2e^{-2n\varepsilon^2}\) at \(n=4309\), \(\varepsilon=.02\), and whether a bound exceeding one is informative.
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\).
For an i.i.d. sample with \(n=4309\), the bound at \(\varepsilon=.02\) is 0.064, giving at least 93.6% simultaneous coverage within .02. The inverted 95% band has half-width .0207. At \(n=100\), the raw bound 1.846 is vacuous; at \(n=1000\), .899 is weak. These are hypothetical sampling guarantees, not certified coverage for Apple’s dependent returns. DKW bounds the same uniform gap used by KS; KDE smooths the ECDF’s point masses. A primary proof states the independent, common-distribution setup.
Estimating a density without a formula: kernel density estimation
The empirical distribution assigns mass to sample points. Gaussian KDE replaces each mass with a kernel of bandwidth \(h\) and sums the kernels:
\[\hat f_h(x) = \frac{1}{nh}\sum_{i=1}^{n} K\!\left(\frac{x - x_i}{h}\right),\]
Small bandwidths follow sample detail; large ones can merge real modes. For a Normal target, Silverman’s rule \(h^\star\approx1.06\hat\sigma n^{-1/5}\) is a starting value derived from an integrated-error criterion. With studentised Apple SD 1 and \(n=4309\), it is about .199. It is not a distribution-free choice. KernelDensity from sklearn.neighbors fits this smoothed density; the cell compares three bandwidths.
Silverman gives .199 on the studentised Apple scale, .00413 on raw returns, and 1.97 for this mixture. The mixture’s between-component separation contributes strongly to its SD, challenging the one-hump reference rule. Sweep \(h\) rather than treating that rule as a verdict. KernelDensity expects a 2-D input, hence mix[:, None]; score_samples returns log-density, so use np.exp before plotting.
Estimating a density with a formula: fit a Student-t
Fit a parametric family by maximum likelihood. scipy.stats.t.fit estimates degrees of freedom \(\nu\), location and scale. As \(\nu\to\infty\) the Student-t tends to Normal; its density tails decay as \(|x|^{-(\nu+1)}\). Its absolute moment of order \(k\) is finite only when \(k<\nu\): variance needs \(\nu>2\) and kurtosis needs \(\nu>4\). Predict the implications of fitting \(\nu\) near 3, while distinguishing the fitted model from the unknown true law.
The fitted values are \(\nu=3.29\), location .008 and scale .663. Scale is not standard deviation: the model SD is \(\text{scale}\sqrt{\nu/(\nu-2)}\) for \(\nu>2\). With \(\nu<4\), the fitted model has no finite fourth moment; this is not proof that Apple’s true fourth moment is infinite. The in-sample log-likelihoods are −5,731.9 (t) and −6,113.7 (Normal), a gain of 382 for one extra parameter, which lowers AIC. Tesla’s 253 observed 2020 returns fit \(\nu=3.31\) and give Shapiro–Wilk \(p=1.8\times10^{-6}\) under its nominal sampling assumptions. These descriptions need stability, dependence and out-of-sample checks. SciPy defines the t family and its location/scale parameters.
yfinance downloads market prices over the network in Colab or local Python. The archived course CSVs remain the reproducible browser examples.
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 =", stats.t.fit(zz)[0], " 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 pairs empirical and reference quantiles. A straight line indicates compatible location/scale and shape over the visible range; heavier sample tails bend away at the ends. stats.probplot selects the reference with dist and passes t degrees of freedom through sparams. Predict the Normal and fitted-t panels before running.
The Normal comparison shows departures at both ends. The fitted t comparison is straighter through much of the sample, but its most extreme right-tail quantiles exceed the observations. A symmetric t does not capture every asymmetry; the sample skew is −0.28. The plot locates departures rather than accepting a law. The following tests must respect parameter estimation and dependence.
Testing for normality: Shapiro–Wilk and Kolmogorov–Smirnov
State the null distribution, statistic and sampling assumptions before interpreting a test. Shapiro–Wilk compares ordered observations with Normal order statistics; KS uses the largest ECDF–reference CDF gap. kstest(returns, "norm") compares raw returns to \(N(0,1)\), whose scale is unrelated to the roughly .02 return SD. Its \(D\approx.47\) is not a pure shape diagnostic.
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 scale mismatch would also produce a large KS gap for perfectly Normal returns with SD .02. Studentising, or supplying fitted location and scale with args=(mean, sd), removes that mismatch. Because those parameters come from the same data, the usual specified-reference KS p-value no longer has its standard calibration.
The descriptive fitted-t distance .015 is smaller than the Normal distance .075. Standardisation and fitting used these same returns, while serial dependence adds another complication. Do not interpret the printed nominal p-values as calibrated evidence of acceptance. Monte Carlo calibration must repeat fitting in every simulated sample; the simulated dependence structure must also be suitable. SciPy describes this refitting procedure.
A valid p-value conditions on \(H_0\) and the test assumptions; it is not \(P(H_0\mid\text{data})\), which requires a prior (§5). Twenty independent true-null tests at 5% produce one false rejection on average (§3.3). Non-rejection does not establish the null. Here the fitted-t nominal p-value is uncalibrated, so it cannot establish non-rejection at its reported level.
A cheap bimodality detector
The bimodality coefficient \(b=(g_1^2+1)/(g_2+3)\) combines sample skewness \(g_1\) and excess kurtosis \(g_2\). The heuristic screen \(b>5/9\) uses the uniform distribution’s value as a reference. Separated mixture components can lower excess kurtosis and increase \(b\), but skewed unimodal variables can also flag. Predict the mixture’s value from \(g_1=.381\) and \(g_2=-1.214\), then compare Apple.
Mixture \(b=.641\) crosses the screen; Apple \(b=.120\) does not. Use a flag to inspect KDEs or mixture models, and heavy-tail evidence to consider a t or an EVT model (§3.5). The coefficient has no p-value; neither result establishes the number of modes or the correct model.
Under i.i.d. sampling, ddof=1 makes sample variance unbiased; standard deviation need not be unbiased. When sampling without replacement from a finite population, its variance target uses denominator \(N-1\). NumPy defaults to ddof=0. The ten-number sample gives SDs 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 under i.i.d. sampling the staircase is within 0.02 of the truth with probability at least 93.6 %; applying that coverage to dependent Apple returns requires further justification. At \(n = 100\) the bound exceeds 1 and says nothing.
Within the fitted t model, moments of order \(k\ge\nu\) are not finite, so \(\nu\approx3.3\) implies no finite fourth moment. This does not prove the true population moment structure. Apple’s observed shares are 76.9% within ±1 sample SD (Normal: 68.3%) and 98.5% within ±3 (Normal: 99.7%); its Q–Q tails depart from the Normal line. The lecture’s archived ticker fits report df from 3 to 6, rather than a universal law for all stocks.
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. Studentising or passing fitted location/scale corrects the scale mismatch, but changes the null calibration. The fitted-t distance is \(D = 0.015\); its nominal printed \(p = 0.27\) is not a calibrated goodness-of-fit probability because parameters were fitted on the same data.
Sampling Uncertainty and Comparisons
| Part | Question |
|---|---|
| 3.2.1 Sampling uncertainty and dependence | Why does the estimate change across samples? |
| 3.2.2 Uncertainty in repeated business measurements | Which observations should we resample together? |
Sampling uncertainty and dependence
Question: Why does the estimate change across samples?
Result: Explain a standard error, confidence interval and dependence assumption.
Sampling variation, bias and standard error
\[\operatorname{Bias}(\widehat\theta)=E[\widehat\theta]-\theta,\qquad \operatorname{MSE}(\widehat\theta)=\operatorname{Var}(\widehat\theta)+\operatorname{Bias}(\widehat\theta)^2.\]
Under i.i.d. sampling, \(\widehat{SE}(\bar X)=s/\sqrt n\). More customer rows within a branch or day are not necessarily independent observations; increasing \(n\) does not remove bias.
Confidence intervals: the meaning of 95%
For an i.i.d. Normal sample with unknown variance, a mean interval is
\[\bar x\ \pm\ t_{n-1,\,0.975}\frac{s}{\sqrt n}.\]
A mean interval is not a prediction interval for an individual customer, delivery or future period. Coverage requires the sampling assumptions; larger samples do not repair selection bias or leakage.
Bootstrap: approximate an estimator’s sampling distribution
Recompute \(\widehat\theta_b^*\) for each of \(B\) resamples; its 2.5th and 97.5th percentiles define the percentile interval. Longer blocks preserve longer local runs but provide fewer distinct blocks.
Assumptions: Enough data and a sufficiently stable process. More resamples reduce simulation noise; they do not create new historical information. BCa adds bias/acceleration adjustments and is developed in the Book.
Why dependence changes a standard error
For a stationary series with autocovariance \(\gamma_k\),
\[\operatorname{Var}(\bar X)=\frac{1}{n}\left[\gamma_0+2\sum_{k=1}^{n-1}(1-k/n)\gamma_k\right].\]
HAC means heteroskedasticity and autocorrelation consistent covariance estimation. In \(r_t=\mu+\varepsilon_t\), the OLS intercept estimates the mean measured outcome; HAC changes its uncertainty. For dependent daily operating measurements or monthly returns, inspect lag/block sensitivity. Chapter 4 develops OLS regression.
Uncertainty in repeated business measurements
Question: Which observations should we resample together?
Result: A block-bootstrap interval checked against block length and HAC uncertainty.
Builds on Sampling uncertainty and dependence.
Business application: operational uncertainty
Estimate mean service time, distinguish SD from SE and resample adjacent operating days. Choose the sampling unit and dependence assumptions before interpreting an interval.
Separate customer variation from uncertainty about the mean
Eight independent service times are observed. Does a 95% mean interval contain 95% of future service times?
Fixed-seed Normal simulation: true mean 20 min, SD 4 min.
SD describes individual service times; SE describes the sample mean across repeated samples. Under iid Normal sampling, 95% refers to repeated interval coverage of the fixed population mean, not individual customer coverage.
Resample consecutive operating days together
Six daily waiting-time deviations are −2, −1, 0, 1, 2 and 1 minutes. How does a two-day block preserve local order?
Invented teaching data.
Blocks preserve adjacency within a sampled run; block boundaries can jump and circular blocks wrap. This tiny fixture demonstrates mechanics, not a calibrated operational interval. Stable time behaviour and enough history remain assumptions; HAC is a different covariance estimate.
The unit of evidence is a realised period
Finance transfer example: The following existing market-data case applies the same general methods. Its sample and limitations remain explicit.
A return estimate needs its frequency, cost basis and sampling unit. Load the 34 realised monthly net returns from Chapter 2. The period end is the outcome date; the separate date column is the signal/rebalance date. Never count a pending next-month outcome as zero.
Resample adjacent months
An ordinary bootstrap samples individual observations independently. Trading returns and IC can have serial dependence, especially with overlapping horizons. A circular block bootstrap instead samples contiguous runs and wraps the last observation back to the first. The function makes one matrix of indices; no table-processing loop is hidden in it.
Each replicate has 34 month indices. Length 1 gives iid draws; length 3 retains three-month runs, wrapping at the last month. Independently sampled blocks assume a sufficiently stable process and cannot recreate unseen regimes or repair selection.
A confidence interval and a sensitivity check
net.to_numpy() copies return values to a NumPy array so a two-dimensional index matrix can form resamples; the ledger retains the month labels.
Estimate uncertainty about the mean monthly net return with a fixed seed and 2,000 replicates. Compare independent resampling with three-month blocks, then change the block length to six. Neither method is required to give a wider interval on this particular sample.
The block-3 interval is narrower here. These 34-month intervals omit strategy search, sample selection and unverified execution assumptions.
An intercept-only model with HAC uncertainty
sm names statsmodels.api; OLS means ordinary least squares. np.ones supplies the single intercept column. Chapter 4 develops regression models.
Reading a fitted constant is enough here: in a regression on a column of ones, the coefficient is the sample mean. HAC changes its estimated uncertainty to allow heteroskedasticity and autocorrelation. A three-lag choice is a teaching setting to challenge, not a validated optimum.
With 34 observations, HAC is approximate; |t| > 2 is not an automatic acceptance rule.
Sample contiguous time blocks. Time blocks retain local ordering within sampled runs, subject to stationarity assumptions.
The estimated covariance and uncertainty. HAC alters the standard error, allowing specified heteroskedasticity/autocorrelation; small-sample limitations remain.
Worked examples: Bootstrap, tests and rolling diagnostics
Scope of these examples. The ordinary bootstrap and permutation examples below illustrate independent sampling or exchangeability. Their nominal intervals and p-values do not automatically apply to serially dependent returns, overlapping holding periods or nonstationary regimes. BCa does not fix missing population moments or make tail-moment inference reliable; the rolling tests are descriptive alarms requiring prospective evaluation, not validated entry signals.
A mean, median, kurtosis, Sharpe ratio, quantile or maximum drawdown changes across samples. Bootstrap uncertainty must match the statistic and sampling process. Heavy-tail moments, dependence and path statistics require additional care; a fitted \(t_{3.3}\) lacks the moments needed for ordinary kurtosis inference. The worked tests below address a mean, a year-to-year difference, marginal distributions and rolling diagnostics.
The same sample, reloaded
Reload the archived AAPL log returns and select 2017: 251 daily observations. The refitted Student-t diagnostics use the full catalogue sample.
2017 mean daily log return: 0.151%; sample SD: 1.109%. Those sample summaries motivate the mean test; they do not validate its independence or stationarity assumptions.
The bootstrap idea
For an i.i.d. sample, the empirical distribution \(\hat F_n\) assigns mass \(1/n\) to each observed value. Draw \(n\) indices with replacement and recompute \(T\); across \(B\) repeats, the SD of \(T^*\) estimates its SE and the percentile endpoints estimate an interval. Glivenko–Cantelli supports empirical-CDF convergence, but that alone does not establish bootstrap validity for every statistic, tail moment or dependent process.
Reinvent the bootstrap SE, then check it against a case you know
The \(B\times n\) index matrix produces \(B\) resamples. Compute the statistic along axis 1. First use 120 illustrative lognormal draws (seed 11) for the median; then 251 AAPL 2017 daily log returns (seed 5) for the mean. The iid reference SE for that mean is \(0.01109/\sqrt{251}\approx0.0007\).
Median = 1.659; bootstrap SE = 0.198; estimated bootstrap bias = +0.053. For the AAPL mean, bootstrap and iid formula SE both round to 0.00070 (ratio 1.000). Agreement checks this calculation; it does not prove validity for other statistics or return dependence.
Percentile interval for Apple’s 2017 mean return
The percentile method reads the 2.5% and 97.5% quantiles of the bootstrap means. For this sample, compare those endpoints with \(\bar x\pm1.96s/\sqrt n\) and zero.
Percentile interval \([0.00014,0.00294]\); Normal approximation \([0.00014,0.00288]\). Both barely exclude zero under their sampling approximations. A corresponding inverted test and interval agree when built from the same model and method; percentile bootstrap and Student-t results are only approximately related. Symmetry is not a proof of correct coverage, especially with tail-sensitive statistics or dependence.
When the percentile interval goes wrong: BCa for the kurtosis of 2017
BCa adjusts the probability levels used to read the bootstrap distribution. Bias correction \(z_0=\Phi^{-1}(\Pr_*(\hat\theta^*<\hat\theta))\) measures median displacement; acceleration \(a\) comes from the standardised third moment of leave-one-out jackknife estimates. These adjust for bias/skewness under regularity conditions and require \(n\) jackknife evaluations; they do not guarantee reliable tail-moment inference.
\[\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 2017 sample excess kurtosis. The fraction of resamples below its observed value determines the sign of \(z_0\); check that fraction before interpreting the correction.
Sample excess kurtosis = 4.580; bootstrap SE = 1.55; 59.35% of replicates lie below the sample value. BCa shifts the percentile interval upward in this draw. It cannot restore nonexistent population moments or validate dependent-return inference. See SciPy bootstrap limitations.
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?
For the pre-specified \(H_0:\mu=0\), \(\hat t=\bar x/(s/\sqrt n)\). It has Student-t\(_{n-1}\) distribution for i.i.d. Normal data. Large-sample approximations need appropriate dependence and moment conditions; \(n\) in the hundreds alone is insufficient. Choose one- or two-sided inference before viewing the sign.
\(\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.
Observed \(\hat t=2.16\). Normal one-sided \(p=0.0154\); Student-t\(_{250}\) two-sided \(p=0.0318\) and greater-sided \(p=0.0159\). A pre-specified one-sided test counts only the chosen tail. These are nominal iid calculations on the 2017 sample; dependence and retrospective year/sidedness selection need separate treatment. They describe that sample and do not forecast 2018.
Permutation test: was 2017 a different year from 2018?
Under equal distributions and exchangeable independent units, relabel pooled values while holding the group sizes fixed. Recompute the mean difference for each shuffle; its two-sided tail share estimates the permutation p-value. No Normal model is needed, but equality of means alone, serial dependence or calendar-year assignment does not establish exchangeability.
AAPL 2017–2018 mean log-return gap = 0.179 percentage points/day; shuffle SE = 0.133 points; permutation \(p=0.1838\) and Welch \(p=0.183\). These nominal calculations fail to establish a difference under their assumptions; they do not establish equality. More shuffles reduce Monte Carlo error, not dependence or selection bias.
A rolling two-sample test as a regime detector
Compare the last ten observed returns with the ten ending twenty rows earlier. The denominator combines their sample variances; equal window sizes make this algebra coincide with the pooled-SE expression.
\[\hat t = \frac{\bar r_{\text{now}} - \bar r_{\text{lag}}}{\sqrt{\left(s^2_{\text{now}} + s^2_{\text{lag}}\right)/10}}.\]
shift(20) moves the comparison window; rolling(10) recomputes its summaries. A \(|\hat t|>2\) flag is a descriptive alarm, not proof of a changed distribution or a validated trading rule.
Twenty of 249 dates are flagged: December 2012 (1), March 2013 (2), May–July (17). The full-sample raw flag share is 5.2%; that denominator includes initial dates without a defined statistic. Short overlapping windows, dependence and repeated testing prevent interpreting this raw rate or clustering as a calibrated regime probability. A shorter window reacts faster but has noisier estimates.
A twenty-versus-twenty KS statistic also reacts to marginal volatility or tail-shape changes. Its usual p-value still requires suitable sampling independence; overlapping daily evaluations add multiplicity.
272 windows yield 33 raw flags at \(p<0.2\); the minimum nominal p-value is 0.034. With 20 observations per side, \(D\) changes in steps of \(1/20\) and power is limited. Treat the loose threshold as a descriptive diagnostic requiring prospective calibration; Chapter 6 develops models that pool more history.
Sampling from \(\hat F_n\) gives each observed value probability \(1/n\), independently and with replacement. Glivenko–Cantelli supports empirical-CDF convergence; validity for the chosen statistic requires further regularity and sampling assumptions. Here SEs were 0.00070 (mean), 0.198 (median) and 1.55 (sample kurtosis 4.58). BCa changed \([1.47,7.31]\) to \([2.23,9.15]\); it does not guarantee reliable tail-moment inference.
Under \(H_0:F_A=F_B\) with exchangeable units, equal-sized relabellings sample the null distribution without a Normal model. Calendar-year serial returns do not automatically satisfy that assumption. The archived gap is 0.179 percentage points/day (1.35 shuffle-SEs); nominal permutation \(p=0.184\) from 5,000 shuffles and Welch \(p=0.183\).
Research Design and Multiple Testing
| Part | Question |
|---|---|
| 3.3.1 Hypotheses and repeated comparisons | How can a search produce an attractive false finding? |
| 3.3.2 Pre-specified questions and a research log | What was chosen before the result was visible? |
Hypotheses and repeated comparisons
Question: How can a search produce an attractive false finding?
Result: Distinguish a p-value, error rates and multiple-testing adjustments.
Hypotheses, p-values and decision errors
State \(H_0\) and \(H_1\) before testing. For example, \(H_0:\mu=0\) against \(H_1:\mu\ne0\) for a fixed strategy and return definition.
Type I rejects a true \(H_0\); Type II fails to reject a specified false \(H_0\). Power is \(1-\beta\) for the stated effect and design; failure to reject does not establish zero effect.
Multiple comparisons: family-wise error and false discoveries
Independent true nulls give FWER \(1-(1-\alpha)^m\): 64% at \(m=20\), \(\alpha=0.05\). Bonferroni uses \(\alpha/m\) and needs no independence. BH finds the largest \(k\) with \(p_{(k)}\le kq/m\) and rejects through \(k\), under suitable dependence. Record the full search; neither method repairs leakage.
Research design: paired evidence and a fixed protocol
Random assignment helps identify causal effects; market backtests are observational. Fix the target, effect, features, split, costs, metrics and stopping rule before seeing results. Adaptive reward maximisation and fixed-rule evaluation answer different questions; the Book develops power, A/B, sequential testing and Bandits.
Pre-specified questions and a research log
Question: What was chosen before the result was visible?
Result: A fixed comparison, complete search history and stopping rule.
Builds on Hypotheses and repeated comparisons.
Business application: process changes and product experiments
Compare the same branches, state the counterfactual limits of before/after data, and account for every product experiment tried. Record the outcome metric, comparison and stopping rule in advance.
Compare a process change on the same branches
The same four branches are measured before and after a process change. What quantity should the paired test analyse?
Invented teaching data.
Pair by branch first; a negative change means fewer waiting minutes. The paired t-test assumes independent branch differences and, for an exact small-sample test, Normal differences. Before/after matching alone does not rule out seasonality or establish causation; a randomised comparison addresses a different design question.
Account for all product experiments tried
Four product changes yield p-values .001, .02, .04 and .20. Does testing only the best result control the full search?
Invented teaching data.
Bonferroni controls the chance of any false rejection; BH controls the expected false-discovery proportion under its dependence conditions. Include all declared tests, retain valid base tests and pre-specify the metric and stopping rule. Peeking requires a sequential design; adaptive Bandits optimise reward rather than automatically supplying fixed-horizon A/B inference.
Fix the research question before reading the answer
Finance transfer example: The following existing market-data case applies the same general methods. Its sample and limitations remain explicit.
Fix the universe, timestamp contract, label horizon, primary metric, costs, model settings and evaluation dates before the result. Retain all candidates tried.
A holdout used repeatedly to revise features becomes development data. Randomisation assumptions from experiments do not transfer automatically to markets.
Two hundred searches, one attractive result
stats names scipy.stats, which supplies the one-sample test. multipletests applies family corrections; its first returned item contains the rejection flags.
This fixed-seed simulation isolates selection risk: 200 independent candidate signals, 150 observations each, ten with a non-zero mean. Run the same family once, then compare an unadjusted screen with Bonferroni and Benjamini–Hochberg. It is a simulation of repeated research, not a measured result from MyStrategy.
These counts describe one seeded draw. Valid base p-values, a complete family and suitable dependence are prerequisites; correction does not repair leakage.
Research log and stopping rule
Log each dated hypothesis, configuration, source hashes, cutoff and all metrics. Apply the stated continue/reject rule; label new exploration separately. A winner crossing a threshold is not a valid retrospective stopping rule.
The Book works through A/B power, peeking, alpha spending and Bandits. Their randomised/stationary simulation assumptions do not establish a trading allocation policy.
Development evidence. Repeated feedback makes the period part of the research process, even if it was originally called OOS.
No; valid inputs and tests are prerequisites. Multiplicity adjustment cannot repair data leakage or invalid base p-values.
Worked examples: Power, sequential testing and Bandits
Scope of these examples. The A/B examples use randomised independent units and the Bandit examples use stationary simulated rewards. A historical backtest does not inherit those conditions. Multiplicity adjustment requires valid base tests; adaptive search and arbitrary cross-feature dependence need additional care.
Experimental design fixes the collection and analysis protocol before outcomes arrive. The worked examples use independent randomised A/B units, a pre-specified sequential null process and stationary simulated Bandit rewards. A market backtest is observational; inspecting and tuning its running results creates an adaptive search, not a randomised trial. All simulations retain their fixed seeds.
Classical A/B set-up: a real effect that reads “not significant”
Randomise independent units, fix \(N\) per arm and the hypothesis in advance, then analyse the mean difference and SE. The Normal approximation rejects when \(|z|>1.96\). Here control outcomes are \(N(0.10,1)\) and treatment outcomes \(N(0.18,1)\): true lift 0.08, 800 units/arm, seed 41.
Seed 41 gives \(\hat\tau=+0.086\), \(SE=0.051\) and \(z=1.70\): this run does not reject. Planned power for lift 0.08 is 0.36, so non-rejection is common even when the effect is real. Power describes repeated experiments under a specified alternative; it does not determine the outcome of this one run.
How big must N be?
For equal arm sizes, plan the detectable absolute effect \(\tau\), significance \(\alpha\) and power \(1-\beta\). The approximation below uses a common variance \(\sigma^2\), \(z_\alpha=1.96\) and \(z_\beta=0.84\).
\[n \;\approx\; \frac{2\sigma^2\,(z_\alpha + z_\beta)^2}{\tau^2},\]
Halving \(\tau\) quadruples \(n\). For Bernoulli outcomes use the baseline variance \(p(1-p)\); at baseline \(p=0.10\), a one-percentage-point lift gives \(2(0.09)(2.8)^2/0.0001\).
The baseline-variance approximation gives 14,112, 3,528 and 565 units/arm. A different effect-size/variance approximation can change these estimates; NormalIndPower with arcsine Cohen \(h\) gives about 14,745 for the one-point lift. Pre-register \(\tau\), \(\alpha\), power and \(N\); a non-significant result alone cannot distinguish low power from a small effect.
- 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.
The stopping and monitoring rule must match the eventual test.
Sequential tests: the cost of peeking
Under the specified true null, \(z_k=\sum_{i\le k}y_i/\sqrt{k}\). Looking once at a fixed \(N\) and stopping at the first crossing among 1,000 looks are different procedures. The simulation records any crossing of \(|z_k|>1.96\) over 2,000 null experiments.
The raw peek-and-stop rate is 0.518 for this seed and horizon. Overlapping looks are correlated, but the fixed-time 1.96 boundary does not control their joint crossing probability. Here \(z_k\) is a standardised running sum, not an ordinary random walk; the law of the iterated logarithm explains eventual crossings under this iid Normal null. A reported fixed-time \(p=0.04\) does not calibrate the whole adaptive stopping rule. Rolling tests likewise need a declared family and dependence-aware calibration.
The fix: spend your alpha
Fix the interim schedule and control \(\alpha\) for the whole crossing event. Bonferroni allocates \(\alpha/K\) per look by the union bound; it needs no independence. Pocock uses a joint-distribution constant boundary (2.413 for these five equally spaced looks). O’Brien–Fleming uses stringent early boundaries that relax later; Lan–DeMets expresses the budget through an alpha-spending function. The code compares five looks at 200, 400, 600, 800 and 1,000 observations.
Continuous peeking (previous seed) gave 0.518; these planned five-look results use seed 43. Bonferroni is conservative because the union bound does not exploit overlapping events. Valid sequential designs can allow early efficacy or futility stopping, with design-specific power/sample-size costs. Always-valid p-values, for example from a mixture sequential probability ratio test, require their own model and stopping-valid construction; they are not ordinary p-values or simply alpha spending in a continuous limit. See Johari and colleagues.
Multi-armed bandits: earning while learning
A Bandit optimises accumulated reward while learning arm quality. \(\varepsilon\)-greedy explores a random arm with probability 0.10; otherwise it uses the largest estimated mean. UCB1 uses \(\hat\mu_k+\sqrt{2\ln(t)/n_k}\), adding optimism for rarely sampled arms. Thompson sampling draws each Bernoulli arm’s rate from its Beta posterior and picks the largest draw. Here stationary true rates are 0.08, 0.10 and 0.12, with 4,000 pulls and seed 51; always selecting arm 2 has expected reward 480.
These policy rewards and pull counts describe one seeded run, not a universal ranking. The printed “regret vs oracle 480” is expected-oracle reward minus realised reward; standard pseudo-regret uses \(\sum_k n_k(p_* - p_k)\) and separates allocation loss from reward noise. Thompson’s Beta updating is developed in Chapter 5; adaptive sampling does not provide an unbiased fixed-rule comparison.
When to use which, and what goes wrong
Choose a fixed-\(N\) A/B design for pre-planned effect estimation, a valid group-sequential design for monitored stopping, or a Bandit for adaptive reward. Protect the declared metric family; for 12 independent null KPIs at 5%, FWER is \(1-0.95^{12}\approx46\%\), not “almost always”. Interference/SUTVA violations may require cluster or switchback randomisation; novelty requires a justified observation horizon. CUPED uses pre-experiment covariates for variance reduction: benefit depends on predictive correlation, coverage and design, not a guaranteed 30–60% saving. See Deng and colleagues’ CUPED analysis.
For \(m\) valid true-null tests at \(\alpha=0.05\), expected false rejections are \(0.05m\), regardless of independence. Bonferroni controls the probability of any false rejection; BH controls \(E[V/\max(R,1)]\) under suitable dependence. The following repeats the complete 200-signal, 150-period simulation (10 real signals; seed 9).
In this draw: naive = 19 discoveries (10 real, 9 false); Bonferroni = 9 (all real), missing the real signal at \(p=0.00067\); BH = 10 (all real). The realised false-discovery fraction need not equal the expected FDR guarantee. Keep every candidate and any changes to metric, subset, sidedness or stopping time; valid base tests and the appropriate dependence conditions remain necessary.
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.
The fixed-time p-value does not calibrate the team’s whole peek-and-stop rule. In these null simulations, continuous monitoring gave 51.8% and five naive looks 13.9%. Pre-plan valid sequential boundaries (Pocock 2.413 gave 4.6% here), use a stopping-valid method under its assumptions, or analyse once at the fixed planned N.
Association and Prediction Evaluation
| Part | Question |
|---|---|
| 3.4.1 Covariance, correlation and dependence | Which kind of association does the statistic measure? |
| 3.4.2 Evaluate predictions within comparable cohorts | Does the score order outcomes within each date? |
Covariance, correlation and dependence
Question: Which kind of association does the statistic measure?
Result: Choose a linear or rank measure and state what zero does not imply.
Covariance, Pearson correlation and scale
\[s_{xy}=\frac{1}{n-1}\sum_i(x_i-\bar x)(y_i-\bar y),\qquad r_{xy}=\frac{s_{xy}}{s_xs_y}.\]
Pearson lies between \(-1\) and \(1\). Constant inputs make it undefined; outliers, pair selection and pooled periods can change the answer.
Explain: A correlation of 0.6 does not mean a 60% probability of a correct trade or a 60% return.
Rank association: Spearman and Kendall
Spearman is Pearson correlation of both rank vectors:
\[\rho_s=\operatorname{Corr}(\operatorname{rank}X,\operatorname{rank}Y).\]
Strictly increasing transformations preserve Spearman ranks. Without ties, \(\tau=(C-D)/\binom n2\); Kendall tau-b adjusts for ties. Both measures need the same valid pairs.
An IC is an association statistic applied to predictions and realised targets. This chapter uses per-date Spearman IC; state that definition whenever reporting IC.
Zero correlation, independence and an IC time series
For symmetric \(X\) and \(Y=X^2\), opposing slopes can cancel covariance. Population distance correlation detects broader dependence under finite-first-moment assumptions; sample coefficients do not prove independence, causation or tradable profit.
For \(T\) date-level IC values:
\[\overline{IC}=T^{-1}\sum_t IC_t,\qquad ICIR=\overline{IC}/s_{IC}.\]
This ICIR is not annualised. Dependence between dates affects uncertainty; inspect coverage, stability and a time-aware interval for mean IC.
Evaluate predictions within comparable cohorts
Question: Does the score order outcomes within each date?
Result: Per-date IC with maturity, coverage and stability reported.
Builds on Covariance, correlation and dependence.
Business application: demand scores and later sales
Inspect linear and rank association, then evaluate stores within each week. Keep predictions with pending sales; exclude them only from evaluation until the outcomes are observed.
Choose an association measure for the question
Do higher demand scores accompany higher unit sales? Would Pearson detect a symmetric U-shaped relationship?
Invented teaching data.
Covariance has score-times-units dimensions. Pearson measures linear association; Spearman compares ranks and Kendall compares pair order. A zero Pearson coefficient can coexist with complete nonlinear dependence. Association alone establishes neither causation nor forecasting value on new data.
Evaluate demand rankings within each week
Stores are ranked before weekly sales arrive. Should a pending sales outcome remove a store from the earlier prediction list?
Invented teaching data.
Keep all seven predictions, but score only the six observed pairs. Compute one rank correlation per week before averaging weeks; pooled correlation and unequal weighting answer other questions. Retain counts, nonconstant vectors and stability. Finance calls this dated rank association IC; ICIR is its mean divided by its dated sample SD, not a success probability.
Define IC on one investable cross-section
Finance transfer example: The following existing market-data case applies the same general methods. Its sample and limitations remain explicit.
Information coefficient (IC) here means the association between a score known at the signal date and its later realised outcome, computed across stocks at that date. Pearson measures linear association; Spearman applies Pearson to ranks. The central MyStrategy task is to obtain one IC per signal date, then study the resulting time series.
Pooled correlation mixes dates. Per-date IC asks whether scores order outcomes in the same cross-section; its dated average does not establish net profit.
Keep prediction eligibility separate from evaluation
Evaluate only prediction rows whose outcomes are mature and finite. Pending rows remain in the candidate/prediction table. For a sector-relative diagnostic, subtract each date-sector mean realised return after the outcome is known; that transformed outcome is evaluation data, never a feature available at the signal date.
Evaluation uses 306 predictions across 34 dates: nine stocks per date, three per sector. Sector-relative outcomes do not establish neutral holdings.
Compute Spearman IC and ICIR by date
Rank both the prediction and each target within the same valid date cross-section, using average ranks for ties. Correlate those ranks within the date. Use the dates as the observations when summarising stability; also retain coverage and investigate constant-score or constant-target dates rather than filling undefined correlations with zero.
Raw and sector-relative targets answer different questions. ICIR uses the dated sample standard deviation and remains unannualised; undefined dates must not be filled with zero.
Stability, redundancy and nonlinear dependence
Inspect annual IC and the full dated series before reporting the average. A favourable year can conceal a later reversal. Feature redundancy should be investigated within development data, and distance correlation can screen non-monotonic dependence; it has no sign and does not by itself show a useful ranking direction.
Inspect dated coverage and annual reversals. Use §3.2’s dependence-aware intervals for mean IC; compare these diagnostics with reconciled net-return evidence.
One IC per signal date. Compute one cross-sectional IC on each date, then analyse that dated series and its dependence.
No; exposures must be checked on actual holdings. An outcome adjustment changes the diagnostic target, not the portfolio accounting or exposures.
Worked examples: Association beyond a straight line
Scope of these examples. These full-sample association examples are descriptive. For research selection, preserve development boundaries and evaluate predictions using their original stock/date columns and horizons. Nonlinear dependence is not a ranking direction, causal effect or guarantee of net profit.
Choose the statistic for the question: Pearson summarises linear co-movement, Spearman/Kendall rank ordering, and distance correlation broader dependence. Population distance correlation characterises independence under finite-first-moment assumptions; its sample estimate needs inference. The worked examples compare these measures and then inspect joint tail events.
Can Pearson see a U-shape?
Pearson’s \(r = \operatorname{Cov}(X,Y)/(\sigma_X \sigma_Y)\) is the cosine between centred vectors. It summarises linear co-movement; computation does not require the true relationship to be linear. Finite variances matter for a population coefficient, and outliers affect a sample estimate. Here \(y = (x - 5)^2 + \varepsilon\) with \(x\) uniform on \([0, 10]\) creates a noisy U-shape. Predict Pearson before running.
\(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\)?
The noisy sample gives \(r = 0.040\), Spearman 0.025 and correlation of \(|x - 5|\) with \(y\) 0.943. Cancellation explains the small first two coefficients; it does not justify the printed shorthand “no relationship”. Folding changes the feature and reveals the U-shape. A non-monotone relationship can still have nonzero Pearson correlation; Pearson alone does not describe its shape.
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 noisy sine. In the archive, tomorrow’s SPY correlations are −0.05 with today’s SPY, −0.06 with S&P 500, −0.03 with Nasdaq and +0.02 with Hang Seng. Same-date US measures are more correlated: S&P/Nasdaq 0.71, S&P/today’s SPY 0.78 and Nasdaq/today’s SPY 0.88. These are descriptive sample coefficients, not proof of efficient markets or absence of predictability. Check each market’s timestamp and information availability before treating it as a forecast feature.
Ranks instead of values: Spearman and Kendall
Spearman’s \(\rho_S\) correlates ranks. A strictly increasing transformation preserves them; deterministic monotonicity without ties gives an extreme rank coefficient. Kendall’s \(\tau\) compares the \(\binom{n}{2}\) observation pairs; without ties its probabilistic interpretation is \(\tau = P(\text{concordant}) - P(\text{discordant})\). Average ranks and tau-b handle ties differently. Ranks limit an observation’s numerical influence but do not eliminate outlier effects. Predict the original log, parabola and cosine examples before running.
The log example gives Pearson 0.921 and Spearman/Kendall 1.000. The parabola gives 0.075/0.059/0.039 and the cosine 0.036/0.033/0.014. Their small rank coefficients reflect cancelling ordering, not independence. A threshold or plateau can create ties, so an increasing-looking relationship need not produce a perfect rank coefficient.
Distance correlation sees any dependence
Distance correlation (Székely, Rizzo and Bakirov, 2007) compares pairwise-distance structure. The population coefficient is zero exactly at independence when the relevant first moments are finite; a finite-sample value alone does not certify independence. For \(a_{ij} = |x_i - x_j|\) and \(b_{ij}\), subtract row and column means and restore the grand mean to get \(A_{ij}\). Averaged products of the two centred distance matrices supply distance covariance; normalisation gives a nonnegative, unsigned correlation. The browser example implements these matrices with SciPy.
dcor supplies distance correlation in Colab. This new illustrative fixture reproduces the original 10,000-point parabola; the earlier NumPy/pandas imports are folded.
Previously introduced setup · 1 line
Copy complete code includes this earlier setup and the visible code below. Run it in Colab or local Python.
Introduced earlier: Can Pearson see a U-shape?
# Previously introduced: Can Pearson see a U-shape?
import numpy as np, pandas as pd!pip install dcor --quiet
import dcor
# Original 10,000-point parabola (illustrative)
x = 0.1 * np.arange(0, 10000)
parabo = pd.DataFrame({"x": x, "y": (x - 500)**2})
dcor.distance_correlation(parabo["x"], parabo["y"]) # ≈ 0.49The implementation below stores \(n \times n\) distance matrices. The Colab package provides optimised implementations for supported inputs; runtime depends on sample size and environment.
The original simulated examples give distance correlation 0.966 for the log, 0.495 for the parabola and 0.366 for the cosine. It detects dependence that the near-zero linear/rank coefficients miss. Its size is not directly comparable with Pearson and it has no direction; inspect the shape rather than treating either coefficient as a fitted model or a profitability test.
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.
In this 505-row historical sample, lag1/target Pearson is −0.214 and distance correlation 0.158; ret5 gives −0.067 and 0.131. These values cannot identify whether dependence is essentially linear, whether a model is sufficient or whether reversal is tradable after costs. Serial dependence affects uncertainty, and selection among many features requires §3.3’s research protocol. Evaluate any resulting signal chronologically with original index values and information cutoffs.
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.
Full-sample Pearson is 0.465 and Spearman 0.459; their similarity does not establish linearity or absence of influential outliers. Lower-tail thresholds are −1.69% for S&P and −5.14% for Tesla. The joint fraction is 41/126 = 32.5% in the lower tail and 30/126 = 23.8% in the upper tail, versus 5% under independence. Pearson within market-tail days is 0.480, higher than 0.465 overall. Conditioning changes the sample and can change correlation; it does not guarantee a decrease. Report conditional co-exceedance with uncertainty and dated thresholds.
Distance correlation is 0.495 for the parabola; Pearson/Spearman/Kendall are 0.075/0.059/0.039. The small directional coefficients reflect cancelling slopes and ordering. Population distance correlation characterises independence under finite-first-moment assumptions; finite-sample values require an inferential check.
On 126 market-tail days, 32.5% are also Tesla-tail days, versus 5% under independence; the upper-tail share is 23.8%. Pearson within the worst days is 0.480 against 0.465 overall, so it rises in this corrected market-tail sample. Conditioning changes the paired sample and can restrict variance, but does not force correlation to fall; the reported conditional tail probability and denominator answer the diversification question directly.
Portfolio Losses and Tail Risk
| Part | Question |
|---|---|
| 3.5.1 Historical losses and path risk | Which loss horizon and denominator belong to the claim? |
| 3.5.2 EVT and block extremes | How does a block-extreme model differ from ordinary returns? |
| 3.5.3 Threshold excesses and daily tail risk | How do excesses above a threshold imply daily tail probabilities? |
| 3.5.4 Tail diagnostics and sensitivity | How sensitive is the extrapolation to sparse tail evidence? |
Historical losses and path risk
Question: Which loss horizon and denominator belong to the claim?
Result: Historical VaR, ES and drawdown on a reconciled account.
Loss quantiles and Expected Shortfall
Define loss as \(L=-R\) for a stated horizon. A large positive value is a large loss.
\[\operatorname{VaR}_p=F_L^{-1}(p),\qquad \operatorname{ES}_p=\frac{1}{1-p}\int_p^1 F_L^{-1}(u)\,du.\]
For a continuous loss distribution, ES equals \(E[L\mid L\ge\operatorname{VaR}_p]\). VaR describes a cutoff; ES describes severity beyond it. Neither is a maximum possible loss.
The classroom code uses an interpolated sample quantile and averages observations at/above it. With finite samples or ties this convention can differ from exact fractional-tail quantile integration; report the actual tail count.
Drawdown is a path statistic
For net portfolio value \(V_t\) including the opening value,
\[H_t=\max_{s\le t}V_s,\qquad D_t=V_t/H_t-1,\qquad D_{\min}=\min_t D_t.\]
Drawdown measures loss from a running peak. We report its negative signed minimum; positive maximum-drawdown magnitude would be \(-D_{\min}\).
Specify horizon, fees, opening NAV and signed versus positive drawdown convention.
Extreme-value theory: model the tail, state assumptions
Block maxima: under suitable conditions, maxima of large blocks are approximated by a generalised extreme-value (GEV) distribution. Apply to losses; worst single-day loss in a month is not the compounded monthly loss.
Peaks over threshold: sufficiently high exceedances \(L-u\mid L>u\) can be approximated by a generalised Pareto distribution (GPD). Threshold choice trades approximation bias against a small tail sample.
The tail-shape parameter \(\xi\) distinguishes heavy unbounded tails (\(\xi>0\)), an exponential limit (\(\xi=0\)) and a bounded upper tail (\(\xi<0\)). Dependence may require declustering or adjusted inference.
Tail fits need stability checks and uncertainty. GEV and POT/GPD fitting, sensitivity and Q–Q diagnostics follow below. Two tail months cannot validate an extreme-risk model.
Measure risk on the reconciled net account
State the return horizon before taking a quantile. Use the Chapter 2 ledger for monthly net portfolio losses. Daily stock returns, monthly portfolio returns and the worst single day inside a 20-day block are different random variables and require different labels.
Fees are already in net_return; do not subtract them twice. Do not use gross returns for a net-risk report. A distributional risk estimate is evidence for a risk review, not a guarantee that losses stop at that level.
Historical VaR and Expected Shortfall
Define loss as minus the simple net return, so positive numbers mean losses. Historical 95% VaR is the empirical 95th loss percentile. The classroom tail-mean ES averages losses at or above that threshold; record the quantile convention, ties and number of tail observations.
At 95%, only two of 34 months enter the tail mean. Changing the percentile exposes the small denominator; a fitted model adds assumptions rather than observed crises.
Maximum drawdown depends on the path
Quantiles ignore the order of returns; drawdown preserves it. Include the opening NAV before the first realised outcome so an initial loss is counted. A drawdown is a decline from the running peak, not a one-period loss quantile or a worst-day-in-block statistic.
The realised teaching path has signed minimum drawdown −0.1833. Review each leg, joint losses, exposure, concentration, borrowing and execution before a risk decision.
Tail co-exceedance and the limits of EVT
Average correlation does not describe how often two assets fall together in their respective tails. Use the original S&P 500/Tesla daily-return example and keep it separate from the monthly ledger. Conditioning on a narrow market-return slice can reduce Pearson mechanically; report the tail event and its denominator.
The archive’s S&P/Tesla headers are reversed; price-series checks identify the columns. The read-time rename below corrects this example while preserving the CSV.
These thresholds use the full sample for description. Prospective monitoring requires prior-only thresholds and uncertainty; block extremes retain a single-day outcome.
A sparse-sample tail summary requiring uncertainty and more evidence. Report the two-of-34 denominator; tail precision is limited by the available observations.
The number of market-tail dates. The conditional fraction counts asset-tail events among market-tail dates.
EVT and block extremes
Question: How does a block-extreme model differ from ordinary returns?
Result: A GEV fit with explicit blocks, tail shape and support.
Builds on Historical losses and path risk.
Why EVT belongs in a trading risk toolkit
The centre of a return distribution describes ordinary days. A risk review also asks how severe the rare losses could be, how frequently a threshold might be exceeded and whether joint stress defeats diversification.
Extreme-value theory (EVT) supplies limiting models for extremes and high-threshold excesses. It supports model-based extrapolation beyond observed quantiles, with assumptions and estimation uncertainty.
The Central Limit Theorem concerns suitably normalised sums or means. EVT concerns suitably normalised maxima or excesses. Averaging and taking a maximum answer different questions.
Check: A Normal-looking sample mean does not imply Normal one-day loss tails.
Two routes from the same daily losses
The worked example uses 20-day blocks; the diagram’s four-day blocks are illustrative. Block maxima and threshold excesses describe different samples and probabilities.
GEV: a distribution for block extremes
For a block maximum loss \(M\) and \(z=(m-\mu)/\sigma\),
\[G(m)=\exp\{-[1+\xi z]^{-1/\xi}\},\qquad 1+\xi z>0,\quad\sigma>0.\]
At \(\xi=0\), the limit is \(G(m)=\exp[-\exp(-z)]\) (Gumbel).
\(\mu\) is location, \(\sigma\) is scale, and \(\xi\) controls tail shape. Location is generally neither the mean nor the mode.
Python convention: scipy.stats.genextreme returns c = -xi. A negative c indicates a positive heavy-tail index.
Tail shape controls support and moment existence
\(\xi>0\): Fréchet; \(\xi=0\): Gumbel; \(\xi<0\): Weibull. These are fitted-family properties, not proof of the population’s tail.
For a positive heavy-tail index, the fitted model has a finite moment of order \(r\) only when \(r\xi<1\). Thus its variance requires \(\xi<1/2\).
An estimated index does not prove the population has or lacks a moment. Support, confidence intervals, dependence and stability all matter.
Predict the tail family
Self-check: genextreme.fit returns c = -0.35. Which interpretation follows from this fitted family?
- A bounded upper loss tail because c is negative
- A heavy unbounded upper loss tail because ξ = 0.35
- A Gaussian tail with standard deviation 0.35
- A 35% probability of losing money
Build daily losses and non-overlapping blocks
Predict: How many complete 20-day blocks come from 1,496 daily returns?
The largest daily log loss is 0.2563, a simple loss of about 22.6%. The same single-day outcome remains an extreme when selected from its 20-day block.
Fit GEV and inspect its support
The fit gives \(\xi\approx0.354\), location 0.0136 and scale 0.0070. Its 95th block-extreme percentile is 0.0506. This is the fitted worst-day threshold exceeded in 5% of 20-day blocks, rather than a compounded 20-day loss.
Starting values seed optimisation; compare alternative starts and likelihoods, support and stability.
Threshold excesses and daily tail risk
Question: How do excesses above a threshold imply daily tail probabilities?
Result: A POT/GPD calculation with the exceedance fraction retained.
Builds on EVT and block extremes.
POT: the generalised Pareto model
Set a high loss threshold \(u\). Conditional on \(L>u\), define the positive excess \(Y=L-u\).
\[H(y)=1-(1+\xi y/\beta)^{-1/\xi},\qquad y\ge0,\quad1+\xi y/\beta>0.\]
\(\beta>0\) is excess scale. At \(\xi=0\), \(H(y)=1-e^{-y/\beta}\). With \(\xi<0\), excesses end at \(-\beta/\xi\).
genpareto uses c = xi, unlike GEV. Fix excess location to zero with floc=0. The threshold \(u\) belongs outside the excess distribution.
Explore the tail shape before fitting
Try: Set ξ to −0.25, 0 and +0.50. Explain the endpoint and the change in the conditional 95th percentile while scale remains fixed.
Fit excesses above a 1.5% daily log loss
The 1.5% threshold gives 91 excesses from 1,496 days: 6.08%, with \(\xi=0.4363\) and scale 0.00643. The whole 1985–1990 fit is descriptive; forecasts require prior-only observations.
From conditional excesses to unconditional daily risk
Let \(p_u=N_u/n\). For \(q>1-p_u\) and \(\xi\ne0\),
\[\operatorname{VaR}_q=u+\frac{\beta}{\xi}\left[\left(\frac{p_u}{1-q}\right)^\xi-1\right].\]
For the fitted continuous GPD tail and \(\xi<1\),
\[\operatorname{ES}_q=\frac{\operatorname{VaR}_q+\beta-\xi u}{1-\xi}.\]
At \(\xi=0\), VaR becomes \(u+\beta\log[p_u/(1-q)]\), and ES becomes VaR \(+\beta\).
At \(\xi\ge1\), the fitted tail has infinite ES. The displayed finite estimates use the original fit with \(\xi<1\).
Compute daily VaR99 and ES99
The fitted daily log-loss VaR99 is 0.0327 and ES99 0.0577. VaR transforms to a simple loss of about 3.21%; transforming the ES average does not give simple-loss ES.
Tail diagnostics and sensitivity
Question: How sensitive is the extrapolation to sparse tail evidence?
Result: A risk statement supported by threshold, Q–Q and tail-shape diagnostics.
Builds on Threshold excesses and daily tail risk.
Threshold selection trades bias against variance
Self-check: Why should we compare several plausible POT thresholds?
- The highest threshold always gives the most precise estimate
- More excesses always mean a better tail approximation
- Threshold choice changes approximation bias, sample size and uncertainty
- We should select whichever threshold gives the smallest reported VaR
Sensitivity is evidence, not a confidence interval
Threshold differences are sensitivity evidence. Check mean-excess plots, parameter stability and dependence-aware uncertainty; the three estimates are not a calibrated confidence interval.
Q–Q diagnostics expose the most fragile observations
The 1987 crash remains the most fragile observation. Clustered exceedances may require declustering or conditional-volatility modelling; record the changed target and denominator.
Return levels extrapolate a probability, not a calendar promise
A 1-in-T-day return level \(x_T\) satisfies \(P(L>x_T)=1/T\) in the fitted stationary marginal model.
\[x_T=u+\frac{\beta}{\xi}\left[(p_uT)^\xi-1\right].\]
The 2,500-day level uses only 1,496 observations. Model estimates extrapolate; interpolated empirical quantiles stay within the sample extremes.
Check: “1-in-1,000” does not mean the next breach arrives exactly 1,000 days later. Clustering and changing regimes further weaken a simple waiting-time interpretation.
Hill estimation provides a second tail-shape diagnostic
For descending positive losses \(L_{(1)}\ge\cdots\ge L_{(k+1)}>0\),
\[\widehat\xi_k=\frac1k\sum_{i=1}^k\log\frac{L_{(i)}}{L_{(k+1)}}.\]
Hill assumes a heavy, regularly varying tail. Plot \(\widehat\xi_k\) against \(k\): a stable region can support a threshold choice, but a plateau is not proof.
Too few extremes give noise; too many admit the body. GEV, GPD and Hill use different finite samples and need not agree exactly.
Continue in the Book: Run the Hill and return-level cells, then compare their assumptions and samples with the two fits above.
Self-check: write a defensible tail-risk statement
Self-check: Which summary correctly reports the POT result?
- The strategy cannot lose more than 3.27% in a month
- A 1985–1990 Dow fit above 1.5% daily log loss uses 91 excesses and estimates daily log-loss VaR99 near 3.27%
- Ninety-one exceedances guarantee a precise future tail probability
- ES99 is a 99% confidence interval for mean return
Explain in your own words: What changes if the threshold rises? Which units differ from the monthly portfolio report? What additional evidence would you need for a prospective risk decision?
Scope of these examples. Returns below are daily log returns; their block minima remain single-day outcomes. They are not monthly simple portfolio returns or cumulative block losses. The GEV/GPD/Hill estimates use sparse, potentially dependent tails and do not justify a hard trading limit. Check threshold/block sensitivity, temporal stability and estimation uncertainty.
Under suitable finite-variance and dependence assumptions, central-limit results concern normalised sums or means. EVT instead concerns normalised maxima or high-threshold excesses under appropriate tail conditions. These historical examples compare original return tails, collect block extremes, fit and diagnose GEV/POT, then inspect Hill estimates and return levels. Neither limit theory establishes a hard loss bound from this finite sample.
Same mean, same sd — same tail?
Draw 4 309 Normal returns with Apple’s observed mean and standard deviation (seed 71). The draw parameters, rather than the realised simulation’s sample moments, match the archive. Predict the 0.1% quantile and counts beyond four observed standard deviations; matching two parameters does not determine tail shape.
The original 0.1% quantiles are about −10.5% for Apple and −5.8% for the Normal simulation. At 1% they are −5.8%/−4.6%, and the conditional tail-return means are −7.8%/−5.1%. Apple has 25 observations beyond four observed sd; the simulation has none (about 0.27 expected under its fitted Normal draw parameters). Dispersion alone does not determine tail losses; these are descriptive archive/simulation results.
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\) and location \(\mu\), scale \(\sigma>0\). The quantile function \(Q(p)\) supplies the fitted extreme quantile. For \(\xi > 0\) (Fréchet), a heavy upper tail behaves as \(1 - F(x) \sim x^{-1/\xi}\); moments of order \(r\) need \(r\xi<1\). \(\xi = 0\) is Gumbel; maxima of Normal, lognormal and exponential parents belong to that domain under appropriate normalisation. \(\xi < 0\) (Weibull) has upper endpoint \(\mu - \sigma/\xi\). Existence of a non-degenerate extreme-value limit is an assumption, not a property of every parent. A Student-\(t_\nu\) parent gives \(\xi = 1/\nu\) (so \(\nu=3.3\) implies about 0.30); the Dow fit is a different sample, not validation of that Apple model.
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.
For \(\mu = 0\) and \(\sigma = 1\), the endpoint \(\mu - \sigma/\xi\) is +3 when \(\xi=-1/3\) and −3 when \(\xi=1/3\). The latter is a lower support boundary; the upper tail remains unbounded. The Gumbel case has no finite endpoint.
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.
The first three block minima are −1.08%, −1.01% and −0.60%, and the median is −1.65%. Every complete block contributes one daily extreme regardless of how many bad days it contains. Non-overlap removes shared observations, not volatility dependence. Check block size, temporal stability and uncertainty.
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?
At 1.5% daily log loss, 91/1,496 days exceed the threshold. The original yearly counts are 1, 16, 29, 14, 10 and 21. These selected days retain stress clusters that one-maximum-per-block compression hides. Declustering changes the sample and event rate; an adjusted fit must state that changed target.
Warm-up: GEV on the block maxima of a Normal
The Normal is in the Gumbel domain of attraction: its suitably normalised block maxima converge to \(\xi = 0\). The original seed-1 demonstration fits 100 maxima from 100 draws each. Finite-block and finite-sample estimates need not equal zero; a materially different estimate calls for checking uncertainty and fitting choices, not automatically declaring a bug.
\(\xi = -0.058\), with location 2.33 and scale 0.40; the sample mean of block maxima is 2.54. That fitted location is neither the sample mean nor a universal “worst case”. This simulation checks the procedure on a known limiting domain; it does not validate tail estimates for market data.
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\) and \(\sigma = 0.0070\) define a fitted Fréchet upper tail. Its finite-moment condition is \(r\xi<1\): this fitted model has finite variance but no third moment; an estimate does not prove the population’s moments. Its lower support starts near −0.0063; negative block loss means all days in that block rose. The sparse histogram reaches 25.6% daily log loss. Comparing this market fit with Normal block maxima requires their different units, block lengths and sampling uncertainty; neither family gives a deterministic future loss 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 reads SciPy’s \(c\) as \(\xi\). The returned c = −0.354 means xi = -c, a heavy unbounded upper tail. The wrong sign instead creates an asserted endpoint near 0.0335, which excludes the observed 0.2563 loss. Check fitted support against every retained observation; “bounded” is a property of the selected fitted family, not approval of a hard loss limit.
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 25.6% crash block lies far above the fitted comparison line at the upper right. Most other points are closer to it; that visual agreement does not validate the far tail. With 74 block observations, the most extreme event has substantial leverage on extrapolation. Report that mismatch, parameter/support stability and dependence-aware uncertainty rather than treating the fit as a guaranteed stress boundary.
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 scores fall below 0.10. The crash score is about 0.0007; the reciprocal is about 1,464 blocks from the unrounded fit. January 1988 and October 1989 score about 0.021. Converting such marginal probabilities into expected recurrence requires stationarity and a specified event process; it does not predict the next event date. Clustering after the crash and estimating the fit on these same observations limit that interpretation.
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.
The 91-excess GPD fit gives \(\xi = 0.436\), against 0.354 from 74 GEV block extremes. Its fitted mean excess \(\beta/(1 - \xi)\) is 1.14%, versus observed 1.19%. Daily log-loss VaR99 is 3.27%, versus empirical 3.05%. At 0.1%, an empirical interpolated quantile is still defined between extreme order statistics; it is sparse and cannot extrapolate beyond the sample maximum. Agreement across fits is a sensitivity check, not calibrated certainty.
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\) is 0.468/0.423/0.441 at \(k\) = 25/50/100 and 0.629 at \(k = 200\), whose threshold is about 0.82%. The last estimate may include more of the body; this change alone does not prove the power-law assumption failed. Compare GEV 0.35, GPD 0.44 and the Hill region under their different samples. Their spread is sensitivity evidence; calibrated uncertainty must handle dependence and threshold selection.
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\), GPD/empirical log-loss levels are 4.86%/4.71%. At \(T = 1000\), they are 8.87%/7.77%, with the empirical interpolation between the second and third worst observations. At \(T = 2500\), they are about 13.2%/15.3%, with fewer than one expected sample observation beyond that probability. The observed maximum is 25.6%; estimates remain uncertain and “return period” is a probability scale, not a scheduled breach.
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 fitted 5% block-extreme return threshold is −5.06%, breached by 3/74 blocks (4.1%). At 1% it is about −9.5%; at 0.1% about −22.3%; the crash breaches both. With the same location/scale and \(\xi = 0\), the 5% threshold is −3.45%, breached by seven blocks. These are fitted worst-single-day quantiles inside 20-day blocks; daily marginal VaR and compounded block loss answer different questions.
Use the fitted extreme quantile as a risk-review diagnostic with support checks, alternative blocks/thresholds, temporal stability and dependence-aware uncertainty. This window ends in 1990 and does not calibrate a present portfolio or justify a hard stop. GPD and Hill agreement is a sensitivity check, not a confidence interval or proof.
In August 2007 quantitative equity strategies suffered substantial losses. Goldman Sachs’ 13 August statement reported overlapping trades, greater volatility/correlation, reduced risk/leverage and a US$3 billion investment in Global Equity Opportunities. The “25-standard-deviation moves” remark is documented in Dowd et al., How Unlucky is 25-Sigma?; Khandani and Lo’s research investigates deleveraging and withdrawal of marketmaking capital. These sources do not turn this Dow tail fit into a model of those portfolios.
A standard Normal one-sided 25σ tail probability is about \(3.1 \times 10^{-138}\) (two-sided \(6.1 \times 10^{-138}\)). A large sigma count tests a distributional model and its estimated scale; it does not identify the probability under a heavy-tail model. The Dow crash’s fitted block score near 0.0007 concerns a different dataset and event definition. Check the relevant portfolio, horizon, leverage and dependence before making a risk claim.
Lesson: standard deviations can describe scale without Gaussianity. Turning a sigma count into a Gaussian tail probability adds a distributional assumption that must be checked; tail fits add their own assumptions and uncertainty.
\(\xi > 0\) gives the fitted Fréchet heavy, unbounded upper tail; a moment of order \(r\) needs \(r\xi<1\). SciPy GEV uses c = −ξ, so c = −0.35 means \(\xi = +0.35\). The opposite reading creates a bounded family inconsistent with the 25.6% observation. genpareto uses c = +ξ.
The fitted model assigns a 5% probability that a 20-day block contains a single-day return below −5.1%; 3/74 blocks (4.1%) breached it in this sample. It is neither compounded monthly loss nor a hard stop. GEV/GPD/Hill estimates are about 0.35/0.44/0.42–0.47 under different samples; their spread is sensitivity, not a calibrated interval.
Chapter Wrap-up
You can now audit a feature distribution, evaluate a score on each signal date, attach uncertainty to realised evidence and describe a portfolio loss with the correct horizon. The running case preserves 639 candidates while evaluating only 306 of 315 predictions on mature outcomes; the account contains 34 realised monthly returns. Those denominators belong beside every conclusion.
The report should distinguish bookkeeping correctness, predictive association, net economic outcomes and risk. A positive IC does not establish a profitable implementation; a reconciled NAV does not establish significance; a tail fit does not establish a hard loss limit. The small current-survivor teaching universe and assumed clocks remain limitations throughout.
Chapter 4 develops predictive modelling, Chapter 6 develops time dependence, and Chapter 7 revisits strategy selection. The Slides and guided notebook use the same handoff and calculations.
Method references
Time-series bootstrap variants: the browser example implements fixed-length circular blocks with NumPy.
statsmodels robust covariance: HAC covariance options and inferential conventions.
statsmodels multiple testing: Bonferroni and FDR procedures.
SciPy goodness-of-fit calibration: refitting estimated parameters in Monte Carlo tests.
SciPy GEV and GPD: support, quantile functions and opposite shape conventions.