First Course Ding: Chapter 23

An econometric perspective on instrumental variables

Chapter 23 is where the design-based IV story becomes an econometrics workflow: controls, overidentification, and alternative representations of the same parameter.

Show code 
from pathlib import Path
import math

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd

import crabbymetrics as cm

np.set_printoptions(precision=4, suppress=True)


def iv_moments(theta, data):
    resid = data["y"] - data["x"] @ theta
    return data["z"] * resid[:, None]


def iv_jacobian(theta, data):
    del theta
    return -(data["z"].T @ data["x"]) / data["x"].shape[0]

1 Card’s Returns-To-Schooling Data

Show code 
card = pd.read_csv(Path("../data/ding") / "card1995.csv")
y = card["lwage"].to_numpy(dtype=float)
d = card["educ"].to_numpy(dtype=float)[:, None]

controls = [
    "exper",
    "expersq",
    "black",
    "south",
    "smsa",
    "reg661",
    "reg662",
    "reg663",
    "reg664",
    "reg665",
    "reg666",
    "reg667",
    "reg668",
    "smsa66",
]
x = card[controls].to_numpy(dtype=float)
z_excl = card[["nearc2", "nearc4"]].to_numpy(dtype=float)

keep = np.isfinite(y) & np.isfinite(d[:, 0]) & np.isfinite(x).all(axis=1) & np.isfinite(z_excl).all(axis=1)
y = y[keep]
d = d[keep]
x = x[keep]
z_excl = z_excl[keep]

2 Closed-Form TwoSLS And Overidentified GMM

Show code 
ols = cm.OLS()
ols.fit(np.column_stack([d[:, 0], x]), y)
ols_hc1 = ols.summary(vcov="hc1")

tsls = cm.TwoSLS()
tsls.fit(d, x, z_excl, y)
tsls_hc1 = tsls.summary(vcov="hc1")

features = np.column_stack([d[:, 0], x])
feature_scale = features.std(axis=0)
instruments = np.column_stack([x, z_excl])
X = np.column_stack([np.ones(len(y)), (features - features.mean(axis=0)) / feature_scale])
Z = np.column_stack([
    np.ones(len(y)),
    (instruments - instruments.mean(axis=0)) / instruments.std(axis=0),
])

gmm = cm.GMM(iv_moments, jacobian_fn=iv_jacobian, max_iterations=200, tolerance=1e-8)
gmm.fit({"x": X, "y": y, "z": Z}, np.zeros(X.shape[1]))
gmm_summary = gmm.summary(vcov="sandwich", omega="iid")

print("OLS return to education:", round(float(ols_hc1["coef"][0]), 4))
print("OLS HC1 SE:", round(float(ols_hc1["coef_se"][0]), 4))
print("TwoSLS return to education:", round(float(tsls_hc1["coef"][0]), 4))
print("TwoSLS HC1 SE:", round(float(tsls_hc1["coef_se"][0]), 4))
print("Two-step GMM return to education:", round(float(gmm_summary["coef"][1] / feature_scale[0]), 4))
print("Two-step GMM SE:", round(float(gmm_summary["se"][1] / feature_scale[0]), 4))
print("GMM J-statistic:", round(float(gmm_summary["j_stat"]), 4), "with df =", gmm_summary["j_df"])
OLS return to education: 0.0747
OLS HC1 SE: 0.0036
TwoSLS return to education: 0.1571
TwoSLS HC1 SE: 0.0526
Two-step GMM return to education: 0.1552
Two-step GMM SE: 0.0522
GMM J-statistic: 1.2988 with df = 1

The GMM calculation centers and scales nonconstant regressors and instruments: squared experience and binary indicators otherwise have very different scales. The education coefficient and SE are transformed back to log-wage per year of schooling. Scaling changes the finite-sample identity-weighted first step, so this two-step estimate is not required to equal TwoSLS exactly. A tighter tolerance on the raw, poorly scaled system can fail the v0.9 convergence guard even with a small score; a small score alone is not enough.

3 Control-Function View

Show code 
first_stage = cm.OLS()
first_stage.fit(np.column_stack([x, z_excl]), d[:, 0])
pi = first_stage.summary()
d_hat = first_stage.predict(np.column_stack([x, z_excl]))
resid = d[:, 0] - d_hat

control_function = cm.OLS()
control_function.fit(np.column_stack([d[:, 0], resid, x]), y)
cf_hc1 = control_function.summary(vcov="hc1")

print("TwoSLS coefficient:", round(float(tsls_hc1["coef"][0]), 4))
print("Control-function coefficient on education:", round(float(cf_hc1["coef"][0]), 4))
TwoSLS coefficient: 0.1571
Control-function coefficient on education: 0.1571

4 Anderson-Rubin Grid

The source R script also inverts a reduced-form test over candidate education returns. For a candidate \(\beta\), regress \(Y - \beta D\) on controls and the excluded instrument. If the excluded instrument is still predictive, that candidate \(\beta\) is inconsistent with the IV moment. This grid uses nearc4, matching the script’s single-instrument Anderson-Rubin calculation.

Show code 
def normal_pvalue(z_stat):
    return math.erfc(abs(float(z_stat)) / math.sqrt(2.0))


beta_grid = np.linspace(0.0, 0.35, 176)
ar_pvalues = []
z_nearc4 = z_excl[:, [1]]
for beta in beta_grid:
    y_beta = y - beta * d[:, 0]
    ar = cm.OLS()
    ar.fit(np.column_stack([z_nearc4[:, 0], x]), y_beta)
    ar_s = ar.summary(vcov="hc1")
    ar_pvalues.append(normal_pvalue(ar_s["coef"][0] / ar_s["coef_se"][0]))

ar_pvalues = np.asarray(ar_pvalues)
accepted = beta_grid[ar_pvalues >= 0.05]

pd.DataFrame(
    {
        "lower": [accepted.min()],
        "upper": [accepted.max()],
        "grid_points_accepted": [accepted.size],
    },
    index=["AR 95% confidence set on grid"],
)
lower upper grid_points_accepted
AR 95% confidence set on grid 0.03 0.28 126
Show code 
fig, ax = plt.subplots(figsize=(6, 4))
ax.plot(beta_grid, ar_pvalues)
ax.axhline(0.05, color="black", linestyle="--", linewidth=1.0)
ax.set_xlabel("Candidate return to schooling")
ax.set_ylabel("Reduced-form p-value")
ax.set_title("Anderson-Rubin inversion with nearc4")
fig.tight_layout()

By Chapter 23, the IV problem has several equivalent front ends:

  • closed-form TwoSLS
  • overidentified linear IV as two-step GMM
  • a control-function representation

In crabbymetrics, those are now all available without leaving the same Rust-backed core.