Regression Discontinuity

48.1 A rule with a cutoff (40 min)

NorthMart customers who spent $1,000 or more last year get Gold status: free shipping and priority support. Does Gold raise next-year spending?

Comparing Gold and non-Gold customers is hopeless: Gold customers were already the big spenders. But customers who spent $990 and $1,010 are almost identical; only one has Gold. If nothing else changes at exactly $1,000, any jump in next-year spending at the cutoff is the effect of Gold. That's a regression discontinuity (RD) design (Thistlethwaite & Campbell, 1960). The variable that decides treatment (last year's spend) is the running variable.

The estimand is narrow: the effect for customers near $1,000, not for a customer who spends $5,000. Keep that in every sentence you write about an RD.

Create causalkit/tiers.py (provided; read only the docstrings):

"""
NorthMart Gold tier (simulated). Provided for Module 5: read the docstrings, not the code.

Customers whose spending last year reached $1,000 are upgraded to Gold (free shipping, priority support).

    from causalkit.tiers import customers, truth

    df = customers()                    # 30,000 customers; sharp rule: Gold if pre_spend >= 1,000
    df = customers(fuzzy=True)          # managers upgrade some keen customers below the line; some above never activate
    df = customers(manipulation=0.3)    # customers just below $1,000 who know the rule buy their way over it
    truth(df)                           # true effect of Gold at the cutoff (for checking only)

Columns an analyst has:
    customer_id, pre_spend (last year, $; the running variable), age, tenure_years, engaged,
    gold (1 = has Gold status), spend (next year, $; the outcome)
Columns starting with '_' are the simulator's hidden truth.
"""
import numpy as np
import polars as pl

CUTOFF = 1000.0


def customers(n: int = 30_000, seed: int = 54, effect: float = 120.0, fuzzy: bool = False,
              manipulation: float = 0.0) -> pl.DataFrame:
    """effect: true effect of Gold on next-year spend for customers AT the cutoff ($); it grows with pre_spend.
    manipulation: share of customers with organic spend in [$850, $1,000) who top up to just over $1,000.
    They are keen customers: their next-year spend is higher even without Gold."""
    rng = np.random.default_rng(seed)
    engaged = rng.binomial(1, 0.5, n)
    tenure = np.round(rng.gamma(2.0, 2.0, n), 1)
    age = np.clip(rng.normal(45, 13, n), 18, 85).round().astype(int)
    intent = rng.normal(0, 1, n)
    pre = np.exp(rng.normal(6.55 + 0.25 * engaged + 0.03 * tenure + 0.15 * intent, 0.55, n))
    keen = np.zeros(n, bool)
    if manipulation > 0:
        near = (pre >= 850) & (pre < CUTOFF)
        keen = near & (rng.random(n) < manipulation)
        pre = np.where(keen, CUTOFF + rng.uniform(0, 40, n), pre)
    pre = np.round(pre, 2)
    x = (pre - CUTOFF) / 1000
    y0 = (900 + 650 * x - 180 * x**2 + 60 * engaged + 8 * tenure + 60 * intent + 150 * keen
          + rng.normal(0, 220, n))
    tau = effect + 60 * x
    above = pre >= CUTOFF
    if fuzzy:
        keen_share = 1 / (1 + np.exp(-1.5 * intent))           # managers upgrade (and keen customers activate)
        p_below = (0.08 + 0.25 * np.exp(-np.abs(x) * 8)) * 2 * keen_share
        gold = np.where(above, rng.random(n) < 0.6 + 0.4 * keen_share, rng.random(n) < p_below)
    else:
        gold = above
    return pl.DataFrame({
        "customer_id": [f"G{i:05d}" for i in range(1, n + 1)],
        "pre_spend": pre, "age": age, "tenure_years": tenure, "engaged": engaged,
        "gold": gold.astype(np.int64), "spend": np.round(y0 + tau * gold, 2),
        "_y0": y0, "_tau": tau, "_keen": keen,
    })


def truth(df: pl.DataFrame) -> dict:
    """The effect at the cutoff, and the average effect for everyone with Gold."""
    at_cutoff = (df["_tau"] - 60 * (df["pre_spend"] - CUTOFF) / 1000).mean()
    return {"effect_at_cutoff": round(float(at_cutoff), 1),
            "ATT": round(df.filter(pl.col("gold") == 1)["_tau"].mean(), 1)}


