nb4.3 — RD with `rdrobust`; CCT bandwidths

Companion to Ch 4.3, “Regression Discontinuity.”

The chapter made a single promise: when a treatment switches on at a known cutoff c in a running variable X, the jump in \mathbb{E}[Y\mid X] at c is the causal effect for the units sitting right at the cutoff. This notebook turns that promise into running code. We follow Sam, who trades on index effects and has fixated on the annual reconstitution of the Russell stock indexes: a stock’s end-of-May market-cap rank decides whether it lands in the Russell 1000 (large-cap) or crosses into the Russell 2000 (small-cap, far more passive money tracking it). The boundary near rank 1{,}000 is the cutoff. We recenter rank so the cutoff is c = 0.

We do seven things:

Simulate a sharp-RD world with a true jump \tau = 0.6 percentage points of passive ownership.
Build the canonical RD plot — binned scatter + two local-linear fits + the visible jump — first by hand, then with rdplot.
Estimate the jump with local-linear regression by hand (triangular kernel) across bandwidths, to see the bias–variance trade-off (h=4 \to \approx 0.56, h=32 \to \approx 0.41).
Estimate with rdrobust using the data-driven CCT bandwidth and robust bias-corrected inference; confirm the CI brackets the true 0.6.
Run a McCrary-style density test for manipulation at the cutoff.
Show the global high-order polynomial pathology (Gelman–Imbens) for contrast.
A fuzzy-RD mini-demo: a jump in treatment probability recovered via the IV/Wald ratio, connecting straight back to Week 3.

Then a Your Turn section where you vary the bandwidth and inject manipulation to watch the density test fire.

Setup

We fix one seeded generator so the notebook reproduces bit-for-bit (the CONVENTIONS.md rule), and force matplotlib’s non-interactive Agg backend so the notebook runs headless — laptop, server, or CI — with no display attached.

The imports to notice are the last two. The production RD toolkit lives in the rdrobust package: from rdrobust import rdrobust, rdplot. rdrobust does the local-polynomial estimation, the data-driven bandwidth, and the robust bias-corrected inference in one call; rdplot draws the canonical RD plot. A companion density test (rddensity, the Cattaneo–Jansson–Ma implementation of McCrary’s idea) ships as a separate package that is not installed in this environment, so in §5 we roll a transparent binned-density version by hand — same logic, fewer dependencies.

import matplotlib
matplotlib.use("Agg")  # headless, non-interactive backend

import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt
from scipy import stats
from rdrobust import rdrobust, rdplot  # production RD: CCT bandwidth + robust bias-corrected inference

rng = np.random.default_rng(20260528)  # one seeded generator for the whole notebook

import rdrobust as _rd
print("numpy", np.__version__, "| pandas", pd.__version__, "| statsmodels", sm.__version__)
print("rdrobust package loaded (rdrobust + rdplot available)")

1. A sharp-RD world, built to order

We reproduce the chapter’s §11 simulation exactly. Sam imagines N = 6{,}000 stocks near the Russell 1000/2000 boundary. The running variable X is each stock’s rank distance from the cutoff, recentered so c = 0: stocks with X < 0 are just inside the Russell 1000, stocks with X \ge 0 have crossed into the Russell 2000. In a sharp design, crossing the cutoff flips treatment with certainty, so the treatment indicator is the deterministic step

D_i = \mathbf{1}\{X_i \ge c\}.

The outcome Y is next-year passive-fund ownership (%). It moves smoothly with rank through the cutoff — a gentle upward trend plus a touch of curvature on the treated side — and then the treatment adds a clean vertical jump of \tau = 0.6 on top. The curvature on the treated side is deliberate: it is exactly what will bias a too-wide local-linear fit in §3, and what a global polynomial will read catastrophically in §6.

The empirical spec, stated CONVENTIONS-§4 style: outcome = next-year passive ownership; treatment = Russell-2000 membership D=\mathbf{1}\{X\ge 0\}; running variable = recentered reconstitution rank; identifying assumption = the potential-outcome mean \mathbb{E}[Y(0)\mid X] is continuous at c=0, so any jump at c is the treatment and nothing else.

N = 6_000
tau_true = 0.6                          # the true jump at the cutoff (passive-ownership pct points)
c = 0.0                                  # cutoff (rank recentered so the boundary is 0)

