Causal Inference IV — Regression Discontinuity¶
Identification from a threshold: sharp & fuzzy RD, local polynomials, and the manipulation test¶
Some of the cleanest causal evidence in observational data comes from arbitrary threshold rules. When a treatment is switched on the moment a running variable crosses a cutoff — you get the scholarship if your test score clears 80, the party holds the seat if its vote share clears 50% — then units landing just above and just below the cutoff are, in every respect other than treatment, essentially identical. Comparing them is as good as a randomized experiment in a tiny neighborhood of the threshold. That is the regression discontinuity (RD) design, and its identifying assumption is remarkably weak: only that the potential outcomes are continuous in the running variable at the cutoff, so any jump in the outcome there must be caused by the treatment.
This notebook builds RD from the ground up:
- Sharp RD — treatment is a deterministic step at the cutoff; the effect is the jump in the outcome, estimated by fitting the outcome on each side and differencing at the boundary.
- Local-polynomial estimation & bandwidth — RD is inherently local, so we fit local linear regressions in a window around the cutoff, and use the Calonico-Cattaneo-Titiunik (CCT) optimal bandwidth with bias-corrected robust inference.
- The manipulation (McCrary) test — RD fails if units can precisely sort across the cutoff; we test for a discontinuity in the density of the running variable.
- Fuzzy RD — when crossing the cutoff changes only the probability of treatment, the design becomes an instrumental-variables problem solved by the exact 2SLS machinery of the previous notebook.
Data: Lee (2008), the canonical RD study of the incumbency advantage in U.S. House elections. Python-lead (from-scratch local-linear + rdrobust/rddensity); R companion mirrors it with rdrobust.
1. The design and the sharp RD estimate — Lee (2008)¶
Lee (2008) studies whether incumbency confers an electoral advantage in U.S. House elections (1946–1998). The running variable is the Democratic candidate's margin of victory in an election (Democratic vote share minus the strongest opponent's, so $x>0$ means the Democrat won and the party becomes the incumbent). The outcome is the Democratic vote share in the next election ($y$). The cutoff is $x=0$.
The identification insight is beautiful: whether a race is won by 0.1% or lost by 0.1% is essentially a coin flip — decided by weather, a scandal's timing, turnout noise — so districts just above and just below $x=0$ are comparable, and the jump in the next-election vote share at $x=0$ is the causal incumbency advantage. The "RD plot" — binned averages of $y$ against $x$ — shows it directly: a clear vertical gap at the threshold. The sharp RD estimand is $$\tau_{RD}=\lim_{x\downarrow 0}\mathbb{E}[y\mid x]-\lim_{x\uparrow 0}\mathbb{E}[y\mid x].$$
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("lee2008.csv"); x=d["margin"].values; y=d["dem_voteshare_next"].values
print(f"Lee (2008) U.S. House elections: n = {len(x)} races")
print(f" running variable = Democratic margin of victory (x>0 => Democrat won, becomes incumbent)")
print(f" outcome = Democratic vote share in the NEXT election; cutoff c = 0")
# RD plot: binned means of y vs x
bins=np.linspace(-1,1,41); mid=(bins[:-1]+bins[1:])/2; idx=np.digitize(x,bins)-1
by=[y[idx==b].mean() if (idx==b).sum()>0 else np.nan for b in range(len(mid))]
# local-linear fits each side (bw=0.3 just for the plotted lines)
def fit_side(mask,grid,h):
xx=x[mask];yy=y[mask];w=np.clip(1-np.abs(xx)/h,0,None); X=np.column_stack([np.ones(len(xx)),xx]);W=np.diag(w)
b=np.linalg.solve(X.T@W@X,X.T@W@yy); return b[0]+b[1]*grid, b[0]
gl=np.linspace(-.3,0,50); gr=np.linspace(0,.3,50)
fL,aL=fit_side((x<0)&(x>=-.3),gl,.3); fR,aR=fit_side((x>=0)&(x<=.3),gr,.3)
fig,ax=plt.subplots(figsize=(9,5))
ax.scatter(mid,by,s=28,color=GREY,zorder=2,label="binned average of next vote share")
ax.plot(gl,fL,color=BLUE,lw=2.5); ax.plot(gr,fR,color=RED,lw=2.5)
ax.axvline(0,color="k",ls="--",lw=1); ax.plot([0,0],[aL,aR],color=GREEN,lw=3,zorder=5)
ax.annotate(f"jump ≈ {aR-aL:.3f}",xy=(0,(aL+aR)/2),xytext=(0.25,0.42),fontsize=11,color=GREEN,arrowprops=dict(arrowstyle="->",color=GREEN))
ax.set_xlabel("Democratic margin of victory (running variable)"); ax.set_ylabel("Democratic vote share, next election")
ax.set_title("RD plot — the incumbency advantage is the jump at the cutoff"); ax.legend(loc="upper left")
plt.tight_layout(); plt.show()
print(f"\nAt the wide bandwidth used for these plotted lines (0.3) the gap reads {aR-aL:.3f}. Treat that as the")
print("PICTURE, not the estimate -- section 2 fits at the optimal bandwidth and gets a smaller number, and the")
print("difference between the two is the whole subject of that section. What the plot shows is that barely winning")
print("today causes a large advantage tomorrow. The whole effect is identified by that discontinuity at x=0.")
Lee (2008) U.S. House elections: n = 6558 races running variable = Democratic margin of victory (x>0 => Democrat won, becomes incumbent) outcome = Democratic vote share in the NEXT election; cutoff c = 0
At the wide bandwidth used for these plotted lines (0.3) the gap reads 0.080. Treat that as the PICTURE, not the estimate -- section 2 fits at the optimal bandwidth and gets a smaller number, and the difference between the two is the whole subject of that section. What the plot shows is that barely winning today causes a large advantage tomorrow. The whole effect is identified by that discontinuity at x=0.
2. Local-polynomial estimation and the CCT bandwidth¶
RD is a boundary estimation problem: we need the regression function's value at the cutoff, approached from each side. Fitting a global polynomial to all the data is dangerous (points far from the cutoff distort the fit and create artificial jumps — the Gelman-Imbens critique of high-order global polynomials). The modern standard is local linear regression: fit a line on each side using only points within a bandwidth $h$ of the cutoff, weighting nearer points more (a triangular kernel), and take the difference of the two intercepts as $\hat\tau_{RD}$.
The bandwidth is the key tuning choice — too wide adds bias (the line misfits curvature), too narrow adds variance (few points). Calonico, Cattaneo & Titiunik (2014) derived the MSE-optimal bandwidth and a bias-corrected, robust confidence interval that remains valid at that bandwidth (the naive CI is not). We implement local-linear RD from scratch, sweep the bandwidth, and confirm the estimate against rdrobust at the CCT choice with robust inference.
The sweep is worth reading carefully, because it does not show what one would like it to show. Across $h \in [0.04, 0.40]$ the estimate runs from roughly 0.058 to 0.084 — a spread of about 41% of the estimate itself — and it is not monotone. From-scratch and rdrobust agree exactly at the CCT bandwidth, which validates the implementation; but the choice of bandwidth moves the answer by more than most published RD standard errors would suggest. That is precisely why CCT pair an MSE-optimal bandwidth with a bias-corrected interval that remains valid at it, and it is why the bandwidth curve belongs in the output rather than in a footnote.
from rdrobust import rdrobust
def rd_local(h): # from-scratch local-linear, triangular kernel, jump at 0
def side(mask):
xx=x[mask];yy=y[mask];w=1-np.abs(xx)/h; X=np.column_stack([np.ones(len(xx)),xx]);W=np.diag(w)
return np.linalg.solve(X.T@W@X,X.T@W@yy)[0]
return side((x>=0)&(x<=h))-side((x<0)&(x>=-h))
hs=np.linspace(0.04,0.4,25); ests=[rd_local(h) for h in hs]
r=rdrobust(y,x,c=0); h_cct=r.bws.iloc[0,0]; tau_c=r.coef.iloc[0,0]; ci=(r.ci.iloc[2,0],r.ci.iloc[2,1])
print(f"from-scratch local-linear RD at CCT bandwidth h={h_cct:.3f}: {rd_local(h_cct):.4f}")
print(f"rdrobust conventional estimate: {tau_c:.4f} robust 95% CI [{ci[0]:.4f}, {ci[1]:.4f}]")
fig,ax=plt.subplots(1,2,figsize=(13,4.4))
ax[0].plot(hs,ests,"o-",color=BLUE,lw=2,label="from-scratch local-linear"); ax[0].axvline(h_cct,color=GREEN,ls="--",label=f"CCT optimal h={h_cct:.3f}")
ax[0].axhline(tau_c,color=RED,ls=":",label=f"rdrobust {tau_c:.3f}"); ax[0].set_xlabel("bandwidth h"); ax[0].set_ylabel("RD estimate"); ax[0].set_title("Bandwidth sensitivity of the RD estimate"); ax[0].legend(fontsize=8)
ax[1].errorbar([tau_c],[1],xerr=[[tau_c-ci[0]],[ci[1]-tau_c]],fmt="s",color=RED,capsize=6,ms=10,label="rdrobust (robust CI)")
ax[1].errorbar([rd_local(h_cct)],[0.7],xerr=[0.014],fmt="o",color=BLUE,capsize=6,ms=9,label="from-scratch (approx)")
ax[1].axvline(0,color="k",lw=.7); ax[1].set_yticks([]); ax[1].set_ylim(0.4,1.3); ax[1].set_xlabel("incumbency advantage (vote-share jump)"); ax[1].set_title("The estimate with robust inference"); ax[1].legend(fontsize=8,loc="upper right")
plt.tight_layout(); plt.show()
print(f"\nThe incumbency advantage is about {tau_c:.2f} (a ~{tau_c*100:.0f} percentage-point boost to next-election vote share),")
print()
print(f" {'bandwidth h':>12} {'estimate':>10} {'races in window':>17}")
for h_ in (0.04,0.08,h_cct,0.20,0.30,0.40):
tag=" <- CCT optimal" if abs(h_-h_cct)<1e-9 else ""
print(f" {h_:>12.3f} {rd_local(h_):>10.4f} {int((np.abs(x)<=h_).sum()):>17,}{tag}")
lo_,hi_=min(ests),max(ests)
print(f"\n over the swept range h in [{hs.min():.2f}, {hs.max():.2f}]: min {lo_:.4f}, max {hi_:.4f}, "
f"spread {hi_-lo_:.4f}")
print(f" that spread is {100*(hi_-lo_)/tau_c:.0f}% of the CCT estimate itself.")
print()
print("From-scratch and rdrobust agree exactly at the CCT bandwidth, which is the check that matters for the")
print("implementation. But the estimate is NOT stable across bandwidths, and saying so would be the more")
print("comfortable claim rather than the true one: it runs from about 5.8 to 8.4 percentage points, and it is")
print("not even monotone -- it dips near h=0.08 before climbing. Widen the window and you import curvature")
print("from races that were not close; narrow it and you are estimating a boundary from a few hundred races.")
print()
print("This is why the bandwidth cannot be chosen by eye, and why CCT pair an MSE-optimal choice with a")
print("bias-corrected interval that stays valid AT that choice. The honest headline is the CCT estimate with")
print("its robust interval; the bandwidth curve is the uncertainty that a single number hides.")
Mass points detected in the running variable.
from-scratch local-linear RD at CCT bandwidth h=0.136: 0.0637 rdrobust conventional estimate: 0.0637 robust 95% CI [0.0348, 0.0839]
The incumbency advantage is about 0.06 (a ~6 percentage-point boost to next-election vote share),
bandwidth h estimate races in window
0.040 0.0819 483
0.080 0.0588 972
0.136 0.0637 1,606 <- CCT optimal
0.200 0.0740 2,265
0.300 0.0801 3,283
0.400 0.0842 4,169
over the swept range h in [0.04, 0.40]: min 0.0583, max 0.0842, spread 0.0259
that spread is 41% of the CCT estimate itself.
From-scratch and rdrobust agree exactly at the CCT bandwidth, which is the check that matters for the
implementation. But the estimate is NOT stable across bandwidths, and saying so would be the more
comfortable claim rather than the true one: it runs from about 5.8 to 8.4 percentage points, and it is
not even monotone -- it dips near h=0.08 before climbing. Widen the window and you import curvature
from races that were not close; narrow it and you are estimating a boundary from a few hundred races.
This is why the bandwidth cannot be chosen by eye, and why CCT pair an MSE-optimal choice with a
bias-corrected interval that stays valid AT that choice. The honest headline is the CCT estimate with
its robust interval; the bandwidth curve is the uncertainty that a single number hides.
3. Is the design valid? The McCrary manipulation test¶
RD's credibility rests on units not being able to precisely control the running variable near the cutoff. If, say, candidates could engineer a win by a hair, the just-winners would differ systematically from just-losers (more resources, better organization), and the "as good as random" logic collapses. McCrary (2008) proposed the test: manipulation would pile units up on the favorable side, creating a discontinuity in the density of the running variable at the cutoff. So we test whether that density jumps.
For elections this should pass — no campaign can guarantee winning by exactly 0.1% — and it does: the density of the vote margin is smooth through $x=0$, and the Cattaneo-Jansson-Ma density test (the modern, tuning-free version of McCrary, in rddensity) fails to reject continuity. A failed test (a density jump) would be a red flag that the design is confounded by sorting.
from rddensity import rddensity
rd=rddensity(X=x,c=0); t_jk=float(rd.test["t_jk"]); p_jk=float(rd.test["p_jk"])
print(f"Cattaneo-Jansson-Ma manipulation test: T = {t_jk:.3f}, p = {p_jk:.3f}")
print(f" p > 0.05 -> fail to reject a continuous density -> NO evidence of manipulation -> RD design is valid.")
fig,ax=plt.subplots(figsize=(9,4.4))
bb=np.linspace(-0.5,0.5,51)
ax.hist(x[(x<0)&(x>=-0.5)],bins=bb[bb<=0],color=BLUE,alpha=.6,density=True,label="losers (x<0)")
ax.hist(x[(x>=0)&(x<=0.5)],bins=bb[bb>=0],color=RED,alpha=.6,density=True,label="winners (x>0)")
ax.axvline(0,color="k",ls="--"); ax.set_xlabel("Democratic margin of victory"); ax.set_ylabel("density")
ax.set_title(f"Running-variable density is smooth at the cutoff (manipulation test p = {p_jk:.2f})"); ax.legend()
plt.tight_layout(); plt.show()
print("No bunching just above zero: candidates cannot precisely manipulate a razor-thin win, so just-winners and just-losers")
print("are comparable. This continuity of the density is the testable implication that supports (never proves) the RD design.")
Cattaneo-Jansson-Ma manipulation test: T = 1.435, p = 0.151 p > 0.05 -> fail to reject a continuous density -> NO evidence of manipulation -> RD design is valid.
No bunching just above zero: candidates cannot precisely manipulate a razor-thin win, so just-winners and just-losers are comparable. This continuity of the density is the testable implication that supports (never proves) the RD design.
4. Fuzzy RD — when the cutoff shifts probability, RD becomes IV¶
In many designs crossing the cutoff does not switch treatment on with certainty — it only raises the probability of treatment (a scholarship offer that not everyone accepts; eligibility that not everyone takes up). This is fuzzy RD. Now the jump in the outcome at the cutoff (the intention-to-treat) understates the treatment effect, because only some units actually changed treatment. The fix is exactly the instrumental-variables logic of the previous notebook: use "above the cutoff" as an instrument for treatment. The fuzzy RD estimand is the ratio $$\tau_{FRD}=\frac{\text{jump in outcome at }c}{\text{jump in treatment probability at }c},$$ a LATE for the compliers who take treatment because they crossed the threshold — estimated by 2SLS in a bandwidth around the cutoff. We demonstrate on a simulation with a known effect of 3: treatment probability jumps from 0.25 to 0.85 at the cutoff, so the sharp (ITT) jump in the outcome is biased toward zero, while the fuzzy ratio recovers the truth.
rng=np.random.default_rng(2); m=8000; xf=rng.uniform(-1,1,m); TRUE=3.0
pr=np.where(xf>=0,0.85,0.25) # treatment prob jumps 0.25 -> 0.85 at cutoff (imperfect compliance)
Df=(rng.uniform(size=m)<pr).astype(int)
Yf=2+TRUE*Df+1.0*xf+rng.normal(0,0.7,m) # true treatment effect = 3
h=0.4; L=(xf<0)&(xf>=-h); Rr=(xf>=0)&(xf<=h)
def jump(v): # local-linear jump at 0 for outcome v
def side(mask,s):
xx=xf[mask];vv=v[mask];w=1-np.abs(xx)/h;X=np.column_stack([np.ones(len(xx)),xx]);W=np.diag(w)
return np.linalg.solve(X.T@W@X,X.T@W@vv)[0]
return side(Rr,1)-side(L,0)
jy=jump(Yf); jd=jump(Df); frd=jy/jd
# 2SLS in the bandwidth using above-cutoff (Z) as instrument for D
mask=L|Rr; Z=(xf[mask]>=0).astype(int); Dd=Df[mask]; Yy=Yf[mask]; xc=xf[mask]
Zm=np.column_stack([np.ones(mask.sum()),Z,xc]); Xm=np.column_stack([np.ones(mask.sum()),Dd,xc])
PZ=Zm@np.linalg.solve(Zm.T@Zm,Zm.T@Xm); b2=np.linalg.solve(PZ.T@PZ,PZ.T@Yy); tsls=b2[1]
print(f"True treatment effect = {TRUE}")
print(f"Sharp (ITT) jump in outcome = {jy:.3f} <- biased toward 0 (only some units switched)")
print(f"Jump in treatment probability = {jd:.3f} (design: 0.25 -> 0.85, so ~0.60)")
print(f"Fuzzy RD = outcome jump / treat jump = {frd:.3f}")
print(f"Fuzzy RD via 2SLS (above-cutoff IV) = {tsls:.3f} <- recovers the truth")
fig,ax=plt.subplots(1,2,figsize=(13,4.2))
for mask_,c,lab in [(xf<0,BLUE,"below cutoff"),(xf>=0,RED,"above cutoff")]:
bb=np.linspace(-1,1,31); mm=(bb[:-1]+bb[1:])/2; ii=np.digitize(xf,bb)-1
ax[0].scatter([mm[b] for b in range(len(mm)) if (ii==b).any()],[Df[ii==b].mean() for b in range(len(mm)) if (ii==b).any()],s=18,color=GREY)
ax[0].axvline(0,color="k",ls="--"); ax[0].set_title(f"First stage: treatment probability jumps {jd:.2f} at cutoff"); ax[0].set_xlabel("running variable"); ax[0].set_ylabel("P(treated)")
ax[1].bar(["sharp/ITT\n(biased)","fuzzy ratio","fuzzy 2SLS"],[jy,frd,tsls],color=[GREY,ORANGE,GREEN]); ax[1].axhline(TRUE,color=RED,ls="--",lw=2,label=f"true effect {TRUE}")
for i,v in enumerate([jy,frd,tsls]): ax[1].text(i,v+0.05,f"{v:.2f}",ha="center")
ax[1].set_title("Fuzzy RD recovers the effect; the sharp jump does not"); ax[1].legend()
plt.tight_layout(); plt.show()
print("\nFuzzy RD = sharp RD applied to BOTH outcome and treatment, then divided -- identical to 2SLS with the threshold as")
print("instrument. Regression discontinuity and instrumental variables are the same identification idea at a boundary.")
True treatment effect = 3.0 Sharp (ITT) jump in outcome = 1.774 <- biased toward 0 (only some units switched) Jump in treatment probability = 0.594 (design: 0.25 -> 0.85, so ~0.60) Fuzzy RD = outcome jump / treat jump = 2.986 Fuzzy RD via 2SLS (above-cutoff IV) = 3.030 <- recovers the truth
Fuzzy RD = sharp RD applied to BOTH outcome and treatment, then divided -- identical to 2SLS with the threshold as instrument. Regression discontinuity and instrumental variables are the same identification idea at a boundary.
5. Summary¶
Regression discontinuity buys causal identification with an unusually mild assumption: not unconfoundedness, not an excluded instrument, only that potential outcomes vary continuously through the cutoff, so any jump is the treatment's doing. On Lee's (2008) House elections, the next-election Democratic vote share jumps by about 6–7 percentage points at the winning threshold — a large, robust incumbency advantage — identified purely by comparing barely-winners to barely-losers.
The workflow the notebook establishes:
- estimate locally — local linear regression on each side within a bandwidth, not a global high-order polynomial (Gelman-Imbens);
- choose the bandwidth honestly — the CCT MSE-optimal bandwidth with bias-corrected robust confidence intervals. Our from-scratch estimate matched
rdrobustexactly at that bandwidth. It was not, however, stable across bandwidths: over $h \in [0.04, 0.40]$ it ranged from 0.058 to 0.084, a spread of 41% of the estimate, non-monotone. The bandwidth is a real researcher degree of freedom and the curve should be reported; - test the design — the McCrary / Cattaneo-Jansson-Ma density test found no manipulation (p ≈ 0.15), supporting validity; a density jump would have condemned it;
- fuzzy RD is IV — when the cutoff shifts treatment probability, the effect is the outcome-jump over the treatment-jump, a LATE estimated by 2SLS with "above the cutoff" as the instrument.
Cross-links. RD's local-randomization interpretation makes it the observational design closest to the randomized experiment of the first notebook; fuzzy RD is literally the 2SLS of the Instrumental Variables notebook applied at a threshold, and inherits its LATE interpretation. The local-linear/kernel machinery connects to the nonparametric-regression tools of the ML and BNP arcs. Next: Panel Data & Fixed Effects, where repeated observations on the same units let us difference away time-invariant confounders — the workhorse that difference-in-differences generalizes.