import numpy as np
import pandas as pd
import matplotlib
matplotlib.use("Agg") # non-interactive backend: render to buffers, never pop a window
import matplotlib.pyplot as plt
from scipy import stats
rng = np.random.default_rng(42) # one seed for the whole notebook -> reproducible
# Sam's world, in PERCENT units (so 0.08 means 0.08% per day).
TRUE_MU_NULL = 0.0 # the skeptic's claim: the rule has no edge
SIGMA = 0.75 # daily return volatility, in percent
N_DAYS = 252 # one trading year
ALPHA = 0.05 # significance level we commit to in advance
print(f"sigma = {SIGMA}%/day, N = {N_DAYS} days, alpha = {ALPHA}")Notebook 1.5 — Power Curves and the t-Test From Scratch
Sam has a trading rule. Each morning, buy yesterday’s biggest losers, hold one day, sell at the close — a bet on short-term reversal. Over 252 trading days Sam measures an average daily return of \bar{x} = 0.08\%, which annualizes to more than 20\%. The returns are noisy, though: some days +1.5\%, some days -1.2\%. The whole question of Chapter 1.5 was: could a rule whose true edge is exactly zero produce +0.08\% over 252 days, just by luck?
In this notebook we build the machine that answers it — the one-sample t-test — by hand, check our arithmetic against scipy, and then do the thing a formula can never do for you: simulate it. We will watch a true null get rejected almost exactly 5\% of the time (that is size), watch a real effect get caught some fraction of the time (that is power), draw the power curves, see a confidence interval fall out as an inverted test, and finally let Sam test twenty worthless rules at once and count the false positives.
1. The t-test by hand
The recipe from the chapter, in three lines of algebra. Given a sample x_1,\dots,x_N:
\bar{x}=\frac{1}{N}\sum_i x_i,\qquad s=\sqrt{\frac{1}{N-1}\sum_i (x_i-\bar{x})^2},\qquad t=\frac{\bar{x}-\mu_0}{s/\sqrt{N}}.
The denominator s/\sqrt{N} is the estimated standard error of the mean. The p-value is then a tail area of Student’s t distribution with N-1 degrees of freedom — which tail, and whether you double it, depends on the alternative. Let us write that as one function and own every step.
def t_test_by_hand(x, mu0=0.0, alternative="greater"):
'''One-sample t-test computed from scratch.
alternative: 'greater' -> H1: mu > mu0 (upper tail; Sam's reversal rule)
'less' -> H1: mu < mu0 (lower tail)
'two-sided'-> H1: mu != mu0 (both tails; Maya's neighborhood gap)
Returns (t_stat, p_value, df).
'''
x = np.asarray(x, dtype=float)
n = x.size
xbar = x.mean()
s = x.std(ddof=1) # ddof=1 -> divide by N-1 (Bessel correction)
se = s / np.sqrt(n) # estimated standard error of the mean
t_stat = (xbar - mu0) / se
df = n - 1
if alternative == "greater":
p = stats.t.sf(t_stat, df) # sf = 1 - cdf = upper-tail area
elif alternative == "less":
p = stats.t.cdf(t_stat, df) # lower-tail area
elif alternative == "two-sided":
p = 2.0 * stats.t.sf(np.abs(t_stat), df) # both tails
else:
raise ValueError("alternative must be 'greater', 'less', or 'two-sided'")
return t_stat, p, dfSam’s numbers
Rather than hard-code \bar{x}=0.08 and s=0.75, let us generate a sample that lands near those values so the rest of the notebook works on a real array of returns. We draw 252 daily returns from a normal with a small true edge and the stated volatility. (Returns are not really normal — they are fat-tailed — but with N=252 the CLT makes the t-test a fine approximation, exactly as the chapter argued.)
# Generate Sam's observed year, then affine-rescale it to land EXACTLY on the chapter's
# headline numbers (x-bar = 0.08%, s = 0.75%) so this notebook matches the prose. This is a
# legitimate teaching move: it fixes the sample's mean and spread while keeping a real,
# noisy array of 252 returns to run the test on.
raw = rng.normal(loc=0.08, scale=SIGMA, size=N_DAYS)
raw = (raw - raw.mean()) / raw.std(ddof=1) # standardize to mean 0, sd 1
sam_returns = 0.08 + 0.75 * raw # now mean=0.08%, s=0.75% exactly
xbar = sam_returns.mean()
s = sam_returns.std(ddof=1)
se = s / np.sqrt(N_DAYS)
print(f"x-bar = {xbar:.4f}% s = {s:.4f}% se = {se:.4f}%")
t_hand, p_hand, df = t_test_by_hand(sam_returns, mu0=0.0, alternative="greater")
print(f"\nBy hand (one-sided, H1: mu > 0):")
print(f" t = {t_hand:.4f} df = {df} p = {p_hand:.4f}")
print(f" Reject H0 at alpha={ALPHA}? {'YES' if p_hand < ALPHA else 'no'}")2. Validate against scipy.stats.ttest_1samp
A hand-rolled test you have not checked is a hand-rolled test you should not trust. scipy reports a two-sided p-value by default; for a one-sided alternative you halve it when the estimate is on the alternative’s side. Newer scipy also accepts alternative= directly, so we check both routes. If our by-hand t and p do not match to machine precision, we have a bug — and we assert exactly that.
res = stats.ttest_1samp(sam_returns, popmean=0.0, alternative="greater")
print(f"scipy (one-sided): t = {res.statistic:.4f} p = {res.pvalue:.4f}")
print(f"by hand : t = {t_hand:.4f} p = {p_hand:.4f}")
# The numbers must agree to ~15 digits, not just 'look close'.
assert np.isclose(t_hand, res.statistic), "t-stat mismatch!"
assert np.isclose(p_hand, res.pvalue), "p-value mismatch!"
# Also confirm the two-sided route: scipy two-sided p == 2 * our one-tail (since xbar>0).
res2 = stats.ttest_1samp(sam_returns, popmean=0.0) # default is two-sided
t_2s, p_2s, _ = t_test_by_hand(sam_returns, mu0=0.0, alternative="two-sided")
assert np.isclose(p_2s, res2.pvalue), "two-sided p mismatch!"
print(f"\nTwo-sided: by hand p = {p_2s:.4f}, scipy p = {res2.pvalue:.4f} (match)")
print("All assertions passed -> our by-hand t-test is correct.")3. Simulating size: does a true null really reject only 5% of the time?
The whole promise of “\alpha = 0.05” is a frequency claim about a rule, not about Sam’s one dataset: if the null is true, the test should wrongly reject at most 5\% of the time. That is checkable. We will conjure a world where the rule is genuinely worthless (\mu = 0), run the test thousands of times on fresh fake years, and count how often it cries victory. If our machine is calibrated, that count should sit right at 5\%. This is the test’s size.
def simulate_rejection_rate(true_mu, n, n_sims, alpha=ALPHA,
sigma=SIGMA, alternative="greater", seed_rng=None):
'''Fraction of simulated samples where the t-test rejects H0: mu=0.
With true_mu = 0 this estimates SIZE; with true_mu != 0 it estimates POWER.
'''
r = seed_rng if seed_rng is not None else np.random.default_rng(0)
rejects = 0
for _ in range(n_sims):
sample = r.normal(loc=true_mu, scale=sigma, size=n)
_, p, _ = t_test_by_hand(sample, mu0=0.0, alternative=alternative)
if p < alpha:
rejects += 1
return rejects / n_sims
N_SIMS = 20_000
size_hat = simulate_rejection_rate(true_mu=0.0, n=N_DAYS, n_sims=N_SIMS,
alternative="greater", seed_rng=rng)
# Monte Carlo standard error of a proportion: sqrt(p(1-p)/n_sims)
mc_se = np.sqrt(size_hat * (1 - size_hat) / N_SIMS)
print(f"Empirical size over {N_SIMS:,} sims: {size_hat:.4f} (+/- {2*mc_se:.4f})")
print(f"Target alpha: {ALPHA} -> empirical size is within ~2 MC std errors of 0.05.")
assert abs(size_hat - ALPHA) < 0.01, "Test is not calibrated!"The empirical size lands on 0.05 to within Monte Carlo noise. The test is calibrated: it keeps its promise. One subtle, beautiful consequence from the chapter — when the null is true, the p-value is uniformly distributed on [0,1]. Every value is equally likely, which is why exactly 5\% of true-null p-values fall below 0.05. Let us see that uniformity directly, because it is the seed of the multiple-testing problem at the end.
# Collect many p-values under the TRUE null and histogram them.
null_pvals = np.empty(N_SIMS)
for i in range(N_SIMS):
sample = rng.normal(loc=0.0, scale=SIGMA, size=N_DAYS)
_, null_pvals[i], _ = t_test_by_hand(sample, mu0=0.0, alternative="greater")
fig, ax = plt.subplots(figsize=(7, 4))
ax.hist(null_pvals, bins=20, range=(0, 1), color="steelblue",
edgecolor="white", density=True)
ax.axhline(1.0, color="crimson", linestyle="--", label="Uniform(0,1) density = 1")
ax.set_xlabel("p-value under the true null")
ax.set_ylabel("density")
ax.set_title("Under H0, p-values are uniform — that's why 5% land below 0.05")
ax.legend()
fig.tight_layout()
print(f"Fraction of null p-values below 0.05: {(null_pvals < 0.05).mean():.4f}")4. Simulating power: catching a real effect
Now flip the world on. Suppose the reversal rule genuinely earns \mu = 0.05\% per day. Power is the probability the test catches it — rejects H_0 — given that H_1 is true. The chapter computed this analytically: the true mean sits 0.05/0.0473 \approx 1.06 standard errors above zero, the rule fires at t>1.65, so power \approx \Pr(Z>1.65-1.06) \approx 0.28. Low power: a real edge this small slips past nearly three times in four. Let us confirm the \approx 0.28 by brute-force simulation, then compare to the formula.
TRUE_EDGE = 0.05 # percent per day -- a small but real edge
power_hat = simulate_rejection_rate(true_mu=TRUE_EDGE, n=N_DAYS, n_sims=N_SIMS,
alternative="greater", seed_rng=rng)
# Analytic (large-N normal) approximation for comparison.
se_true = SIGMA / np.sqrt(N_DAYS)
ncp = TRUE_EDGE / se_true # how many SEs the true mean sits above 0
z_crit = stats.norm.ppf(1 - ALPHA) # one-sided critical value ~ 1.645
power_analytic = stats.norm.sf(z_crit - ncp) # P(Z > z_crit - ncp)
print(f"Simulated power at mu={TRUE_EDGE}%: {power_hat:.4f}")
print(f"Analytic power (normal approx): {power_analytic:.4f}")
print(f"True mean sits {ncp:.2f} standard errors above zero.")5. The power curve vs. effect size
Power is not one number — it is a curve. Hold N = 252 fixed and sweep the true edge from 0 (where power must equal \alpha, since that is the null) up to a clearly detectable size. At \mu = 0 the curve passes through 0.05; it then climbs toward 1 as the effect grows. We overlay the analytic curve so you can see the simulation tracing it.
effect_grid = np.linspace(0.0, 0.20, 11) # true edges from 0% to 0.20%/day
sim_power = np.array([
simulate_rejection_rate(true_mu=mu, n=N_DAYS, n_sims=5_000,
alternative="greater", seed_rng=rng)
for mu in effect_grid
])
ana_power = stats.norm.sf(z_crit - effect_grid / se_true)
fig, ax = plt.subplots(figsize=(7, 4.5))
ax.plot(effect_grid, sim_power, "o-", color="darkorange", label="simulated")
ax.plot(effect_grid, ana_power, "--", color="navy", label="analytic (normal approx)")
ax.axhline(ALPHA, color="grey", linestyle=":", label=f"alpha = {ALPHA}")
ax.axhline(0.80, color="green", linestyle=":", label="80% power")
ax.set_xlabel("true daily edge mu (percent)")
ax.set_ylabel("power = P(reject H0)")
ax.set_title(f"Power curve vs. effect size (N={N_DAYS}, one-sided, alpha={ALPHA})")
ax.set_ylim(0, 1.02)
ax.legend()
fig.tight_layout()
print("At mu=0 the curve sits at alpha; it climbs toward 1 as the edge grows.")6. The power curve vs. sample size
The lever Sam actually controls is N. Fix the true edge at 0.05\%/day and sweep the number of trading days. Because the standard error shrinks like 1/\sqrt{N}, power rises — but slowly, and that slowness is the whole lesson. To reach 80\% power for this tiny edge, Sam needs years of data. We mark the 80\% line and read off roughly where it is crossed.
N_grid = np.array([63, 126, 252, 504, 756, 1008, 1512, 2016, 2520])
sim_power_N = np.array([
simulate_rejection_rate(true_mu=TRUE_EDGE, n=n, n_sims=5_000,
alternative="greater", seed_rng=rng)
for n in N_grid
])
ana_power_N = stats.norm.sf(z_crit - TRUE_EDGE / (SIGMA / np.sqrt(N_grid)))
fig, ax = plt.subplots(figsize=(7, 4.5))
ax.plot(N_grid, sim_power_N, "o-", color="darkorange", label="simulated")
ax.plot(N_grid, ana_power_N, "--", color="navy", label="analytic")
ax.axhline(0.80, color="green", linestyle=":", label="80% power")
ax.set_xlabel("sample size N (trading days)")
ax.set_ylabel("power")
ax.set_title(f"Power vs. N (true edge = {TRUE_EDGE}%/day, one-sided, alpha={ALPHA})")
ax.set_ylim(0, 1.02)
ax.legend()
fig.tight_layout()
print("Power climbs with N, but reaching 80% for this tiny edge takes thousands of days.")7. A confidence interval is an inverted test
The chapter’s cleanest idea: a 95\% confidence interval is exactly the set of null values \mu_0 you could not reject with a two-sided test. So we can build the CI two ways and watch them agree. (1) The textbook formula \bar{x}\pm t^\*\,\widehat{\operatorname{se}}. (2) The inversion: scan candidate \mu_0 values and keep the ones the two-sided test fails to reject. The endpoints of the surviving set are the CI endpoints.
# (1) Formula CI.
t_star = stats.t.ppf(1 - ALPHA / 2, df=N_DAYS - 1) # two-sided critical value
ci_lo = xbar - t_star * se
ci_hi = xbar + t_star * se
print(f"Formula 95% CI for mu: [{ci_lo:.4f}%, {ci_hi:.4f}%]")
# (2) Inversion CI: which mu0 survive a two-sided test at alpha?
candidates = np.linspace(-0.30, 0.40, 14001)
survives = np.array([
t_test_by_hand(sam_returns, mu0=m, alternative="two-sided")[1] >= ALPHA
for m in candidates
])
kept = candidates[survives]
print(f"Inversion 95% CI: [{kept.min():.4f}%, {kept.max():.4f}%]")
# They should match to grid resolution.
assert abs(kept.min() - ci_lo) < 0.01 and abs(kept.max() - ci_hi) < 0.01
zero_in = ci_lo <= 0.0 <= ci_hi
print(f"\nIs 0 inside the 95% CI? {zero_in} "
f"-> two-sided test {'does NOT reject' if zero_in else 'rejects'} at alpha={ALPHA}.")
print("The two constructions agree: a CI is a test, inverted.")8. The multiple-testing peek
Change one detail of Sam’s story: suppose Sam tested not one rule but twenty — reversal, momentum, day-of-week, moon phase, the lot — and reported only the best. If all twenty are truly worthless, what is the chance at least one looks “significant” at \alpha = 0.05 by pure luck? The chapter’s arithmetic: 1 - 0.95^{20} \approx 0.64. Better than a coin flip. We simulate it: run 20 independent worthless rules, a thousand experiments over, and count how often at least one false positive appears — then compare to the formula.
N_RULES = 20
N_EXPERIMENTS = 5_000
false_positive_counts = np.empty(N_EXPERIMENTS, dtype=int)
for e in range(N_EXPERIMENTS):
# 20 worthless rules, each a fresh year of pure-noise returns (true mu = 0).
n_signif = 0
for _ in range(N_RULES):
sample = rng.normal(loc=0.0, scale=SIGMA, size=N_DAYS)
_, p, _ = t_test_by_hand(sample, mu0=0.0, alternative="greater")
if p < ALPHA:
n_signif += 1
false_positive_counts[e] = n_signif
any_fp_rate = (false_positive_counts >= 1).mean()
formula = 1 - (1 - ALPHA) ** N_RULES
print(f"P(at least one 'significant' rule), simulated: {any_fp_rate:.4f}")
print(f"P(at least one), formula 1 - 0.95^20: {formula:.4f}")
print(f"Mean number of false 'discoveries' per experiment: "
f"{false_positive_counts.mean():.3f} (expected ~{N_RULES*ALPHA})")
print("\nTest 20 dead rules and cherry-pick the winner -> a spurious 'edge' is the norm.")Your Turn
You now have the whole machine. Three extensions, in rising difficulty:
Find the N for 80% power. Write a function
n_for_power(target_power, true_mu, sigma, alpha)that uses the analytic normal approximation to solve for the smallest N reachingtarget_power(hint: invert \text{power}=\Phi\!\big(\tfrac{\mu\sqrt{N}}{\sigma}-z_\alpha\big) for N), then verify it by simulation. For Sam’s \mu=0.05\%, \sigma=0.75\%, \alpha=0.05, you should land near N\approx 1{,}390 days — about five and a half years (this uses the one-sided z_\alpha=1.645, matching Section 4; with a two-sided z_{\alpha/2}=1.96 the requirement rises to N\approx 1{,}760 days, about seven years).Maya, two-sided. Maya compares loan-approval rates in two neighborhoods; her null is a difference of zero and either direction is a finding. Redraw the power-vs-effect curve with
alternative="two-sided"and explain why, at every nonzero effect, two-sided power is a bit lower than the one-sided power you plotted in Section 5.Bonferroni rescue. In the Section 8 simulation, re-run with each of the 20 tests judged at \alpha/20 = 0.0025 instead of 0.05. Confirm that the chance of any false positive drops back near 5\% — the Bonferroni correction restoring the family-wide error rate.