import matplotlib
matplotlib.use("Agg") # headless: render figures to buffers, never to a screen
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
from linearmodels.iv import IV2SLS
rng = np.random.default_rng(42) # one pinned generator 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__, "| numpy", np.__version__)
import linearmodels
print("linearmodels", linearmodels.__version__)nb3.5 — The Weak-Instrument Pathology, Reproduced
Week 3, Chapter 3.5. Companion notebook.
Devon read an IV paper whose 2SLS estimate was three times the OLS estimate, with a tight standard error and dense stars. Chapter 3.5 argued that a tight standard error on a 2SLS coefficient tells you nothing about whether the instrument was strong enough to trust it. This notebook makes that argument by simulation, where we know the true effect — we set it to exactly zero — so any nonzero estimate is pure bias.
The single most important habit this notebook drills: a single draw can mislead. One simulated dataset might, by luck, give a 2SLS estimate near the truth and a confidence interval that covers it. The pathology only becomes undeniable when we average over many draws. So almost everything here is a Monte Carlo: we redraw the world 1,000+ times and look at averages, coverage rates, and distributions, never at one lucky (or unlucky) sample.
What we will establish, all by averaging over draws:
- With a weak instrument, the average 2SLS estimate is biased toward OLS — the very bias IV was supposed to cure.
- The conventional 2SLS 95% confidence interval under-covers badly: it contains the true zero far less than 95% of the time.
- The first-stage F is usually small but occasionally spuriously large — significance is not strength.
- Anderson–Rubin confidence intervals, built by inverting a test over a grid, recover ≈95% coverage even when the instrument is weak — and they are sometimes unbounded, which is them telling the truth.
- The Olea–Pflueger effective F flags the weak case.
- With a strong instrument, every method behaves and agrees.
The “Your Turn” closes with many-instruments bias: padding one weak instrument with junk makes the toward-OLS bias worse.
1. The data-generating process: a world where we know the truth
We build the world from Chapter 3.5.8, deliberately rigged so that everything that can go wrong is something we chose and can measure. The structural equation is
y = \beta\, x + u, \qquad \beta = 0 \;\;(\text{the true effect is exactly zero}).
The regressor x is endogenous: it shares the shock u with the outcome, so OLS will be biased upward. The instrument z is valid (it never enters y except through x, so exclusion holds), but weak: its first-stage coefficient \pi is small.
x = \pi\, z + c\, u + v, \qquad y = \beta\, x + u.
Here c controls the strength of the endogeneity (how badly OLS is biased), \pi controls the first-stage strength (how strong the instrument is), and z, u, v are independent standard normals. Because \beta = 0, the outcome is literally y = u — so z touches y only through x, exactly as exclusion demands. We pick a small \pi = 0.07 and a strong endogeneity c = 2.0 to make the pathology vivid; later we flip \pi to 0.8 for the strong-instrument contrast.
def simulate(N, pi, beta_true=0.0, endog=2.0, rng=None):
'''Draw one dataset: endogenous x, valid-but-(maybe)-weak instrument z, true effect beta_true.'''
z = rng.normal(size=N) # valid instrument
u = rng.normal(size=N) # structural error (the endogeneity)
v = rng.normal(size=N) # idiosyncratic first-stage noise
x = pi * z + endog * u + v # weak first stage (pi) + endogeneity (endog * u)
y = beta_true * x + u # exclusion holds: z absent from y given x
return z, x, y
# Settings for the whole notebook
N = 400 # sample size
PI_WEAK = 0.07 # weak first stage
PI_STRONG = 0.8 # strong first stage
ENDOG = 2.0 # endogeneity strength -> OLS badly biased upward
BETA_TRUE = 0.0 # the truth; any nonzero estimate is pure bias
z, x, y = simulate(N, PI_WEAK, rng=rng)
print(f"Drew one dataset: N={N}, pi={PI_WEAK}, true beta={BETA_TRUE}")
print(f"corr(x, structural shock) is built in via endog={ENDOG} -> OLS biased up")2. Hand-rolled 2SLS (and a check against linearmodels)
Two-stage least squares with one instrument and a constant is just OLS run twice, or equivalently the projection formula
\hat{\boldsymbol\beta}_{\text{2SLS}} = (\mathbf{X}'\mathbf{P}_Z\mathbf{X})^{-1}\mathbf{X}'\mathbf{P}_Z\mathbf{y}, \qquad \mathbf{P}_Z = \mathbf{Z}(\mathbf{Z}'\mathbf{Z})^{-1}\mathbf{Z}',
where \mathbf{X} = [\mathbf{1}, x] and \mathbf{Z} = [\mathbf{1}, z]. The conventional (homoskedastic) standard error comes from \hat\sigma^2 (\mathbf{X}'\mathbf{P}_Z\mathbf{X})^{-1}. We write it by hand so nothing is hidden, then confirm it reproduces linearmodels.IV2SLS to machine precision. We also grab the first-stage F: with one instrument, F = t^2 from regressing x on z.
def fit_2sls(z, x, y):
'''Hand-rolled just-identified 2SLS: returns (beta_hat, conventional SE).'''
N = len(y)
Z = np.column_stack([np.ones(N), z])
X = np.column_stack([np.ones(N), x])
ZtZ_inv = np.linalg.inv(Z.T @ Z)
Pz_X = Z @ (ZtZ_inv @ (Z.T @ X)) # P_Z X
XtPzX = X.T @ Pz_X
beta = np.linalg.solve(XtPzX, Pz_X.T @ y) # (X'PzX)^-1 X'Pz y
resid = y - X @ beta # structural residual (uses x, NOT x-hat)
sigma2 = (resid @ resid) / (N - X.shape[1])
cov = sigma2 * np.linalg.inv(XtPzX)
return beta[1], np.sqrt(cov[1, 1])
def ols_slope(x, y):
N = len(y)
X = np.column_stack([np.ones(N), x])
return np.linalg.solve(X.T @ X, X.T @ y)[1]
def first_stage_F(z, x):
'''One-instrument first-stage F = t^2 from regressing x on [1, z].'''
N = len(x)
Z = np.column_stack([np.ones(N), z])
ZtZ_inv = np.linalg.inv(Z.T @ Z)
b = ZtZ_inv @ (Z.T @ x)
resid = x - Z @ b
sigma2 = (resid @ resid) / (N - 2)
se1 = np.sqrt(sigma2 * ZtZ_inv[1, 1])
return (b[1] / se1) ** 2
# Cross-check against linearmodels on the single dataset we drew
b_hand, se_hand = fit_2sls(z, x, y)
df = pd.DataFrame({"y": y, "x": x, "z": z, "const": 1.0})
lm = IV2SLS(df["y"], df[["const"]], df[["x"]], df[["z"]]).fit(cov_type="unadjusted")
print(f"hand-rolled 2SLS beta = {b_hand:+.6f}")
print(f"linearmodels 2SLS beta = {lm.params['x']:+.6f}")
print(f"max abs difference = {abs(b_hand - lm.params['x']):.2e}")
print(f"\nOLS beta = {ols_slope(x, y):+.3f} (true 0.000; biased UP by endogeneity)")
print(f"2SLS beta = {b_hand:+.3f} SE = {se_hand:.3f} first-stage F = {first_stage_F(z, x):.2f}")3. Why one draw is not enough
Look at the numbers above. On this particular draw the 2SLS estimate landed somewhere, the SE looks respectable, the first-stage F is whatever it is. A student could eyeball this single result and draw the wrong conclusion — maybe 2SLS happened to land near zero, maybe the F happened to be large. A single draw is an anecdote, not evidence. The whole point of the chapter is a statement about the sampling distribution of these estimators, and you cannot see a distribution from one point.
So from here on we do the honest thing: we redraw the world many times and study the average behavior. The function below runs one Monte Carlo experiment — many draws at a fixed \pi — and records, for each draw, the OLS estimate, the 2SLS estimate, whether the conventional 95% interval covered the truth, and the first-stage F.
4. Anderson–Rubin by hand, and the Monte Carlo engine
The cure from Chapter 3.5.4 is Anderson–Rubin (AR): never divide by the weak first stage. To test H_0:\beta=\beta_0, form the adjusted outcome y - \beta_0 x and regress it on the instrument z. If \beta_0 is the truth, z has no relationship with the adjusted outcome (exclusion), so the coefficient on z should be zero — test it with an F-test. To get a confidence interval, sweep \beta_0 over a grid and keep every value the test fails to reject at 5%.
We vectorize the grid sweep with a small algebra trick. Regressing (y - \beta_0 x) on \mathbf{Z}, the coefficient and the residual are both linear in \beta_0:
\text{coef}_z(\beta_0) = b^{(y)}_z - \beta_0\, b^{(x)}_z, \qquad \hat r(\beta_0) = r^{(y)} - \beta_0\, r^{(x)},
where b^{(y)}, b^{(x)} are the OLS coefficients of y and x on \mathbf{Z} and r^{(y)}, r^{(x)} their residuals. So we compute the four pieces once and evaluate the AR F-statistic cheaply at every grid point. An AR interval that touches a grid edge we flag as unbounded — the honest signature of a weak instrument.
def ar_interval(z, x, y, grid):
'''Anderson-Rubin 95% CI by grid inversion. Returns (lo, hi, accept_mask, empty, unbounded).'''
N = len(y)
Z = np.column_stack([np.ones(N), z])
ZtZ_inv = np.linalg.inv(Z.T @ Z)
crit = stats.f.ppf(0.95, 1, N - 2) # 1 instrument -> F(1, N-2)
by = ZtZ_inv @ (Z.T @ y) # coef of y on [1, z]
bx = ZtZ_inv @ (Z.T @ x) # coef of x on [1, z]
ry = y - Z @ by # residuals, linear pieces
rx = x - Z @ bx
var_coef1 = ZtZ_inv[1, 1]
accept = np.zeros(len(grid), dtype=bool)
for i, b0 in enumerate(grid):
coef1 = by[1] - b0 * bx[1] # coef on z for this candidate
resid = ry - b0 * rx
sigma2 = (resid @ resid) / (N - 2)
F_ar = coef1 ** 2 / (sigma2 * var_coef1)
accept[i] = F_ar < crit # fail to reject -> keep b0
if not accept.any():
return np.nan, np.nan, accept, True, False # empty
idx = np.where(accept)[0]
lo, hi = grid[idx[0]], grid[idx[-1]]
unbounded = bool(accept[0] or accept[-1]) # touches a grid edge
return lo, hi, accept, False, unbounded
GRID = np.linspace(-10, 10, 1001) # candidate beta values for AR inversion
def run_mc(N, pi, reps, grid, rng, beta_true=0.0, endog=2.0):
'''One Monte Carlo experiment at fixed pi. Returns a dict of arrays / coverage rates.'''
ols_e = np.empty(reps); iv_e = np.empty(reps); Fs = np.empty(reps)
conv_cov = 0; ar_cov = 0; ar_unb = 0; ar_empty = 0
for r in range(reps):
z, x, y = simulate(N, pi, beta_true=beta_true, endog=endog, rng=rng)
b_iv, se_iv = fit_2sls(z, x, y)
iv_e[r] = b_iv
ols_e[r] = ols_slope(x, y)
Fs[r] = first_stage_F(z, x)
# conventional 95% CI coverage of the truth
if b_iv - 1.96 * se_iv <= beta_true <= b_iv + 1.96 * se_iv:
conv_cov += 1
# AR 95% CI coverage of the truth
lo, hi, _, empty, unb = ar_interval(z, x, y, grid)
if empty:
ar_empty += 1
else:
if unb: ar_unb += 1
if lo <= beta_true <= hi: ar_cov += 1
return {
"ols": ols_e, "iv": iv_e, "F": Fs,
"conv_coverage": conv_cov / reps,
"ar_coverage": ar_cov / reps,
"ar_unbounded": ar_unb / reps,
"ar_empty": ar_empty / reps,
}
print("AR + Monte Carlo engine ready.")5. The weak instrument: averaging over 1,200 draws
Now the main event. We run 1,200 draws of the weak-instrument world and read the averages. Watch four things: (a) the average OLS estimate sits far above the true zero — endogeneity at work; (b) the average 2SLS estimate sits between OLS and zero — it did not purge the bias, it shrank it partway back toward OLS, exactly the 1/(F+1) pull; (c) the conventional coverage is far below the promised 95%; (d) the AR coverage is back near 95%, and AR is often unbounded.
Because 2SLS is a ratio with a near-zero denominator, its sampling distribution is fat-tailed and the mean can be thrown around by a few wild draws. The honest summary of “where the estimate typically lands” is the median, so we report both but lean on the median for the central-tendency claim.
REPS = 1200
weak = run_mc(N, PI_WEAK, REPS, GRID, rng, beta_true=BETA_TRUE, endog=ENDOG)
print(f"WEAK instrument (pi={PI_WEAK}), {REPS} draws, true beta = {BETA_TRUE}\n")
print(f" median OLS estimate = {np.median(weak['ols']):+.3f} (biased up by endogeneity)")
print(f" median 2SLS estimate = {np.median(weak['iv']):+.3f} (sits BETWEEN OLS and 0 -> toward-OLS bias)")
print(f" mean 2SLS estimate = {np.mean(weak['iv']):+.3f} (fat-tailed; mean is unstable)")
print()
print(f" conventional 95% coverage = {weak['conv_coverage']:.3f} <-- should be 0.95, badly UNDER")
print(f" Anderson-Rubin coverage = {weak['ar_coverage']:.3f} <-- back near 0.95, HONEST")
print(f" fraction of AR CIs unbounded = {weak['ar_unbounded']:.3f}")
print()
print(f" mean first-stage F = {np.mean(weak['F']):.2f}")
print(f" fraction of draws F < 10 = {np.mean(weak['F'] < 10):.3f} (klaxon almost always sounding)")
print(f" but MAX first-stage F = {np.max(weak['F']):.1f} <-- sometimes spuriously LARGE")
# How far is 2SLS pulled back toward OLS? (fraction of OLS bias retained)
pull = np.median(weak['iv']) / np.median(weak['ols'])
print(f"\n 2SLS retains ~{pull:.0%} of OLS's bias (a strong instrument would retain ~0%).")The averages tell the whole story and no single draw could have. The median 2SLS estimate is positive and large — it landed between the biased OLS number and the truth of zero, carrying a big chunk of the very bias IV was meant to remove. The conventional 95% interval covered the truth far less than 95% of the time: a researcher reporting “estimate ± 1.96 SE” from a single weak-IV dataset is wrong about their own uncertainty more than one time in five. Anderson–Rubin, by contrast, covers near 95% — at the price of frequently being unbounded, which is simply AR refusing to pretend it knows an answer it does not.
We assert the three load-bearing claims so the notebook fails loudly if the pathology ever stops reproducing.
# Load-bearing claims of the chapter, asserted on the Monte Carlo output.
assert weak["conv_coverage"] < 0.90, "conventional coverage should be well below 0.95"
assert weak["ar_coverage"] > 0.93, "AR coverage should be near 0.95"
assert 0.0 < np.median(weak["iv"]) < np.median(weak["ols"]), \
"median 2SLS must sit strictly between the truth (0) and OLS (biased toward OLS)"
print("All weak-IV pathology assertions passed:")
print(f" conventional coverage {weak['conv_coverage']:.3f} < 0.95 (under-covers)")
print(f" AR coverage {weak['ar_coverage']:.3f} >= 0.93 (honest)")
print(f" 0 < median 2SLS {np.median(weak['iv']):.3f} < median OLS {np.median(weak['ols']):.3f} (toward OLS)")6. Significance is not strength: the first-stage F distribution
Bound–Jaeger–Baker’s lesson was that a statistically significant first stage can still be far too weak to support inference. We see it directly: across the 1,200 weak-IV draws the first-stage F is usually tiny, but its distribution has a long right tail — every so often a draw produces an F well above 10 purely by chance. A researcher who happened to draw one of those lucky samples would declare the instrument “strong” and march confidently into the pathology. The F > 10 rule is a guard against the typical draw, not a guarantee about your draw.
fig, ax = plt.subplots()
ax.hist(weak["F"], bins=60, color="#c0392b", alpha=0.8, edgecolor="white")
ax.axvline(10, color="black", linestyle="--", linewidth=2, label="F = 10 rule of thumb")
ax.axvline(np.mean(weak["F"]), color="#2c3e50", linestyle=":", linewidth=2,
label=f"mean F = {np.mean(weak['F']):.1f}")
ax.set_xlabel("first-stage F-statistic")
ax.set_ylabel("count across draws")
ax.set_title(f"Weak instrument: first-stage F is usually small, occasionally spuriously large\n"
f"(fraction of draws with F > 10: {np.mean(weak['F'] > 10):.1%})")
ax.legend()
fig.tight_layout()
print(f"Of {REPS} draws, {(weak['F'] > 10).sum()} produced a spurious F > 10 despite a truly weak instrument.")7. The contrast: a strong instrument, where everything works
The machinery is not broken in general — it is broken by weakness. Set \pi = 0.8 (a strong first stage) and re-run the identical Monte Carlo. Now 2SLS lands on the truth, the conventional interval recovers its 95% coverage, AR agrees and is never unbounded, and the first-stage F is enormous on every draw.
strong = run_mc(N, PI_STRONG, REPS, GRID, rng, beta_true=BETA_TRUE, endog=ENDOG)
print(f"STRONG instrument (pi={PI_STRONG}), {REPS} draws, true beta = {BETA_TRUE}\n")
print(f" median OLS estimate = {np.median(strong['ols']):+.3f} (still biased; OLS never fixed)")
print(f" median 2SLS estimate = {np.median(strong['iv']):+.3f} (on the truth -> bias purged)")
print(f" conventional 95% coverage = {strong['conv_coverage']:.3f} (recovered)")
print(f" Anderson-Rubin coverage = {strong['ar_coverage']:.3f} (agrees)")
print(f" fraction AR unbounded = {strong['ar_unbounded']:.3f} (never)")
print(f" mean first-stage F = {np.mean(strong['F']):.0f} (off the charts)")
assert strong["conv_coverage"] > 0.93, "strong-IV conventional coverage should recover to ~0.95"
assert abs(np.median(strong["iv"])) < 0.10, "strong-IV 2SLS should land on the truth"
print("\nStrong-IV assertions passed: conventional coverage recovers and 2SLS hits the truth.")8. The picture: coverage, conventional vs. Anderson–Rubin
One bar chart makes the central result unmissable. For the weak instrument, the conventional interval sits far below the 95% line while AR clears it; for the strong instrument, both clear it. The conventional method is not “approximately right” with a weak instrument — it is systematically overconfident.
fig, ax = plt.subplots()
labels = ["Weak\nconventional", "Weak\nAnderson-Rubin", "Strong\nconventional", "Strong\nAnderson-Rubin"]
vals = [weak["conv_coverage"], weak["ar_coverage"], strong["conv_coverage"], strong["ar_coverage"]]
colors = ["#c0392b", "#27ae60", "#e67e22", "#16a085"]
bars = ax.bar(labels, vals, color=colors, alpha=0.85, edgecolor="black")
ax.axhline(0.95, color="black", linestyle="--", linewidth=2, label="nominal 95%")
ax.set_ylabel("empirical coverage of the true value")
ax.set_ylim(0, 1.05)
ax.set_title("Coverage of the truth: conventional 2SLS CIs under-cover when the instrument is weak")
for b, v in zip(bars, vals):
ax.text(b.get_x() + b.get_width() / 2, v + 0.02, f"{v:.2f}", ha="center", fontweight="bold")
ax.legend(loc="lower right")
fig.tight_layout()
print("Conventional coverage collapses only under weakness; AR holds near 95% in both worlds.")9. The picture: distribution of 2SLS estimates, weak vs. strong
Overlay the sampling distributions of the 2SLS estimate from the two experiments against the true value (red line) and the OLS bias (orange line). The strong-instrument distribution is a tight bell centered on the truth. The weak-instrument distribution is shifted toward OLS and far wider — and it has fat tails we clip for the plot. This is the toward-OLS bias and the fat-tailed ratio, both visible at once.
fig, ax = plt.subplots()
clip = (-2, 2) # clip fat tails so the shapes are visible
w_clip = np.clip(weak["iv"], *clip)
s_clip = np.clip(strong["iv"], *clip)
ax.hist(w_clip, bins=60, density=True, alpha=0.55, color="#c0392b", label="weak 2SLS")
ax.hist(s_clip, bins=60, density=True, alpha=0.55, color="#16a085", label="strong 2SLS")
ax.axvline(BETA_TRUE, color="red", linewidth=2.5, label=f"truth = {BETA_TRUE}")
ax.axvline(np.median(weak["ols"]), color="#e67e22", linewidth=2.5, linestyle="--",
label=f"OLS bias ≈ {np.median(weak['ols']):.2f}")
ax.axvline(np.median(weak["iv"]), color="#7d1010", linewidth=2, linestyle=":",
label=f"weak 2SLS median ≈ {np.median(weak['iv']):.2f}")
ax.set_xlabel("2SLS estimate (tails clipped to [-2, 2] for display)")
ax.set_ylabel("density across draws")
ax.set_title("Sampling distribution of 2SLS: weak is shifted toward OLS and far wider")
ax.legend(fontsize=9)
fig.tight_layout()
print("The weak-IV distribution is dragged off the truth toward OLS; the strong-IV one sits on the truth.")10. The diagnostic that flags it: Olea–Pflueger effective F
Finally, the diagnostic. linearmodels reports the Olea–Pflueger effective F in its first-stage diagnostics (robust to heteroskedasticity, and the right object in many-instrument settings). On a weak draw it sits far below the conventional threshold of ~10; on a strong draw it is huge. The effective F is the number to report and to demand from any IV paper — it does not reward you for adding junk instruments, as we will see in the “Your Turn”.
def effective_F(z, x, y):
'''Pull the Olea-Pflueger effective F from linearmodels' first-stage diagnostics.'''
d = pd.DataFrame({"y": y, "x": x, "z": z, "const": 1.0})
res = IV2SLS(d["y"], d[["const"]], d[["x"]], d[["z"]]).fit(cov_type="robust")
return float(res.first_stage.diagnostics.loc["x", "f.stat"])
# Average the effective F over a modest set of fresh draws in each world
rng_F = np.random.default_rng(7)
eff_weak = np.array([effective_F(*simulate(N, PI_WEAK, endog=ENDOG, rng=rng_F)) for _ in range(100)])
eff_strong = np.array([effective_F(*simulate(N, PI_STRONG, endog=ENDOG, rng=rng_F)) for _ in range(100)])
print("Olea-Pflueger effective F (averaged over 100 fresh draws each):")
print(f" weak instrument: mean effective F = {eff_weak.mean():6.2f} (below ~10 -> FLAGGED weak)")
print(f" strong instrument: mean effective F = {eff_strong.mean():6.1f} (far above -> strong)")
assert eff_weak.mean() < 10 < eff_strong.mean(), "effective F must flag weak and clear strong"
print("\nThe effective F correctly flags the weak case and clears the strong case.")Your Turn — Many-instruments bias: when more is worse
Sam has one weak instrument and reasons: more instruments must mean a stronger first stage. So he pads his single weak instrument with extra ones. If those extras are junk (irrelevant to x), watch what happens to the bias.
The mechanism, from Chapter 3.5.7: 2SLS replaces x with its first-stage fitted value \hat x, the projection of x onto all the instruments. Every extra junk instrument lets that projection chase a little more of the idiosyncratic noise in x — and that noise is exactly the endogenous part. So \hat x creeps back toward x itself, and 2SLS creeps back toward OLS. In the extreme (instruments approaching the sample size), 2SLS collapses onto OLS — the worst case of the toward-OLS bias.
Run the experiment below: keep one weak relevant instrument, pad it with K-1 pure-noise instruments, and track the median 2SLS estimate across draws as K grows. The bias should march monotonically from its just-identified value up toward the OLS bias.
Extensions to try yourself: (1) Replace the junk instruments with weakly relevant ones (give each a tiny \pi) and confirm the bias still worsens. (2) Add the effective F at each K and watch it stay low even as the naive first-stage F on many instruments climbs — the trap Sam falls into. (3) Push \pi in the main experiment down to 0.02 and confirm the AR interval is unbounded on essentially every draw.
def fit_2sls_many(Zmat, x, y):
'''Hand-rolled 2SLS slope with a matrix of instruments Zmat (N x K).'''
N = len(y)
Z = np.column_stack([np.ones(N), Zmat])
X = np.column_stack([np.ones(N), x])
ZtZ_inv = np.linalg.inv(Z.T @ Z)
Pz_X = Z @ (ZtZ_inv @ (Z.T @ X))
XtPzX = X.T @ Pz_X
return np.linalg.solve(XtPzX, Pz_X.T @ y)[1]
K_values = [1, 3, 10, 30, 60]
median_iv = []
median_ols = []
for K in K_values:
rng_K = np.random.default_rng(123) # same seed each K -> compare like with like
ests = np.empty(800); olss = np.empty(800)
for r in range(800):
z1 = rng_K.normal(size=N)
u = rng_K.normal(size=N); v = rng_K.normal(size=N)
x = PI_WEAK * z1 + ENDOG * u + v # ONE weak relevant instrument
y = BETA_TRUE * x + u # truth still zero
if K == 1:
Zmat = z1.reshape(-1, 1)
else:
junk = rng_K.normal(size=(N, K - 1)) # K-1 PURE-NOISE instruments
Zmat = np.column_stack([z1, junk])
ests[r] = fit_2sls_many(Zmat, x, y)
olss[r] = ols_slope(x, y)
median_iv.append(np.median(ests))
median_ols.append(np.median(olss))
print(f" K = {K:2d} instruments (1 relevant + {K-1:2d} junk): median 2SLS = {np.median(ests):.3f}")
ols_ref = np.median(median_ols)
print(f"\n OLS bias reference (median) = {ols_ref:.3f}")
assert median_iv[-1] > median_iv[0] + 0.05, "many junk instruments should push 2SLS toward OLS"
print(f" many-instruments bias: median 2SLS rose from {median_iv[0]:.3f} (K=1) "
f"to {median_iv[-1]:.3f} (K={K_values[-1]}), marching toward OLS = {ols_ref:.3f}.")fig, ax = plt.subplots()
ax.plot(K_values, median_iv, "o-", color="#8e44ad", linewidth=2, markersize=9, label="median 2SLS")
ax.axhline(ols_ref, color="#e67e22", linestyle="--", linewidth=2, label=f"OLS bias = {ols_ref:.2f}")
ax.axhline(BETA_TRUE, color="red", linewidth=2, label=f"truth = {BETA_TRUE}")
ax.set_xlabel("number of instruments K (1 relevant + K-1 junk)")
ax.set_ylabel("median 2SLS estimate")
ax.set_title("Many-instruments bias: padding a weak instrument with junk drags 2SLS toward OLS")
ax.legend()
fig.tight_layout()
print("Takeaway: a few strong instruments beat many weak ones, always.")
print("Adding junk does not strengthen identification — it reintroduces the endogeneity that IV was")
print("supposed to purge. The right fix for a weak instrument is a BETTER instrument, not MORE of them.")
print("\nThe discipline, from the chapter's ledger: before believing any 2SLS number — yours or a")
print("published paper's — check the first-stage F, the Olea-Pflueger effective F, and the")
print("Anderson-Rubin interval. A tight conventional SE proves nothing about a weak instrument.")