2SLS in two regressions
Two-stage least squares generalises the Wald estimator to continuous treatments, several instruments and covariates:
- First stage: regress D on Z (and covariates); keep the fitted values D̂, the part of D driven by Z.
- Second stage: regress Y on D̂ (and covariates).
Never run the second stage by hand for inference: its standard errors use the wrong residuals. Use a
dedicated IV routine such as linearmodels.IV2SLS.
!pip install -q numpy pandas statsmodels linearmodels matplotlib causaldata
import numpy as np, pandas as pd
import statsmodels.api as sm
from linearmodels.iv import IV2SLS
import matplotlib.pyplot as plt
rng = np.random.default_rng(101)
The ability-bias simulation
Unobserved ability raises both schooling and earnings. Living close to a college lowers the cost of schooling, and (we assume for now) has no other effect on earnings. True return: $3,000 per year.
N = 5_000
ability = rng.normal(0, 1, N)
near_college = rng.binomial(1, 0.5, N)
education = 12 + 1.5 * near_college + 2.0 * ability + rng.normal(0, 1, N)
earnings = 25_000 + 3_000 * education + 5_000 * ability + rng.normal(0, 4_000, N)
df = pd.DataFrame({"earnings": earnings, "education": education, "near_college": near_college, "const": 1.0})
ols = sm.OLS(df.earnings, df[["const", "education"]]).fit()
# 2SLS by hand (point estimate only!)
d_hat = sm.OLS(df.education, df[["const", "near_college"]]).fit().fittedvalues
manual = sm.OLS(df.earnings, sm.add_constant(d_hat.rename("education_hat"))).fit()
iv = IV2SLS(dependent=df.earnings, exog=df.const, endog=df.education, instruments=df.near_college).fit(cov_type="robust")
print(f"OLS: {ols.params['education']:,.0f}")
print(f"Manual 2SLS: {manual.params['education_hat']:,.0f} (SE {manual.bse['education_hat']:,.0f}: WRONG)")
print(f"IV2SLS: {iv.params['education']:,.0f} (SE {iv.std_errors['education']:,.0f}: correct)")
Real data: distance to college
Card (1995) used whether a man grew up near a four-year college (nearc4) as an instrument for schooling,
using the US National Longitudinal Survey of Young Men.
from causaldata import close_college
card = close_college.load_pandas().data.dropna()
controls = "1 + exper + I(exper**2) + black + south + smsa + married"
ols_card = IV2SLS.from_formula(f"lwage ~ {controls} + educ", card).fit(cov_type="robust")
iv_card = IV2SLS.from_formula(f"lwage ~ {controls} + [educ ~ nearc4]", card).fit(cov_type="robust")
print(f"OLS return: {100 * ols_card.params['educ']:.1f}% per year (SE {100 * ols_card.std_errors['educ']:.1f})")
print(f"2SLS return: {100 * iv_card.params['educ']:.1f}% per year (SE {100 * iv_card.std_errors['educ']:.1f})")
print(iv_card.first_stage.diagnostics.round(4))
The IV estimate is larger than OLS, the opposite of what simple ability bias predicts. A common reading: compliers are people for whom college costs, not ability, were the barrier, and their returns are high (a LATE, Lesson 9.1). Note the much larger standard error: IV uses only the variation driven by the instrument.
Weak instruments
When the first stage is weak, 2SLS is biased toward OLS and its confidence intervals are misleading. The classic rule of thumb (Staiger and Stock) asks for a first-stage F above 10. Recent work (Lee, McCrary, Moreira and Porter, 2022) shows that valid 5% t-tests need F above about 104, or a robust procedure such as the Anderson–Rubin test.
def simulate_iv(strength, n=1_000, reps=300):
est, fstats = [], []
for _ in range(reps):
a = rng.normal(size=n); z = rng.binomial(1, 0.5, n)
educ = 12 + strength * z + 2 * a + rng.normal(size=n)
y = 25_000 + 3_000 * educ + 5_000 * a + rng.normal(0, 4_000, n)
d = pd.DataFrame({"y": y, "educ": educ, "z": z, "const": 1.0})
r = IV2SLS(d.y, d.const, d.educ, d.z).fit(cov_type="robust")
est.append(r.params["educ"]); fstats.append(r.first_stage.diagnostics["f.stat"].iloc[0])
return np.median(est), np.median(fstats)
rows = [(s, *simulate_iv(s)) for s in [1.5, 0.5, 0.25, 0.1, 0.05]]
print(pd.DataFrame(rows, columns=["first-stage coefficient", "median 2SLS", "median F"]).round(2))
As the first stage weakens, the median estimate drifts from 3,000 toward the OLS value (about 5,000 here).
Always report the first stage: its coefficient, its F-statistic, and a plot of D against Z. If F is small, say so, and report weak-instrument-robust inference (Anderson–Rubin) rather than the usual t-test.
For the Card data, compute the reduced form (lwage on nearc4 and controls) and the first stage (educ on nearc4 and controls) with statsmodels. Check that their ratio equals the 2SLS coefficient.
Show a solution
import statsmodels.formula.api as smf
ctrl = "exper + I(exper**2) + black + south + smsa + married"
rf = smf.ols(f"lwage ~ nearc4 + {ctrl}", data=card).fit().params["nearc4"]
fs = smf.ols(f"educ ~ nearc4 + {ctrl}", data=card).fit().params["nearc4"]
print(f"Reduced form {rf:.4f} / first stage {fs:.4f} = {rf / fs:.4f} | 2SLS {iv_card.params['educ']:.4f}")
Checklist and recap
- I use an IV routine for inference, never a manual second stage.
- I can estimate and interpret Card's college-proximity IV.
- I report first-stage strength and know why weak instruments push 2SLS toward OLS.