# Running variable X: rank distance from the Russell 1000/2000 cutoff. X >= 0 -> Russell 2000.
X = rng.uniform(-50, 50, size=N)
D = (X >= c).astype(float)              # SHARP RD: treatment is a deterministic step at c

# Smooth, continuous-through-c mean: gentle slope + a little curvature on the treated side only.
mu = 5.0 + 0.02 * X + 0.0020 * np.where(X > 0, X**2, 0.0)
Y = mu + tau_true * D + rng.normal(0.0, 1.0, size=N)

df = pd.DataFrame({"Y": Y, "X": X, "D": D})
print(df.describe().round(3).loc[["mean", "min", "max"]])
print(f"\nshare treated (X >= 0): {df['D'].mean():.3f}  (sharp design -> exactly the share with X>=0)")

2. The RD plot: see the jump before you trust a regression

RD is the most visual design in the book, and a good RD plot is half the argument. Before any regression we plot the data, because the human eye is an excellent jump-detector. The canonical plot layers three things (chapter §7):

Binned scatter. Chop X into bins, average Y within each bin, plot one dot per bin. Raw scatters are too noisy; binning averages away idiosyncratic noise and reveals the shape of \mathbb{E}[Y\mid X].
Two fitted curves. A separate local-linear fit on each side of c, fit only within a window around the cutoff and extrapolated to the cutoff but never across it — so the reader sees the window the estimate rests on. These are the lines whose boundary heights we difference.
The jump. The vertical gap between the two fitted curves at c is \hat\tau — and it should be visible to the naked eye.

We draw it by hand first so nothing is hidden, then let rdplot do it the production way.

def binned_means(x, y, c=0.0, binwidth=2.0, lo=-50, hi=50):
    '''Mean of y within equal-width bins of x; returns bin centers and bin means.'''
    edges = np.arange(lo, hi + binwidth, binwidth)
    idx = np.digitize(x, edges) - 1
    idx = np.clip(idx, 0, len(edges) - 2)
    centers = (edges[:-1] + edges[1:]) / 2
    means = np.array([y[idx == k].mean() if np.any(idx == k) else np.nan for k in range(len(centers))])
    return centers, means

def side_line(d, cutoff, treated_side, h_plot=12.0):
    '''Local OLS line within a window of half-width h_plot on one side of the cutoff.

    We fit *locally* (not on all data) so each line's height at the boundary is read from
    points near c, exactly as RD intends -- a global line would be bent by far-away curvature.
    '''
    if treated_side:
        sub = d[(d["X"] >= cutoff) & (d["X"] <= cutoff + h_plot)]
    else:
        sub = d[(d["X"] < cutoff) & (d["X"] >= cutoff - h_plot)]
    Xc = (sub["X"] - cutoff).to_numpy()
    M = sm.add_constant(Xc)
    r = sm.OLS(sub["Y"].to_numpy(), M).fit()
    return r

h_plot = 12.0
cen, bm = binned_means(df["X"].to_numpy(), df["Y"].to_numpy(), c=c, binwidth=2.0)
rl = side_line(df, c, treated_side=False, h_plot=h_plot)
rr = side_line(df, c, treated_side=True, h_plot=h_plot)

xl = np.linspace(-h_plot, 0, 50); xr = np.linspace(0, h_plot, 50)
yl = rl.params[0] + rl.params[1] * (xl - c)
yr = rr.params[0] + rr.params[1] * (xr - c)
jump_at_c = rr.params[0] - rl.params[0]   # difference of the two intercepts (boundary heights) at c

fig, ax = plt.subplots(figsize=(8, 5))
ax.scatter(cen[cen < 0], bm[cen < 0], s=22, color="#3b6", label="binned mean (control, R1000)")
ax.scatter(cen[cen >= 0], bm[cen >= 0], s=22, color="#c43", label="binned mean (treated, R2000)")
ax.plot(xl, yl, color="#3b6", lw=2)
ax.plot(xr, yr, color="#c43", lw=2)
ax.axvline(c, color="0.4", ls="--", lw=1)
ax.annotate(f"jump $\\hat\\tau \\approx$ {jump_at_c:.2f}", xy=(0, (rl.params[0]+rr.params[0])/2),
            xytext=(8, rl.params[0]-0.8),
            arrowprops=dict(arrowstyle="->", color="0.2"))