def public(df: pl.DataFrame) -> pl.DataFrame:
    return df.select([c for c in df.columns if not c.startswith("_")])

48.2 The RD plot (40 min)

Create day48/rd_plot.py:

# day48/rd_plot.py   (python -m day48.rd_plot)
"""
Gold status starts at $1,000 of spending. Customers at $990 and $1,010 are nearly identical, except one
has Gold. Plot next-year spend against last year's spend in $25 bins: a jump at $1,000 is the effect.
"""
import matplotlib.pyplot as plt
import polars as pl

from causalkit.tiers import CUTOFF, customers, public

df = public(customers())
print(df.head(3))
print(f"{df.height:,} customers; {df['gold'].mean():.1%} have Gold")

bins = (df.filter(pl.col("pre_spend").is_between(400, 1600))
        .with_columns(bin=(pl.col("pre_spend") // 25) * 25 + 12.5)
        .group_by("bin").agg(pl.col("spend").mean(), pl.col("gold").mean(), n=pl.len()).sort("bin"))
fig, axes = plt.subplots(1, 2, figsize=(11, 4))
for ax, col, label in [(axes[0], "spend", "next-year spend ($)"), (axes[1], "gold", "share with Gold")]:
    left, right = bins.filter(pl.col("bin") < CUTOFF), bins.filter(pl.col("bin") >= CUTOFF)
    ax.scatter(left["bin"], left[col], s=12, color="#64748b")
    ax.scatter(right["bin"], right[col], s=12, color="#d97706")
    ax.axvline(CUTOFF, color="#dc2626", ls="--", lw=0.8)
    ax.set_xlabel("last year's spend ($)")
    ax.set_ylabel(label)
fig.tight_layout()
fig.savefig("figures/day48_rd_plot.png", dpi=120)
print("Saved figures/day48_rd_plot.png")

python -m day48.rd_plot

│ customer_id ┆ pre_spend ┆ age ┆ tenure_years ┆ engaged ┆ gold ┆ spend   │
│ G00001      ┆ 1304.93   ┆ 46  ┆ 11.6         ┆ 0       ┆ 1    ┆ 1492.84 │
│ G00002      ┆ 782.06    ┆ 61  ┆ 14.1         ┆ 1       ┆ 0    ┆ 826.5   │
│ G00003      ┆ 391.97    ┆ 30  ┆ 10.5         ┆ 0       ┆ 0    ┆ 574.76  │
30,000 customers; 42.1% have Gold
Saved figures/day48_rd_plot.png

Open the figure. Left: next-year spend rises smoothly with last year's spend, then jumps by roughly $120 at $1,000. Right: the share with Gold goes from 0 to 1 at the cutoff, so this is a sharp RD. Always show the RD plot: if you can't see a jump, a regression shouldn't convince anyone there is one.

48.3 Estimating the jump (65 min)

Concept. Fit a line on each side of the cutoff, using only customers close to it, and measure the gap between the two lines at $1,000. Two choices:

• Bandwidth h: how close is "close"? Narrow = less bias (lines fit better near the cutoff), more noise.
• Kernel: a triangular kernel weights customers more the closer they are to the cutoff.

rdrobust (Calonico, Cattaneo & Titiunik, 2014) picks the bandwidth that minimises mean squared error, and reports a robust bias-corrected confidence interval: the MSE-optimal bandwidth leaves some bias, and the robust CI accounts for it. It's the standard RD tool.

Create day48/local_linear.py:

# day48/local_linear.py   (python -m day48.local_linear)
"""
Estimating the jump. Three ways:
  a global polynomial (common, fragile),
  a local linear regression by hand: only customers within $h of the cutoff, nearer ones weighted more,
  rdrobust (Calonico, Cattaneo & Titiunik): data-driven bandwidth and a bias-corrected robust CI.
"""
import numpy as np
import statsmodels.formula.api as smf
from rdrobust import rdrobust

from causalkit.tiers import CUTOFF, customers, public, truth

full = customers()
df = public(full).to_pandas()
df["x"] = (df.pre_spend - CUTOFF) / 1000              # distance from the cutoff, $ thousands
df["above"] = (df.x >= 0).astype(int)


def local_linear(d, h):
    """Triangular-kernel local linear regression with separate slopes on each side."""
    near = d[d.x.abs() <= h / 1000].copy()
    near["w"] = 1 - near.x.abs() / (h / 1000)
    fit = smf.wls("spend ~ above * x", near, weights=near.w).fit(cov_type="HC1")
    return fit.params["above"], fit.bse["above"], len(near)


quartic = smf.ols("spend ~ above * (x + I(x**2) + I(x**3) + I(x**4))", df).fit(cov_type="HC1")
print(f"global quartic polynomial, all customers:  {quartic.params['above']:6.1f}  (SE {quartic.bse['above']:.1f})")
est, se, n = local_linear(df, 200)
print(f"local linear, h = $200 ({n:,} customers):   {est:6.1f}  (SE {se:.1f})")

rd = rdrobust(df.spend, df.pre_spend, c=CUTOFF)
h = rd.bws.iloc[0, 0]
print(f"rdrobust: bandwidth ${h:.0f}, {sum(rd.N_h):,} customers")
print(f"   conventional estimate {rd.coef.iloc[0, 0]:6.1f}  95% CI [{rd.ci.iloc[0, 0]:.1f}, {rd.ci.iloc[0, 1]:.1f}]")
print(f"   robust bias-corrected  {rd.coef.iloc[1, 0]:6.1f}  95% CI [{rd.ci.iloc[2, 0]:.1f}, {rd.ci.iloc[2, 1]:.1f}]")
t = truth(full)
print(f"Truth: effect at the cutoff {t['effect_at_cutoff']:.1f}; average for all Gold customers {t['ATT']:.1f}")

python -m day48.local_linear

global quartic polynomial, all customers:   112.0  (SE 9.6)
local linear, h = $200 (8,050 customers):    124.7  (SE 11.3)
rdrobust: bandwidth $302, 12,065 customers
   conventional estimate  123.9  95% CI [106.0, 141.9]
   robust bias-corrected   123.4  95% CI [101.9, 145.0]
Truth: effect at the cutoff 120.0; average for all Gold customers 158.9

Local linear and rdrobust are within $5 of the truth, and the robust CI contains it. Note the last line: Gold is worth $159 on average for all Gold customers (bigger spenders gain more), but the RD estimates the $120 effect at the cutoff. An RD can't tell you the effect of Gold for a $3,000 customer.

48.4 Bandwidths and polynomials (40 min)

Try it. Estimate the jump by local linear regression with h = $50, $100, $200, $400, $800, and with global polynomials of order 1, 2, 4 and 6 fitted to all customers. Create day48/bandwidths.py.

# day48/bandwidths.py   (python -m day48.bandwidths)
import statsmodels.formula.api as smf

from causalkit.tiers import CUTOFF, customers, public

df = public(customers()).to_pandas()
df["x"] = (df.pre_spend - CUTOFF) / 1000
df["above"] = (df.x >= 0).astype(int)
for h in [50, 100, 200, 400, 800]:
    near = df[df.x.abs() <= h / 1000].copy()
    near["w"] = 1 - near.x.abs() / (h / 1000)
    fit = smf.wls("spend ~ above * x", near, weights=near.w).fit(cov_type="HC1")
    print(f"local linear, h = ${h:<4} {fit.params['above']:6.1f}  (SE {fit.bse['above']:5.1f}, n = {len(near):,})")
for k in [1, 2, 4, 6]:
    terms = " + ".join(f"I(x**{p})" for p in range(1, k + 1))
    fit = smf.ols(f"spend ~ above * ({terms})", df).fit(cov_type="HC1")
    print(f"global polynomial, order {k}:  {fit.params['above']:6.1f}  (SE {fit.bse['above']:5.1f})")

local linear, h = $50    129.6  (SE  22.4, n = 2,001)
local linear, h = $100   131.4  (SE  15.9, n = 4,051)
local linear, h = $200   124.7  (SE  11.3, n = 8,050)
local linear, h = $400   124.9  (SE   8.1, n = 15,756)
local linear, h = $800   124.6  (SE   6.0, n = 26,255)
global polynomial, order 1:   304.3  (SE  13.2)
global polynomial, order 2:   133.1  (SE   6.2)
global polynomial, order 4:   112.0  (SE   9.6)
global polynomial, order 6:   121.9  (SE  13.2)

48.5 Day 48 mini-project: a fuzzy discontinuity (65 min)

In reality, store managers upgrade some keen customers below $1,000, and some customers above it never activate Gold. Crossing $1,000 still makes Gold much more likely, so it's an instrument (Day 37), local to the cutoff. Create day48/fuzzy_rd.py:

# day48/fuzzy_rd.py   (python -m day48.fuzzy_rd)
"""
Mini-project: a FUZZY discontinuity. Store managers upgrade some keen customers below $1,000, and some
customers above it never activate Gold. Crossing $1,000 still makes Gold much more likely: an instrument,
exactly like Day 37's invitation. Effect = jump in spend / jump in Gold share (a LATE at the cutoff).
"""
import numpy as np
import statsmodels.formula.api as smf
from rdrobust import rdrobust

from causalkit.tiers import CUTOFF, customers, public, truth


def near_cutoff(df, h=200):
    df = df.assign(x=(df.pre_spend - CUTOFF) / 1000, above=(df.pre_spend >= CUTOFF).astype(int))
    near = df[df.x.abs() <= h / 1000].copy()
    near["w"] = 1 - near.x.abs() / (h / 1000)
    return near


def wald_and_naive(df):
    near = near_cutoff(df)
    first = smf.wls("gold ~ above * x", near, weights=near.w).fit().params["above"]
    reduced = smf.wls("spend ~ above * x", near, weights=near.w).fit().params["above"]
    naive = smf.wls("spend ~ gold + above:x + x", near, weights=near.w).fit().params["gold"]
    return first, reduced, naive


full = customers(fuzzy=True)
df = public(full).to_pandas()
first, reduced, naive = wald_and_naive(df)
print(f"Gold share jumps by {first:.2f} at $1,000; spend jumps by {reduced:.1f}")
print(f"Wald ratio (h = $200):                        {reduced / first:6.1f}")
print(f"Regress spend on gold near the cutoff:        {naive:6.1f}   <- who gets upgraded is not random")
rd = rdrobust(df.spend, df.pre_spend, c=CUTOFF, fuzzy=df.gold)
print(f"rdrobust fuzzy: {rd.coef.iloc[0, 0]:.1f}, robust 95% CI [{rd.ci.iloc[2, 0]:.1f}, {rd.ci.iloc[2, 1]:.1f}]")
print(f"Truth at the cutoff: {truth(full)['effect_at_cutoff']:.1f}")

waldS, naiveS = [], []
for s in range(20):
    f, r, n = wald_and_naive(public(customers(seed=900 + s, fuzzy=True)).to_pandas())
    waldS.append(r / f)
    naiveS.append(n)
print(f"\nOver 20 datasets: Wald ratio {np.mean(waldS):.1f} (SD {np.std(waldS):.1f}); "
      f"naive regression {np.mean(naiveS):.1f} (SD {np.std(naiveS):.1f})")

python -m day48.fuzzy_rd

Gold share jumps by 0.51 at $1,000; spend jumps by 66.7
Wald ratio (h = $200):                         129.9
Regress spend on gold near the cutoff:         153.3   <- who gets upgraded is not random
rdrobust fuzzy: 131.8, robust 95% CI [75.5, 190.1]
Truth at the cutoff: 120.0

Over 20 datasets: Wald ratio 124.3 (SD 17.4); naive regression 154.6 (SD 7.2)

• The fuzzy RD is the Wald estimator from Day 37: reduced-form jump ÷ first-stage jump. It's unbiased on average ($124 over 20 datasets), but noisy: dividing by a jump of 0.5 doubles the noise.
• Regressing spend on Gold near the cutoff is biased by +$35 on average, consistently. Managers upgrade keen customers, who would spend more anyway: selection bias, even within $200 of the cutoff.
• The CI is three times wider than in the sharp RD. Fuzziness costs precision, like a weak first stage.

48.6 Save your work with Git (20 min)

git add .
git commit -m "Day 48: Gold tier simulator, RD plot, local linear, rdrobust, fuzzy RD"
git push

Day 48 checklist

• I can explain the RD design, the running variable, and what "effect at the cutoff" means
• I always show an RD plot with binned means
• I can estimate the jump by local linear regression and with rdrobust's robust bias-corrected CI
• I avoid global polynomials, and show stability across bandwidths
• I can estimate a fuzzy RD as a Wald ratio, and explain why it's noisier