import matplotlib
matplotlib.use("Agg") # headless: render to buffers, never to a screen
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.formula.api as smf
from statsmodels.stats.diagnostic import linear_reset
rng = np.random.default_rng(42) # one generator, pinned, drives the whole notebook
plt.rcParams["figure.figsize"] = (9, 5)
plt.rcParams["axes.grid"] = True
plt.rcParams["grid.alpha"] = 0.3
print("pandas", pd.__version__)
import statsmodels
print("statsmodels", statsmodels.__version__)nb2.5 — The Biased-Estimator Lab
Week 2, Chapter 2.5. Companion notebook.
Chapters 2.1–2.4 taught you to compute \hat{\beta} correctly and to quantify its wobble with good standard errors. But a standard error tells you how much your estimate bounces if you redrew the sample — it says nothing about whether the number you are bouncing around is the right number. This notebook is about being precisely wrong. We simulate the three classic threats from the chapter and watch each one push \hat{\beta} off target in a way that no amount of extra data can fix:
- Omitted variable bias (OVB) — verify \hat{\tilde\beta}_1 \xrightarrow{p} \beta_1 + \beta_2\delta_1, show the bias does not shrink with N, and confirm the two-sign rule by flipping correlations.
- Measurement error / attenuation — corrupt a regressor and watch \hat\beta \to \lambda\beta with reliability ratio \lambda = \sigma^2_{x^\*}/(\sigma^2_{x^\*}+\sigma^2_m); plot it.
- Mismeasured y — show classical noise in the outcome inflates standard errors but leaves the slope unbiased (the crucial asymmetry).
- Functional form / RESET — fit a line to a quadratic truth, catch it with Ramsey’s RESET test, fix it with a polynomial term.
Then a “Your Turn”: make the omitted variable orthogonal to the regressor and watch the bias vanish. Everything is a controlled experiment: we know the true parameters because we wrote them down, so we can check the estimator against the truth to three decimals.
1. Omitted Variable Bias: \hat{\tilde\beta}_1 \xrightarrow{p} \beta_1 + \beta_2\,\delta_1
The result in one sentence. Leave out a variable that belongs in the model and is correlated with a regressor you kept, and your coefficient absorbs part of the omitted variable’s effect — the contamination is exactly \beta_2\,\delta_1: the omitted variable’s true effect \beta_2 times the slope \delta_1 from regressing the omitted variable on the included one.
We build Maya’s fair-lending world from the chapter, but as a clean simulation. The true (long) model is
y = \beta_0 + \beta_1 x + \beta_2 z + \varepsilon,
with a known causal effect \beta_1. We then run the short regression that drops z, and check its slope against the formula. To control the correlation between x and z exactly, we generate z from x through an auxiliary regression z = \delta_0 + \delta_1 x + \text{noise} — so we know \delta_1 by construction. Here x is the demographic indicator’s continuous analog (think a standardized risk-relevant trait), z is creditworthiness, and \varepsilon is everything else.
# Known truth -- we write the parameters down, then check the estimator against them.
N = 100_000
beta0, beta1, beta2 = 1.0, 0.0, 0.30 # beta1 = 0: ZERO true causal effect of x on y
delta0, delta1 = 0.0, -0.50 # z correlated with x: slope of z on x is -0.50
def simulate_ovb(n, beta1, beta2, delta1, rng):
"""Generate (x, z, y) from the long model y = beta0 + beta1*x + beta2*z + eps,
with z = delta0 + delta1*x + noise so Cov(x,z)/Var(x) = delta1 exactly (in plim)."""
x = rng.normal(0.0, 1.0, size=n)
z = delta0 + delta1 * x + rng.normal(0.0, 1.0, size=n)
eps = rng.normal(0.0, 1.0, size=n)
y = beta0 + beta1 * x + beta2 * z + eps
return pd.DataFrame({"y": y, "x": x, "z": z})
df = simulate_ovb(N, beta1, beta2, delta1, rng)
df.head()1.1 Long regression (clean) vs. short regression (contaminated)
The long regression includes z and should recover the true \beta_1 = 0. The short regression drops z and should land on \beta_1 + \beta_2\delta_1 = 0 + (0.30)(-0.50) = -0.15 — a 15-point penalty against group x in a world with literally zero discrimination. The entire “effect” is borrowed from creditworthiness through the correlation.
long_fit = smf.ols("y ~ x + z", data=df).fit()
short_fit = smf.ols("y ~ x", data=df).fit()
b1_long = long_fit.params["x"]
b1_short = short_fit.params["x"]
# delta1 measured directly from the auxiliary regression z ~ x
aux_fit = smf.ols("z ~ x", data=df).fit()
delta1_hat = aux_fit.params["x"]
predicted_short = beta1 + beta2 * delta1 # the OVB formula, true parameters
predicted_hat = b1_long + long_fit.params["z"] * delta1_hat # formula from estimates
print(f"True beta1 (causal) : {beta1:+.4f}")
print(f"Long-regression beta1_hat (x|z) : {b1_long:+.4f} <- recovers the truth")
print(f"Short-regression beta1_tilde (x) : {b1_short:+.4f} <- contaminated")
print()
print(f"OVB formula beta1 + beta2*delta1 : {predicted_short:+.4f}")
print(f"Short estimate : {b1_short:+.4f}")
print(f"Match (abs diff) : {abs(b1_short - predicted_short):.4f}")# The short slope lands on the formula, not on the truth.
assert abs(b1_short - predicted_short) < 0.02, "short slope should match beta1 + beta2*delta1"
assert abs(b1_long - beta1) < 0.02, "long slope should recover the true beta1"
print("CHECK PASSED: simulated OVB matches beta2 * delta1 to within Monte Carlo noise.")1.2 More data does not help: precisely wrong
The single most important thing about bias: it lives in the probability limit, not in the sampling noise. So as N grows, the short slope does not drift toward the truth — it homes in tighter on the wrong number -0.15. We sweep N from a few hundred to a million and watch the short estimate clamp onto the biased target while its standard error shrinks. A biased estimator with a million observations is a precisely wrong number.
Ns = [200, 1_000, 5_000, 25_000, 100_000, 1_000_000]
rows = []
for n in Ns:
d = simulate_ovb(n, beta1, beta2, delta1, rng)
f = smf.ols("y ~ x", data=d).fit()
rows.append({
"N": n,
"beta1_short": f.params["x"],
"SE": f.bse["x"],
"target (beta1+beta2*delta1)": predicted_short,
})
sweep = pd.DataFrame(rows)
print(sweep.to_string(index=False, float_format=lambda v: f"{v:.4f}"))fig, ax = plt.subplots()
ax.errorbar(sweep["N"], sweep["beta1_short"], yerr=1.96 * sweep["SE"],
fmt="o-", capsize=4, label=r"short $\hat{\tilde\beta}_1$ (95% CI)")
ax.axhline(predicted_short, color="crimson", ls="--",
label=r"biased target $\beta_1+\beta_2\delta_1=-0.15$")
ax.axhline(beta1, color="green", ls=":", label=r"true $\beta_1=0$")
ax.set_xscale("log")
ax.set_xlabel("sample size N (log scale)")
ax.set_ylabel(r"short-regression slope on $x$")
ax.set_title("More data -> precisely wrong: the CI tightens onto the biased target, not the truth")
ax.legend()
fig.tight_layout()
print("Estimate converges to -0.15, SE -> 0. The bias never shrinks.")1.3 The two-sign rule: flip the signs, flip the bias
The bias \beta_2\delta_1 is a product of two signs. “Same signs push up, opposite signs push down.” We verify all four cells of the chapter’s table by flipping the sign of \beta_2 (effect of the omitted variable on y) and \delta_1 (correlation of omitted with included), holding |\beta_2|=0.30 and |\delta_1|=0.50 fixed so the magnitude of the bias is always 0.15.
signs = []
for b2 in (+0.30, -0.30):
for d1 in (+0.50, -0.50):
d = simulate_ovb(200_000, beta1, b2, d1, rng)
f = smf.ols("y ~ x", data=d).fit()
bias = b2 * d1
direction = "up (too +)" if bias > 0 else "down (too -)"
signs.append({
"sign(beta2)": "+" if b2 > 0 else "-",
"sign(delta1)": "+" if d1 > 0 else "-",
"bias = beta2*delta1": bias,
"beta1_short": f.params["x"],
"pushes": direction,
})
signtab = pd.DataFrame(signs)
print(signtab.to_string(index=False, float_format=lambda v: f"{v:+.4f}"))
# Verify each simulated short slope matches its predicted bias (beta1 = 0 here).
for r in signs:
assert abs(r["beta1_short"] - r["bias = beta2*delta1"]) < 0.02
print("\nCHECK PASSED: all four sign combinations match the two-sign rule.")2. Measurement Error: Attenuation Toward Zero
The result in one sentence. Regress on a noisy copy of the true variable — classical errors-in-variables, x = x^\* + m with random static m — and the OLS slope is pulled toward zero by the reliability ratio \lambda = \sigma^2_{x^\*}/(\sigma^2_{x^\*}+\sigma^2_m), the share of the regressor’s variance that is genuine signal: \hat\beta_1 \xrightarrow{p} \lambda\,\beta_1.
Now there is no omitted variable — the model is perfectly specified in the true regressor x^\*. The only sin is that we measure x^\* with independent noise. We set a known \beta_1 = 2, fix \sigma^2_{x^\*}=1, and crank up the noise SD \sigma_m. Because the noise leaks into the composite error u = \varepsilon - \beta_1 m and sits inside the regressor x = x^\* + m, regressor and error become mechanically correlated, and the slope attenuates.
beta1_true = 2.0
N_att = 200_000
x_star = rng.normal(0.0, 1.0, size=N_att) # Var(x*) = 1
eps = rng.normal(0.0, 1.0, size=N_att)
y_att = 1.0 + beta1_true * x_star + eps # well-specified in x*
rows = []
for sigma_m in [0.0, 0.5, 1.0, 1.5, 2.0, 3.0]:
m = rng.normal(0.0, sigma_m, size=N_att)
x_obs = x_star + m # classical EIV
fit = smf.ols("y ~ x", data=pd.DataFrame({"y": y_att, "x": x_obs})).fit()
lam = 1.0 / (1.0 + sigma_m**2) # Var(x*)/(Var(x*)+Var(m)), since Var(x*)=1
rows.append({
"sigma_m": sigma_m,
"lambda (reliability)": lam,
"beta_hat": fit.params["x"],
"predicted lambda*beta": lam * beta1_true,
})
att = pd.DataFrame(rows)
print(att.to_string(index=False, float_format=lambda v: f"{v:.4f}"))# Each estimate should sit on lambda * beta1_true.
err = (att["beta_hat"] - att["predicted lambda*beta"]).abs().max()
print(f"\nMax abs deviation from lambda*beta : {err:.4f}")
assert err < 0.02, "attenuated slope must match lambda * beta1"
print("CHECK PASSED: attenuated slope matches lambda * beta to within Monte Carlo noise.")fig, ax = plt.subplots()
ax.plot(att["lambda (reliability)"], att["beta_hat"], "o-",
label=r"OLS $\hat\beta_1$ vs reliability")
# straight line from origin (lambda=0 -> 0) to (lambda=1 -> beta1_true)
lam_grid = np.linspace(0, 1, 50)
ax.plot(lam_grid, beta1_true * lam_grid, "--", color="crimson",
label=r"prediction $\hat\beta_1=\lambda\beta_1$")
ax.axhline(beta1_true, color="green", ls=":", label=r"true $\beta_1=2$")
ax.set_xlabel(r"reliability ratio $\lambda=\sigma^2_{x^*}/(\sigma^2_{x^*}+\sigma^2_m)$")
ax.set_ylabel(r"estimated slope $\hat\beta_1$")
ax.set_title("Attenuation: more noise (lower lambda) shrinks the slope linearly toward zero")
ax.legend()
fig.tight_layout()
print("As noise rises, lambda falls, and the slope slides from 2 toward 0 along the predicted line.")3. Mismeasured y Inflates SEs but Does Not Bias the Slope
The crucial asymmetry: classical noise in the regressor biases the slope (Section 2), but classical noise in the outcome does not. If y = y^\* + v with v random static uncorrelated with x, then y = \beta_0 + \beta_1 x + (\varepsilon + v) — the noise just folds into the equation error. Since v is uncorrelated with x, the zero-conditional-mean assumption survives and \hat\beta_1 stays consistent. The only cost is a larger error variance, which widens the standard error. We add rising noise to y (the regressor x is now measured perfectly) and watch the slope hold steady while the SE balloons.
N_y = 50_000
x_clean = rng.normal(0.0, 1.0, size=N_y) # regressor measured perfectly
eps_y = rng.normal(0.0, 1.0, size=N_y)
y_star = 1.0 + beta1_true * x_clean + eps_y # true outcome
rows = []
for sigma_v in [0.0, 1.0, 2.0, 4.0]:
v = rng.normal(0.0, sigma_v, size=N_y)
y_obs = y_star + v # classical noise in the OUTCOME
fit = smf.ols("y ~ x", data=pd.DataFrame({"y": y_obs, "x": x_clean})).fit()
rows.append({
"sigma_v (noise in y)": sigma_v,
"beta_hat": fit.params["x"],
"SE(beta_hat)": fit.bse["x"],
})
ytab = pd.DataFrame(rows)
print(ytab.to_string(index=False, float_format=lambda v: f"{v:.4f}"))
# Slope stays at ~2 regardless of outcome noise; SE grows.
assert (ytab["beta_hat"] - beta1_true).abs().max() < 0.05, "slope must stay unbiased"
assert ytab["SE(beta_hat)"].is_monotonic_increasing, "SE must grow with outcome noise"
print("\nCHECK PASSED: slope ~ 2.0 throughout (unbiased); SE rises monotonically (precision cost only).")4. Functional Form and the RESET Test
The result in one sentence. Fit a straight line to a curved truth, and the neglected curvature lands in the error term, correlates with x, and the reported slope describes a line that exists nowhere. Ramsey’s RESET test catches it by checking whether powers of the fitted values (\hat y^2, \hat y^3) add explanatory power; if they do, your shape is wrong.
We build Priya’s convex insurance world: the true model is quadratic, y = \beta_0 + \beta_1 x + \beta_2 x^2 + \varepsilon with \beta_2 > 0. We fit a (wrong) line, run RESET, then add the x^2 term and confirm RESET stops rejecting.
Confirmed statsmodels 0.14.6 signature: linear_reset(res, power=3, test_type='fitted', use_f=False, cov_type='nonrobust', cov_kwargs=None). The first positional argument is a fitted RegressionResults object; test_type='fitted' (the default) regresses on powers of the fitted values, which is exactly Ramsey’s procedure; use_f=True gives the F-form. We use these defaults explicitly.
x_ff = rng.uniform(0, 5, 4000)
y_ff = 1 + 0.5 * x_ff + 0.4 * x_ff**2 + rng.normal(0, 1, 4000) # true beta2 = 0.4 > 0
ff = pd.DataFrame({"y": y_ff, "x": x_ff})
linear_fit = smf.ols("y ~ x", data=ff).fit()
reset_lin = linear_reset(linear_fit, power=3, test_type="fitted", use_f=True) # adds yhat^2, yhat^3
print(f"Linear fit: RESET F = {reset_lin.fvalue:.1f}, p = {reset_lin.pvalue:.2e}")
print("Tiny p-value -> REJECT linearity: there is curvature the straight line is missing.")quad_fit = smf.ols("y ~ x + I(x**2)", data=ff).fit()
reset_quad = linear_reset(quad_fit, power=3, test_type="fitted", use_f=True)
print(f"Quadratic fit: RESET F = {reset_quad.fvalue:.2f}, p = {reset_quad.pvalue:.2f}")
print("Large p-value -> FAIL to reject: the shape now looks right.\n")
print(f"True (beta1, beta2) = (0.50, 0.40)")
print(f"Quadratic estimates (x, x^2)= ({quad_fit.params['x']:.3f}, {quad_fit.params['I(x ** 2)']:.3f})")
assert reset_lin.pvalue < 0.01, "linear model should fail RESET"
assert reset_quad.pvalue > 0.05, "quadratic model should pass RESET"
print("\nCHECK PASSED: RESET rejects the line, passes the quadratic; estimates recover the truth.")order = np.argsort(x_ff)
fig, ax = plt.subplots()
ax.scatter(x_ff, y_ff, s=4, alpha=0.15, color="gray", label="data (convex truth)")
ax.plot(x_ff[order], linear_fit.fittedvalues.to_numpy()[order], color="crimson", lw=2,
label="linear fit (wrong shape)")
ax.plot(x_ff[order], quad_fit.fittedvalues.to_numpy()[order], color="green", lw=2,
label="quadratic fit (right shape)")
ax.set_xlabel("emissions intensity x")
ax.set_ylabel("insurance premium y")
ax.set_title("A straight line through a convex world: understates high-x, overstates low-x")
ax.legend()
fig.tight_layout()
print("The line cuts across the curve; the quadratic tracks it. RESET told us which to trust.")Your Turn — Make the Omitted Variable Orthogonal and Watch the Bias Vanish
The OVB formula \text{bias} = \beta_2\,\delta_1 is a product: it is poison only when both factors are nonzero — the omitted variable must matter for y (\beta_2\neq 0) and be correlated with your regressor (\delta_1\neq 0). Below we keep the omitted variable just as relevant as before (\beta_2 = 0.30, it still drives y) but make it orthogonal to x by setting \delta_1 = 0. Now z is generated independently of x. The short regression should recover the true \beta_1 = 0 — the bias term collapses to \beta_2\cdot 0 = 0 — and the only remaining cost of omitting z is a larger error variance (a precision hit), not bias.
This is the same lesson as the bias–consistency ledger: a relevant variable uncorrelated with x lives in the SE column, not the bias column.
# delta1 = 0 -> z is orthogonal to x (independent draw), but still affects y (beta2 = 0.30).
df_orth = simulate_ovb(200_000, beta1=0.0, beta2=0.30, delta1=0.0, rng=rng)
short_orth = smf.ols("y ~ x", data=df_orth).fit()
long_orth = smf.ols("y ~ x + z", data=df_orth).fit()
print(f"True beta1 : {0.0:+.4f}")
print(f"Short slope on x (z omitted, delta1=0): {short_orth.params['x']:+.4f} <- bias has vanished")
print(f"Long slope on x (z included) : {long_orth.params['x']:+.4f}")
print()
print(f"Short-model residual SD (z omitted) : {np.std(short_orth.resid):.4f} <- larger: lost precision")
print(f"Long-model residual SD (z included) : {np.std(long_orth.resid):.4f}")
assert abs(short_orth.params["x"] - 0.0) < 0.02, "with delta1=0 the bias should vanish"
print("\nCHECK PASSED: orthogonal omitted variable -> zero bias (only a precision cost).")
# TRY NEXT: set delta1 = 0 but beta2 = 0 too (irrelevant variable) -- still no bias.
# Then set BOTH nonzero again and confirm the bias returns. Bias needs BOTH slopes.