Regression discontinuity in Python: rdrobust with robust bias correction

Many programmes are assigned by a rule on a score. Students who score at least 60 on an exam get a scholarship, candidates with more than half the vote win the seat, firms below a size threshold are exempt from a regulation. Units just above and just below the threshold are alike in almost every respect except treatment. A regression discontinuity (RD) design compares them. It estimates the effect from the jump in the outcome at the cutoff.

What regression discontinuity estimates

Let x be the running variable (the score) and c the cutoff, with treatment for x ≥ c. In a sharp design treatment switches on exactly at the cutoff, and the estimand is the jump in the conditional mean of the outcome:

τ = lim_{x↓c} E[y | x] − lim_{x↑c} E[y | x]

It is the average treatment effect for units at the cutoff. The identifying assumption is that the conditional means of the potential outcomes are continuous at c: nothing else changes discontinuously there, and units cannot sort precisely onto one side of the threshold. A heap of observations just above the cutoff is a warning sign; McCrary (2008) proposes testing for a discontinuity in the density of the running variable.

In a fuzzy design crossing the cutoff only changes the probability of treatment. The ratio of the jump in the outcome to the jump in treatment take-up is then a local average treatment effect for compliers at the cutoff, an instrumental-variables estimate with the cutoff as instrument.

The standard estimator fits a local linear regression on each side of the cutoff, using only observations within a bandwidth h and weighting those closest to c most (Imbens and Lemieux 2008). The bandwidth trades bias against variance. A wider window brings in more data but lets curvature away from the cutoff bias the jump. When h is chosen to minimize the mean squared error, conventional confidence intervals are invalid, because the bias is of the same order as the standard error. The robust bias-corrected interval of Calonico, Cattaneo and Titiunik (2014) subtracts an estimate of the bias, obtained from a higher-order fit with a pilot bandwidth b, and widens the standard error to account for the extra estimation. Calonico, Cattaneo and Farrell (2020) derive a smaller bandwidth that minimizes the coverage error of that interval.

In OpenEconometrics

oe.rdrobust implements this approach. The call shape is oe.rdrobust(data=df, y="gpa", running="score", cutoff=60); observations with running >= cutoff are on the treated side. It corresponds to the community Stata command rdrobust, here rdrobust gpa score, c(60), and stores the selected bandwidths, which Stata's rdbwselect reports, in metrics.

• The result has three rows. Conventional is the local polynomial estimate with its usual standard error. Bias-corrected subtracts the estimated bias but keeps the conventional standard error. Robust has the bias-corrected estimate with the robust standard error; its confidence interval is the one to report.
• p (default 1, local linear) and q (default p + 1) set the polynomial orders; kernel is "triangular" (default), "epanechnikov" or "uniform".
• bwselect="mserd" (default) picks one MSE-optimal bandwidth for both sides; "msetwo" picks one per side and "cerrd" the coverage-error-optimal one. h= and b= fix the bandwidths, as a number or [left, right].
• vce="nn" (default, nnmatch=3 neighbours) or "hc0" to "hc3" chooses the residuals of the variance; covariance="cluster" with cluster= clusters it.
• fuzzy="received" names the treatment received for a fuzzy design; extra["first_stage"] and extra["reduced_form"] hold the two jumps. covariates= adds a common-slope adjustment, deriv=1 estimates a kink, and weights= accepts analytic weights.

oe.rdplot returns the binned means and global polynomial fits of an RD plot. oe.rddensity tests for manipulation of the running variable at bandwidths you supply. For the local-randomization approach, oe.rdrandinf and oe.rdwinselect are covered in the local randomization guide. The treatment-effects family guide gives every rdrobust formula and bandwidth rule.

Supported options

Options of oe.rdrobust

oe.rdplot

Stata command: rdplot.

oe.rddensity

Stata command: rddensity.

Generated from oe.capabilities() of openecon 0.3.20a1. Stata names identify the corresponding workflow, not numerical parity.

Example

The script simulates 3,000 students with exam scores drawn from a normal distribution with mean 55 and standard deviation 10. Those scoring 60 or more receive a scholarship. First-year GPA is a smooth quadratic function of the score, plus a jump of 0.25 for scholarship recipients and normal noise with standard deviation 0.5. The true effect at the cutoff is 0.25. The script fits the default MSE-optimal design and then repeats the fit with the coverage-error-optimal bandwidth.

import numpy as np
import pandas as pd
import openecon as oe


