Endogeneity, Instrumental Variables, and Two-Stage Least Squares

1. Endogeneity

A regressor \(X_j\) in the structural equation \(Y=X'\beta+u\) is endogenous if \(E[X_ju]\neq0\). Here \(\beta\) is the causal or structural parameter, not the projection coefficient. OLS then converges to \(\beta+Q_{XX}^{-1}E[Xu]\neq\beta\). The classic sources are:

• omitted variables (Lecture 3): ability in a wage equation;
• simultaneity (Lecture 6): price and quantity determined jointly by supply and demand;
• measurement error (Lecture 3): observed schooling equals true schooling plus noise;
• selection on gains: those who expect larger returns choose more treatment.

2. The instrumental-variables estimator

Partition \(X=(X_1',X_2')'\) into exogenous regressors \(X_1\) (including the constant) and endogenous \(X_2\). An instrument vector \(Z=(X_1',Z_2')'\), with \(\ell\ge k\) elements, must satisfy:

Exogeneity has two parts: the instrument is as good as randomly assigned with respect to \(u\), and it affects \(Y\) only through \(X\) (the exclusion restriction). With \(\ell=k\) (just identification), premultiplying \(Y=X'\beta+u\) by \(Z\) and taking expectations gives \(E[ZY]=E[ZX']\beta\), so

One regressor, one instrument. \(\beta_1=\operatorname{Cov}(Z,Y)/\operatorname{Cov}(Z,X)\): the ratio of the reduced form to the first stage. With a binary instrument this is the Wald estimator

3. Two-stage least squares

