Causal Inference: Potential Outcomes, DiD, and Regression Discontinuity

1. Potential outcomes

Each unit \(i\) has potential outcomes \(Y_i(1)\) (treated) and \(Y_i(0)\) (untreated). We observe \(Y_i=D_iY_i(1)+(1-D_i)Y_i(0)\). The individual effect \(Y_i(1)-Y_i(0)\) is never observed: the fundamental problem of causal inference. The estimands are averages:

The notation embeds SUTVA: one unit's treatment does not affect another's outcomes (no spillovers), and there is a single version of treatment.

Selection bias. The naive difference in means decomposes as

Randomization makes \(D\) independent of \((Y(1),Y(0))\), so the selection bias is zero and ATT = ATE. In a randomized experiment, regressing \(Y\) on \(D\) and pre-treatment covariates improves precision without introducing bias asymptotically. The fully interacted specification (Lin, 2013) is robust to heterogeneous effects.

2. Selection on observables

Under these assumptions, \(ATE=E\big[\mu_1(X)-\mu_0(X)\big]\), with \(\mu_d(x)=E[Y\mid D=d,X=x]\). Three estimation approaches follow.

• Regression adjustment: fit \(\widehat\mu_1\) and \(\widehat\mu_0\) separately, then average \(\widehat\mu_1(X_i)-\widehat\mu_0(X_i)\). A single regression of \(Y\) on \(D\) and \(X\) instead estimates a variance-weighted average of conditional effects, which differs from the ATE when effects are heterogeneous (Angrist, 1998; Słoczyński, 2022).
• Inverse probability weighting (IPW): \(\widehat{ATE}=n^{-1}\sum_i\big[D_iY_i/\widehat e(X_i)-(1-D_i)Y_i/(1-\widehat e(X_i))\big]\). This is unstable when \(\widehat e\) is near 0 or 1; use normalized (Hájek) weights.
• Doubly robust (AIPW):

It is consistent if either the outcome model or the propensity model is correct. It is efficient if both are. Combined with cross-fitting, it allows flexible machine-learning estimates of the nuisance functions while preserving \(\sqrt n\) inference (double/debiased machine learning; Chernozhukov et al., 2018).

Matching on covariates or on the propensity score is an alternative. Use bias-corrected matching estimators (Abadie and Imbens, 2011); bootstrap standard errors are invalid for nearest-neighbor matching.

Overlap in practice. Plot \(\widehat e(X)\) by treatment group. Trim units with extreme scores (Crump et al., 2009), and report how the estimand changes. Sensitivity: unconfoundedness is untestable. Report how strong unobserved confounding would have to be to overturn the result (Imbens, 2003; Cinelli and Hazlett, 2020).

3. Difference-in-differences

The 2×2 design. Two groups (treated \(G=1\), control \(G=0\)) and two periods (pre, post). The DiD estimand is

The regression \(Y_{it}=\alpha+\beta G_i+\gamma Post_t+\tau(G_i\times Post_t)+\varepsilon_{it}\) gives \(\widehat\tau=\widehat\tau_{DiD}\). Parallel trends is an assumption about untreated potential outcomes, so it is untestable in the post period and scale-dependent: it cannot hold in both levels and logs unless trends are zero (Roth and Sant'Anna, 2023). Choose the scale on economic grounds. Cluster at the level of treatment assignment (Lecture 4).

Event studies. With multiple periods, estimate dynamic effects relative to the period before treatment:

where \(E_i\) is the treatment date. Pre-period coefficients near zero are consistent with parallel trends, but pre-trend tests have low power, and conditioning on passing them can distort inference (Roth, 2022). The approach of Rambachan and Roth (2023) gives inference that is robust to bounded violations of parallel trends.

Staggered adoption and the TWFE problem. When units adopt treatment at different times, the static TWFE regression \(Y_{it}=\alpha_i+\lambda_t+\tau D_{it}+\varepsilon_{it}\) is a weighted average of all 2×2 DiDs in the data. These include comparisons that use already-treated units as controls (Goodman-Bacon, 2021). If effects grow over time, those "forbidden comparisons" subtract treatment effects, and some implicit weights on group-time effects are negative. \(\widehat\tau\) can then have the wrong sign even if every unit's effect is positive (de Chaisemartin and D'Haultfœuille, 2020). Dynamic TWFE event studies suffer from the same contamination (Sun and Abraham, 2021). Heterogeneity-robust alternatives:

• Callaway and Sant'Anna (2021): estimate group-time ATTs, \(ATT(g,t)\), using never-treated or not-yet-treated units as controls, then aggregate with explicit, non-negative weights;
• Sun and Abraham (2021): an interaction-weighted event-study estimator;
• Borusyak, Jaravel and Spiess (2024): imputation — fit unit and time effects on untreated observations only, impute \(Y(0)\) for the treated, and average the differences.

4. Synthetic control

With a single (or few) treated aggregate unit(s) — a country, a state — and a long pre-period, Abadie, Diamond and Hainmueller (2010) construct a synthetic control: a convex combination of donor units, with weights \(w_j\ge0\) and \(\sum w_j=1\), chosen to match the treated unit's pre-treatment outcomes and predictors. The effect estimate is the post-treatment gap between the treated unit and its synthetic counterpart. The convexity constraint prevents extrapolation and makes weights interpretable. Inference uses placebo (permutation) tests: apply the method to each donor and compare the treated unit's post/pre RMSPE ratio with the placebo distribution. Extensions include synthetic DiD (Arkhangelsky et al., 2021), which combines unit and time weights, and conformal inference (Chernozhukov, Wüthrich and Zhu, 2021).

5. Regression discontinuity

Treatment is determined, fully or partly, by whether a running variable \(X\) crosses a cutoff \(c\).

Sharp RD: \(D=1\{X\ge c\}\). If \(E[Y(0)\mid X=x]\) and \(E[Y(1)\mid X=x]\) are continuous at \(c\), then

Fuzzy RD: the probability of treatment jumps at \(c\) but by less than one. The estimand is the ratio of the outcome jump to the treatment jump, a LATE for compliers at the cutoff. It is IV with \(1\{X\ge c\}\) as the instrument.

Estimation. Use local linear regression on each side within a bandwidth \(h\), with a triangular kernel. Avoid global high-order polynomials, which give noisy estimates and misleading inference (Gelman and Imbens, 2019). Choose \(h\) to minimize MSE (Imbens and Kalyanaraman, 2012). Because the MSE-optimal bandwidth leaves a bias of the same order as the standard error, conventional confidence intervals undercover. Use robust bias-corrected inference (Calonico, Cattaneo and Titiunik, 2014), implemented in rdrobust.

Validation.

• Manipulation: if agents can precisely sort around the cutoff, continuity fails. Test for a discontinuity in the density of \(X\) at \(c\) (McCrary, 2008; Cattaneo, Jansson and Ma, 2020).
• Covariate balance: pre-determined covariates should not jump at \(c\). Run the RD with each covariate as the outcome.
• Placebo cutoffs and bandwidth sensitivity: no effects at fake cutoffs, and stable estimates across a range of bandwidths.
• External validity: the estimate is local to \(X=c\). Extrapolation away from the cutoff requires further assumptions.

import numpy as np
rng = np.random.default_rng(13)
N, T = 300, 10
cohort = rng.choice([4, 7, 99], size=N)                  # adoption period (99 = never treated)
alpha = rng.normal(size=N)[:, None]
lam = np.linspace(0, 1, T)[None, :]
t = np.arange(T)[None, :]
D = (t >= cohort[:, None]).astype(float)
effect = D * 0.5 * (t - cohort[:, None] + 1)              # effects grow with time since adoption
Y = alpha + lam + effect + rng.normal(0, 0.5, (N, T))

def twfe(Y, D):
    dm = lambda M: M - M.mean(1, keepdims=True) - M.mean(0, keepdims=True) + M.mean()
    y, d = dm(Y), dm(D)
    return (y * d).sum() / (d * d).sum()

true_att = effect[D == 1].mean()
print("true ATT (average over treated cells):", round(true_att, 3))
print("static TWFE:", round(twfe(Y, D), 3))

# imputation estimator: fit unit + time effects on untreated cells only, impute Y(0)
mask = D == 0
a, l = np.zeros(N), np.zeros(T)
for _ in range(200):                                      # alternating least squares
    a = np.nansum(np.where(mask, Y - l, np.nan), 1) / mask.sum(1)
    l = np.nansum(np.where(mask, Y - a[:, None], np.nan), 0) / mask.sum(0)
tau_hat = (Y - a[:, None] - l[None, :])[D == 1].mean()
print("imputation estimator:", round(tau_hat, 3))

import numpy as np, wooldridge as wd, statsmodels.formula.api as smf
exp_, obs = wd.data("jtrain2"), wd.data("jtrain3")
X = "age + agesq + educ + black + hisp + married + re74 + re75 + unem74 + unem75"
# TODO: experimental ATE; then on obs: ps = smf.logit(f"train ~ {X}", obs).fit().predict()