ax.set_xlabel("running variable $X$ = rank distance from cutoff  ($c=0$)")
ax.set_ylabel("next-year passive ownership (%)")
ax.set_title("RD plot (by hand): binned means + separate fits + the jump")
ax.legend(loc="upper left", fontsize=9)
fig.tight_layout()
print(f"hand-drawn boundary jump (local OLS within |X|<={h_plot:.0f} each side) = "
      f"{jump_at_c:.3f}  (true tau = {tau_true})")

The same plot, the production way: `rdplot`

rdplot automates the binning (it picks a sensible number of bins via its binselect rule) and overlays a polynomial fit on each side. Note the API defaults: rdplot(y, x, c=0, p=4, kernel='uni'). We pass hide=True so it returns the plot object without forcing its own show(), then we read off its chosen bin count J. The shape should match our hand-drawn version: smooth trends on each side and a clear step at c=0.

rp = rdplot(y=df["Y"].to_numpy(), x=df["X"].to_numpy(), c=0,
            title="rdplot: binned scatter + local-polynomial fits",
            x_label="X = rank distance from cutoff", y_label="passive ownership (%)",
            hide=True)
print("rdplot chosen number of bins (left, right):", rp.N_h if hasattr(rp, "N_h") else "n/a")
print("rdplot bin-selection method:", rp.binselect)
print("rdplot J (bins each side):", rp.J)

3. Local-linear estimation by hand: the bias–variance trade-off

Now we estimate the jump properly — not with a full-data OLS line, but with local-linear regression inside a window of half-width h (the bandwidth) around the cutoff, weighting points by a triangular kernel that gives weight 1 at c and declines linearly to 0 at the window edge. This is one weighted least squares (WLS) regression of Y on a constant, the treatment dummy D, the recentered running variable (X-c), and their interaction; the coefficient on D is the estimated jump (chapter §5).

The bandwidth is the single most consequential choice in RD, and it is a bias–variance trade-off (chapter §6). Narrow h: few points, so the linear approximation is nearly exact (low bias) but the fit jitters (high variance, wide SE). Wide h: many points (low variance, tight SE) but distant data with curvature bend the line and mis-read its height at the boundary (high bias). With our true jump of 0.6, watch the estimate sit near the truth at h=4 with a wide SE and drift downward (biased) as h grows, even as the SE shrinks.

def local_linear_rd(df, c=0.0, h=8.0):
    '''Triangular-kernel local-linear RD: one WLS with a cutoff dummy and separate slopes.

    The coefficient on D is the estimated jump tau at the cutoff.
    '''
    d = df.copy()
    d["Xc"] = d["X"] - c
    d = d[d["Xc"].abs() <= h].copy()                  # keep only the window |X - c| <= h
    w = (1.0 - (d["Xc"].abs() / h)).clip(lower=0)     # triangular kernel: 1 at c, 0 at the edge
    M = pd.DataFrame({"const": 1.0, "D": d["D"], "Xc": d["Xc"], "Xc_D": d["Xc"] * d["D"]})
    res = sm.WLS(d["Y"], M, weights=w).fit(cov_type="HC1")
    return res

rows = []
for h in (4, 8, 16, 32):
    r = local_linear_rd(df, c=c, h=h)
    rows.append({"h": h, "tau_hat": r.params["D"], "se": r.bse["D"], "n_in_window": int(r.nobs)})
ll_table = pd.DataFrame(rows)
print(ll_table.round(3).to_string(index=False))
print(f"\ntrue tau = {tau_true}.  Note: h=4 -> ~0.56 (near truth, wide SE); "
      f"h=32 -> ~0.41 (biased low, tight SE).")

A picture makes the trade-off unmistakable. We plot \hat\tau against h with \pm 1.96\,\text{SE} error bars. As h grows the error bars shrink (variance falls) but the point estimate slides away from the true 0.6 (bias grows). The “best” bandwidth is the one that balances these — which is exactly what the data-driven CCT rule in §4 computes for us, so we are not eyeballing it (or, worse, bandwidth-shopping until the stars appear).

hs = np.array([2, 4, 6, 8, 12, 16, 24, 32, 40])
taus, ses = [], []
for h in hs:
    r = local_linear_rd(df, c=c, h=float(h))
    taus.append(r.params["D"]); ses.append(r.bse["D"])
taus = np.array(taus); ses = np.array(ses)

