import matplotlib
matplotlib.use("Agg") # headless, non-interactive backend
import numpy as np
import pandas as pd
import statsmodels.formula.api as smf
rng = np.random.default_rng(42)
# Panel dimensions, in the spirit of the CPS state-year panels BDM used.
N_STATES = 50 # ~50 states
N_YEARS = 20 # ~20 years -> a LONG panel, where serial correlation bites
ALPHA = 0.05 # nominal test size we promised to honor
print(f"Panel: {N_STATES} states x {N_YEARS} years = {N_STATES * N_YEARS} cells")
print(f"Nominal test size we are auditing: {ALPHA:.0%}")nb5.5 — The BDM Placebo-Law Simulation
Companion to Ch 5.5, the reader’s guide to Bertrand, Duflo & Mullainathan (2004), “How Much Should We Trust Differences-in-Differences Estimates?” (QJE 119(1): 249–275).
The chapter promised you a trick so clean you could run it on your own laptop. Here it is. BDM caught a whole literature reporting standard errors that were too small — not because anyone cheated, but because the conventional difference-in-differences (DiD) formula assumes the outcome is independent across years within a state, and slow-moving economic outcomes flatly violate that. Too-small standard errors mean too-large t-statistics, which mean you reject the null far too often, which means you “discover” effects that are not there.
We are going to measure that over-rejection. The plan, lifted straight from §2 and §7 of the chapter:
- Build a state-year panel whose outcome is serially correlated (an AR(1) process) and has no real treatment effect whatsoever.
- Sprinkle a fake “law” onto a random subset of states starting in a random year, run the standard two-way fixed-effects (TWFE) DiD with naive (classical) standard errors, and check whether the placebo coefficient comes back “significant” at 5%.
- Repeat hundreds of times. Because the true effect is zero by construction, the fraction of “significant” runs is the true size of the test. If the machinery were honest it would be 5%. We will watch it come back far higher.
- Apply the fixes — cluster by state, and collapse to pre/post — and watch the rejection rate fall back toward 5%.
- Turn the dial on the AR(1) persistence and confirm the over-rejection grows with serial correlation.
This is the Week-1 idea of the size of a test (the Type I error rate, the \alpha you promised to honor) wearing the Week-4 clothes of DiD and clustering. By the end you will have rebuilt BDM’s headline result from scratch.
Setup
We fix one seeded generator, rng = np.random.default_rng(42), so the notebook reproduces bit-for-bit (the CONVENTIONS.md rule), and force matplotlib’s non-interactive Agg backend so it runs headless — laptop, server, or CI — with no display attached.
We lean on statsmodels for OLS with both classical and cluster-robust standard errors. Everything else is numpy and pandas.
1. A state-year panel with serial correlation and NO real effect
We simulate one outcome Y_{st} for state s in year t. Three pieces, all labeled:
Y_{st} = \underbrace{\mu_s}_{\text{state level}} + \underbrace{\lambda_t}_{\text{year shock}} + \underbrace{u_{st}}_{\text{AR(1) within-state error}}, \qquad u_{st} = \rho\, u_{s,t-1} + \varepsilon_{st}.
The state means \mu_s and common year shocks \lambda_t are exactly what the state and year fixed effects in a DiD will absorb, so they are harmless. The villain is u_{st}: a first-order autoregressive — AR(1) — error, where this year’s shock is \rho times last year’s plus fresh noise. When \rho is large, a state that runs high one year tends to keep running high: its residuals form long, gently drifting runs. That is serial correlation, the within-unit autocorrelation you first met as a disease in Chapter 2.4, and it is precisely the genuine time-series structure of the economy that the conventional DiD standard error ignores.
Crucially, there is no treatment term in this equation. The true policy effect is zero, by construction. Any “significant” placebo law we find later is, with no ambiguity, a false positive.
We draw the innovations with the AR(1) stationary variance so that \operatorname{Var}(u_{st}) does not depend on \rho — that way, when we later sweep \rho, we are changing the correlation structure and not secretly also changing how noisy the data are.
def simulate_panel(rng, n_states=N_STATES, n_years=N_YEARS, rho=0.8,
sigma_u=1.0, sigma_state=1.0, sigma_year=0.5):
'''One state-year panel with AR(1) within-state errors and NO treatment effect.
Returns a tidy long DataFrame with columns: state, year, Y.
rho controls serial correlation; the innovation variance is scaled so that the
stationary Var(u) = sigma_u**2 regardless of rho (clean dial for the rho sweep).
'''
mu_s = rng.normal(0.0, sigma_state, size=n_states) # state fixed levels
lam_t = rng.normal(0.0, sigma_year, size=n_years) # common year shocks
# AR(1) errors, one independent series per state, drawn at the stationary variance.
eps_sd = sigma_u * np.sqrt(1.0 - rho**2) if rho < 1.0 else sigma_u
u = np.empty((n_states, n_years))
# initialise t=0 from the stationary distribution so there is no burn-in transient
u[:, 0] = rng.normal(0.0, sigma_u, size=n_states)
for t in range(1, n_years):
u[:, t] = rho * u[:, t-1] + rng.normal(0.0, eps_sd, size=n_states)
Y = mu_s[:, None] + lam_t[None, :] + u # (n_states, n_years)
states = np.repeat(np.arange(n_states), n_years)
years = np.tile(np.arange(n_years), n_states)
df = pd.DataFrame({"state": states, "year": years, "Y": Y.reshape(-1)})
return df
demo = simulate_panel(rng, rho=0.8)
print(demo.head(3).to_string(index=False))
print(f"\nshape: {demo.shape[0]} rows ({demo['state'].nunique()} states x {demo['year'].nunique()} years)")Seeing the serial correlation before we estimate
Two pictures make the disease concrete. On the left, a few states’ outcome paths over time at \rho = 0.8: notice how each line drifts — high years cluster next to high years, low next to low — rather than jittering independently. On the right, the lag-1 autocorrelation of the within-state residuals (after removing each state’s mean) for two worlds: \rho = 0 (independent errors, the assumption the naive SE makes) and \rho = 0.8 (the reality). The independent world piles up near zero; the persistent world sits up near 0.8. The naive DiD standard error treats the right-hand world as if it were the left-hand one — and that is the whole bug.
import matplotlib.pyplot as plt
def within_state_lag1_acf(df):
'''Lag-1 autocorrelation of residuals after removing each state's own mean.'''
g = df.copy()
g["resid"] = g["Y"] - g.groupby("state")["Y"].transform("mean")
acfs = []
for _, blk in g.groupby("state"):
r = blk.sort_values("year")["resid"].to_numpy()
if r[:-1].std() > 0 and r[1:].std() > 0:
acfs.append(np.corrcoef(r[:-1], r[1:])[0, 1])
return np.array(acfs)
viz_rng = np.random.default_rng(7)
panel_hi = simulate_panel(viz_rng, rho=0.8)
panel_lo = simulate_panel(viz_rng, rho=0.0)
fig, axes = plt.subplots(1, 2, figsize=(13, 4.2))
ax = axes[0]
for s in range(6):
blk = panel_hi[panel_hi["state"] == s].sort_values("year")
ax.plot(blk["year"], blk["Y"], marker="o", ms=3, lw=1.3, label=f"state {s}")
ax.set_title(r"A. Outcome paths, $\rho=0.8$: each state drifts (serial correlation)")
ax.set_xlabel("year"); ax.set_ylabel(r"$Y_{st}$"); ax.legend(ncol=2, fontsize=8)
ax = axes[1]
ax.hist(within_state_lag1_acf(panel_lo), bins=15, alpha=0.6,
color="#999999", label=r"$\rho=0$ (independent)")
ax.hist(within_state_lag1_acf(panel_hi), bins=15, alpha=0.6,
color="#c0392b", label=r"$\rho=0.8$ (persistent)")
ax.axvline(0.0, color="k", lw=0.8, ls=":")
ax.set_title("B. Within-state lag-1 autocorrelation of residuals")
ax.set_xlabel("lag-1 autocorrelation"); ax.set_ylabel("# states"); ax.legend()
fig.tight_layout(); fig.savefig("nb5.5-serial-correlation.png", dpi=110)
plt.close(fig)
print("saved nb5.5-serial-correlation.png")
print(f"mean within-state ACF: rho=0 -> {within_state_lag1_acf(panel_lo).mean():+.3f}",
f"| rho=0.8 -> {within_state_lag1_acf(panel_hi).mean():+.3f}")2. Assign a fake law and run one placebo DiD
Now the mischief. We pick a random half of the states to be “treated,” and for each treated state we pick a random year at which its fake law “switches on” (and stays on). This is the persistent treatment dummy D_{st} — off for years, then on for years — that the chapter (via the Moulton intuition) flagged as the worst case when paired with serially correlated errors. The untreated states never switch on.
We then estimate the exact TWFE DiD specification you carried through all of Week 4,
Y_{st} = \alpha_s + \lambda_t + \beta\, D_{st} + \varepsilon_{st},
with state fixed effects \alpha_s (C(state)), year fixed effects \lambda_t (C(year)), and the placebo treatment D_{st} (D). We read off \hat\beta and its classical (iid) standard error, form t = \hat\beta / \mathrm{SE}, and ask the Week-1 question: is |t| > 1.96? Remember — the honest answer should be “yes” only about 5% of the time, because \beta truly is zero.
Empirical-spec box. Outcome: Y_{st} (simulated, serially correlated). Treatment: placebo law D_{st} (random states, random onset year; true effect = 0). Controls: none. Fixed effects: state + year. Clustering: none yet (classical SE — the thing on trial). Sample: full simulated panel. Identifying assumption (deliberately one we will violate): residuals iid across years within a state.
def assign_placebo(df, rng, treat_frac=0.5):
'''Add a placebo treatment dummy D: a random treat_frac of states switch on
at a random year (and stay on). True effect is zero -- nothing is added to Y.'''
out = df.copy()
states = out["state"].unique()
n_treat = max(1, int(round(treat_frac * len(states))))
treated = set(rng.choice(states, size=n_treat, replace=False))
years = np.sort(out["year"].unique())
# one random onset year per treated state (not the very last year, so D varies)
onset = {}
for s in treated:
onset[s] = rng.choice(years[:-1]) if len(years) > 1 else years[0]
D = np.zeros(len(out), dtype=float)
st = out["state"].to_numpy(); yr = out["year"].to_numpy()
for i in range(len(out)):
if st[i] in treated and yr[i] >= onset[st[i]]:
D[i] = 1.0
out["D"] = D
out["ever_treated"] = out["state"].isin(treated).astype(float)
# a single common split year (median of treated onsets) for the collapse-to-2-periods fix
out.attrs["split_year"] = int(np.median(list(onset.values())))
return out
# One illustrative placebo draw with classical SEs.
one = assign_placebo(simulate_panel(rng, rho=0.8), rng)
fit_naive = smf.ols("Y ~ D + C(state) + C(year)", data=one).fit() # classical/iid SE
b = fit_naive.params["D"]
se = fit_naive.bse["D"]
t = b / se
print(f"placebo beta_hat = {b:+.4f}")
print(f"classical SE = {se:.4f}")
print(f"t-stat = {t:+.3f} -> |t| > 1.96 ? {abs(t) > 1.96}")
print("(true effect is exactly 0, so a 'significant' result here is a FALSE POSITIVE)")One draw proves nothing — a single coin flip cannot tell you whether a coin is fair. The size of a test is a statement about the long run: rerun the placebo experiment many times and count how often it cries wolf. That is exactly the Week-1 coin-flip lab, where you built a universe with a genuinely fair coin and measured how often your test wrongly called it biased. Here the “fair coin” is the genuinely null law.
We package one placebo replication into a function that returns the t-statistic under each standard-error flavor we care about, so the Monte Carlo loop can simply count rejections.
def one_replication(rng, rho, n_states=N_STATES, n_years=N_YEARS, treat_frac=0.5):
'''Simulate a null panel, assign a placebo law, fit TWFE DiD, and return the
placebo t-stat under three inference methods:
'naive' : classical/iid SE (the procedure on trial)
'cluster' : cluster-robust SE by state (BDM fix 1)
'collapse' : collapse each state to one pre + one post mean, two-period DiD (BDM fix 2)
Returns a dict of t-stats.
'''
df = simulate_panel(rng, n_states=n_states, n_years=n_years, rho=rho)
df = assign_placebo(df, rng, treat_frac=treat_frac)
out = {}
# --- TWFE fit, reused for naive and clustered (same point estimate, different SE) ---
fit = smf.ols("Y ~ D + C(state) + C(year)", data=df).fit()
out["naive"] = fit.params["D"] / fit.bse["D"]
fit_cl = smf.ols("Y ~ D + C(state) + C(year)", data=df).fit(
cov_type="cluster", cov_kwds={"groups": df["state"]})
out["cluster"] = fit_cl.params["D"] / fit_cl.bse["D"]
# --- Collapse to pre/post: reduce each state to ONE before- and ONE after-mean,
# split at a single COMMON year, then run a 2x2 DiD with treated AND control states.
# Keeping the controls is what lets the common year shocks difference out; the
# interaction treat:post is the collapsed placebo effect. With effectively two
# periods per state there is no within-period serial correlation left to misread.
g = df.copy()
split = g.attrs["split_year"]
g["post"] = (g["year"] >= split).astype(float)
coll = (g.groupby(["state", "post", "ever_treated"], observed=True)["Y"]
.mean().reset_index())
# each state appears once with post=0 and once with post=1 (a 2-period panel)
if coll["post"].nunique() == 2 and coll["ever_treated"].nunique() == 2:
fit_c = smf.ols("Y ~ ever_treated * post + C(state)", data=coll).fit()
key = "ever_treated:post"
out["collapse"] = fit_c.params[key] / fit_c.bse[key]
else:
out["collapse"] = np.nan
return out
print(one_replication(rng, rho=0.8))3. The headline: the naive rejection rate is far above 5%
Now we industrialize. We run N_REPS placebo replications at a strongly persistent \rho = 0.8, and for the naive classical standard error we count the fraction of runs with |t| > 1.96. That fraction is a direct, brute-force estimate of the test’s true size.
If the standard errors were honest it would land near 0.05. Watch what actually happens. (We use 600 reps — modest, so the notebook runs in well under a minute, but plenty to make the over-rejection unmistakable.)
N_REPS = 600
mc_rng = np.random.default_rng(42) # fresh, seeded generator for the Monte Carlo
methods = ["naive", "cluster", "collapse"]
tstats = {m: np.empty(N_REPS) for m in methods}
for r in range(N_REPS):
res = one_replication(mc_rng, rho=0.8)
for m in methods:
tstats[m][r] = res[m]
def rejection_rate(tarr, crit=1.96):
tarr = tarr[~np.isnan(tarr)]
return float(np.mean(np.abs(tarr) > crit))
rej_naive = rejection_rate(tstats["naive"])
print(f"NAIVE (classical/iid) SE, rho=0.8, {N_REPS} placebo laws")
print(f" empirical rejection rate at 5% = {rej_naive:.3f}")
print(f" nominal (honest) rate = {ALPHA:.3f}")
print(f" --> the advertised 5% test fires {rej_naive/ALPHA:.0f}x too often on PURE NOISE")There it is — BDM’s headline table, rebuilt on your laptop. The conventional standard error, applied to serially correlated data, turns an advertised 5% false-positive rate into something many times larger. The abstract worry (“the SEs might be a little off”) has become a number you measured (“the test that claims 5% actually fires on noise a third of the time or so”). And a number changes behavior in a way a caveat never does.
Why does this happen? The classical SE divides as if you have N \times T independent observations. But when the within-state errors are persistent, the effective number of independent observations is far smaller — each state’s twenty years carry much less than twenty observations’ worth of independent information. The formula over-counts information, so it reports a standard error that is too small, a t-stat that is too big, and a verdict that is too confident.
4. The fixes: cluster by state, and collapse to pre/post
The same Monte Carlo loop already recorded the t-stats under the two cures from §4 of the chapter, so we can read their rejection rates straight off.
- Cluster by state. The cluster-robust standard error of Chapter 2.4, with the cluster set to the unit at which treatment is assigned — the state. It lets each state’s residuals be arbitrarily correlated across all its years (the block-diagonal \boldsymbol\Omega, one block per state), which is exactly the serial-correlation pattern the disease consists of. With ~50 states we are in the safe many-clusters zone.
- Collapse to pre/post. Throw away the time dimension’s internal structure: split at one common year and average each state into a single “before” and “after” number, keeping both treated and control states so the common year shocks still difference out. This reduces the panel to effectively two periods and runs a 2\times2 DiD (the
ever_treated:postinteraction). With only two periods there is no within-period serial correlation left to misread. The cost — which §6 of the chapter is upfront about — is power: you discarded information, so true effects get harder to detect.
We expect both to drag the rejection rate back toward the honest 5%.
rej = {m: rejection_rate(tstats[m]) for m in methods}
summary = pd.DataFrame({
"method": ["Naive (classical SE)", "Cluster by state", "Collapse to pre/post"],
"rejection_rate_5pct": [rej["naive"], rej["cluster"], rej["collapse"]],
})
summary["x_nominal"] = (summary["rejection_rate_5pct"] / ALPHA).round(1)
print(f"rho = 0.8, {N_REPS} placebo replications, nominal size = {ALPHA:.0%}\n")
print(summary.to_string(index=False))Bar chart: rejection rate by method
The single picture that tells the whole story. The dashed line is the honest 5% we were promised. The naive bar towers over it; the two fixes sit back down near the line.
fig, ax = plt.subplots(figsize=(8.2, 4.6))
labels = summary["method"].tolist()
rates = summary["rejection_rate_5pct"].to_numpy()
colors = ["#c0392b", "#2e86c1", "#27ae60"]
bars = ax.bar(labels, rates, color=colors, width=0.6)
ax.axhline(ALPHA, color="k", ls="--", lw=1.2, label=f"nominal 5% (honest size)")
for b, v in zip(bars, rates):
ax.text(b.get_x() + b.get_width()/2, v + 0.005, f"{v:.2f}",
ha="center", va="bottom", fontsize=11, fontweight="bold")
ax.set_ylabel("placebo rejection rate at 5%")
ax.set_title(r"BDM placebo experiment: false-positive rate by method ($\rho=0.8$)")
ax.set_ylim(0, max(rates.max() * 1.25, 0.1))
ax.legend()
fig.tight_layout(); fig.savefig("nb5.5-rejection-by-method.png", dpi=110)
plt.close(fig)
print("saved nb5.5-rejection-by-method.png")5. Turning the dial: over-rejection grows with serial correlation
The chapter’s §7 Exercise 2 asks you to prove the culprit is serial correlation, not chance, by sweeping the AR(1) parameter \rho from 0 (independent errors) up toward 1 (highly persistent). The prediction is sharp: at \rho = 0 the naive test should sit near the honest 5%, because the iid assumption is actually true; as \rho climbs, the naive rejection rate should climb with it, while the clustered rate stays anchored near 5%.
We use fewer reps per \rho (the pattern is robust and we are sweeping several values) and a fresh seeded generator so the sweep is reproducible.
rhos = [0.0, 0.2, 0.4, 0.6, 0.8, 0.9]
N_REPS_SWEEP = 300
sweep_rng = np.random.default_rng(42)
rows = []
for rho in rhos:
tn = np.empty(N_REPS_SWEEP); tc = np.empty(N_REPS_SWEEP)
for r in range(N_REPS_SWEEP):
res = one_replication(sweep_rng, rho=rho)
tn[r] = res["naive"]; tc[r] = res["cluster"]
rows.append({"rho": rho,
"naive": rejection_rate(tn),
"cluster": rejection_rate(tc)})
sweep = pd.DataFrame(rows)
print(sweep.to_string(index=False))The mechanism, in one plot
The naive curve climbs steeply with \rho — at \rho = 0 it hugs the honest 5% line (the iid assumption is true there, so the naive SE is fine), and it lifts off as persistence grows, reaching tens of percent at \rho = 0.8–0.9. The clustered curve stays flat near 5% across the whole range. That contrast is the proof: the over-rejection is caused by serial correlation, and clustering by state defuses it. This is BDM’s mechanism evidence, rebuilt.
fig, ax = plt.subplots(figsize=(8.2, 4.8))
ax.plot(sweep["rho"], sweep["naive"], "o-", color="#c0392b", lw=2,
label="naive (classical SE)")
ax.plot(sweep["rho"], sweep["cluster"], "s-", color="#2e86c1", lw=2,
label="cluster by state")
ax.axhline(ALPHA, color="k", ls="--", lw=1.2, label="nominal 5%")
ax.set_xlabel(r"AR(1) serial-correlation parameter $\rho$")
ax.set_ylabel("placebo rejection rate at 5%")
ax.set_title("Over-rejection grows with serial correlation")
ax.set_ylim(0, max(sweep["naive"].max() * 1.2, 0.1))
ax.legend()
fig.tight_layout(); fig.savefig("nb5.5-rho-sweep.png", dpi=110)
plt.close(fig)
print("saved nb5.5-rho-sweep.png")What to carry forward
Three things from this notebook will keep paying off.
The placebo design measures the size of a test directly. Because the fake law has a true effect of exactly zero, every rejection is a false positive with no ambiguity. Counting them over many draws turns the abstract Week-1 idea of “size” into a number you read off your own output. You did not assume the naive SE was broken — you caught it red-handed.
Serial correlation, ignored, makes DiD over-confident — and the longer/more-persistent the panel, the worse. At \rho = 0.8 the advertised 5% test fired many times too often on pure noise; the \rho-sweep showed the damage scaling smoothly with persistence and vanishing at \rho = 0. The classical formula over-counts independent information, shrinks the SE, inflates the t-stat.
The fixes work, and they are the ones the field adopted. Clustering by state — the unit at which treatment is assigned — let each state’s residuals be arbitrarily correlated across its years and pulled the rejection rate back to roughly 5% (this is why Priya clustered by state in Week 4). Collapsing to pre/post did the same by destroying the within-period time structure, at a cost in power. Both are one line of code away.
Your Turn
You have rebuilt BDM’s indictment and seen the cures work. Now drive it yourself.
Sweep the panel length T. BDM show the over-rejection gets worse with more years. Hold \rho = 0.8 fixed and rerun the naive vs clustered rejection rates for
n_yearsin{5, 10, 20, 40}(pass it throughone_replication). Plot rejection rate against T. Does the naive rate climb with the panel length? Explain why a longer panel gives serial correlation more room to do damage.Rediscover the few-clusters caveat. The chapter’s §6 warns that clustering itself over-rejects when there are too few clusters (the 30–50 rule of thumb; recall Priya’s single treated state). Shrink
n_statesto{6, 10, 20, 50}at \rho = 0.8 and watch the clustered rejection rate. Below ~20 states, does the clustered fix start creeping back above 5%? (Stretch: read up on the Cameron–Gelbach–Miller (2008) wild cluster bootstrap, the fix for this fix.)Add a block bootstrap. Implement BDM’s third cure: resample whole states (each state’s entire time series, kept intact) with replacement, refit the TWFE DiD on each resample, and build the bootstrap distribution of \hat\beta to get a bootstrap p-value. Confirm it, too, restores honest size — because resampling whole states preserves each state’s serial-correlation structure instead of shattering it.
Put a real effect back in. Add a genuine treatment effect (e.g.
Y += 0.3 * D) and compare the power — the fraction of runs that correctly detect it — of the clustered TWFE versus the collapsed two-period design. You should see the collapse buy honest size at a real cost in power, exactly the size/power trade-off from Week 1 made concrete.