import matplotlib
matplotlib.use("Agg") # headless: render to buffers, never to a screen
import numpy as np
import matplotlib.pyplot as plt
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.3nb2.2 — Gauss–Markov by Monte Carlo: Seeing “BLUE”
Week 2, Chapter 2.2. Companion notebook.
Chapter 2.2 proved two things about Maya’s loan regression — that OLS is unbiased (it aims at the truth) and that its variance is exactly \sigma^2(\mathbf{X}'\mathbf{X})^{-1} — and then asserted a third: the Gauss–Markov theorem, which says OLS is BLUE, the Best Linear Unbiased Estimator. “Best” means smallest variance among all linear unbiased recipes. Proofs are convincing, but there is a different kind of conviction that comes from watching a theorem hold. That is what this notebook is for.
We fake Maya’s universe — a world where we secretly know the true \boldsymbol{\beta} — and then play God: we draw thousands of alternate datasets and, on each one, re-run the estimators. The histogram of the estimates across those alternate universes is the sampling distribution the theorems are statements about. We will:
- Verify unbiasedness. Average \hat{\boldsymbol{\beta}} over many samples; watch it land on the truth.
- Verify the variance formula. Compare the Monte Carlo spread of \hat{\beta}_1 to \sigma^2(\mathbf{X}'\mathbf{X})^{-1}.
- See efficiency (BLUE). Pit OLS against a rival linear unbiased estimator and watch OLS win on variance.
- Break homoskedasticity. Watch OLS stay unbiased but lose its efficiency crown, and watch the classical SE go wrong — teeing up Ch 2.4.
Then a Your Turn that breaks the assumption that matters most.
1. Maya’s universe, by decree
In the real world Maya never sees the true \beta_1 — that is the whole point of Chapter 1.3. But we are running the simulation, so we get to set it. We decree the true model
\text{rate}_i = \beta_0 + \beta_1\,\text{income}_i + \varepsilon_i, \qquad \beta_0 = 12,\ \beta_1 = -0.42,
where income is measured in tens of thousands of dollars, the rate is an APR in percentage points, and the errors \varepsilon_i are mean-zero with standard deviation \sigma = 1.5 points. So each extra $10,000 of income truly lowers the rate by 0.42 points — the number Maya estimated in the chapter, now promoted to ground truth.
Two design decisions matter. First, we fix the regressor matrix \mathbf{X} once and reuse it in every simulated dataset. This is exactly the “condition on \mathbf{X}” move from the proofs: the only thing we let vary across alternate universes is the random error \boldsymbol{\varepsilon}. Second, we draw incomes with healthy spread — that spread is what gives the slope something to grip, as the \sigma^2/\sum(x_i-\bar{x})^2 formula warned.
beta_true = np.array([12.0, -0.42]) # [intercept, slope]
sigma = 1.5 # error std, in percentage points
N = 400 # loans per simulated dataset
# Fix the regressors ONCE (we condition on X). Incomes spread around 6 (=$60k).
income = rng.normal(6.0, 3.0, size=N)
X = np.column_stack([np.ones(N), income]) # N x 2 design matrix, full rank
# Pre-compute the pieces that never change across universes.
XtX_inv = np.linalg.inv(X.T @ X)
print("X shape:", X.shape)
print("income spread (sd):", round(income.std(ddof=1), 3))
print("(X'X)^{-1} =\n", XtX_inv)2. One draw, then ten thousand
Here is the engine. draw_dataset makes one fresh error vector and returns the simulated y; ols solves the normal equations \mathbf{X}'\mathbf{X}\,\hat{\boldsymbol\beta} = \mathbf{X}'\mathbf{y} without ever forming the inverse explicitly (numerically the right habit). We use np.linalg.solve rather than inverting because it is both faster and more stable.
Run it once to see a single \hat{\boldsymbol\beta} — it will be close to (12, -0.42) but not exact, because this is one noisy sample. The magic is in the averaging.
def draw_dataset(rng):
eps = rng.normal(0.0, sigma, size=N) # homoskedastic, mean-zero, independent
return X @ beta_true + eps # the TRUE model generates y
def ols(y):
return np.linalg.solve(X.T @ X, X.T @ y) # (X'X)^{-1} X'y, solved not inverted
y0 = draw_dataset(rng)
print("one sample's beta-hat:", np.round(ols(y0), 4))
print("true beta: ", beta_true)Close, but off by a little — that “little” is sampling noise. Now we draw 20,000 alternate universes, fit OLS on each, and stack the 20,000 estimates into an array. Each row is one \hat{\boldsymbol\beta} from one parallel Maya in one parallel world.
n_sims = 20_000
betas_ols = np.array([ols(draw_dataset(rng)) for _ in range(n_sims)])
print("collected betas shape:", betas_ols.shape) # (20000, 2)3. Unbiasedness, seen
The proof said \mathbb{E}[\hat{\boldsymbol\beta}] = \boldsymbol\beta. The Monte Carlo version of “expected value” is “average over the 20,000 draws.” If the proof is right, the column means of our array should sit on top of (12, -0.42).
mean_hat = betas_ols.mean(axis=0)
print("true beta ........", beta_true)
print("mean of beta-hats.", np.round(mean_hat, 4))
print("difference .......", np.round(mean_hat - beta_true, 4))
assert np.allclose(mean_hat, beta_true, atol=0.02), "unbiasedness check failed"
print("\nUNBIASED: the average estimate lands on the truth (within Monte Carlo noise).")There it is — the average slope estimate is -0.42 to two decimals, and the intercept is 12. The recipe aims at the truth. Notice we never used homoskedasticity or normality to get this; the proof in §2 of the chapter leaned only on \mathbb{E}[\boldsymbol\varepsilon\mid\mathbf{X}]=\mathbf{0}, which our draw_dataset honors by construction (the errors are mean-zero and independent of X).
A picture makes it visceral. Below is the sampling distribution of the slope \hat\beta_1 — the histogram of all 20,000 estimates — with the truth marked. It is centered on -0.42, not shifted off it.
fig, ax = plt.subplots()
ax.hist(betas_ols[:, 1], bins=60, color="steelblue", alpha=0.8, edgecolor="white", lw=0.3)
ax.axvline(beta_true[1], color="crimson", ls="--", lw=2, label=r"true $\beta_1 = -0.42$")
ax.axvline(betas_ols[:, 1].mean(), color="black", ls=":", lw=2, label="mean of estimates")
ax.set_xlabel(r"$\hat{\beta}_1$ (slope estimate)")
ax.set_ylabel("count across 20,000 simulated samples")
ax.set_title("Sampling distribution of the OLS slope — centered on the truth")
ax.legend()
fig.tight_layout()
fig.savefig("/tmp/nb22_unbiased.png", dpi=90)
print("slope histogram centered at", round(betas_ols[:,1].mean(), 4))4. The variance formula, checked to the decimal
The chapter’s boxed result is \operatorname{Var}(\hat{\boldsymbol\beta}\mid\mathbf{X}) = \sigma^2(\mathbf{X}'\mathbf{X})^{-1}. The Monte Carlo version of “variance” is the empirical variance of our 20,000 estimates. If the formula is right, the empirical variance of the slope should match the bottom-right entry of \sigma^2(\mathbf{X}'\mathbf{X})^{-1}.
This is the same trick as the histogram for the sample mean in Week 1: the analytic formula and the empirical spread are two views of one object.
theory_cov = sigma**2 * XtX_inv # the boxed formula
emp_var_slope = betas_ols[:, 1].var() # Monte Carlo variance of the slope
th_var_slope = theory_cov[1, 1]
print(f"Empirical Var(slope) ...... {emp_var_slope:.6e}")
print(f"Theory Var(slope) ....... {th_var_slope:.6e}")
print(f"ratio (should be ~1) ...... {emp_var_slope / th_var_slope:.4f}")
# also check the intercept variance and the full covariance matrix
print("\nEmpirical cov matrix:\n", np.round(np.cov(betas_ols.T), 5))
print("Theory cov matrix:\n", np.round(theory_cov, 5))
assert abs(emp_var_slope / th_var_slope - 1) < 0.05, "variance formula check failed"
print("\nVARIANCE FORMULA CONFIRMED: empirical spread matches sigma^2 (X'X)^{-1}.")The ratio sits within a few percent of 1.0, and the whole 2\times2 covariance matrices agree. The formula is not a fudge — it is the literal spread of the estimates across repeated samples.
5. Efficiency: OLS versus a rival
Now the heart of the matter. Gauss–Markov claims OLS is Best — minimum variance — among all linear unbiased estimators. To see “best,” we need a competitor that is also linear and also unbiased, so that the only thing separating them is variance.
Our rival is a “split-sample” estimator: it throws away half the data and runs OLS on a fixed half (say, the first N/2 loans). Is it linear in \mathbf{y}? Yes — it is still \mathbf{C}\mathbf{y} for a matrix \mathbf{C} built from that half of \mathbf{X}. Is it unbiased? Yes — the §2 proof works on any fixed subsample, since \mathbb{E}[\boldsymbol\varepsilon \mid\mathbf{X}]=\mathbf{0} holds for those rows too. So it is a card-carrying member of the linear-unbiased club. It is just wasteful: it ignores half the spread in income, and spread is precision. Gauss–Markov predicts it must have larger variance than OLS. Let us watch.
half = N // 2
X_half = X[:half].copy() # fixed subsample; .copy() avoids chained-indexing surprises
XhtXh_inv = np.linalg.inv(X_half.T @ X_half)
def rival_split(y):
yh = y[:half]
return np.linalg.solve(X_half.T @ X_half, X_half.T @ yh)
betas_rival = np.array([rival_split(draw_dataset(rng)) for _ in range(n_sims)])
print("OLS mean beta-hat:", np.round(betas_ols.mean(axis=0), 4))
print("Rival mean beta-hat:", np.round(betas_rival.mean(axis=0), 4), " <- also unbiased")
print()
print(f"Var(slope), OLS .... {betas_ols[:,1].var():.6e}")
print(f"Var(slope), Rival .... {betas_rival[:,1].var():.6e}")
print(f"Rival/OLS variance ... {betas_rival[:,1].var()/betas_ols[:,1].var():.2f}x (>1: OLS wins)")Both estimators are unbiased — their slope distributions both center on -0.42 — but the rival’s variance is roughly twice the OLS variance (throwing away half the data costs you, just as the \sigma^2/\sum(x_i-\bar x)^2 formula predicts: halving the effective sample roughly doubles the variance). The two histograms below tell the whole story at a glance: same center, different width. That difference in width is what “BLUE” means. OLS is the tighter cluster.
fig, ax = plt.subplots()
ax.hist(betas_rival[:, 1], bins=70, color="darkorange", alpha=0.55,
edgecolor="white", lw=0.2, label="rival (half-sample)")
ax.hist(betas_ols[:, 1], bins=70, color="steelblue", alpha=0.65,
edgecolor="white", lw=0.2, label="OLS (full sample)")
ax.axvline(beta_true[1], color="crimson", ls="--", lw=2, label=r"true $\beta_1$")
ax.set_xlabel(r"$\hat{\beta}_1$ (slope estimate)")
ax.set_ylabel("count across 20,000 samples")
ax.set_title("Both unbiased (same center); OLS is narrower (smaller variance) — BLUE")
ax.legend()
fig.tight_layout()
fig.savefig("/tmp/nb22_efficiency.png", dpi=90)
print("OLS std:", round(betas_ols[:,1].std(), 4), " Rival std:", round(betas_rival[:,1].std(), 4))6. Break homoskedasticity: OLS keeps its aim, loses its crown
Gauss–Markov’s “Best” rests on Assumption 4, homoskedasticity: every error has the same variance \sigma^2. In Maya’s real data that is false — the scatter of rates is wider for low-income borrowers (more risk-based pricing discretion) and tighter for high earners. Let us inject exactly that: make the error standard deviation shrink with income, \operatorname{sd}(\varepsilon_i) \propto 1/\sqrt{w_i} where w_i grows with income. Two things should happen, per the chapter:
- OLS stays unbiased. Heteroskedasticity touches the variance of \boldsymbol\varepsilon, not its mean, and the unbiasedness proof never used Assumption 4. So \hat{\boldsymbol\beta} still centers on the truth.
- OLS is no longer BLUE. A weighted least squares (WLS) estimator that down-weights the noisy low-income observations now achieves smaller variance than OLS. OLS is dethroned.
We build the heteroskedastic world, then compare OLS to the WLS estimator that uses the (here known) weights w_i.
# Heteroskedastic errors: variance LARGER for low income. Pick weights w_i that rise with income.
w = 0.05 + 0.95 * (income - income.min()) / np.ptp(income) # in (0.05, 1.0], rises with income
sd_i = sigma / np.sqrt(w) # sd is large where w is small (low income)
def draw_hetero(rng):
eps = rng.normal(0.0, sd_i) # per-observation sd -> heteroskedastic
return X @ beta_true + eps
# WLS: minimize sum w_i (y_i - x_i' b)^2 == OLS on sqrt(w)-scaled rows.
sw = np.sqrt(w)
Xw = X * sw[:, None]
XwtXw_inv = np.linalg.inv(Xw.T @ Xw)
def wls(y):
yw = y * sw
return np.linalg.solve(Xw.T @ Xw, Xw.T @ yw)
betas_ols_h = np.array([ols(draw_hetero(rng)) for _ in range(n_sims)])
# redraw the SAME error realizations for a fair comparison by reusing a child generator:
rng_w = np.random.default_rng(123)
betas_wls_h = np.array([wls(draw_hetero(rng_w)) for _ in range(n_sims)])
print("True slope ...........", beta_true[1])
print("OLS mean slope .......", round(betas_ols_h[:,1].mean(), 4), "(still unbiased)")
print("WLS mean slope .......", round(betas_wls_h[:,1].mean(), 4), "(still unbiased)")
print()
print(f"Var(slope) OLS ...... {betas_ols_h[:,1].var():.6e}")
print(f"Var(slope) WLS ...... {betas_wls_h[:,1].var():.6e}")
print(f"WLS is {betas_ols_h[:,1].var()/betas_wls_h[:,1].var():.2f}x tighter -> OLS no longer BLUE")Both still aim true, but now WLS is the tighter cluster — under heteroskedasticity a weighted estimator beats OLS, so OLS has lost the “B.” That is the chapter’s claim, made flesh.
There is a second, sneakier casualty: the classical standard-error formula breaks. The formula \sigma^2(\mathbf{X}'\mathbf{X})^{-1} assumed \operatorname{Var}(\boldsymbol\varepsilon\mid\mathbf{X}) = \sigma^2\mathbf{I}, which is now false. So the SE that comes off a naive regression printout no longer matches the actual spread of \hat\beta_1. Below we compare the classical-formula SE to the true Monte Carlo SE under heteroskedasticity — they disagree, which means every t-stat and confidence interval built from the classical SE is untrustworthy. This is exactly the problem Chapter 2.4 fixes with robust (HC1/HC2/HC3) standard errors.
# What a naive analyst would report: classical SE using a single pooled sigma-hat estimate.
# Use the average estimated residual variance across sims to mimic "what the formula says".
classical_se_slope = np.sqrt((sigma**2 * XtX_inv)[1, 1]) # the formula, assuming homoskedasticity
true_mc_se_slope = betas_ols_h[:, 1].std() # the HONEST spread under heteroskedasticity
print(f"Classical-formula SE(slope) ... {classical_se_slope:.5f} (assumes homoskedasticity)")
print(f"True Monte Carlo SE(slope) .... {true_mc_se_slope:.5f} (what actually happens)")
print(f"ratio (classical / true) ...... {classical_se_slope/true_mc_se_slope:.3f}")
print()
print("The classical SE MISREPORTS the real uncertainty -> t-stats & CIs are wrong.")
print("Chapter 2.4 repairs this with heteroskedasticity-robust standard errors.")Your Turn — break the assumption that matters
You just saw that breaking homoskedasticity (Assumption 4) costs OLS its efficiency crown and wrecks its standard errors — annoying, but survivable: \hat\beta_1 still aimed at the truth, and Ch 2.4 patches the SEs. Now break the assumption whose failure is existential: zero conditional mean, \mathbb{E}[\boldsymbol\varepsilon\mid\mathbf{X}]=\mathbf{0} (Assumption 3).
Edit draw_biased below so the error is correlated with income — for instance, let part of the error grow with income, \varepsilon_i = (\text{mean-zero noise}) + \gamma\,(\text{income}_i - \bar{\text{income}}), with \gamma \neq 0. This mimics an omitted variable (creditworthiness) that rides along with income. Then run the cell and answer:
- Is \hat\beta_1 still centered on the true -0.42? By how much is it off, and which way?
- Crank N up to 40,000 and rerun. Does the bias shrink the way variance did, or does the estimate converge confidently to the wrong number? (This is the difference between unbiased-and-consistent versus consistent-for-the-wrong-target from §6 — more data makes you precisely wrong.)
- Which exact line of the §2 unbiasedness proof breaks when \mathbb{E}[\boldsymbol\varepsilon\mid \mathbf{X}]\neq\mathbf{0}? This is the engine of omitted-variable bias — the subject of Chapter 2.5.
# YOUR TURN: set gamma != 0 to make the error depend on income (violating E[eps|X]=0).
gamma = 0.0 # <-- change me, e.g. gamma = 0.30, then try larger
income_centered = income - income.mean()
def draw_biased(rng):
noise = rng.normal(0.0, sigma, size=N)
eps = noise + gamma * income_centered # gamma != 0 => E[eps | X] != 0
return X @ beta_true + eps
betas_biased = np.array([ols(draw_biased(rng)) for _ in range(5_000)])
print(f"gamma = {gamma}")
print("true slope .........", beta_true[1])
print("mean OLS slope .....", round(betas_biased[:,1].mean(), 4))
print("bias ...............", round(betas_biased[:,1].mean() - beta_true[1], 4))
print()
print("With gamma=0 the bias is ~0 (Assumption 3 holds). Set gamma != 0 and watch it appear,")
print("then raise N and watch the WRONG number get more precise, not more correct. -> Ch 2.5")