fig, ax = plt.subplots(figsize=(8, 5))
ax.errorbar(hs, taus, yerr=1.96 * ses, fmt="o-", capsize=3, color="#36c", label=r"$\hat\tau \pm 1.96\,$SE")
ax.axhline(tau_true, color="#c43", ls="--", lw=1.5, label=f"true $\\tau = {tau_true}$")
ax.set_xlabel("bandwidth $h$ (narrow $\\to$ wide)")
ax.set_ylabel(r"local-linear $\hat\tau$")
ax.set_title("Bias–variance trade-off: SEs shrink but bias grows as $h$ widens")
ax.legend()
fig.tight_layout()
print("As h grows: SE falls", f"({ses[0]:.3f} -> {ses[-1]:.3f})",
      "but |bias| grows", f"(|{taus[0]-tau_true:.3f}| -> |{taus[-1]-tau_true:.3f}|).")

4. The production estimate: `rdrobust` with the CCT bandwidth

Doing this by hand taught us the moving parts; in practice you call rdrobust, which implements the modern standard:

Calonico–Cattaneo–Titiunik (2014) data-driven bandwidth (bwselect='mserd', the MSE-optimal choice — the default).
Local-linear fit (p=1) with a triangular kernel (kernel='tri', the default).
Robust bias-corrected inference: the MSE-optimal bandwidth is deliberately wide enough that a leftover bias remains, so the naive interval \hat\tau \pm 1.96\,\widehat{\text{se}} undercovers. CCT estimate the bias, subtract it, and inflate the SE for the extra noise the correction introduces. The result is the Robust row, which achieves nominal 95% coverage even at the optimal bandwidth.

rdrobust returns three rows — Conventional, Bias-Corrected, and Robust — as labeled pandas objects (.coef, .ci, .se, .bws, .pv). Report the Robust row. The chosen bandwidth lives in .bws (h is the estimation bandwidth, b the bias-correction bandwidth).

API note (numpy arrays). We pass .to_numpy() arrays rather than pandas Series. The current rdrobust mishandles a Series carrying a non-default index in some code paths (notably fuzzy=), so numpy arrays are the safe, portable choice.

out = rdrobust(y=df["Y"].to_numpy(), x=df["X"].to_numpy(), c=0)  # CCT bandwidth + robust inference (defaults)

h_cct = out.bws.loc["h", "left"]
tau_robust = float(out.coef.loc["Robust"].iloc[0])
ci_robust = out.ci.loc["Robust"]
print(f"CCT estimation bandwidth h = {h_cct:.2f}  (bias-correction bandwidth b = {out.bws.loc['b','left']:.2f})")
print(f"effective sample in window (left, right): {out.N_h}")
print()
print("rdrobust coefficient table:")
print(out.coef.round(4).to_string())
print()
print("rdrobust 95% CIs:")
print(out.ci.round(4).to_string())
print()
brackets = (ci_robust.iloc[0] <= tau_true <= ci_robust.iloc[1])
print(f"Robust point estimate = {tau_robust:.3f}, 95% CI = "
      f"[{ci_robust.iloc[0]:.3f}, {ci_robust.iloc[1]:.3f}]")
print(f"Does the robust CI bracket the true tau = {tau_true}?  {brackets}")
assert brackets, "robust CI failed to bracket the true jump"
print(f"Compare: hand-rolled local-linear at h~{h_cct:.0f} gave "
      f"{local_linear_rd(df, c=c, h=h_cct).params['D']:.3f} — same neighborhood, as expected.")

5. The McCrary density test: did anyone game the cutoff?

The local-randomization story has one fatal enemy: manipulation. If units can precisely control which side of c they land on, the just-above and just-below groups stop being exchangeable. The fingerprint is a discontinuity in the density of the running variable at c — too few observations just below, too many just above (chapter §8).

The McCrary (2008) test estimates the density of X separately on each side of c and tests whether the two agree at the boundary. The null is continuity (no sorting) — what you hope to see; rejecting it is evidence of manipulation. The modern Cattaneo–Jansson–Ma implementation ships in the rddensity package (a separate install — from rddensity import rddensity), which is not available in this environment. So we build a transparent, dependency-free version with the same logic: bin X, normalize to a density, fit a line to the binned density on each side, and test whether the two intercepts at c differ.

For the Russell design this should pass — a firm cannot finely control its market-cap rank to the dollar on the ranking day, which is much of why the cutoff is clean. We verify that, then in Your Turn we inject manipulation and watch the test fire.