def simulate(n=3000, seed=8):
    """Students scoring 60 or more on an entrance exam get a scholarship. True jump in GPA = 0.25."""
    rng = np.random.default_rng(seed)
    score = rng.normal(55, 10, size=n)
    treated = (score >= 60).astype(float)
    s = score - 60
    gpa = 2.0 + 0.03 * s - 0.0008 * s**2 + 0.25 * treated + rng.normal(scale=0.5, size=n)
    return pd.DataFrame({"score": score, "gpa": gpa})


if __name__ == "__main__":
    df = simulate()
    rd = oe.rdrobust(data=df, y="gpa", running="score", cutoff=60)
    print(rd.summary())
    cer = oe.rdrobust(data=df, y="gpa", running="score", cutoff=60, bwselect="cerrd")
    robust = cer.coefficients[2]
    print(f"\nCER-optimal bandwidth h = {cer.metrics['h_left']:.3f}: robust estimate {robust.estimate:.4f}, "
          f"95% CI [{robust.ci_low:.4f}, {robust.ci_high:.4f}]")

With openecon 0.3.20a1 from PyPI the script prints:

Regression discontinuity: sharp local polynomial — gpa
Observations: 3000  |  Covariance: robust  |  Confidence: 95%

Term            Estimate  Std. error        z     P>|stat|   CI lower  CI upper
--------------  --------  ----------  -------  -----------  ---------  --------
Conventional    0.224784    0.067084  3.35078  0.000805846  0.0933014  0.356266
Bias-corrected   0.20686    0.067084   3.0836   0.00204511  0.0753781  0.338342
Robust           0.20686   0.0789103  2.62146   0.00875541  0.0521989  0.361522

h_left: 5.36683  |  h_right: 5.36683  |  b_left: 8.97242  |  b_right: 8.97242
n_left: 2080  |  n_right: 920  |  n_h_left: 595  |  n_h_right: 487
n_b_left: 1023  |  n_b_right: 662  |  p: 1  |  q: 2
cutoff: 60

CER-optimal bandwidth h = 3.596: robust estimate 0.1947, 95% CI [0.0199, 0.3695]

Reading the output

• The estimate is local. The MSE-optimal bandwidth is 5.37 points on each side, so only 595 students below and 487 above the cutoff carry weight, out of 3,000. The bias bandwidth b is wider, 8.97 points.
• Report the robust row. The conventional estimate is 0.2248; the bias-corrected one is 0.2069. The robust standard error, 0.0789, is about 18% larger than the conventional 0.0671 because it includes the uncertainty of the bias estimate. The robust 95% interval runs from 0.052 to 0.362 and contains the true 0.25.
• The bias-corrected row mixes conventions. It pairs the corrected estimate with the conventional standard error and is shown for comparison only; its interval ignores the variability of the bias estimate.
• The CER bandwidth is narrower. At h = 3.596 the robust estimate is 0.1947 with an interval from 0.020 to 0.370. Fewer observations widen the interval; similar estimates across reasonable bandwidths are reassuring.
• The design is sharp because treatment is a deterministic function of the score. With imperfect take-up, pass fuzzy= and read the first-stage jump in extra["first_stage"] before the effect.

Limitations

• The common oe.predict and oe.margins adapters do not cover oe.rdrobust; use its own result tables and family helpers.
• Cluster-robust covariance for oe.rdrobust allows at most 1 clustering dimension.
• The three rows are alternative estimates of one parameter; their covariance matrix is diagonal, and the covariances between rows are not reported.
• Covariate adjustment, derivative (kink) estimation and weights require an in-memory data frame; the oe.Dataset replay rejects them.
• oe.rddensity needs explicit bandwidths: automatic density bandwidth selection, restricted fits and the plug-in density variance are outside its validated domain.
• With bwrestrict=True (the default) every bandwidth is capped at the distance from the cutoff to the farthest observation on that side.
• No numerical parity with Stata is claimed: the inventory's stata_parity_validated flag is False.

References

• Calonico, S., Cattaneo, M. D. and Farrell, M. H. (2020). Optimal bandwidth choice for robust bias-corrected inference in regression discontinuity designs. Econometrics Journal, 23(2), 192–210. doi:10.1093/ectj/utz022
• Calonico, S., Cattaneo, M. D. and Titiunik, R. (2014). Robust nonparametric confidence intervals for regression-discontinuity designs. Econometrica, 82(6), 2295–2326. doi:10.3982/ECTA11757
• Imbens, G. W. and Lemieux, T. (2008). Regression discontinuity designs: A guide to practice. Journal of Econometrics, 142(2), 615–635. doi:10.1016/j.jeconom.2007.05.001
• McCrary, J. (2008). Manipulation of the running variable in the regression discontinuity design: A density test. Journal of Econometrics, 142(2), 698–714. doi:10.1016/j.jeconom.2007.05.005