Regression discontinuity in Python: getting the effect at the cutoff right

A regression-discontinuity design estimates the outcome discontinuity at a treatment threshold. Units just above a cutoff are treated; units just below are not; under the design's identifying assumptions, the discontinuity identifies the local treatment effect. In Python, estimating that discontinuity takes only a few lines of regression code. That ease creates a risk: treating the printed coefficient as sufficient evidence without evaluating the specification and identification assumptions.

To see why, we can plant a known effect in simulated data and ask different estimators to recover it. The true jump at the cutoff is set to 0.75. The outcome varies smoothly with the running variable in a way that is not exactly polynomial, and the treatment effect increases away from the cutoff. That last feature matters: the effect at the cutoff (0.75) is not the same as the average effect among the treated (about 1.5). RDD is designed to recover the local effect at the cutoff, not the average effect elsewhere.

Problem 1: Functional form drives the estimate

Start with a tempting shortcut: regress the outcome on the running variable, its square, and a treatment indicator across the full sample, then read the treatment coefficient as the effect. This deliberately overconstrained specification forces the slope and curvature to be identical on both sides of the cutoff. It returns 1.82, against a true value of 0.75.

import numpy as np, statsmodels.api as sm

rng = np.random.default_rng(20260704)
X = rng.uniform(-1, 1, 2000)              # running variable
D = (X >= 0).astype(float)                # treated above the cutoff
Y = np.sin(3*X) + 0.5*np.exp(X) \
    + (0.75 + 1.5*X)*D + rng.normal(0, 0.7, 2000)   # true jump at cutoff = 0.75

shortcut = sm.add_constant(np.column_stack([D, X, X**2]))
print(sm.OLS(Y, shortcut).fit().params[1])  # -> 1.82, not 0.75

# If a global quadratic is used as a diagnostic, allow each side its own
# slope and curvature. At X=0, the coefficient on D is the fitted jump.
side_specific = sm.add_constant(np.column_stack([
    D, X, D*X, X**2, D*(X**2)
]))
global_quadratic = sm.OLS(Y, side_specific).fit(cov_type="HC3")

The shortcut runs without warning, but it is not a correctly specified two-sided RD curve. By omitting interactions between treatment status and the running-variable terms, it makes one curve serve both sides and can charge misspecified shape to the treatment indicator. The side-specific diagnostic above removes that particular restriction, although a high-order global fit remains a poor default for boundary estimation.

Increasing the shortcut's polynomial order exposes rather than solves the problem. The following values record a sensitivity exercise under the same common-coefficient restriction. They are not three equally recommended RD estimators.

Overconstrained shortcut order	Printed treatment coefficient (true jump 0.75)
Quadratic	1.82
Cubic	0.81
Fifth-order	0.74

None of these values validates the shortcut. High-order global polynomials are especially unstable at a boundary and should not be selected because one estimate happens to be near the planted truth. The table shows how the common-coefficient restriction changes the printed coefficient across polynomial orders; it does not establish that the three specifications are equally acceptable.

The fix: estimate locally

A local approach avoids guessing the global shape. Instead of fitting one curve to the entire sample, we fit simple models on each side of the cutoff using only nearby observations, and weight points closer to the cutoff more heavily. The reason is that any smooth function is approximately linear in a small neighborhood, so a local linear fit works without committing to a global polynomial.

def local_linear(X, D, Y, h):
    m = np.abs(X) <= h                    # keep points within bandwidth h of the cutoff
    w = 1 - np.abs(X[m]) / h              # triangular weights, heavier near the cutoff
    Z = sm.add_constant(np.column_stack([D[m], X[m], D[m]*X[m]]))
    return sm.WLS(Y[m], Z, weights=w).fit(cov_type="HC3")

fit = local_linear(X, D, Y, h=0.20)
print(fit.params[1])                      # -> 0.87 on this draw

On the reference draw we get 0.87, with a 95% confidence interval that covers the true value of 0.75. Across 200 draws from this data-generating process, the mean local estimate is 0.75. The repeated draws measure sampling behavior for this specified simulation; they do not validate every bandwidth rule or outcome function.

What still requires judgment

Local estimation does not eliminate all choices. First, the bandwidth matters. A narrow window uses observations very close to the cutoff and targets the local effect, but with higher variance. A wider window increases precision but pulls in observations where the treatment effect is larger, shifting the estimate upward toward the treated average of 1.5. As the bandwidth approaches the full sample, the local estimator becomes global again.

Second, recovery is approximate, not exact. Even in clean simulations, estimates vary across samples and bandwidths, and that variation is part of the result, not a nuisance to ignore. In practice, bandwidth selection should follow established data-driven methods, and results should be reported with sensitivity to reasonable alternatives.

Problem 2: The cutoff may not isolate treatment

Even a well-executed local design can answer the wrong question. RDD identifies a treatment effect only if crossing the cutoff changes one thing: the treatment of interest. If multiple programs or policies activate at the same threshold, the observed jump combines all of them.

Consider a concrete case: an income cutoff that determines eligibility for both a health program and a housing subsidy. An RDD at that threshold estimates the combined effect of both interventions, and no choice of bandwidth or functional form can separate them within the design. A density-continuity diagnostic can test for evidence of sorting in the running variable, but failure to reject a discontinuity does not confirm the absence of sorting and says nothing about whether another treatment changes at the threshold. In this simulation, the running variable is uniform by construction. The package reports a simpler descriptive screen, not a McCrary density estimate: within 0.1 of the cutoff it counts 102 observations below and 110 above, producing a binomial-balance z statistic of 0.55. When a second 0.50 jump is added at the same threshold, the local fit estimates about 1.18 rather than the target program's 0.75 effect. The estimator runs as specified, but its discontinuity combines both interventions.

Reference compound-treatment simulation and descriptive window-balance screen

Quantity	Value
Target program effect at cutoff	0.75
Second program jump at cutoff	0.50
Local estimate of combined discontinuity	1.18
Observations in ±0.1 window (below / above)	102 / 110
Binomial-balance z statistic	0.55

This is the more serious threat. Functional-form issues bias the estimate; identification failures redefine it.

What to prioritize, before the estimate

Three checks matter more than the coefficient itself, in order.

1. Does the cutoff isolate the treatment? If several interventions change at the threshold, the estimate bundles them. This has to be argued from institutional knowledge, not inferred from the data alone. Clarify whether the design is sharp or fuzzy, and what estimand follows.
2. Is the estimate robust to defensible choices? Vary bandwidths, functional forms, and specifications. If the estimate moves across them, report the range. Stability is evidence; a single preferred number is not.
3. What exactly is being estimated? RDD recovers a local effect at the cutoff. When treatment effects vary with the running variable, that differs from the average effect among participants.

The value the code returns is the output, read alongside its uncertainty and its identification argument; the finding is whether the design and the estimate hold up.

Reproduction

The full simulation framework, including the data generation, global and local estimators, bandwidth variation, compound-treatment case, frozen JSON output, and assertion-based recovery check, is in the pinned public package. The documented commands run from a clean checkout and regenerate the fixed-seed results. Modifying the outcome shape or treatment effect shows how the estimators behave under a different data-generating process; it does not validate identification in an empirical application.