def mccrary_density_test(x, c=0.0, binwidth=2.0, h=20.0, lo=-50, hi=50):
    '''Binned local-linear density test (McCrary-style): compare one-sided density intercepts at c.

    Returns the left/right density at c, their difference, a z-stat and a two-sided p-value.
    Null: density is continuous at c (no manipulation).
    '''
    edges = np.arange(lo, hi + binwidth, binwidth)
    counts, _ = np.histogram(x, bins=edges)
    centers = (edges[:-1] + edges[1:]) / 2
    dens = counts / (len(x) * binwidth)               # normalized so it integrates to ~1
    m = np.abs(centers - c) <= h
    cc, dd = centers[m], dens[m]

    def side(mask):                                    # one-sided line; intercept = density at c
        xs = cc[mask] - c
        r = sm.OLS(dd[mask], sm.add_constant(xs)).fit(cov_type="HC1")
        return r.params[0], r.bse[0]

    f_left, se_left = side(cc < c)
    f_right, se_right = side(cc >= c)
    diff = f_right - f_left
    se = np.sqrt(se_left**2 + se_right**2)
    z = diff / se
    pv = 2 * (1 - stats.norm.cdf(abs(z)))
    return dict(dens_left=f_left, dens_right=f_right, diff=diff, se=se, z=z, pv=pv), centers, dens

res, cen_d, dens_d = mccrary_density_test(df["X"].to_numpy(), c=c)
print("McCrary-style density test (clean Russell design):")
for k, v in res.items():
    print(f"  {k:>10s} = {v: .4f}")
print(f"\n  -> p = {res['pv']:.3f}. We FAIL to reject continuity (p > 0.05): no evidence of "
      f"manipulation. Good.")

fig, ax = plt.subplots(figsize=(8, 4.5))
ax.bar(cen_d[cen_d < 0], dens_d[cen_d < 0], width=1.8, color="#3b6", alpha=0.7, label="control side")
ax.bar(cen_d[cen_d >= 0], dens_d[cen_d >= 0], width=1.8, color="#c43", alpha=0.7, label="treated side")
ax.axvline(c, color="0.3", ls="--", lw=1)
ax.set_xlabel("running variable $X$"); ax.set_ylabel("density of $X$")
ax.set_title(f"Density of the running variable (clean): no cliff at $c$  (test $p={res['pv']:.2f}$)")
ax.legend()
fig.tight_layout()
print("A smooth density through c is what 'no manipulation' looks like.")

6. The global high-order polynomial pathology (Gelman–Imbens)

For years RD papers fit a single global high-order polynomial to all the data — Y on D plus X, X^2, X^3, X^4 over the whole range — and read the coefficient on D as the jump. Gelman and Imbens (2019) showed why this is a trap: polynomials are erratic near boundaries, distant data dictate the curvature at the cutoff, and the answer swings with the arbitrary polynomial order (chapter §5).

We make the pathology concrete. Our DGP put curvature only on the treated side. A global quartic fit to all 6{,}000 points lets that distant curvature — and the far-left large-cap points — twist the fit at the boundary, producing a tight-looking SE around a badly biased estimate. We then sweep the order p = 1, 2, 3, 4, 5 and watch the “jump” lurch around with no principled way to choose. Contrast with the local-linear / CCT estimate, which stays near the truth.

def global_poly_rd(df, c=0.0, order=4):
    '''Global polynomial RD on ALL data: Y on D and X-powers (separate slopes via D interactions).'''
    d = df.copy()
    d["Xc"] = d["X"] - c
    cols = {"const": 1.0, "D": d["D"]}
    for k in range(1, order + 1):
        cols[f"X{k}"] = d["Xc"] ** k
        cols[f"X{k}_D"] = (d["Xc"] ** k) * d["D"]      # let slopes differ on each side
    M = pd.DataFrame(cols)
    return sm.OLS(d["Y"], M).fit(cov_type="HC1")

print("Global polynomial RD (fit to ALL 6,000 points), tau_hat by order:")
for p in (1, 2, 3, 4, 5):
    r = global_poly_rd(df, c=c, order=p)
    print(f"  order {p}:  tau_hat = {r.params['D']: .3f}   se = {r.bse['D']:.3f}")
print(f"\ntrue tau = {tau_true}; CCT robust = {tau_robust:.3f}.")
print("The global estimate is order-sensitive and biased — a deceptively precise number in the wrong "
      "place.\nLesson: stay LOCAL, stay LOW-ORDER.  Trust the local-linear / CCT answer.")