import numpy as np, wooldridge as wd, statsmodels.formula.api as smf
exp_, obs = wd.data("jtrain2"), wd.data("jtrain3")
exp_["agesq"] = exp_["age"]**2
X = "age + agesq + educ + black + hisp + married + re74 + re75 + unem74 + unem75"
e1 = smf.ols("re78 ~ train", data=exp_).fit(cov_type="HC1"); e2 = smf.ols(f"re78 ~ train + {X}", data=exp_).fit(cov_type="HC1")
print(f"EXPERIMENT   diff in means {e1.params['train']:.3f} (se {e1.bse['train']:.3f});  with covariates {e2.params['train']:.3f}")

def aipw(d):
    ps = smf.logit(f"train ~ {X}", data=d).fit(disp=0).predict(d)
    m1 = smf.ols(f"re78 ~ {X}", data=d[d.train == 1]).fit().predict(d)
    m0 = smf.ols(f"re78 ~ {X}", data=d[d.train == 0]).fit().predict(d)
    D, Y = d["train"], d["re78"]
    ra_att = (Y - m0)[D == 1].mean()                                     # regression adjustment, ATT
    w = ps / (1 - ps)                                                    # doubly robust ATT
    dr_att = ((D * (Y - m0)).sum() - ((1 - D) * w * (Y - m0)).sum()) / D.sum()
    return ps, ra_att, dr_att

naive = obs.loc[obs.train == 1, "re78"].mean() - obs.loc[obs.train == 0, "re78"].mean()
ps, ra, dr = aipw(obs)
print(f"PSID CONTROLS naive {naive:.3f};  regression adjustment ATT {ra:.3f};  doubly robust ATT {dr:.3f}")
print(f"share of PSID controls with propensity < 0.01: {(ps[obs.train == 0] < 0.01).mean():.3f}")
trim = obs[(ps >= 0.01) & (ps <= 0.99)].copy()
_, ra_t, dr_t = aipw(trim)
print(f"TRIMMED (n = {len(trim)}, treated {int(trim.train.sum())}): regression adjustment {ra_t:.3f}; doubly robust {dr_t:.3f}")

EXPERIMENT   diff in means 1.794 (se 0.671);  with covariates 1.670
PSID CONTROLS naive -15.205;  regression adjustment ATT 0.894;  doubly robust ATT 3.654
share of PSID controls with propensity < 0.01: 0.814
TRIMMED (n = 646, treated 182): regression adjustment 1.785; doubly robust 2.497

import numpy as np, wooldridge as wd, statsmodels.formula.api as smf
df = wd.data("kielmc")
means = df.groupby(["nearinc", "y81"])["rprice"].mean().unstack()
# TODO: DiD from means; regression with interaction; add controls

import numpy as np, wooldridge as wd, statsmodels.formula.api as smf
df = wd.data("kielmc")
mu = df.groupby(["nearinc", "y81"])["rprice"].mean().unstack()
print(mu.round(0).rename(index={0: "far", 1: "near"}, columns={0: "1978", 1: "1981"}))
did = (mu.loc[1, 1] - mu.loc[1, 0]) - (mu.loc[0, 1] - mu.loc[0, 0])
print(f"\nDiD from means: {did:,.0f} dollars")
r1 = smf.ols("rprice ~ nearinc * y81", data=df).fit(cov_type="HC1")
print(f"regression:       {r1.params['nearinc:y81']:,.0f} (robust se {r1.bse['nearinc:y81']:,.0f})")
ctrl = "age + I(age**2) + intst + land + area + rooms + baths"
r2 = smf.ols(f"rprice ~ nearinc * y81 + {ctrl}", data=df).fit(cov_type="HC1")
print(f"with controls:    {r2.params['nearinc:y81']:,.0f} (robust se {r2.bse['nearinc:y81']:,.0f})")
r3 = smf.ols(f"np.log(rprice) ~ nearinc * y81 + {ctrl}", data=df).fit(cov_type="HC1")
b = r3.params["nearinc:y81"]
print(f"log specification with controls: {b:.3f} -> {100*(np.exp(b)-1):.1f}% (robust se {r3.bse['nearinc:y81']:.3f})")

y81         1978      1981
nearinc                   
far      82517.0  101308.0
near     63693.0   70619.0

DiD from means: -11,864 dollars
regression:       -11,864 (robust se 8,636)
with controls:    -14,178 (robust se 6,512)
log specification with controls: -0.093 -> -8.8% (robust se 0.061)