With \(\ell>k\) (over-identification), \(E[Z(Y-X'b)]=0\) is \(\ell\) equations in \(k\) unknowns. 2SLS uses the combination that projects \(X\) onto the instruments:

Two-stage interpretation. Stage 1 regresses each column of \(X\) on \(Z\) to get \(\widehat X=P_ZX\). Stage 2 regresses \(y\) on \(\widehat X\). The point estimates coincide, because \(\widehat X'\widehat X=X'P_ZX\) and \(\widehat X'y=X'P_Zy\). The second-stage OLS standard errors are wrong, because they use residuals \(y-\widehat X\widehat\beta\) instead of \(y-X\widehat\beta\). Always use a dedicated IV routine.

Never do "forbidden regressions". Plugging in fitted values from a nonlinear first stage (probit, squared fitted values) does not give a consistent estimator in general. Linear projection of \(X\) onto \(Z\) is what makes the first-stage residuals orthogonal to everything in the second stage.

2SLS is the efficient GMM estimator under homoskedasticity. With heteroskedasticity, two-step GMM with weight matrix \(\widehat\Omega^{-1}\) is more efficient (Lecture 8).

IV is less precise than OLS. In the just-identified single-regressor case under homoskedasticity, \(\operatorname{avar}(\widehat\beta_{IV})=\operatorname{avar}(\widehat\beta_{OLS})/\rho^2_{XZ}\), where \(\rho_{XZ}\) is the correlation between \(X\) and \(Z\). With \(\rho_{XZ}=0.2\), the variance is 25 times larger. IV trades bias for variance.

4. Weak instruments

When the first stage is weak, the asymptotic normal approximation of Theorem 5.1 is poor even in large samples.

• Bias toward OLS. With many weak instruments, 2SLS is biased toward the OLS probability limit. Approximately, \(E[\widehat\beta_{2SLS}]-\beta\approx(\operatorname{plim}\widehat\beta_{OLS}-\beta)/(F+1)\), where \(F\) is the population first-stage \(F\) statistic (Staiger and Stock, 1997).
• Size distortion. Nominal 5% Wald tests can reject far more often than 5%, and confidence intervals undercover.
• Bound et al. (1995) showed that quarter-of-birth estimates of the return to schooling with many weak instruments could be reproduced with randomly generated instruments.

Diagnosis. With one endogenous regressor, report the first-stage \(F\) for the excluded instruments, computed with a variance estimator matching the error structure. Under heteroskedasticity or clustering, use the effective \(F\) of Montiel Olea and Pflueger (2013). The rule of thumb \(F>10\) (Stock and Yogo, 2005) limits the 2SLS bias to about 10% of the OLS bias. Lee, McCrary, Moreira and Porter (2022) show that a 5% \(t\)-test with \(|t|>1.96\) requires a first-stage \(F\) of about 104.7 for correct size, or else adjusted critical values (the tF procedure).

Robust inference. The Anderson–Rubin (AR) test of \(H_0:\beta=\beta_0\) regresses \(Y-X_2'\beta_0\) on \(Z\) and tests the coefficients on the excluded instruments \(Z_2\). Under \(H_0\) this is a valid test however weak the instruments are, because \(Y-X_2'\beta_0=X_1'\gamma+u\) has no endogenous regressor. Inverting the AR test over a grid of \(\beta_0\) yields a confidence set with correct coverage. The set may be unbounded when instruments are weak: that is an honest statement that the data are uninformative. In the just-identified case, AR is efficient among robust tests. With over-identification, the conditional likelihood ratio (CLR) test of Moreira (2003) has better power.

5. Testing

Over-identification (Sargan–Hansen \(J\)). With \(\ell>k\), test whether all moment conditions can hold simultaneously:

Under homoskedasticity, \(J=nR^2\) from regressing the 2SLS residuals on \(Z\) (Sargan). Limitations: \(J\) tests coherence among instruments, not validity. If all instruments are invalid in the same direction, \(J\) has no power. And with heterogeneous effects, different valid instruments identify different LATEs, so \(J\) can reject even when all instruments are valid.

Endogeneity (Durbin–Wu–Hausman). Regress the endogenous regressor on \(Z\), save the first-stage residual \(\widehat v\), and add it to the structural equation (the control function regression). The OLS coefficients on \(X\) in this augmented regression equal 2SLS. A robust \(t\) test on \(\widehat v\) tests \(H_0:E[X_2u]=0\). Failing to reject is not evidence of exogeneity when the instruments are weak, because the test then has little power.

6. Heterogeneous effects and LATE

Let \(D\in\{0,1\}\) be a treatment with potential outcomes \(Y(1),Y(0)\), and \(Z\in\{0,1\}\) an instrument with potential treatments \(D(1),D(0)\). The population divides into compliers (\(D(1)=1,D(0)=0\)), always-takers, never-takers and defiers.

The Wald/IV estimand is therefore the average effect for compliers, the people whose treatment the instrument changes. It is not, in general, the ATE, nor the effect on the treated. Different instruments move different people and identify different LATEs. External validity has to be argued, not assumed. The characteristics of compliers can be described from the data (Abadie, 2003), which helps that argument.

import numpy as np
rng = np.random.default_rng(5)

def draw(n=500, pi=0.05, rho=0.8, beta=1.0):
    z = rng.normal(size=n)
    e = rng.multivariate_normal([0, 0], [[1, rho], [rho, 1]], size=n)
    x = pi * z + e[:, 1]                       # first stage
    y = beta * x + e[:, 0]                     # structural equation, Corr(u, v) = rho
    b_ols = (x @ y) / (x @ x)
    b_iv = (z @ y) / (z @ x)
    u = y - b_iv * x
    se_iv = np.sqrt((z**2 * u**2).sum()) / abs(z @ x)  # robust just-identified IV s.e.
    wald_rej = abs((b_iv - beta) / se_iv) > 1.96
    r = y - beta * x                           # AR test of H0: beta = 1
    g = (z @ r) / (z @ z); er = r - g * z
    ar_rej = abs(g / (np.sqrt((z**2 * er**2).sum()) / (z @ z))) > 1.96
    F = ((z @ x)**2 / (z @ z)) / (((x - (z @ x) / (z @ z) * z)**2).sum() / (n - 1))
    return b_ols, b_iv, wald_rej, ar_rej, F

for pi in (0.05, 0.15, 0.5):
    R = np.array([draw(pi=pi) for _ in range(3000)])
    print(f"pi={pi}: median F={np.median(R[:,4]):6.1f}  median OLS={np.median(R[:,0]):.2f}  "
          f"median IV={np.median(R[:,1]):.2f}  Wald size={R[:,2].mean():.3f}  AR size={R[:,3].mean():.3f}")

import numpy as np, wooldridge as wd, statsmodels.formula.api as smf
from linearmodels.iv import IV2SLS
df = wd.data("card")
ctrl = "exper + expersq + black + smsa + south + smsa66 + " + " + ".join(f"reg66{i}" for i in range(2, 10))
iv = IV2SLS.from_formula(f"lwage ~ 1 + {ctrl} + [educ ~ nearc4]", data=df).fit(cov_type="robust")
# TODO: OLS, first stage F, AR test: regress (lwage - b0*educ) on nearc4 + controls for a grid of b0

import numpy as np, wooldridge as wd, statsmodels.formula.api as smf
from linearmodels.iv import IV2SLS
df = wd.data("card")
ctrl = "exper + expersq + black + smsa + south + smsa66 + " + " + ".join(f"reg66{i}" for i in range(2, 10))
ols = smf.ols(f"lwage ~ educ + {ctrl}", data=df).fit(cov_type="HC1")
fs = smf.ols(f"educ ~ nearc4 + {ctrl}", data=df).fit(cov_type="HC1")
iv = IV2SLS.from_formula(f"lwage ~ 1 + {ctrl} + [educ ~ nearc4]", data=df).fit(cov_type="robust")
print(f"OLS  return to educ: {ols.params['educ']:.4f} (se {ols.bse['educ']:.4f})")
print(f"first stage: nearc4 -> educ = {fs.params['nearc4']:.4f}, robust F = {fs.tvalues['nearc4']**2:.1f}")
print(f"2SLS return to educ: {iv.params['educ']:.4f} (se {iv.std_errors['educ']:.4f})")

acc = []
for b0 in np.arange(-0.10, 0.60, 0.0025):
    df["ytil"] = df["lwage"] - b0 * df["educ"]
    t = smf.ols(f"ytil ~ nearc4 + {ctrl}", data=df).fit(cov_type="HC1").tvalues["nearc4"]
    if abs(t) < 1.96:
        acc.append(b0)
print(f"Anderson-Rubin 95% set: [{min(acc):.3f}, {max(acc):.3f}]  vs Wald: [{iv.params['educ']-1.96*iv.std_errors['educ']:.3f}, {iv.params['educ']+1.96*iv.std_errors['educ']:.3f}]")

OLS  return to educ: 0.0747 (se 0.0036)
first stage: nearc4 -> educ = 0.3199, robust F = 14.1
2SLS return to educ: 0.1315 (se 0.0540)
Anderson-Rubin 95% set: [0.030, 0.280]  vs Wald: [0.026, 0.237]

import wooldridge as wd, statsmodels.formula.api as smf
from linearmodels.iv import IV2SLS
df = wd.data("mroz").query("inlf == 1").copy()
iv = IV2SLS.from_formula("lwage ~ 1 + exper + expersq + [educ ~ fatheduc + motheduc]", data=df).fit(cov_type="robust")
# TODO: iv.first_stage, iv.wooldridge_overid, control function: add first-stage residual to OLS

import wooldridge as wd, statsmodels.formula.api as smf
from linearmodels.iv import IV2SLS
df = wd.data("mroz").query("inlf == 1").copy()
ols = smf.ols("lwage ~ educ + exper + expersq", data=df).fit(cov_type="HC1")
iv = IV2SLS.from_formula("lwage ~ 1 + exper + expersq + [educ ~ fatheduc + motheduc]", data=df).fit(cov_type="robust")
print(f"OLS  educ: {ols.params['educ']:.4f} (se {ols.bse['educ']:.4f})")
print(f"2SLS educ: {iv.params['educ']:.4f} (se {iv.std_errors['educ']:.4f})")

fs = smf.ols("educ ~ fatheduc + motheduc + exper + expersq", data=df).fit(cov_type="HC1")
print(f"first-stage robust F (2 instruments): {fs.wald_test('fatheduc = 0, motheduc = 0', scalar=True).statistic/2:.1f}")
o = iv.wooldridge_overid
print(f"over-identification (Wooldridge score): stat = {o.stat:.3f}, p = {o.pval:.3f}")

df["vhat"] = smf.ols("educ ~ fatheduc + motheduc + exper + expersq", data=df).fit().resid
cf = smf.ols("lwage ~ educ + exper + expersq + vhat", data=df).fit(cov_type="HC1")
print(f"control function: educ = {cf.params['educ']:.4f} (equals 2SLS), vhat t = {cf.tvalues['vhat']:.2f}, p = {cf.pvalues['vhat']:.3f}")

OLS  educ: 0.1075 (se 0.0132)
2SLS educ: 0.0614 (se 0.0332)
first-stage robust F (2 instruments): 49.5
over-identification (Wooldridge score): stat = 0.443, p = 0.505
control function: educ = 0.0614 (equals 2SLS), vhat t = 1.60, p = 0.110