7. Fuzzy RD: when the cutoff only nudges treatment (IV in disguise)

The real Russell design is fuzzy (chapter §4, §10). Crossing the rank cutoff does not force a fixed amount of passive money into a stock; it raises the probability and intensity of Russell-2000-tracking ownership — index weights are float-adjusted and banding rules add slack. So treatment T does not jump from 0 to 1 at c; its probability jumps.

This is instrumental variables with the above-cutoff indicator Z = \mathbf{1}\{X \ge c\} as the instrument. The effect is the Wald ratio from Week 3:

\tau_{\text{FRD}} = \frac{\text{jump in outcome at } c}{\text{jump in treatment at } c} = \frac{\text{reduced form}}{\text{first stage}}.

We simulate a strong, crisp first stage — treatment probability jumps from 0.15 just inside the R1000 to 0.85 just inside the R2000 (a first-stage jump of \approx 0.70) — and a true complier effect of \tau = 0.6. We use a larger N = 12{,}000 so the ratio is well-estimated, then call rdrobust(..., fuzzy=T), which does both jumps and the ratio in one call with standard errors that respect the ratio (never hand-divide two regressions and trust the naive SE — the Week 3 warning).

N_f = 12_000
Xf = rng.uniform(-50, 50, size=N_f)
mu_f = 5.0 + 0.02 * Xf                                  # symmetric smooth trend (clean fuzzy LATE)
p_treat = np.where(Xf >= c, 0.85, 0.15)                 # treatment PROBABILITY jumps 0.15 -> 0.85
T = (rng.uniform(size=N_f) < p_treat).astype(float)     # realized (fuzzy) treatment
Yf = mu_f + tau_true * T + rng.normal(0.0, 1.0, size=N_f)

# First stage and reduced form, each as its own sharp RD, then the ratio:
fs = rdrobust(y=T, x=Xf, c=0)
rf = rdrobust(y=Yf, x=Xf, c=0)
fs_jump = float(fs.coef.loc["Conventional"].iloc[0])
rf_jump = float(rf.coef.loc["Conventional"].iloc[0])
print(f"first-stage jump in P(treat) at c : {fs_jump:.3f}   (designed ~0.70)")
print(f"reduced-form jump in Y at c       : {rf_jump:.3f}")
print(f"naive Wald ratio (RF / FS)        : {rf_jump / fs_jump:.3f}\n")

# The right way: one fuzzy-RD call (numpy arrays; propagates ratio SEs correctly).
outf = rdrobust(y=Yf, x=Xf, c=0, fuzzy=T)
tau_f = float(outf.coef.loc["Robust"].iloc[0])
ci_f = outf.ci.loc["Robust"]
print(outf.rdmodel)
print(f"fuzzy-RD (IV/Wald) tau = {tau_f:.3f}, 95% robust CI = [{ci_f.iloc[0]:.3f}, {ci_f.iloc[1]:.3f}]")
brackets_f = (ci_f.iloc[0] <= tau_true <= ci_f.iloc[1])
print(f"Does it bracket the true complier effect tau = {tau_true}?  {brackets_f}")
assert brackets_f, "fuzzy-RD CI failed to bracket the true jump"
print("This is the Week-3 LATE for *compliers at the cutoff*: stocks that take on heavy passive\n"
      "ownership BECAUSE they crossed into the Russell 2000.")

Watch precision collapse as the first stage weakens

The fuzzy-RD denominator is the first-stage jump — the IV relevance condition. If the cutoff barely nudges the treatment probability, you are dividing by something near zero and the weak-instrument pathology of Chapter 3.5 returns in RD clothing: the point estimate gets unstable and the confidence interval explodes. We hold the true effect at 0.6 and shrink the first-stage jump from 0.70 to a sliver, watching the robust CI width balloon.

print("first-stage jump   tau_hat    robust 95% CI            CI width")
print("-" * 64)
for lo_p, hi_p in [(0.15, 0.85), (0.30, 0.70), (0.40, 0.60), (0.45, 0.55), (0.48, 0.52)]:
    pt = np.where(Xf >= c, hi_p, lo_p)
    Tw = (rng.uniform(size=N_f) < pt).astype(float)
    Yw = mu_f + tau_true * Tw + rng.normal(0.0, 1.0, size=N_f)
    ow = rdrobust(y=Yw, x=Xf, c=0, fuzzy=Tw)
    tw = float(ow.coef.loc["Robust"].iloc[0])
    ci = ow.ci.loc["Robust"]
    width = ci.iloc[1] - ci.iloc[0]
    print(f"     {hi_p - lo_p:.2f}        {tw: .3f}   [{ci.iloc[0]: .2f}, {ci.iloc[1]: .2f}]   {width:8.2f}")
