import numpy as np, pandas as pd, matplotlib.pyplot as plt, warnings
warnings.filterwarnings("ignore")
BLUE="#2b6cb0"; RED="#c53030"; GREEN="#2f855a"; ORANGE="#dd6b20"; GREY="#a0aec0"; PURP="#6b46c1"
d=pd.read_csv("card.csv")
ctrl=["exper","expersq","black","south","smsa","reg661","reg662","reg663","reg664","reg665","reg666","reg667","reg668","smsa66"]
d=d.dropna(subset=["lwage","educ","nearc4"]+ctrl).reset_index(drop=True)
n=len(d); Y=d["lwage"].values; D=d["educ"].values; Z=d["nearc4"].values
C=np.column_stack([np.ones(n)]+[d[c].values for c in ctrl])          # controls incl. intercept
print(f"Card (1995): n = {n} men; outcome = log wage, treatment = years of education, instrument = grew up near a 4-yr college")
print(f"  {Z.mean()*100:.0f}% grew up near a 4-year college;  mean education {D.mean():.1f} yrs;  mean log wage {Y.mean():.2f}")
# raw first-stage contrast and OLS of wage on schooling (unadjusted, for intuition)
print(f"\nRELEVANCE (raw): near college -> {D[Z==1].mean():.2f} yrs educ vs {D[Z==0].mean():.2f} for those not near "
      f"(+{D[Z==1].mean()-D[Z==0].mean():.2f} yrs)")
fig,ax=plt.subplots(1,2,figsize=(13,4.2))
ax[0].hist(D[Z==1],bins=np.arange(2,19)-.5,alpha=.6,color=GREEN,density=True,label="near a 4-yr college")
ax[0].hist(D[Z==0],bins=np.arange(2,19)-.5,alpha=.5,color=GREY,density=True,label="not near")
ax[0].axvline(D[Z==1].mean(),color=GREEN,lw=2); ax[0].axvline(D[Z==0].mean(),color=GREY,lw=2)
ax[0].set_xlabel("years of education"); ax[0].set_ylabel("density"); ax[0].set_title("First stage (relevance): the instrument shifts schooling"); ax[0].legend(fontsize=8)
ax[1].scatter(D+np.random.uniform(-.25,.25,n),Y,s=4,alpha=.12,color=BLUE); ax[1].set_xlabel("years of education"); ax[1].set_ylabel("log wage")
ax[1].set_title("Wage rises with schooling — but how much is ability?");
plt.tight_layout(); plt.show()
print("The instrument clearly moves education (relevance holds). Whether it is excludable -- affects wages ONLY via education")
print("-- is an assumption we must argue, not test. That untestable exclusion restriction is the price of IV.")

Card (1995): n = 3010 men; outcome = log wage, treatment = years of education, instrument = grew up near a 4-yr college
  68% grew up near a 4-year college;  mean education 13.3 yrs;  mean log wage 6.26

RELEVANCE (raw): near college -> 13.53 yrs educ vs 12.70 for those not near (+0.83 yrs)

The instrument clearly moves education (relevance holds). Whether it is excludable -- affects wages ONLY via education
-- is an assumption we must argue, not test. That untestable exclusion restriction is the price of IV.

def ols(X,y):
    b=np.linalg.lstsq(X,y,rcond=None)[0]; r=y-X@b; s2=r@r/(len(y)-X.shape[1])
    se=np.sqrt(np.diag(s2*np.linalg.inv(X.T@X))); return b,se
