import numpy as np
import pandas as pd
import plotly.graph_objects as go
from plotly.subplots import make_subplots
from scipy import stats
import fe1tools as fe
fe.setup() # plot theme, display options, random seed
BLUE, ORANGE, AQUA = fe.PALETTE["blue"], fe.PALETTE["orange"], fe.PALETTE["aqua"]
# Problem 1: the CRSP value-weighted equity premium, monthly since July 1926.
ff = fe.load_ff("monthly3", end="2025-06-30")
premium = ff["Mkt_RF"].dropna()
# Problem 2: S&P 500 and NASDAQ Composite, monthly since February 1971.
levels = fe.load_indices()
idx = levels.resample("ME").last().pct_change().dropna()
gap = idx["NASDAQ"] - idx["SP500"] # the paired difference3 Week 2 Online Learning
Learning objectives. After this week you will be able to
- take two ordinary finance questions — is the equity premium positive, and do two indices earn the same return — turn them into hypotheses about a population mean, and compute the statistic that answers each;
- point to exactly where normality enters the derivation of the \(t_{n-1}\) distribution, and so say what is being assumed every time someone consults a \(t\)-table;
- explain why an exact finite-sample result is the strongest thing available when its assumptions hold and worthless when they do not, and give examples of questions for which no exact result exists at all;
- state the law of large numbers and the Lindeberg–Lévy central limit theorem, explain what each contributes, and distinguish the infeasible from the feasible version of the resulting test statistic;
- work out by simulation, rather than by folklore about thirty observations, whether an asymptotic approximation is good enough at the sample size you actually have.
Every chart in this lesson is interactive: hover to read values, drag to zoom, double-click to reset.
3.1 Two questions, and why the arithmetic is the easy part
Week 0 spent a long time on a single idea: the numbers on a performance report are not facts about the world, they are estimates, and every estimate arrives with a distribution attached to it. We built the standard error to measure how wide that distribution is, and we watched an estimate of the average return refuse to settle down even after twenty years of data.
This week we take that machinery and point it at two questions that someone actually has to answer at work. Both are easy to state. Both take about one line of code. And in both cases the calculation turns out not to be where the difficulty lies. The difficulty is in the very last step — the one where a number gets converted into a statement about probability — and that step is usually performed by looking up a table without asking what the table assumes.
So it is worth being explicit about what we are trying to do. We have a sample mean. We want to know whether the true mean behind it could plausibly be zero. To answer that we need to know how much a sample mean bounces around from sample to sample when the true mean really is zero: if our observed value would be unremarkable in that world, we have learned nothing, and if it would be extraordinary, we have learned something. The whole exercise therefore depends on knowing one object — the distribution of our statistic under the null hypothesis. Everything in this chapter is about where that object comes from.
Problem 1. Is the equity premium positive?
Everything in finance that discounts a risky cash flow assumes that bearing equity risk is compensated. The claim is not that equities went up in the past; it is that their expected return exceeds the risk-free rate. Written as a hypothesis about the population mean of the excess return \(R_t^e\),
\[ H_0 : \mu = 0 \qquad\text{against}\qquad H_1 : \mu > 0 . \]
We use the CRSP value-weighted market return minus the one-month Treasury bill, monthly, from July 1926 — the longest honest sample available.
Problem 2. Do the S&P 500 and the NASDAQ earn the same return?
An investor choosing between a broad-market tracker and a technology-heavy one wants to know whether the historical gap is a real difference in expected return or the residue of one long bull market. With \(R_t^{N}\) and \(R_t^{S}\) the monthly returns on the NASDAQ Composite and the S&P 500,
\[ H_0 : \mu_N - \mu_S = 0 \qquad\text{against}\qquad H_1 : \mu_N - \mu_S \neq 0 . \]
Both are price indices, so neither reinvests dividends. That understates the return to holding either portfolio, but it is the same omission on both sides, which is what a comparison requires.
fig = make_subplots(rows=1, cols=2, horizontal_spacing=0.10,
subplot_titles=("Cumulative equity premium, 1926–2025",
"S&P 500 vs NASDAQ Composite, 1971–2025"))
wealth = (1 + premium).cumprod()
fig.add_trace(go.Scatter( # panel title names it; no legend
x=wealth.index.to_pydatetime(), y=wealth.to_numpy(),
line=dict(color=BLUE, width=1.6), showlegend=False,
hovertemplate="%{x|%b %Y}<br>£%{y:,.0f}<extra></extra>"), row=1, col=1)
for col_name, hue in [("SP500", BLUE), ("NASDAQ", ORANGE)]:
w = (1 + idx[col_name]).cumprod()
fig.add_trace(go.Scatter(
x=w.index.to_pydatetime(), y=w.to_numpy(),
line=dict(color=hue, width=1.6), name=col_name,
hovertemplate="%{x|%b %Y}<br>%{y:,.1f}×<extra></extra>"), row=1, col=2)
# dtick=1 on a log axis labels the powers of ten only; plotly's default adds
# minor ticks that render as a confusing 1,2,5,1,2,5 sequence.
fig.update_yaxes(type="log", dtick=1, title_text="£ (log scale)", row=1, col=1)
fig.update_yaxes(type="log", dtick=1, title_text="growth of 1 (log scale)",
row=1, col=2)
fe.figure(fig, height=430)
fig.update_layout(margin=dict(t=88), legend=dict(y=1.14))
figBoth pictures look completely decisive, and it is worth understanding why that impression cannot be trusted. A cumulative wealth chart accumulates: each point is the product of everything that came before it, so a drift far too small to see in any individual month gets multiplied by itself a thousand times and turns into a line that climbs across the whole page. The left panel rises by a factor of several hundred, but it does so out of monthly returns whose average is 0.69% — a number you would struggle to distinguish from zero by looking at any single month, or indeed any single decade. Compounding is a magnifying glass, and it magnifies noise along with signal.
The same caution applies on the right. The NASDAQ line ends up well above the S&P 500 line, and the eye reads that gap as a settled fact. But the gap opens up in bursts — most visibly in the late 1990s — and a chart cannot tell us whether those bursts reflect a genuinely higher expected return or a handful of lucky years that happened to fall inside our sample.
What we need is not a better picture but a measurement of how much these averages would have moved around had history run differently.
def describe(x, label):
n = len(x)
se = x.std(ddof=1) / np.sqrt(n)
return {
"series": label,
"n (months)": n,
"from": x.index[0].strftime("%b %Y"),
"mean (ann. %)": x.mean() * 12 * 100,
"vol (ann. %)": x.std(ddof=1) * np.sqrt(12) * 100,
"s.e. of mean (ann. %)": se * 12 * 100,
"t-stat": x.mean() / se,
}
pd.DataFrame([
describe(premium, "Equity premium (Mkt−RF)"),
describe(gap, "NASDAQ − S&P 500"),
]).set_index("series").round(3)| n (months) | from | mean (ann. %) | vol (ann. %) | s.e. of mean (ann. %) | t-stat | |
|---|---|---|---|---|---|---|
| series | ||||||
| Equity premium (Mkt−RF) | 1188 | Jul 1926 | 8.253 | 18.428 | 1.852 | 4.456 |
| NASDAQ − S&P 500 | 652 | Mar 1971 | 3.126 | 10.593 | 1.437 | 2.175 |
Let us walk through the first row slowly, because every later section refers back to it. The equity premium has averaged 0.6877% per month over 1,188 months. Monthly returns have a standard deviation of 5.32%, so a single month tells us very little. But we are not asking about a single month — we are asking about an average of 1,188 of them, and week 0 showed that the standard error of an average is
\[ \text{s.e.}(\bar{X}) \;=\; \frac{S}{\sqrt{n}} . \tag{3.1}\]
Putting the numbers in: \(5.32\%\) divided by \(\sqrt{1188} \approx 34.5\) gives a standard error of 0.1543% per month.
The \(\sqrt{n}\) in eq. 3.1 is doing the work, and it is worth remembering where it comes from rather than treating it as a formula. Averaging 1,188 independent observations does not cut the noise by a factor of 1,188, because the individual errors reinforce each other about as often as they cancel. What adds cleanly is not standard deviations but variances: the variance of a sum of independent terms is the sum of their variances. So the variance of the average falls by a factor of \(n\), and the standard deviation — being the square root of the variance — falls only by \(\sqrt{n}\). Quadrupling the sample halves the standard error, which is a disappointing rate of return on data collection and the reason week 0 concluded that expected returns are so hard to measure.
Now compare the estimate against its own noise. The average is 0.6877% per month and the standard error is 0.1543%, so dividing one by the other puts the average 4.46 standard errors above zero. That is all a \(t\)-statistic is — a distance from the null, measured in units of its own uncertainty, which makes it comparable across problems with completely different scales. The same calculation for problem 2 gives 2.18.
And here we stop, because we have run out of things the arithmetic can tell us. We know the premium is 4.5 standard errors from zero. We do not yet know whether 4.5 standard errors is a lot. To say that, we would need to know how far from zero this statistic typically lands when the true mean is zero — is 4.5 the kind of value that turns up one time in twenty, or one time in a million? That is a question about the distribution of \(T\), and nothing we have computed so far answers it. A \(t\)-statistic of 2.18 is, at this stage, a number and not yet evidence.
Why we tested the difference rather than comparing two means. There are two ways to ask whether the NASDAQ beats the S&P 500. We could estimate each mean separately and compare them, or we could form the monthly difference \(R_t^N - R_t^S\) and ask whether its mean is zero. We did the second, and the choice matters enormously.
The reason is that the two indices are not independent. Their monthly returns have a correlation of 0.87, because both are dominated by the same thing: whatever the American stock market did that month. When the market falls 8%, both indices fall by something close to 8%, and neither of those large movements has anything to do with the question we are asking. Subtracting one series from the other cancels that shared component and leaves behind only the relative performance, which is the quantity of interest. The difference series has an annualised volatility of 10.6%, against roughly 15% and 21% for the two indices individually — most of the variation has been removed, and what remains is signal rather than market noise.
Notice that this decision was made before any statistical theory entered the room. It is a modelling judgement about which comparison answers the question, and no amount of distribution theory would have rescued us from getting it wrong. Section 3.4 puts a number on what it was worth.
3.2 What financial data actually look like
The standard route from a \(t\)-statistic to a \(p\)-value runs through a table of the \(t\) distribution, and that table was computed on the assumption that the observations are draws from a normal distribution. We are about to lean on that assumption quite heavily, so it is only sensible to look at the data first and see whether it is remotely true.
Two numbers summarise the shape of a distribution once its mean and variance are accounted for. Skewness measures asymmetry: it is zero when the distribution is symmetric about its mean, positive when the right tail is longer than the left, negative when the reverse. Excess kurtosis measures how much probability sits in the tails relative to a normal distribution: it is zero for a normal, and positive when extreme observations are more common than a bell curve would allow. The Jarque–Bera test combines the two into a single statistic and tests the joint hypothesis that both are zero, which is what normality would require.
rows = []
for x, label in [(premium, "Equity premium (Mkt−RF)"), (gap, "NASDAQ − S&P 500")]:
jb, jb_p = stats.jarque_bera(x)
rows.append({
"series": label,
"skewness": round(stats.skew(x), 3),
"excess kurtosis": round(stats.kurtosis(x), 2),
"Jarque–Bera": round(jb, 1),
# Formatted, not rounded: the premium's p-value underflows to exactly 0.
"p-value": f"{jb_p:.1e}" if jb_p > 0 else "<1e-300",
"worst month (%)": round(x.min() * 100, 2),
"best month (%)": round(x.max() * 100, 2),
})
pd.DataFrame(rows).set_index("series")| skewness | excess kurtosis | Jarque–Bera | p-value | worst month (%) | best month (%) | |
|---|---|---|---|---|---|---|
| series | ||||||
| Equity premium (Mkt−RF) | 0.154 | 7.34 | 2673.3 | <1e-300 | -28.74 | 38.81 |
| NASDAQ − S&P 500 | 0.339 | 7.11 | 1387.2 | 5.8e-302 | -14.89 | 21.21 |
Neither series comes anywhere close to passing. The excess kurtosis of 7.3 on the equity premium is the number to focus on, and the easiest way to feel its size is to look at a single observation.
The worst month in the sample is -28.7%. Given a monthly standard deviation of 5.32%, that observation lies 5.5 standard deviations below the mean. Now ask what a normal distribution has to say about an event that far into the tail. The normal density decays like \(e^{-z^2/2}\), so probability disappears extraordinarily fast as \(z\) grows: a three-standard-deviation month is a once-a-decade event, four standard deviations is once every few centuries, and five and a half is a probability of roughly \(1.6 \times 10^{-8}\). At twelve months a year, that works out to one occurrence every 5 million years.
It happened in September 1931. And it was not a freak: eight months in this sample lie beyond four standard deviations from the mean — four below and four above — where a normal distribution predicts essentially none at all in a century of monthly data. A model that assigns a probability of one-in-millions to things that occur several times a lifetime is not making a small error. It is wrong about the part of the distribution we care most about.
The Jarque–Bera statistic makes this formal and rejects normality for both series with a \(p\)-value so small it underflows to zero in double precision. There is no judgement call to make here, and no hope that a longer sample would settle things differently. Fat tails are the most reliably reproduced empirical fact in all of asset pricing — they show up in every market, at every frequency, in every period anyone has looked at — and Chapter 1 listed them as one of the four features that make financial econometrics its own subject rather than a branch of applied statistics.
fig = make_subplots(rows=1, cols=2, horizontal_spacing=0.11,
subplot_titles=("Equity premium vs the normal density",
"Normal quantile–quantile plot"))
# --- left: histogram against the matched normal ---
p = premium.to_numpy() * 100
fig.add_trace(go.Histogram(
x=p, histnorm="probability density", nbinsx=90,
marker=dict(color=fe.band_colour("blue", 0.55), line=dict(width=0)),
name="monthly premia",
hovertemplate="%{x:.1f}%<extra></extra>"), row=1, col=1)
grid = np.linspace(p.min(), p.max(), 400)
fig.add_trace(go.Scatter(
x=grid, y=stats.norm.pdf(grid, p.mean(), p.std(ddof=1)),
line=dict(color=ORANGE, width=2), name="matched normal",
hovertemplate="%{x:.1f}%<extra></extra>"), row=1, col=1)
# --- right: QQ plot of both series, standardised ---
for x, hue, label in [(premium, BLUE, "Equity premium"),
(gap, ORANGE, "NASDAQ − S&P 500")]:
z = ((x - x.mean()) / x.std(ddof=1)).to_numpy()
z.sort()
theo = stats.norm.ppf((np.arange(1, len(z) + 1) - 0.5) / len(z))
fig.add_trace(go.Scatter(
x=theo, y=z, mode="markers",
marker=dict(color=hue, size=3.5, opacity=0.65), name=label,
hovertemplate="theory %{x:.2f}<br>data %{y:.2f}<extra></extra>"),
row=1, col=2)
lim = 5.0
fig.add_trace(go.Scatter(
x=[-lim, lim], y=[-lim, lim], mode="lines",
line=dict(color=fe.INK["muted"], width=1, dash="dash"),
showlegend=False, hoverinfo="skip"), row=1, col=2)
fig.update_xaxes(title_text="monthly excess return (%)", row=1, col=1)
fig.update_yaxes(title_text="density", row=1, col=1)
fig.update_xaxes(title_text="normal quantile", range=[-lim, lim], row=1, col=2)
fig.update_yaxes(title_text="standardised return", range=[-8.5, 8.5], row=1, col=2)
fe.figure(fig, height=430)
fig.update_layout(margin=dict(t=88), legend=dict(y=1.14), bargap=0.02)
figThe left panel shows the problem in the shape everyone recognises. The histogram is too tall in the middle and too heavy at the edges, while the fitted normal curve is too wide in the middle and falls away far too quickly at the sides. Quiet months are more common than a normal distribution expects, ordinary months are less common, and disastrous months are very much more common. This is exactly what excess kurtosis of 7 looks like.
The right panel is worth learning to read properly, because it is the diagnostic you will use most often. A quantile–quantile plot sorts the observations from smallest to largest and plots each one against the value a normal distribution would have put in that position. If the data really were normal, the two would agree all the way along and every point would sit on the dashed 45-degree line. Here they agree in the middle, where the points lie neatly on the line, and then peel away at both ends. Reading the far right: the largest standardised observation in the equity premium is about \(+7\), whereas a normal sample of this size would have produced a largest value of roughly \(+3.3\). The data’s extremes are more than twice as extreme as they should be. The same thing happens on the left, which is why the plot forms a stretched S rather than bending in only one direction.
That two-sided departure is kurtosis, and it is overwhelmingly the dominant feature of both series.
Skewness turns out to be a much smaller effect here — and, curiously, it runs the opposite way to the story everyone tells about markets. The received wisdom is that equity returns are negatively skewed because crashes are violent and recoveries are gradual. Measured over the full CRSP sample, though, the skewness of the equity premium is +0.15, which is slightly positive. The explanation is visible in Figure 3.1: the single largest monthly moves in a century of data are not the crashes of 1929 and 1931 but the violent rebounds that followed them. April 1933 gained 39%, which is a larger move than any month in the sample lost.
This is not a general fact about equities, and it would be a mistake to generalise from it — it is a fact about a sample that happens to include the Great Depression. Over shorter, post-war windows the familiar negative skew reappears: the S&P 500’s monthly returns since 1971 have a skewness of -0.44. The reason to care about the sign at all is that Section 3.4 will show it determines the direction in which a finite-sample test misleads you, so it is not a detail we can leave unexamined.
3.3 Where the exact answer comes from, and what it costs
Suppose for the moment that we ignore everything in the previous section and carry on regardless. The statistic we computed above, written out in full, is
\[ T \;=\; \frac{\bar{X} - \mu_0}{S/\sqrt{n}}, \qquad S^2 = \frac{1}{n-1}\sum_{i=1}^{n}(X_i - \bar{X})^2 , \tag{3.2}\]
and every statistical package in existence will tell you that under the null hypothesis \(T\) follows a \(t\) distribution with \(n-1\) degrees of freedom. That claim is exact. It is not an approximation that improves with sample size; it is an equality that holds at \(n = 4\) as precisely as at \(n = 4{,}000\). Results of that quality are rare, and it is worth understanding exactly what we had to assume to get one.
The derivation is in Wackerly, Mendenhall, and Scheaffer (2008), and it comes in three steps. Watch for the assumption of normality, because it is used in all three, in a different way each time.
Step 1: the numerator is exactly normal. Suppose \(X_1,\dots,X_n\) are independent draws from \(N(\mu, \sigma^2)\). The sample mean is a weighted sum of them, \(\bar{X} = \frac{1}{n}\sum X_i\), and the normal family has a special property: any linear combination of independent normal variables is itself normal, exactly, with no approximation. So \(\bar{X}\) is normal, with mean \(\mu\) and variance \(\sigma^2/n\), and standardising gives
\[ Z \;=\; \frac{\bar{X} - \mu}{\sigma/\sqrt{n}} \;\sim\; N(0,1). \tag{3.3}\]
This closure property is unusual. Average a set of uniform variables and you do not get a uniform; average a set of exponentials and you do not get an exponential. For any parent other than the normal, \(\bar{X}\) has some other distribution at finite \(n\), and in general it is one without a convenient name. As a rough guide, the average inherits the family traits of the parent: the mean of skewed data is itself skewed, and the mean of heavy-tailed data is itself heavy-tailed, though both effects weaken as \(n\) grows. That weakening is the subject of Section 3.5, but at this stage of the argument it is no help at all, because we are after an exact statement.
Step 2: the denominator is exactly chi-square. Still assuming a normal sample,
\[ \frac{(n-1)S^2}{\sigma^2} \;\sim\; \chi^2_{n-1}. \tag{3.4}\]
The intuition is that \(S^2\) is built from a sum of squared deviations, and a sum of squares of independent standard normals is by definition a chi-square. The degrees of freedom are \(n-1\) rather than \(n\) because the deviations are taken around \(\bar{X}\) rather than around the unknown \(\mu\): once you know \(n-1\) of the deviations, the last one is determined, since they must sum to zero. Week 0 met this result when it explained why the sampling distribution of a variance leans to the right rather than being symmetric. Like step 1, it is a statement about normal populations and no others.
Step 3: the numerator and the denominator are independent. This is the step that textbooks tend to state and move past, and it is by far the most demanding of the three. Recall how a \(t\)-distributed variable is defined: it is the ratio \(Z/\sqrt{W/\nu}\) where \(Z\) is standard normal, \(W\) is chi-square with \(\nu\) degrees of freedom, and the two are independent of each other. Without independence the ratio has some other distribution entirely, because the way the denominator behaves would then carry information about the numerator.
So to reach \(t_{n-1}\) we need \(\bar{X}\) and \(S^2\) — the estimated mean and the estimated variance, computed from the very same data — to be statistically independent. That should strike you as a strong requirement, and it is.
Why should it hold at all? Intuitively, in a normal sample the location of the data and its spread carry no information about each other: shifting every observation up by a constant changes \(\bar{X}\) and leaves \(S^2\) untouched, and the symmetry of the normal means that a sample which happens to land high is no more or less spread out than one which lands low. That intuition is correct, and it is also considerably more fragile than it looks.
Independence of \(\bar{X}\) and \(S^2\) is not a convenience. It is normality. For i.i.d. observations with finite variance, the sample mean and the sample variance are independent if and only if the population is normal (Lukacs 1942). The implication runs in both directions: the property does not merely follow from normality, it characterises it. If you found a non-normal population whose sample mean and sample variance were independent, you would have disproved a theorem.
This closes off the escape route one would naturally reach for. With most assumptions there is some hope of weakening them — replacing “normal” with “symmetric”, or with “finite fourth moment”, or with some regularity condition that most real data satisfy. Here there is no such hope, because the property we need is logically equivalent to the assumption we are trying to avoid. Any departure from normality, however slight, breaks the independence in step 3, and once step 3 goes the \(t_{n-1}\) result goes with it.
Which is to say: the distribution tabulated in the back of every statistics textbook is the correct distribution for exactly one family of data-generating processes. Section 3.2 established that financial returns are not in that family.
Put the three steps together. Substituting eq. 3.3 and eq. 3.4 into the definition of a \(t\) variable, the unknown \(\sigma\) cancels between numerator and denominator — which is the whole point of the construction, since \(\sigma\) is not something we know:
\[ T \;=\; \frac{Z}{\sqrt{\dfrac{(n-1)S^2/\sigma^2}{n-1}}} \;=\; \frac{(\bar{X} - \mu)/(\sigma/\sqrt{n})}{S/\sigma} \;=\; \frac{\bar{X} - \mu}{S/\sqrt{n}} \;\sim\; t_{n-1}. \tag{3.5}\]
Step back and look at the structure of what just happened, because every exact finite-sample result in statistics has this same shape. A complete specification of how the data were generated goes in at the top; a named, tabulated distribution for the statistic comes out at the bottom. The conclusion is precise and unimprovable — and it is precise about the model you fed in, not about the data on your desk.
Change the assumed distribution of \(X_i\) and eq. 3.5 does not degrade gracefully into something approximately right. It simply becomes false, and is replaced by some other distribution which depends on the new assumption in a complicated way, usually has no closed form, and generally has to be found by simulation — when it can be found at all. That is the subject of the next two sections: first how wrong things go, and then what we do instead.
3.4 Fragility: exact for the assumed model, wrong for yours
An exact result carries no approximation error at all, which is what makes it the strongest kind of statement statistics can offer. The very same property is what makes it brittle. Because the derivation used normality three times and in three different ways, it tells us what happens under normality and says nothing whatsoever about what happens otherwise. There is no continuity clause in eq. 3.5 — no theorem promising that data which are nearly normal produce a distribution which is nearly \(t_{n-1}\).
That the guarantee vanishes is a matter of logic. Whether it matters in practice is an empirical question, and the honest way to settle it is to simulate.
The experiment is simple enough to describe in a sentence. Pick a population whose distribution we know. Draw a sample of size \(n\) from it, compute the \(t\)-statistic in eq. 3.2, and write the answer down. Do that several hundred thousand times, and the collection of answers traces out the true sampling distribution of \(T\) for that population — no theory required, just brute force. Then lay the \(t_{n-1}\) density that the textbook prescribes on top and see whether they match.
To make the comparison fair we use three populations standardised to have mean 0 and variance 1, so that the only thing that differs between them is the shape. Any discrepancy we find can then be blamed on skewness or kurtosis alone, rather than on one population simply being more spread out than another.
rng = np.random.default_rng(fe.SEED)
DGPS = {
"Normal": dict(hue=BLUE, note="skew 0, excess kurtosis 0"),
"Student's t, ν=5": dict(hue=ORANGE, note="skew 0, excess kurtosis 6"),
"Chi-square, ν=2": dict(hue=AQUA, note="skew 2, excess kurtosis 6"),
}
def draw(kind, size):
"""Draw from one of the three populations, standardised to mean 0, var 1."""
if kind == "Normal":
return rng.standard_normal(size)
if kind == "Student's t, ν=5":
return rng.standard_t(5, size) / np.sqrt(5 / 3) # var of t_5 is 5/3
if kind == "Chi-square, ν=2":
return (rng.chisquare(2, size) - 2) / 2 # mean 2, var 4
raise ValueError(kind)
def t_stats(kind, n, reps, chunk=20_000):
"""Distribution of the t-statistic under the null, by simulation."""
out = []
done = 0
while done < reps:
m = min(chunk, reps - done)
X = draw(kind, (m, n))
out.append(X.mean(1) / (X.std(1, ddof=1) / np.sqrt(n)))
done += m
return np.concatenate(out)The three populations are chosen to separate the two effects we care about. The normal is the control: here the textbook result is exactly right, so anything other than a perfect match would mean the simulation itself is broken. The Student’s \(t\) with 5 degrees of freedom is symmetric like the normal but has much heavier tails, with an excess kurtosis of 6 — close to what we measured on the real data. The chi-square with 2 degrees of freedom, recentred and rescaled, has the same excess kurtosis of 6 but is strongly asymmetric, with a skewness of 2.
That last pairing is the useful one. Because the \(t_5\) and the chi-square share a kurtosis and differ only in skewness, any difference in how they behave has to be caused by the asymmetry rather than by the tails. This is the kind of controlled comparison that a simulation makes easy and algebra makes painful.
One detail in the code deserves a note. t_stats draws the samples in chunks of 20,000 rather than all at once, because a matrix of 400,000 replications by 1,000 observations would need several gigabytes of memory. The results are identical; it simply keeps the machine alive.
N_SMALL, REPS = 20, 400_000
fig = go.Figure()
edges = np.linspace(-6, 4, 101) # 0.1-wide bins: fine enough, not noisy
mid = (edges[:-1] + edges[1:]) / 2
for kind, spec in DGPS.items():
T = t_stats(kind, N_SMALL, REPS)
dens, _ = np.histogram(T, bins=edges, density=True)
fig.add_trace(go.Scatter(
x=mid, y=dens, line=dict(color=spec["hue"], width=2),
name=f"{kind} ({spec['note']})",
hovertemplate="t = %{x:.2f}<br>density %{y:.3f}<extra></extra>"))
# Drawn last, and dark, so it stays visible where the normal line sits on it.
fig.add_trace(go.Scatter(
x=mid, y=stats.t.pdf(mid, N_SMALL - 1),
line=dict(color=fe.INK["primary"], width=1.4, dash="6px,5px"),
name=f"t_{N_SMALL-1} (what the table assumes)",
hovertemplate="t = %{x:.2f}<extra></extra>"))
fig.update_xaxes(title_text="t-statistic", range=[-6, 4])
fig.update_yaxes(title_text="density")
fe.figure(fig, height=430)
fig.update_layout(margin=dict(t=96), legend=dict(y=1.10))
figStart with the normal line, which is the sanity check. It lies on the dashed reference so precisely that you have to look for the dashes to see them, which is eq. 3.5 doing exactly what it promised at a sample size of only 20.
The heavy-tailed line is the first surprise. The \(t_5\) population has an excess kurtosis of 6, meaning its extreme observations are dramatically more common than a normal’s — and yet its \(t\)-statistic distribution is very nearly the textbook one. If you had been asked to guess in advance which of our two violations of normality would do the most damage, heavy tails would have been the natural answer, and it would have been wrong.
The skewed line is the second surprise, and there are two things happening to it. It is shifted, so its centre no longer sits at zero, and it is tilted, with a long tail stretching out to the left and a compressed, foreshortened right side.
The direction of that tilt is genuinely counter-intuitive and repays a careful explanation. The chi-square population is skewed to the right — its rare extreme observations are large and positive. Yet the distribution of its \(t\)-statistic is skewed to the left. The sign has flipped.
Here is why. In right-skewed data, the sample mean and the sample standard deviation are positively related to each other, and it is the same observations that drive both. Suppose a particular sample happens to contain two or three of the rare large values. Those observations pull the sample mean \(\bar{X}\) upward, which enlarges the numerator of eq. 3.2. But being far from the rest of the data, they also inflate \(S\), which enlarges the denominator. The two effects partly offset, so \(T\) does not get as large as the numerator alone would suggest; the statistic is pulled back toward zero exactly when it would otherwise have been most extreme.
Now consider a sample that happens to miss all the large values. Its mean is below average, giving a negative numerator, but its spread is also small, giving a small denominator — and a small denominator makes the ratio bigger in magnitude. So on the downside the two effects reinforce rather than offset, and \(T\) is free to run a long way negative.
Large positive values get damped; large negative values get amplified. That asymmetry in the statistic is precisely the left tail we see in the figure, and it is the mirror image of the asymmetry in the data. Remember the direction: it is what determines which way a test goes wrong, and Table 3.3 is about to make the consequences concrete.
rows = []
for n in [20, 60, 250, 1000]:
row = {"n": n}
for kind in DGPS:
T = t_stats(kind, n, 200_000)
row[f"{kind} — upper"] = np.mean(T > stats.t.ppf(0.95, n - 1)) * 100
row[f"{kind} — 2-sided"] = np.mean(np.abs(T) > stats.t.ppf(0.975, n - 1)) * 100
rows.append(row)
pd.DataFrame(rows).set_index("n").round(2)| Normal — upper | Normal — 2-sided | Student's t, ν=5 — upper | Student's t, ν=5 — 2-sided | Chi-square, ν=2 — upper | Chi-square, ν=2 — 2-sided | |
|---|---|---|---|---|---|---|
| n | ||||||
| 20 | 5.06 | 5.02 | 5.01 | 4.67 | 1.89 | 8.10 |
| 60 | 4.91 | 4.94 | 5.05 | 4.94 | 2.78 | 6.28 |
| 250 | 5.01 | 4.98 | 5.01 | 4.93 | 3.73 | 5.30 |
| 1000 | 4.96 | 4.95 | 5.00 | 4.98 | 4.36 | 5.09 |
Before reading the numbers, be clear about what they measure. Each entry is the proportion of simulated samples in which a test designed to reject 5% of the time actually rejected — with the null hypothesis true by construction, so every rejection is a mistake. This quantity is called the size of the test, and a correctly calibrated test reports 5.00. Anything above 5 means the test finds significance where there is none more often than advertised; anything below means it is unduly reluctant, and will miss real effects when they exist.
Read one column pair at a time. Under normality, every entry is 5.00 at every sample size, which is what “exact” means: the result does not improve as \(n\) grows, because it was never wrong to begin with.
Under heavy tails, the entries are also close to 5, and — the point worth noticing — they are close even at \(n = 20\). Kurtosis of 6 costs us essentially nothing.
Under skewness, the picture falls apart. At \(n = 20\) a one-sided test that should reject 5% of the time rejects less than 2%, while a two-sided test that should reject 5% rejects more than 8%. Both are badly wrong, and note that they are wrong in opposite directions — which follows directly from the tilt we diagnosed in Figure 3.3. Because the \(t\)-statistic’s distribution has a thin right tail, the upper-tail critical value from the table sits too far out and almost nothing exceeds it: too few rejections. Because it simultaneously has a fat left tail, a two-sided test picks up a flood of extreme negative values that the table did not anticipate: too many rejections.
Skewness is what bites, not heavy tails. This is the most useful practical lesson in the chapter, and it contradicts what most people worry about. Everyone knows financial returns have fat tails, and everyone assumes fat tails are what break the standard inference. The simulation says otherwise.
The reason is that the \(t\)-statistic is self-normalising. It does not compare \(\bar{X}\) against a fixed yardstick; it compares \(\bar{X}\) against \(S\), a yardstick computed from the same data. When a sample happens to contain unusually extreme observations, they enlarge the numerator — but they enlarge the denominator too, and the ratio is left roughly where it was. Fat tails inflate both halves of the fraction at once, and the fraction barely notices.
Skewness gets no such cancellation. Asymmetry moves the numerator up and down by different amounts depending on direction, while the denominator, which depends on squared deviations, cannot tell the difference between an observation that is far above the mean and one equally far below. The two halves of the statistic therefore stop moving together, and the ratio is distorted.
Chapter 1 predicted this from the algebra. Expanding the distribution of a standardised mean in powers of \(n\), the leading correction term is proportional to the parent’s skewness and shrinks like \(1/\sqrt{n}\); the kurtosis correction appears only at order \(1/n\) and dies much faster. What we have just simulated is that expansion made visible, and it explains both the direction of the error and the slow rate at which Figure 3.4 will show it disappearing.
3.4.1 Problems where no exact answer exists at all
Everything so far has been about an exact result that exists but breaks when its assumptions fail. That is the milder of the two difficulties with finite-sample theory. The harder one is that for a great many entirely ordinary questions, there is no exact result to break in the first place — not because nobody has been clever enough, but because none exists. And this is true even when the data are perfectly normal.
The Behrens–Fisher problem. Take problem 2 in its unpaired form. We have two independent samples, both normal, and we want to test whether their means are equal. The complication is that their variances are different and both unknown.
This is about as favourable a setting as statistics ever offers: everything is normal, everything is independent, there are only four parameters in the whole problem. And yet no test statistic is known whose exact null distribution is free of the unknown ratio of the two variances. That qualifier is the crux. One can certainly write down statistics whose distribution is calculable — but the answer depends on \(\sigma_1^2/\sigma_2^2\), which we do not know, so the “exact” distribution cannot actually be used to produce a critical value. A nuisance parameter has contaminated the answer.
Nearly a century of effort has gone into this, and what it has produced is approximations. Welch’s (Welch 1947) is the one your software uses by default; it works by choosing a fractional degrees-of-freedom parameter so that the \(t\) distribution matches the true one as closely as possible. It is very good. It is not exact, and it is not claimed to be.
The same difficulty appears in finance whenever a quantity is a ratio of estimates. The Sharpe ratio is the obvious case: its numerator and denominator are both estimated, both noisy, and correlated with each other. Its exact distribution is unknown outside artificial special cases, and every standard error you have ever seen attached to a Sharpe ratio came from an asymptotic approximation (Lo 2002).
The unpaired form is not an academic curiosity here, because the two indices have visibly different variances — annualised volatility of 15.2% for the S&P 500 against 20.7% for the NASDAQ. It is exactly the Behrens–Fisher setting.
n = len(gap)
se_paired = gap.std(ddof=1) / np.sqrt(n)
se_unpaired = np.sqrt(idx["NASDAQ"].var(ddof=1) / n + idx["SP500"].var(ddof=1) / n)
pd.DataFrame([
{"approach": "Paired (mean of the monthly difference)",
"s.e. of mean (ann. %)": se_paired * 12 * 100,
"t-stat": gap.mean() / se_paired,
"exact distribution known?": "only if returns are normal"},
{"approach": "Unpaired (two independent samples)",
"s.e. of mean (ann. %)": se_unpaired * 12 * 100,
"t-stat": gap.mean() / se_unpaired,
"exact distribution known?": "no — Behrens–Fisher"},
]).set_index("approach").round(3)| s.e. of mean (ann. %) | t-stat | exact distribution known? | |
|---|---|---|---|
| approach | |||
| Paired (mean of the monthly difference) | 1.437 | 2.175 | only if returns are normal |
| Unpaired (two independent samples) | 3.488 | 0.896 | no — Behrens–Fisher |
The two approaches do not merely disagree about precision — they disagree about the answer. Pairing gives a \(t\)-statistic of 2.18, which is on the edge of conventional significance. Treating the samples as independent gives 0.90, which nobody would call evidence of anything at all. Same data, same point estimate of the mean difference, opposite conclusions.
The arithmetic behind the gap is worth seeing. The unpaired standard error adds the two variances, as though the indices moved independently: \(\sqrt{(\sigma_N^2 + \sigma_S^2)/n}\). The paired standard error uses the variance of the difference, which for correlated series is \(\sigma_N^2 + \sigma_S^2 - 2\rho\,\sigma_N \sigma_S\). With \(\rho = 0\.87\) that final term removes most of the total, which is why the paired standard error comes out 2.4 times smaller. The correlation is not a nuisance to be averaged away; it is the single most valuable feature of the data, and the unpaired calculation throws it in the bin.
3.5 The asymptotic alternative
At this point we have a choice. We can keep trying to derive the exact distribution of our statistic under increasingly baroque assumptions about returns, or we can give up on exactness and ask for something weaker that we can actually obtain. Asymptotic analysis takes the second route.
The idea is to stop asking what the distribution of a statistic is at our particular sample size, and ask instead what it converges to as the sample size grows without bound. That limiting distribution turns out to be obtainable under assumptions so mild that real data satisfy them, and we then use it as an approximation at whatever \(n\) we happen to have. Two theorems do the work: one tells us the estimate goes to the right place, the other tells us how it is scattered around it on the way.
3.5.1 The law of large numbers
For an i.i.d. sequence \(\{X_i\}\) with \(\mathbb{E}[|X_i|] < \infty\) and \(\mathbb{E}[X_i] = \mu\),
\[ \bar{X}_n \;\xrightarrow{\ a.s.\ }\; \mu . \tag{3.6}\]
In words: as we collect more data, the sample mean settles onto the true mean and stays there. The arrow with “a.s.” above it means almost sure convergence, which is the strongest of the standard modes — it says that the sequence of sample means converges to \(\mu\) for essentially every possible realisation of the data, not merely that it is usually close.
The condition is remarkably weak. All that is required is that \(\mathbb{E}[|X_i|]\) be finite, which is to say that the population mean exists at all. Nothing is assumed about symmetry, tails, or shape. This is what justifies computing a sample average in the first place, and it is why we could report a point estimate in Table 3.1 without any of the machinery that follows.
What the law of large numbers does not do is tell us anything about precision. It says \(\bar{X}_n\) ends up at \(\mu\) eventually; it says nothing about how far away we are after 1,188 observations. For that we need the second theorem.
3.5.2 The central limit theorem
In the Lindeberg–Lévy form, for i.i.d. \(\{X_i\}\) with \(\mathbb{E}[X_i] = \mu\) and \(0 < \operatorname{Var}(X_i) = \sigma^2 < \infty\),
\[ \sqrt{n}\,\frac{\bar{X}_n - \mu}{\sigma} \;\xrightarrow{\ d\ }\; N(0,1). \tag{3.7}\]
Three features of this statement deserve unpacking.
First, the \(\sqrt{n}\) is not decoration. By the law of large numbers, \(\bar{X}_n - \mu\) is heading to zero, so on its own it would converge to a degenerate distribution sitting entirely at a point — true but useless. Dividing by \(\sigma/\sqrt{n}\), its own standard deviation, rescales the shrinking quantity so that it neither collapses nor explodes. \(\sqrt{n}\) is precisely the rate at which the mean approaches \(\mu\), and blowing the picture up at exactly that rate is what leaves something non-trivial to look at.
Second, the convergence is in distribution, which is a weaker notion than the one in eq. 3.6. It does not claim that the standardised mean gets close to any particular normal random variable. It claims only that its cumulative distribution function approaches the normal cdf: the probability of landing below any given value converges to what a normal would give. That is exactly the kind of statement we need for computing \(p\)-values, and no more.
Third, look at what is absent from the assumptions. The exact result in eq. 3.5 required the entire distribution of \(X_i\) to be normal. eq. 3.7 requires a finite mean and a finite, non-zero variance. Skewness does not appear. Kurtosis does not appear. The data may be discrete, bounded, wildly asymmetric, or all three. Whatever shape the parent distribution has, the standardised mean forgets it.
Why should averaging destroy the shape? Loosely, because a sum of many independent contributions is dominated by the accumulation of typical fluctuations rather than by the character of any individual draw. The idiosyncratic features of the parent — the long tail on one side, the lump near zero — get diluted as more and more terms are added, and only the mean and the variance survive the process. Figure 1.2 in week 0 showed this happening.
That is the trade being made, and it should be stated plainly: we accept a weaker conclusion, one that is formally true only in a limit we never reach, in exchange for assumptions that our data can actually satisfy. Whether that is a good bargain depends entirely on how close the limit is at \(n = 1{,}188\), which is a question we will answer by simulation in Section 3.6 rather than by assertion.
3.5.3 From infeasible to feasible
There is still something wrong with eq. 3.7 as a practical tool: it contains \(\sigma\), the population standard deviation, which is not a number we possess. A result we cannot evaluate is no result at all, and for that reason eq. 3.7 is called the infeasible central limit theorem.
The obvious move is to replace \(\sigma\) with the sample standard deviation \(S\), giving the feasible statistic — which is, at last, the thing we actually computed at the start of the chapter:
\[ T_n \;=\; \sqrt{n}\,\frac{\bar{X}_n - \mu}{S} \;\xrightarrow{\ d\ }\; N(0,1). \tag{3.8}\]
The obvious move is not automatically a legitimate one, though. We have replaced a constant with a random variable, and \(T_n\) now has two sources of noise where the infeasible version had one. Why does that not change the limiting distribution? The justification is a chain of three results, each doing a specific job (Sheppard 2021):
- The law of large numbers, applied to squared deviations rather than to the data themselves, gives \(S^2 \xrightarrow{p} \sigma^2\). The sample variance converges in probability to the true variance.
- The continuous mapping theorem says that convergence survives passage through a continuous function. Since the square root is continuous at \(\sigma^2 > 0\), we get \(S \xrightarrow{p} \sigma\) directly.
- Slutsky’s theorem finishes the job. It states that if \(A_n\) converges in distribution to some \(A\), and \(B_n\) converges in probability to a non-zero constant \(c\), then \(A_n/B_n\) converges in distribution to \(A/c\). Setting \(A_n\) to the infeasible statistic and \(B_n\) to \(S/\sigma\), which converges in probability to 1, the limiting distribution is unchanged.
The word “constant” in step 3 is what makes the whole thing work. \(S\) is random at any finite \(n\), but in the limit it stops being random and settles on a fixed number, and a quantity that has stopped fluctuating cannot contribute any uncertainty. The practical upshot is worth stating in plain terms: asymptotically, estimating the standard error is free. Notice how different this is from the finite-sample story, where the randomness of \(S\) was the entire reason we needed a \(t\) distribution with \(n-1\) degrees of freedom instead of a normal.
A distinction people routinely blur. eq. 3.8 does not say that \(T_n\) follows a \(t_{n-1}\) distribution. It says \(T_n\) is approximately standard normal. These are different claims with different justifications: the \(t\)-table is a finite-sample device belonging to eq. 3.5 and its normality assumption, whereas the asymptotic argument produces a normal and never mentions \(t\) at all.
In practice the difference is invisible here. Because \(t_{n-1} \to N(0,1)\) as \(n\) grows, the two critical values at our sample size are 1.960 and 1.962 — a discrepancy in the third decimal place. Reporting either is fine. Knowing which one you are entitled to, and why, is what stops you from also believing the normality assumption that came bundled with the \(t\).
There is one more thing to say before we test any of this, and it is the caveat that the rest of the chapter turns on. eq. 3.8 is a statement about a limit, and a limit says nothing whatsoever about any particular finite \(n\). A sequence can approach its limit immediately or take a million terms to get close, and the theorem does not distinguish the two cases. So the practical question is never “is the central limit theorem true?” — it is — but “how close to the limit are we at the sample size I actually have?” That is an empirical question, it has a different answer for every data-generating process, and it can only be settled by looking.
GRID = [10, 15, 20, 30, 50, 75, 100, 150, 250, 400, 650, 1000]
sizes = {kind: {"upper": [], "two": []} for kind in DGPS}
for n in GRID:
cu, ct = stats.t.ppf(0.95, n - 1), stats.t.ppf(0.975, n - 1)
for kind in DGPS:
T = t_stats(kind, n, 100_000)
sizes[kind]["upper"].append(np.mean(T > cu) * 100)
sizes[kind]["two"].append(np.mean(np.abs(T) > ct) * 100)
fig = make_subplots(rows=1, cols=2, horizontal_spacing=0.09, shared_yaxes=True,
subplot_titles=("One-sided 5% test (problem 1)",
"Two-sided 5% test (problem 2)"))
for col, key in [(1, "upper"), (2, "two")]:
for kind, spec in DGPS.items():
fig.add_trace(go.Scatter(
x=GRID, y=sizes[kind][key], mode="lines+markers",
line=dict(color=spec["hue"], width=2), marker=dict(size=5),
name=kind, showlegend=(col == 1),
hovertemplate="n = %{x}<br>true size %{y:.2f}%<extra></extra>"),
row=1, col=col)
fig.add_hline(y=5, line=dict(color=fe.INK["muted"], width=1, dash="dash"),
row=1, col=col)
TICKS = [10, 20, 50, 100, 250, 1000]
for c in (1, 2):
fig.update_xaxes(type="log", title_text="sample size n (log scale)",
tickvals=TICKS, ticktext=[f"{v:,}" for v in TICKS],
row=1, col=c)
fig.update_yaxes(title_text="true rejection rate (%)", range=[0, 11], row=1, col=1)
fe.figure(fig, height=430)
fig.update_layout(margin=dict(t=96), legend=dict(y=1.10))
figThis figure answers the question the theorem refuses to. Each line traces how quickly a test’s true rejection rate settles onto its advertised 5% as the sample grows, for one of our three populations.
The normal line is flat at 5 across the whole range, from \(n = 10\) to \(n = 1{,}000\). That is what an exact result buys you: sample size is simply irrelevant, because there was never any approximation error to shrink.
The heavy-tailed line is close to 5 essentially from the start. Even at \(n = 10\) — a sample so small that most people would refuse to do inference at all — a population with excess kurtosis of 6 produces a test that is very nearly correctly sized. Self-normalisation does its work immediately and does not need a large sample to do it.
The skewed line is the one that matters. It does converge, exactly as the central limit theorem promises, but it takes its time. In the left panel the one-sided test is still rejecting around 3% instead of 5% at \(n = 100\), and it has not fully arrived even at \(n = 1{,}000\). In the right panel the two-sided test converges faster, reaching acceptable calibration by \(n \approx 250\). Notice again that the two panels are wrong in opposite directions — below 5 on the left, above 5 on the right — which is the tilt from Figure 3.3 showing up in the rejection rates.
Two lessons follow, and they are the practical content of this whole chapter.
The first is that there is no universal minimum sample size. The folklore rule that thirty observations are enough for the central limit theorem is visible here as nonsense in both directions: for the heavy-tailed population, ten is plenty, while for the skewed population, three hundred is marginal. The number that matters is not \(n\) by itself but \(n\) in combination with how asymmetric the data are.
The second is that the one-sided and two-sided tests are not equally trustworthy, and which one you are running is not a detail. Problem 1 asks whether the premium is positive, which is a one-sided question and therefore the more demanding case. Problem 2 asks whether two returns differ, which is two-sided and more forgiving. Had our data been strongly skewed, that distinction alone could have decided which of the two answers we were entitled to believe.
3.6 Does the approximation hold at our sample sizes?
The three populations in the last section were built to isolate mechanisms, not to imitate markets. A chi-square with skewness 2 tells us what asymmetry does in principle; it does not tell us what happens to our two test statistics on our two datasets. For that we can do something considerably better, and it is the most useful practical technique in the chapter.
The idea is to stop guessing at the population and use the sample as a stand-in for it. We have 1,188 observed monthly premia. Treat that collection as if it were the population — it carries the real kurtosis of 7, the real mild skewness, the real everything — and then simulate as before, drawing samples from the data itself. This is the bootstrap, and the procedure has three steps:
- Impose the null. Subtract the sample mean from every observation, so that the collection we are drawing from has a mean of exactly zero. This matters: we want the distribution of \(T\) when the null is true, and the observed data have a mean that is not zero. Without this step we would be simulating under the alternative and would learn nothing about critical values.
- Resample with replacement. Draw \(n\) observations at random from that recentred collection, allowing repeats, to build one artificial sample of the same size as the real one.
- Compute and repeat. Calculate the \(t\)-statistic on the artificial sample, store it, and go back to step 2 a hundred thousand times.
The resulting collection of \(t\)-statistics is an estimate of the true sampling distribution under the null — using the actual shape of monthly returns rather than any assumption about it. If the asymptotic argument is working at our sample size, this distribution should look like a standard normal. If it is not, we will see the discrepancy directly, and we will know not to trust the \(p\)-values in Table 3.6.
def bootstrap_t(x, reps=100_000, chunk=20_000):
"""Distribution of T when H0 is true, resampling the observed returns."""
v = x.to_numpy()
v = v - v.mean() # impose the null
k, out, done = len(v), [], 0
while done < reps:
m = min(chunk, reps - done)
S = v[rng.integers(0, k, size=(m, k))]
out.append(S.mean(1) / (S.std(1, ddof=1) / np.sqrt(k)))
done += m
return np.concatenate(out)
fig = make_subplots(rows=1, cols=2, horizontal_spacing=0.10,
subplot_titles=(f"Equity premium (n = {len(premium)})",
f"NASDAQ − S&P 500 (n = {len(gap)})"))
edges = np.linspace(-4.5, 4.5, 91)
mid = (edges[:-1] + edges[1:]) / 2
checks = {}
for col, (x, label) in enumerate([(premium, "Equity premium"),
(gap, "NASDAQ − S&P 500")], start=1):
T = bootstrap_t(x)
checks[label] = T
dens, _ = np.histogram(T, bins=edges, density=True)
fig.add_trace(go.Scatter(
x=mid, y=dens, line=dict(color=BLUE, width=2),
name="resampled returns", showlegend=(col == 1),
hovertemplate="t = %{x:.2f}<extra></extra>"), row=1, col=col)
fig.add_trace(go.Scatter(
x=mid, y=stats.norm.pdf(mid),
line=dict(color=ORANGE, width=1.6, dash="dash"),
name="N(0,1)", showlegend=(col == 1),
hovertemplate="t = %{x:.2f}<extra></extra>"), row=1, col=col)
fig.update_xaxes(title_text="t-statistic under the null", row=1, col=1)
fig.update_xaxes(title_text="t-statistic under the null", row=1, col=2)
fig.update_yaxes(title_text="density", row=1, col=1)
fe.figure(fig, height=430)
fig.update_layout(margin=dict(t=88), legend=dict(y=1.14))
figrows = []
for label, T in checks.items():
rows.append({
"series": label,
"n": len(premium) if label == "Equity premium" else len(gap),
"true 2.5th pct": np.percentile(T, 2.5),
"true 97.5th pct": np.percentile(T, 97.5),
"true size of nominal 5% (2-sided)": np.mean(np.abs(T) > 1.96) * 100,
"true size of nominal 5% (upper)": np.mean(T > 1.645) * 100,
})
pd.DataFrame(rows).set_index("series").round(3)| n | true 2.5th pct | true 97.5th pct | true size of nominal 5% (2-sided) | true size of nominal 5% (upper) | |
|---|---|---|---|---|---|
| series | |||||
| Equity premium | 1188 | -1.975 | 1.955 | 5.059 | 4.881 |
| NASDAQ − S&P 500 | 652 | -1.973 | 1.931 | 4.921 | 4.754 |
The verdict is about as good as it could be. Despite excess kurtosis above 7 in both series, the true critical values sit within roughly 0.03 of the \(\pm 1.960\) that the normal approximation claims, and a test with a nominal size of 5% really does reject about 5% of the time — 5.0% two-sided and a shade under 5% in the upper tail, the small shortfall being the mild positive skewness we identified in Section 3.2 behaving exactly as Table 3.3 said it would.
It is worth appreciating how much has been achieved and how little was assumed. We have a distribution for our test statistic that is accurate to two decimal places, and we obtained it without ever claiming that returns are normal — which would have been false — and without needing an exact finite-sample result, which in the unpaired version of problem 2 does not even exist. At \(n\) in the hundreds, with skewness this mild, the asymptotic approximation is not a compromise we are tolerating. It is effectively the exact answer.
That said, keep the scope of the claim in view. What has been validated is the normal approximation for these two series at these two sample sizes. It is not a general licence, and Section 3.8 lists the ways the same procedure could mislead on a different dataset.
With the distribution in hand, we can finally convert our two statistics into answers.
def verdict(x, label, one_sided):
n = len(x)
se = x.std(ddof=1) / np.sqrt(n)
t = x.mean() / se
p = 1 - stats.norm.cdf(t) if one_sided else 2 * (1 - stats.norm.cdf(abs(t)))
return {
"problem": label,
"mean (ann. %)": round(x.mean() * 12 * 100, 3),
"t-stat": round(t, 3),
"p-value": f"{p:.2e}", # 1's p-value rounds to 0 otherwise
"95% CI lower (ann. %)": round((x.mean() - 1.96 * se) * 12 * 100, 3),
"95% CI upper (ann. %)": round((x.mean() + 1.96 * se) * 12 * 100, 3),
}
pd.DataFrame([
verdict(premium, "1. Equity premium > 0?", one_sided=True),
verdict(gap, "2. NASDAQ − S&P 500 ≠ 0?", one_sided=False),
]).set_index("problem")| mean (ann. %) | t-stat | p-value | 95% CI lower (ann. %) | 95% CI upper (ann. %) | |
|---|---|---|---|---|---|
| problem | |||||
| 1. Equity premium > 0? | 8.253 | 4.456 | 4.18e-06 | 4.623 | 11.883 |
| 2. NASDAQ − S&P 500 ≠ 0? | 3.126 | 2.175 | 2.96e-02 | 0.309 | 5.942 |
3.7 Reading the output
The two problems have arrived at very different places, and the difference is instructive.
One is settled beyond argument; the other is genuinely marginal. The equity premium’s \(t\)-statistic of 4.46 corresponds to a \(p\)-value of roughly 4.2 in a million. Ask what it would take to overturn that. The largest distortion anywhere in Table 3.3 shifts a rejection rate by three percentage points, at a sample size sixty times smaller than ours and a skewness thirteen times larger. Errors of that order do not turn a \(p\)-value of \(4\times 10^{-6}\) into a non-result; they would have to be wrong by several orders of magnitude. When evidence is this strong, the choice of distribution theory genuinely does not matter, and an hour spent worrying about it is an hour wasted.
The index comparison is the opposite situation. Its \(p\)-value of 3.0% sits just below the conventional 5% line. A distortion of even two percentage points in the size of the test would move it across that line, so here the verdict does depend on the approximation being accurate — and that is exactly why Figure 3.5 was worth an hour of computer time instead of an assumption. Had this been a sixty-month sample rather than a 652-month one, the honest report would have been that we cannot tell.
The general lesson: how much care the distribution theory deserves depends on how close the answer is to the decision boundary. Check the approximation when it could change your mind, and do not agonise over it when it cannot.
A confidence interval says more than a verdict does. The premium’s interval runs from about 4.6% to 11.9% a year. We have established that it is positive. We have emphatically not established its size: the top of that interval is 2.6 times the bottom, and that is after ninety-nine years of data. If you are discounting a pension liability, a 4.6% equity premium and an 11.9% one imply completely different answers, and a century of history is not enough to choose between them. This is week 0’s conclusion restated at maximum sample length — volatility is estimable, expected return barely is — and it is worth taking personally, because it means the input that most valuation models are most sensitive to is the one we know least about.
Reject the null you tested, not the one you had in mind. Problem 2 rejects equality of average monthly price returns between 1971 and 2025. Read carefully, that is a much narrower statement than “the NASDAQ is the better investment”, and three things stand between the two. The NASDAQ carried 36% more volatility to earn its extra return, so a risk-adjusted comparison is a different exercise entirely. The sample opens in 1971, near the start of a technology expansion that may not repeat. And neither index includes dividends — an omission that systematically favours the higher-yielding S&P 500, which means our estimate of the gap is, if anything, too generous to the NASDAQ.
3.8 What can go wrong
The independence assumption has been doing unadvertised work throughout. Both eq. 3.6 and eq. 3.7 were stated for i.i.d. data, and the bootstrap in Figure 3.5 drew months independently at random — which quietly assumes that the order of the months carries no information.
For the mean of monthly equity returns this is not too damaging, because successive returns are close to uncorrelated. For the variance it is plainly false: volatility clusters, so turbulent months arrive in groups, and a century of data covers regimes whose volatility differs by a factor of several. The consequence is that the true standard error is larger than the one we computed, because dependent observations carry less information than independent ones — a sample of 1,188 correlated months is worth fewer than 1,188 independent draws.
The reassuring part is that the framework survives this. Central limit theorems exist for dependent and heterogeneous data, and that generality is precisely why asymptotic theory is usable in finance at all, where nothing is ever i.i.d. What changes is not the normal limit but the formula for the variance in it. Week 0’s block bootstrap, which resamples contiguous stretches rather than individual months, and the Newey–West estimator are the two standard responses; later weeks make them the default rather than the refinement.
Asymptotic does not mean assumption-free. eq. 3.7 needs a finite variance, and that is a real restriction rather than a formality. Note that it is an assumption about the population, not an observation about the sample: any finite collection of numbers has a finite variance, so no amount of staring at the data can confirm it. Some researchers have argued that returns are drawn from distributions with infinite variance, in which case every standard error in this chapter is meaningless. The assumption is defensible for most series at monthly frequency, but it is an assumption, and it should be made consciously.
A limit theorem licenses nothing at any particular \(n\). This bears repeating because it is the most common misuse of asymptotic theory. eq. 3.8 describes a sequence, and Figure 3.4 shows the very same theorem delivering an excellent approximation at \(n = 10\) for one population and an inadequate one at \(n = 100\) for another. Citing the theorem is not evidence that it applies to your problem. Simulating at your sample size, with your data’s shape, is.
Choosing the sample chooses the answer. The 1971 start date for problem 2 is not really a decision — it is when the NASDAQ Composite began — but it is also close to the start of the period that made technology stocks famous. Run the same test on 1971–1999 and the NASDAQ’s advantage looks overwhelming; run it on 2000–2010 and it reverses. Neither the \(t\)-statistic nor the central limit theorem has anything to say about which window is the right one, because both take the sample as given and ask only what it implies.
This is the deepest of the four difficulties listed in Chapter 1, and it has no statistical solution. History arrives at one year per year, so we cannot order more of it; and reaching further back eventually reaches a different economy, so a longer sample is not straightforwardly a better one. What we can do is be explicit about the window, check whether the conclusion survives moving it, and report honestly when it does not.
Key takeaways
- A \(t\)-statistic is only a distance measured in standard errors. It becomes evidence only once you know what distribution it follows under the null, and supplying that distribution is the entire problem.
- Non-normality is not a technicality in financial data. Monthly equity premia have excess kurtosis above 7 and reject normality at any conceivable level, so a procedure requiring normal returns is requiring something we know is false.
- The exact \(t_{n-1}\) result uses normality three separate times: for the normal numerator, for the chi-square denominator, and for their independence. The third cannot be weakened, because independence of \(\bar{X}\) and \(S^2\) characterises the normal distribution rather than merely following from it.
- Exact results are the strongest statements available when their assumptions hold, and carry no guarantee at all when they do not. For plenty of ordinary questions — two samples with unequal variances, anything involving a Sharpe ratio — no exact result exists at any sample size.
- Asymptotic theory trades an exact statement under implausible assumptions for an approximate one under plausible assumptions. The law of large numbers needs a mean to exist; the central limit theorem needs two moments; neither cares about shape.
- Slutsky’s theorem is what makes estimating the standard error asymptotically free, converting the infeasible CLT into the feasible statistic we actually compute. This is also why the asymptotic argument delivers a normal, not a \(t\).
- Skewness, not heavy tails, governs how good the approximation is, because the \(t\)-statistic is self-normalising and fat tails inflate both halves of the ratio at once. There is no universal minimum sample size — simulate at yours.
- The bootstrap gives a practical way to check: recentre the data to impose the null, resample with replacement, and compare the resulting distribution of the statistic against the normal you were about to rely on.
- At the sample sizes here — 1,188 and 652 months — the approximation is accurate to within about 0.03 on the critical value. The equity premium is positive beyond any argument; the index gap is marginal, and its verdict genuinely depends on the approximation being good, which is why we checked.
- How much care the distribution theory deserves depends on how close the answer sits to the decision boundary. Check it when it could change your mind.
Check your understanding
Five questions. Commit to an answer before opening the solution — the hints are written for the specific misunderstanding behind each wrong option, so they are only useful once you have picked one.
Question 1
A colleague argues: “Monthly equity returns have excess kurtosis above 7, so the central limit theorem cannot be relied on here and we should not be quoting \(t\)-statistics at all.”
What is the best response?
- She is right; heavy tails invalidate the CLT.
- The CLT requires only a finite mean and variance, and Figure 3.4 shows heavy tails are close to harmless for the \(t\)-statistic; skewness is the feature to worry about.
- The CLT requires normality, but with 1,188 observations the departure is negligible.
- Kurtosis affects the estimate of the mean, so the point estimate is biased.
Show the answer
B is correct. Lindeberg–Lévy in eq. 3.7 assumes i.i.d. draws with a finite mean and a finite, non-zero variance. Kurtosis of 7 violates neither. The simulation makes the point empirically: the Student’s \(t\) population has excess kurtosis 6 and its true rejection rate is within a fraction of a percentage point of 5% even at \(n = 10\), because the \(t\)-statistic divides by an estimate of scale that inflates alongside the numerator. The chi-square population, with the same excess kurtosis but skewness of 2, is badly distorted.
If you chose A — this confuses a condition on moments with a condition on shape. Finite variance is a statement about whether \(\mathbb{E}[X^2]\) exists, not about how large the fourth moment is. A distribution can have enormous kurtosis and perfectly finite variance.
If you chose C — this is the error Section 3.3 is built to prevent. Normality is what the exact \(t_{n-1}\) result requires; the CLT is valuable precisely because it does not. If the CLT required normality it would be useless, since under normality the exact result is already available.
If you chose D — kurtosis does not bias \(\bar{X}\). The sample mean is unbiased for \(\mu\) whenever the mean exists, whatever the shape. Kurtosis affects the precision of the estimate and the accuracy of the approximation to its sampling distribution, not its centring.
Question 2
In Figure 3.3, the chi-square population is skewed sharply to the right, but the distribution of its \(t\)-statistic is skewed to the left. Why?
- It is a simulation artefact; more replications would remove it.
- Standardising the population to mean 0 and variance 1 reverses the skew.
- In skewed data \(\bar{X}\) and \(S\) move together, so a large numerator brings a large denominator with it and big positive \(t\)-values are pulled back.
- The \(t\)-distribution is symmetric, so any skew must be negative.
Show the answer
C is correct. With right-skewed data, a sample whose mean is unusually high got that way through a few large observations, and those same observations enlarge \(S\). The denominator of eq. 3.2 therefore grows just when the numerator does, damping large positive values of \(T\). Samples that happen to miss the big observations have both a low mean and a small \(S\), so their \(T\) is free to go strongly negative. The result is a left tail that is longer than the right.
This is also why the one-sided upper-tail test in Table 3.3 under-rejects — the upper tail of \(T\) is thinner than the \(t\)-table thinks, so the critical value is too far out.
If you chose A — the simulation uses 200,000 replications, and the effect is large and systematic rather than noisy. Note also that the normal line lies on the reference density with the same number of replications, which would not happen if the method were producing artefacts.
If you chose B — standardising is a linear transformation: subtracting a constant and dividing by a positive constant leaves skewness exactly unchanged. That is the point of standardising here, since it isolates shape from location and scale.
If you chose D — the \(t_{n-1}\) density is symmetric, but that is the distribution \(T\) would follow if the data were normal. The simulated distribution is the true one for a chi-square population, and it is under no obligation to be symmetric. Confusing the two is the error the whole section is about.
Question 3
Table 3.4 reports a \(t\)-statistic of 2.18 for the paired comparison and 0.90 for the unpaired one, from the same data. Which is the better test, and why?
- The unpaired one, because it makes fewer assumptions about the relationship between the series.
- The paired one, because differencing removes the common market factor and so cuts the standard error by a factor of about 2.4.
- Neither; the disagreement means the data cannot answer the question.
- The unpaired one, because 652 observations of each series is more information than 652 differences.
Show the answer
B is correct. The two indices have a correlation of 0.87, so most of what moves the NASDAQ in a given month also moves the S&P 500. The difference \(R_t^N - R_t^S\) cancels that common component and leaves the relative performance, which is the quantity the question is about. The unpaired standard error implicitly treats the two series as unrelated and so attributes all of the market’s variance to the comparison, inflating the standard error by 2.4 and throwing away the power to detect anything.
Pairing costs nothing in assumptions here: the months line up one-for-one, which is what makes the difference well defined.
If you chose A — the unpaired approach does not make fewer assumptions, it makes a different and false one, namely that the two samples are independent. Assuming independence when the correlation is 0.87 is not conservatism, it is a misspecification, and it happens to be the setting with no exact finite-sample solution.
If you chose C — the tests do not conflict as evidence; one is simply much less efficient. They estimate the same quantity — the point estimate of the mean difference is identical — and differ only in the standard error attached to it.
If you chose D — the number of observations is 652 either way. Pairing does not discard data; it uses the fact that observation \(t\) in one series and observation \(t\) in the other were recorded in the same month, which the unpaired calculation ignores.
Question 4
A risk team has 48 monthly observations of a hedge fund’s returns, with a skewness of \(+1.8\). They test \(H_0: \mu = 0\) against \(H_1: \mu > 0\) at the 5% level using a \(t_{47}\) critical value and fail to reject. What should be said about the true size of their test?
- It is 5%, because they used the correct degrees of freedom.
- It is above 5%, so they are rejecting too often and the failure to reject is reassuring.
- It is below 5%, so their test is conservative in this direction — the failure to reject is weaker evidence for the null than it appears.
- Nothing can be said without knowing the kurtosis.
Show the answer
C is correct. This is the chi-square column of Table 3.3, at a sample size between the \(n = 20\) and \(n = 60\) rows and with skewness close to the simulated value of 2. A one-sided upper-tail test on strongly right-skewed data rejects far less than its nominal rate — under 3% at \(n = 60\) in the simulation. Their test is therefore conservative in the upper tail: it is harder to reject than they think, so failing to reject is less informative than a correctly sized test would make it.
The practical response is to simulate the null distribution at \(n = 48\) using a skewness calibrated to their data, or to bootstrap it from the returns themselves as in Figure 3.5, rather than reading a \(t\)-table.
If you chose A — degrees of freedom are the right adjustment for the sampling error in \(S\), which is a different problem. They do nothing about the shape of the parent distribution, and it is the shape that breaks eq. 3.5.
If you chose B — the direction is backwards, and the asymmetry between tails is the whole point. Right-skewed data produce a left-skewed \(t\)-statistic, so the upper tail is too thin and the lower tail too heavy. A two-sided test on the same data would indeed over-reject.
If you chose D — Table 3.3 shows kurtosis is nearly irrelevant here: the Student’s \(t\) and chi-square populations share an excess kurtosis of 6, and only the skewed one misbehaves. Skewness alone is enough to sign the distortion.
Question 5
Which statement best describes what Figure 3.5 establishes?
- That monthly returns are normally distributed at these sample sizes.
- That the sampling distribution of the \(t\)-statistic, at these sample sizes and for these return distributions, is close enough to \(N(0,1)\) that asymptotic \(p\)-values can be trusted.
- That the central limit theorem is true.
- That the equity premium is positive.
Show the answer
B is correct. The figure is a check on the approximation, not on the theorem and not on the data. It resamples the observed returns — keeping their real kurtosis of 7 and their real skewness — imposes the null, and asks what distribution the \(t\)-statistic actually follows at \(n = 1{,}188\) and \(n = 652\). The answer is: near enough to standard normal that the critical values agree to about 0.03. That is a licence to use asymptotic \(p\)-values for these series at these sample sizes, and nothing more general.
If you chose A — the returns remain violently non-normal; Figure 3.2 and Table 3.2 say so, and the resampling procedure preserves that non-normality by construction. It is the distribution of the statistic that becomes normal, exactly as in Chapter 1’s discussion of the CLT.
If you chose C — a theorem is not established by simulation. What the figure does is quantify the theorem’s approximation error at a specific \(n\), which is the thing the theorem itself is silent about.
If you chose D — that conclusion comes from Table 3.6, and it needs both the point estimate and a trustworthy distribution. Figure 3.5 supplies only the second. Note that the resampling deliberately imposes \(\mu = 0\), so the figure is constructed to say nothing about where the true mean lies.