print("\nAs the first-stage jump shrinks, the CI width explodes — the weak-instrument lesson, in RD form.")

Your Turn

Two extensions, both runnable below.

(A) Vary the bandwidth. You already saw the bias–variance trade-off in §3. Now compare your honest, data-driven choice against a hand-picked one. Call rdrobust with a forced wide bandwidth, e.g. rdrobust(y=..., x=..., c=0, h=40), and compare its estimate and CI to the CCT-chosen result from §4. Which is closer to the true 0.6? Is the forced-wide CI tighter but mis-located? (That tightness is exactly the trap from Question 2 of the chapter’s “Your Turn.”)

(B) Inject manipulation and watch the density test fire. Take the clean running variable and move a fraction of the just-below-cutoff stocks to just-above (as if a firm gamed its rank), then rerun mccrary_density_test. The clean design gave p \approx 0.57 (fail to reject — good). After manipulation the density should show a cliff at c and the test should reject. Try different fractions and see how much manipulation it takes before the test fires.

The cell below does both. Edit forced_h and manip_fraction and rerun to explore.

# (A) data-driven CCT vs. a forced-wide bandwidth -----------------------------------------
forced_h = 40.0
out_wide = rdrobust(y=df["Y"].to_numpy(), x=df["X"].to_numpy(), c=0, h=forced_h)
tw = float(out_wide.coef.loc["Robust"].iloc[0]); ciw = out_wide.ci.loc["Robust"]
print("(A) Bandwidth comparison (true tau = 0.6):")
print(f"    CCT data-driven (h={h_cct:.1f}): tau = {tau_robust:.3f}, "
      f"CI = [{ci_robust.iloc[0]:.3f}, {ci_robust.iloc[1]:.3f}], width = {ci_robust.iloc[1]-ci_robust.iloc[0]:.3f}")
print(f"    forced wide   (h={forced_h:.1f}): tau = {tw:.3f}, "
      f"CI = [{ciw.iloc[0]:.3f}, {ciw.iloc[1]:.3f}], width = {ciw.iloc[1]-ciw.iloc[0]:.3f}")
print("    -> the wide CI is often tighter but mis-located: precision is not accuracy.\n")

# (B) inject manipulation and rerun the density test --------------------------------------
manip_fraction = 0.6                       # share of just-below stocks that 'jump' the cutoff
X_manip = df["X"].to_numpy().copy()
victims = (X_manip >= -4) & (X_manip < 0)  # stocks just inside the Russell 1000
move = victims & (rng.uniform(size=len(X_manip)) < manip_fraction)
X_manip[move] = X_manip[move] + 4.0        # nudge them across the cutoff

res_m, cen_m, dens_m = mccrary_density_test(X_manip, c=c)
print(f"(B) Density test after manipulation (moved {move.sum()} stocks across c):")
print(f"    clean      : z = {res['z']: .2f}, p = {res['pv']:.3f}  (fail to reject -> no manipulation)")
print(f"    manipulated: z = {res_m['z']: .2f}, p = {res_m['pv']:.3f}  "
      f"({'REJECT -> manipulation detected' if res_m['pv'] < 0.05 else 'still not detected'})")

fig, ax = plt.subplots(figsize=(8, 4.5))
ax.bar(cen_m[cen_m < 0], dens_m[cen_m < 0], width=1.8, color="#3b6", alpha=0.7, label="control side")
ax.bar(cen_m[cen_m >= 0], dens_m[cen_m >= 0], width=1.8, color="#c43", alpha=0.7, label="treated side")
ax.axvline(c, color="0.3", ls="--", lw=1)
ax.set_xlabel("running variable $X$ (manipulated)"); ax.set_ylabel("density of $X$")
ax.set_title(f"Manipulated density: a cliff at $c$  (test $p={res_m['pv']:.3f}$)")
ax.legend()
fig.tight_layout()
print("\nThe cliff at c is the fingerprint of sorting — exactly what the McCrary test is built to catch.")