# OLS: wage ~ educ + controls
Xo=np.column_stack([D,C]); bo,so=ols(Xo,Y); ols_educ,ols_se=bo[0],so[0]
# 2SLS from scratch
Zmat=np.column_stack([Z,C])                       # instruments = nearc4 + controls
Xmat=np.column_stack([D,C])                       # regressors  = educ    + controls
PZ=Zmat@np.linalg.solve(Zmat.T@Zmat, Zmat.T@Xmat)  # X-hat = projection of X on Z
b2=np.linalg.solve(PZ.T@PZ, PZ.T@Y); iv_educ=b2[0]
# proper 2SLS SE (use actual D, not D-hat, in residuals)
res=Y-Xmat@b2; s2=res@res/(n-Xmat.shape[1]); V=s2*np.linalg.inv(PZ.T@PZ); iv_se=np.sqrt(V[0,0])
# first stage F on the instrument
Dhat=PZ[:,0]; bf,_=ols(Zmat,D); resf=D-Zmat@bf
XtXinv=np.linalg.inv(Zmat.T@Zmat); s2f=resf@resf/(n-Zmat.shape[1]); seZ=np.sqrt(s2f*XtXinv[0,0]); Fstat=(bf[0]/seZ)**2
# Wald (no controls)
wald=(Y[Z==1].mean()-Y[Z==0].mean())/(D[Z==1].mean()-D[Z==0].mean())
print(f"OLS  return to schooling = {ols_educ:.4f}  (SE {ols_se:.4f})   ~{ols_educ*100:.1f}% per year")
print(f"2SLS return to schooling = {iv_educ:.4f}  (SE {iv_se:.4f})   ~{iv_educ*100:.1f}% per year   [from scratch]")
print(f"Wald (no controls)       = {wald:.4f}")
print(f"first-stage F on nearc4  = {Fstat:.2f}")
fig,ax=plt.subplots(figsize=(7.5,3.4))
ax.errorbar([ols_educ],[1],xerr=[1.96*ols_se],fmt="o",color=BLUE,capsize=6,ms=9,label=f"OLS  {ols_educ:.3f}")
ax.errorbar([iv_educ],[0.7],xerr=[1.96*iv_se],fmt="s",color=RED,capsize=6,ms=9,label=f"2SLS {iv_educ:.3f}")
ax.set_yticks([]); ax.set_ylim(0.4,1.3); ax.set_xlabel("estimated return to a year of schooling (log points)")
ax.set_title("OLS vs 2SLS — Card returns to schooling"); ax.legend()
plt.tight_layout(); plt.show()
print("\n2SLS EXCEEDS OLS here -- surprising if you expect ability to bias OLS upward. The LATE framework (next) explains why:")
print("2SLS estimates the return for a particular subgroup (compliers), not the population average.")

OLS  return to schooling = 0.0747  (SE 0.0035)   ~7.5% per year
2SLS return to schooling = 0.1315  (SE 0.0550)   ~13.2% per year   [from scratch]
Wald (no controls)       = 0.1881
first-stage F on nearc4  = 13.26

2SLS EXCEEDS OLS here -- surprising if you expect ability to bias OLS upward. The LATE framework (next) explains why:
2SLS estimates the return for a particular subgroup (compliers), not the population average.

# complier share = first-stage jump in P(treatment) — here treat 'some college' (educ>=13) as the margin the instrument moves
T=(D>=13).astype(int)                                   # crossed into college
comp_share=T[Z==1].mean()-T[Z==0].mean()               # Wald first stage on the binary margin
always=T[Z==0].mean(); never=1-T[Z==1].mean()
print(f"Decomposition at the 'entered college (educ>=13)' margin the instrument moves:")
print(f"  P(educ>=13 | near a college)          = {T[Z==1].mean():.3f}")
print(f"  P(educ>=13 | not near)                = {T[Z==0].mean():.3f}")
print(f"  complier share  (the difference)      = {comp_share:.3f}")
print(f"  always-takers   (college regardless)  = {always:.3f}")
print(f"  never-takers    (no college regardless)= {never:.3f}")
fig,ax=plt.subplots(1,2,figsize=(13,4.2))
ax[0].bar(["compliers","always-takers","never-takers"],[comp_share,always,never],color=[GREEN,BLUE,GREY])
for i,v in enumerate([comp_share,always,never]): ax[0].text(i,v+0.01,f"{v:.2f}",ha="center")
ax[0].set_ylabel("share of sample"); ax[0].set_title("Who does the instrument move? (LATE is the complier effect)")
# schematic: LATE vs ATE vs OLS
ax[1].bar(["OLS\n(all, biased)","2SLS = LATE\n(compliers)"],[ols_educ,iv_educ],color=[BLUE,RED])
for i,v in enumerate([ols_educ,iv_educ]): ax[1].text(i,v+0.003,f"{v:.3f}",ha="center")
ax[1].set_ylabel("return to schooling"); ax[1].set_title("2SLS estimates the LATE, not the population ATE")
plt.tight_layout(); plt.show()
print("\nIV's estimand is local: the return for the compliers the instrument shifts. For proximity-to-college those are")
print("disproportionately credit-constrained students whose return is high -- why the Card LATE lands above the OLS average.")

