Thresholds everywhere
Governments love cutoffs: an income ceiling for welfare, an age limit for a pension, a test score for a scholarship, a pollution level that triggers inspection. Units just above and just below the cutoff are nearly identical, except that one group is treated. Comparing them near the threshold is like a small randomised experiment, which is why the regression discontinuity design (RDD) is considered among the most credible quasi-experimental methods.
In a sharp RDD, treatment switches on exactly at the cutoff c of a running variable X: D = 1 if X ≥ c. The effect at the cutoff is the jump in E[Y | X] at c.
Continuity assumption: without the treatment, average outcomes would change smoothly through the cutoff. It fails if units can precisely manipulate their score (firms fudging their distress score to qualify), or if another policy switches on at the same threshold. Lesson 10.3 shows how to probe both.
!pip install -q numpy pandas statsmodels matplotlib seaborn
import numpy as np, pandas as pd
import statsmodels.formula.api as smf
import matplotlib.pyplot as plt, seaborn as sns
rng = np.random.default_rng(2026)
sns.set_theme(style="whitegrid")
A municipal grant at a distress score of 50
A city gives a grant to small firms whose "distress score" is 50 or more. Does the grant raise revenue one year later? Distressed firms naturally earn less, so a naive comparison of grant recipients with non-recipients is badly biased. The true effect at the cutoff is +$12,000.
N, CUTOFF = 3_000, 50
score = rng.uniform(20, 80, N)
grant = (score >= CUTOFF).astype(int)
revenue = 100_000 - 500 * score + 4 * (score - 50) ** 2 + 12_000 * grant + rng.normal(0, 6_000, N)
df = pd.DataFrame({"score": score, "grant": grant, "revenue": revenue})
df["x"] = df.score - CUTOFF # centre the running variable at the cutoff
naive = df.revenue[df.grant == 1].mean() - df.revenue[df.grant == 0].mean()
print(f"Naive difference: {naive:,.0f} (true effect at the cutoff: 12,000)")
Local linear estimation
Fit a separate line on each side of the cutoff, using only observations within a bandwidth h of it.
With the running variable centred, the coefficient on grant in
revenue ~ grant + x + grant:x
is exactly the gap between the two lines at x = 0. A triangular kernel gives more weight to observations closer to the cutoff.
def local_linear(data, h, kernel="triangular"):
sub = data[data.x.abs() <= h].copy()
w = (1 - sub.x.abs() / h) if kernel == "triangular" else np.ones(len(sub))
fit = smf.wls("revenue ~ grant * x", data=sub, weights=w).fit(cov_type="HC1")
return fit.params["grant"], fit.bse["grant"], len(sub)
for h in [5, 10, 15, 30]:
est, se, n = local_linear(df, h)
print(f"h = {h:>2}: effect = {est:>8,.0f} (SE {se:,.0f}), n = {n:,}")
Narrow bandwidths have little bias but few observations (noisy); wide ones are precise but let curvature leak into the estimate. Lesson 10.2 picks h in a principled, data-driven way.
The RD plot
The standard RD figure: average the outcome within bins of the running variable, plot the bin means, and overlay the fitted lines from each side. The reader should see the jump before trusting any number.
h = 15
est, se, _ = local_linear(df, h)
bins = pd.cut(df.score, np.arange(20, 81, 2))
binned = df.groupby(bins, observed=True).agg(score=("score", "mean"), revenue=("revenue", "mean"))
plt.figure(figsize=(10, 5))
plt.scatter(binned.score, binned.revenue, color="#334155", zorder=3, label="Bin means")
for side in (0, 1):
sub = df[(df.grant == side) & (df.x.abs() <= h)]
line = smf.ols("revenue ~ x", data=sub).fit()
xs = np.linspace(sub.x.min(), sub.x.max(), 50)
plt.plot(xs + CUTOFF, line.params["Intercept"] + line.params["x"] * xs, lw=3,
color="#f97316" if side else "#2563eb", label="Grant" if side else "No grant")
plt.axvline(CUTOFF, ls="--", color="red")
plt.title(f"Sharp RDD: estimated jump = ${est:,.0f} (SE {se:,.0f})")
plt.xlabel("Distress score"); plt.ylabel("Revenue one year later ($)"); plt.legend(); plt.show()
An RD estimate is local: it is the effect for firms with a distress score near 50. It may differ for very distressed firms. That is the price of credibility, and it should be stated plainly in any memo.
Exercise: the danger of global polynomials
A common older practice fits a high-order polynomial to the whole sample instead of local lines. Estimate
the jump with global polynomials of order 1 to 5 (interacted with grant) using all observations. How
stable are the estimates? (Gelman and Imbens, 2019, argue against orders above 2.)
Show a solution
for order in range(1, 6):
terms = " + ".join([f"I(x**{k})" for k in range(1, order + 1)])
fit = smf.ols(f"revenue ~ grant * ({terms})", data=df).fit(cov_type="HC1")
print(f"order {order}: {fit.params['grant']:>8,.0f} (SE {fit.bse['grant']:,.0f})")
Checklist and recap
- I can explain why units near a cutoff form a local experiment, and state continuity.
- I can estimate a sharp RDD with local linear regression and a kernel.
- I always show a binned RD plot alongside the estimate.