Decomposition at the 'entered college (educ>=13)' margin the instrument moves:
  P(educ>=13 | near a college)          = 0.544
  P(educ>=13 | not near)                = 0.422
  complier share  (the difference)      = 0.122
  always-takers   (college regardless)  = 0.422
  never-takers    (no college regardless)= 0.456

IV's estimand is local: the return for the compliers the instrument shifts. For proximity-to-college those are
disproportionately credit-constrained students whose return is high -- why the Card LATE lands above the OLS average.

rng=np.random.default_rng(0); beta=0.5; reps=300
strengths=np.array([0.02,0.05,0.1,0.15,0.25,0.4,0.6])
def twosls(y,dd,z):
    z1=np.column_stack([np.ones(len(z)),z]); x1=np.column_stack([np.ones(len(dd)),dd])
    ph=z1@np.linalg.solve(z1.T@z1,z1.T@x1); b=np.linalg.solve(ph.T@ph,ph.T@y); return b[1]
def olsslope(y,dd):
    x=np.column_stack([np.ones(len(dd)),dd]); return np.linalg.solve(x.T@x,x.T@y)[1]
ivm=[];olm=[];Fm=[]
for pi in strengths:
    ivs=[];ols_=[];Fs=[]
    for _ in range(reps):
        m=800; U=rng.normal(0,1,m); Z_=rng.normal(0,1,m)
        Dd=pi*Z_+U+rng.normal(0,1,m)                    # U (unobserved) makes D endogenous
        Yy=beta*Dd+U+rng.normal(0,1,m)                  # U also in outcome; Z excluded
        ivs.append(twosls(Yy,Dd,Z_)); ols_.append(olsslope(Yy,Dd))
        z1=np.column_stack([np.ones(m),Z_]); bf=np.linalg.solve(z1.T@z1,z1.T@Dd); rf=Dd-z1@bf
        se=np.sqrt((rf@rf/(m-2))*np.linalg.inv(z1.T@z1)[1,1]); Fs.append((bf[1]/se)**2)
    ivm.append(np.mean(ivs)); olm.append(np.mean(ols_)); Fm.append(np.mean(Fs))
ivm=np.array(ivm); Fm=np.array(Fm)
fig,ax=plt.subplots(1,2,figsize=(13,4.4))
ax[0].plot(Fm,ivm,"o-",color=RED,lw=2,label="mean 2SLS estimate"); ax[0].axhline(beta,color=GREEN,ls="--",lw=2,label=f"true effect = {beta}")
ax[0].axhline(np.mean(olm),color=BLUE,ls=":",lw=2,label=f"OLS (confounded) ≈ {np.mean(olm):.2f}"); ax[0].axvline(10,color="k",ls=":",label="F=10 rule")
ax[0].set_xscale("log"); ax[0].set_xlabel("first-stage F (instrument strength)"); ax[0].set_ylabel("2SLS estimate"); ax[0].set_title("Weak instruments bias 2SLS toward OLS"); ax[0].legend(fontsize=8)
ax[1].plot(strengths,Fm,"s-",color=PURP,lw=2); ax[1].axhline(10,color="k",ls=":",label="F=10 (Staiger-Stock)")
ax[1].set_xlabel("true instrument strength π"); ax[1].set_ylabel("first-stage F"); ax[1].set_title("First-stage F falls with instrument strength"); ax[1].legend(fontsize=8)
plt.tight_layout(); plt.show()
print(f"Card's real first-stage F ≈ {Fstat:.1f} (adequate). In the simulation, when F drops below ~10 the 2SLS estimate")
print(f"collapses from the truth ({beta}) toward the confounded OLS value -- the Bound-Jaeger-Baker weak-instrument trap.")
print("ALWAYS report the first-stage F: a weak instrument is worse than none, because it looks like an answer.")

Card's real first-stage F ≈ 13.3 (adequate). In the simulation, when F drops below ~10 the 2SLS estimate
collapses from the truth (0.5) toward the confounded OLS value -- the Bound-Jaeger-Baker weak-instrument trap.
ALWAYS report the first-stage F: a weak instrument is worse than none, because it looks like an answer.

Causal Inference III — Instrumental Variables¶

When the confounder is unobserved: 2SLS, the LATE framework, and weak instruments¶

1. The endogeneity problem and the instrument¶

2. Two-stage least squares — from scratch¶

3. What IV actually estimates — the LATE framework (Imbens & Angrist, 1994)¶

4. Weak instruments — the diagnostic that made IV honest¶

5. Summary¶