nb4.1 — Event Studies & Parallel Trends: Playing God with Priya’s Climate-Insurance Shock

Week 4, Chapter 4.1. Companion notebook.

In Chapter 4.1 Priya had a clean natural experiment: on January 1 of one year, one state changed how home insurers may price wildfire risk, while every other state left its rules alone. She wants the regulator’s number — did the rule change premiums, and by how much? — and the chapter showed that the honest answer comes from difference-in-differences (DiD): take the treated state’s before-and-after change, subtract the control states’ before-and-after change, and what survives is the regulation’s effect, under the parallel-trends assumption.

This notebook makes that machinery concrete and, because it is a simulation, lets us play God. We will generate state-year premiums with a known true treatment effect and a tunable pre-existing differential trend, so we can always compare what DiD recovers against the truth it is chasing. We will:

  1. Build the panel — one treated state, several controls, unit fixed effects (\alpha_i), a common time shock (\lambda_t), and a known effect baked in. Then compute the four 2\times2 cell means.
  2. Estimate the 2×2 two ways — the double-difference of means by hand, and the TWFE regression Y_{it}=\alpha_i+\lambda_t+\beta D_{it}+\varepsilon_{it} — and confirm they agree to floating point.
  3. Run the event study — replace the single dummy with leads-and-lags, normalize k=-1 to zero, and make the canonical event-study plot of \hat\beta_k with confidence bands.
  4. Break parallel trends on purpose — inject a differential pre-trend in the treated state, watch the leads tilt away from zero and the estimated effect inflate, and discuss why the pre-trend test has low power.
  5. Stress the inference — compare classical, HC1, and state-clustered standard errors, and confront the few-treated-clusters caveat of Bertrand–Duflo–Mullainathan (2004): Priya has exactly one treated state.

We reproduce the chapter’s illustrative ATT of $70 throughout.

import matplotlib
matplotlib.use("Agg")  # headless: render to buffers, never to a screen

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.formula.api as smf
import pyfixest as pf

rng = np.random.default_rng(42)  # one generator, pinned, drives the whole notebook

plt.rcParams["figure.figsize"] = (9, 5)
plt.rcParams["axes.grid"] = True
plt.rcParams["grid.alpha"] = 0.3

print("numpy      ", np.__version__)
print("pandas     ", pd.__version__)
print("pyfixest   ", pf.__version__)

1. Build the panel: one treated state, several controls, a known effect

We generate a balanced panel of states observed once per year. Each state i carries a permanent level \alpha_i — its fixed fire geography, building stock, baseline price — which is exactly what the unit fixed effect will absorb. Every year carries a common shock \lambda_t (nationwide reinsurance costs, inflation) that hits all states equally; the time fixed effect absorbs that. The treatment switches on for the one treated state in the post period, adding a known effect.

The data-generating process (DGP) is the TWFE model written forward:

Y_{it} = \underbrace{\alpha_i}_{\text{state level}} + \underbrace{\lambda_t}_{\text{common trend}} + \underbrace{\beta\,D_{it}}_{\text{true effect}} + \underbrace{\delta\,(t-t^*)\cdot\text{treated}_i}_{\text{differential pre-trend, OFF for now}} + \varepsilon_{it}.

We set \beta = \$70 (the chapter’s ATT), the treated state’s level to $1,500 and the controls’ around $1,150 (so the treated state is permanently more expensive — the fixed difference DiD must remove), and a common upward drift of $30/year. The pre-trend slope \delta is 0 for now; we turn it on in Section 4. Noise is small so the single treated state’s wiggles do not drown the signal.

Why one treated state? Because that is Priya’s actual design, and Card–Krueger’s: one treated jurisdiction (New Jersey), the rest controls. It is wonderfully transparent — and, as Section 5 shows, genuinely hard to do inference on. We keep it honest.

YEARS   = list(range(2014, 2024))   # 10 years
T_STAR  = 2019                       # treatment switches on in 2019 (k = 0)
ATT_TRUE = 70.0                      # the chapter's true treatment effect, in dollars
N_CTRL  = 8                          # number of control states (state 0 is the single treated one)
SIGMA   = 2.5                        # idiosyncratic noise sd (small: one treated state is noisy enough)

def build_panel(pretrend_slope, sigma=SIGMA, rho=0.0, n_ctrl=N_CTRL, seed=42):
    """Balanced state-year panel generated from the TWFE DGP.

    State 0 is the SINGLE treated state (permanently pricier, level 1500).
    Controls have staggered fixed levels around 1150 (the alpha_i the unit FE will absorb).
    A common +30/yr drift is the lambda_t the time FE will absorb.
    pretrend_slope = delta: a differential LINEAR trend present ONLY in the treated state,
                     active in every period (so it shows up in leads AND lags). 0 = clean.
    rho = AR(1) coefficient on each state's error: e_t = rho*e_{t-1} + innovation. rho>0 makes
          residuals SERIALLY CORRELATED within a state -- the disease BDM (2004) warned about.
    Returns a fresh DataFrame; D = 1 only for the treated state in post years (k >= 0).
    """
    gen = np.random.default_rng(seed)            # local generator -> reproducible per call
    rows = []
    for s in range(1 + n_ctrl):                  # state 0 treated; 1..n_ctrl are controls
        treated = (s == 0)
        level = 1500.0 if treated else (1150.0 + 60.0 * (s % 5))   # fixed unit level alpha_i
        e = 0.0
        errors = []                              # build the state's AR(1) error path
        for _ in YEARS:
            e = rho * e + gen.normal(0.0, sigma)
            errors.append(e)
        for j, y in enumerate(YEARS):
            k = y - T_STAR                       # event time relative to treatment
            D = int(treated and k >= 0)          # treatment indicator
            common = 30.0 * (y - YEARS[0])       # lambda_t: common upward drift
            pretrend = pretrend_slope * k if treated else 0.0      # differential trend (treated only)
            premium = level + common + ATT_TRUE * D + pretrend + errors[j]
            rows.append({"state": s, "year": y, "treated": int(treated),
                         "post": int(k >= 0), "D": D, "k": k, "premium": premium})
    return pd.DataFrame(rows)

panel = build_panel(pretrend_slope=0.0)   # clean panel: no differential pre-trend
print(f"{panel.shape[0]} rows = {1 + N_CTRL} states x {len(YEARS)} years")
print(f"Treated state: 0   |   control states: 1..{N_CTRL}   |   t* = {T_STAR} (k=0)")
print(panel.head(3).to_string(index=False))

2. The 2×2: four cell means and the double difference

Collapse the panel to the chapter’s 2\times2 skeleton: two groups (treated / control) by two periods (before / after). In each cell we take the average premium. The four numbers are all we need for the double difference

\widehat{\text{DiD}} =\big(\bar Y_T^{\,\text{post}}-\bar Y_T^{\,\text{pre}}\big)-\big(\bar Y_C^{\,\text{post}}-\bar Y_C^{\,\text{pre}}\big) =\big(\bar Y_T^{\,\text{post}}-\bar Y_C^{\,\text{post}}\big)-\big(\bar Y_T^{\,\text{pre}}-\bar Y_C^{\,\text{pre}}\big),

which the chapter showed can be read two equivalent ways: difference over time, then across groups or difference across groups, then over time. We compute both orders and confirm they land on the same number — and that it is close to the true $70.

def cell_mean(df, tr, po):
    """Average premium in one 2x2 cell. tr = treated flag, po = post flag."""
    sub = df[(df["treated"] == tr) & (df["post"] == po)]
    return sub["premium"].mean()

YT_pre  = cell_mean(panel, tr=1, po=0)
YT_post = cell_mean(panel, tr=1, po=1)
YC_pre  = cell_mean(panel, tr=0, po=0)
YC_post = cell_mean(panel, tr=0, po=1)

print("                 Before        After")
print(f"  Treated     {YT_pre:9.2f}   {YT_post:9.2f}")
print(f"  Control     {YC_pre:9.2f}   {YC_post:9.2f}")

# Order 1: difference over time first (within group), then across groups.
treated_change = YT_post - YT_pre
control_change = YC_post - YC_pre
did_order1 = treated_change - control_change

# Order 2: difference across groups first (within period), then over time.
gap_after  = YT_post - YC_post
gap_before = YT_pre  - YC_pre
did_order2 = gap_after - gap_before

print("-" * 40)
print(f"  treated change (post-pre)   : {treated_change:7.2f}")
print(f"  control change (post-pre)   : {control_change:7.2f}")
print(f"  DiD  (order 1: time, groups): {did_order1:7.2f}")
print(f"  DiD  (order 2: groups, time): {did_order2:7.2f}")
print(f"  true ATT                    : {ATT_TRUE:7.2f}")

assert abs(did_order1 - did_order2) < 1e-9, "the double difference must be symmetric in its two orders"
print("\nCHECK PASSED: both orderings of the double difference agree exactly.")

2.1 The same number as an imputed counterfactual

The double difference is not just arithmetic; it is DiD’s guess for the missing counterfactual. Recall from §4.1.2 that DiD imputes the treated state’s untreated post outcome by starting from its own pre level and adding the control group’s change:

\widehat{\mathbb{E}[Y_{it}(0)\mid \text{treated, post}]} = \bar Y_T^{\,\text{pre}} + \big(\bar Y_C^{\,\text{post}} - \bar Y_C^{\,\text{pre}}\big), \qquad \widehat{\text{ATT}} = \bar Y_T^{\,\text{post}} - \widehat{\mathbb{E}[Y_{it}(0)\mid \text{treated, post}]}.

We compute that imputed counterfactual and confirm the ATT it produces equals the double difference — the chapter’s two costumes for one calculation.

imputed_Y0_post = YT_pre + control_change          # treated pre level lifted by control's trend
att_from_cf     = YT_post - imputed_Y0_post        # actual minus counterfactual

print(f"Treated actual post premium             : {YT_post:8.2f}")
print(f"Imputed counterfactual  E[Y(0)|T,post]  : {imputed_Y0_post:8.2f}  (pre level + control trend)")
print(f"ATT = actual - counterfactual           : {att_from_cf:8.2f}")
print(f"Double difference from Section 2        : {did_order1:8.2f}")

assert abs(att_from_cf - did_order1) < 1e-9, "imputed-counterfactual ATT must equal the double difference"
print("\nCHECK PASSED: the counterfactual story and the double difference are the same number.")

3. The TWFE regression: DiD as one coefficient \beta

Now the regression form the rest of Week 4 builds on. The chapter’s boxed equation is

Y_{it} = \alpha_i + \lambda_t + \beta\,D_{it} + \varepsilon_{it},

with \alpha_i a unit (state) fixed effect and \lambda_t a time (year) fixed effect. In the 2\times2 design the coefficient \beta on the treatment indicator D_{it} is exactly the double difference. We estimate it two ways and confirm both match the by-hand $70:

  • pyfixest.feols — the modern fixed-effects workhorse. We write premium ~ D | state + year: everything after the | is swept out as fixed effects (the regression encoding of the two differencing operations). We cluster by state for honest inference (more on that in Section 5).
  • statsmodels with dummies — premium ~ D + C(state) + C(year), the same model with explicit category dummies. Slower and clumsier with many units, but it makes the FE machinery visible.

pyfixest vs. statsmodels. Both return the identical \hat\beta here — fixed effects are fixed effects. pyfixest is far more convenient for panels: | state + year absorbs the dummies (no giant design matrix), vcov={"CRV1": "state"} clusters in one argument, and .tidy() / .confint() hand back coefficients and CIs ready to plot. statsmodels is the transparent fallback when you want to see the dummies or when pyfixest is unavailable.

# (a) pyfixest: absorb state + year FE; cluster SEs by state.
m_pf = pf.feols("premium ~ D | state + year", data=panel, vcov={"CRV1": "state"})
beta_pf = float(m_pf.coef().iloc[0])

# (b) statsmodels: the same model with explicit dummies (C(...) = categorical).
m_sm = smf.ols("premium ~ D + C(state) + C(year)", data=panel).fit()
beta_sm = float(m_sm.params["D"])

print(f"by-hand double difference : {did_order1:8.4f}")
print(f"pyfixest    beta_hat(D)   : {beta_pf:8.4f}")
print(f"statsmodels beta_hat(D)   : {beta_sm:8.4f}")
print(f"true ATT                  : {ATT_TRUE:8.4f}")

# The 2x2 DiD == TWFE beta is an algebraic identity on a BALANCED panel: equal to float precision.
assert abs(beta_pf - did_order1) < 1e-7, "pyfixest TWFE beta must equal the double difference"
assert abs(beta_sm - did_order1) < 1e-7, "statsmodels TWFE beta must equal the double difference"
assert abs(beta_pf - beta_sm)   < 1e-7, "pyfixest and statsmodels must agree on beta"
print("\nCHECK PASSED: double difference == pyfixest beta == statsmodels beta (to 1e-7).")

4. The event study: leads, lags, and the canonical plot

The 2\times2 blurs all of “before” into one number and all of “after” into another. The event study un-blurs them by replacing the single D_{it} with one dummy per event-time period k=t-t^*:

Y_{it} = \alpha_i + \lambda_t + \sum_{k\neq -1}\beta_k\,\mathbb{1}\{t - t_i^* = k\} + \varepsilon_{it}.

Two rules make this work and make the plot readable:

  1. Build the event-time dummies only for the treated state. A control state has no t^*, so all its event-time indicators are 0 — it sits in the baseline, supplying the common-trend comparison.
  2. Omit k=-1 (the period just before treatment). One event-time dummy must be dropped or the model is collinear with the fixed effects (the dummy trap); by convention we drop k=-1, which pins \beta_{-1}=0 and makes every other \hat\beta_k read as “how far the treated-control gap has moved away from where it stood the instant before treatment.”

The leads (k<0) are the pre-trend test; the lags (k\ge 0) are the dynamic treatment effect.

def event_study(df, drop_k=-1, cluster="state"):
    """Estimate the leads-and-lags TWFE. Returns a tidy DataFrame indexed by event time k.

    Builds one indicator per event time (treated state only), omitting drop_k as the baseline.
    Uses pyfixest with FE state + year and clustered SEs; reads off Estimate and the 95% CI.
    """
    d = df.copy()                                   # never mutate the caller's frame
    ks = sorted(kk for kk in d["k"].unique() if kk != drop_k)
    cols = []
    for kk in ks:
        col = f"e_{'m' if kk < 0 else 'p'}{abs(kk)}"   # e_m3 = k=-3, e_p2 = k=+2
        d[col] = ((d["treated"] == 1) & (d["k"] == kk)).astype(int)
        cols.append((kk, col))
    fml = "premium ~ " + " + ".join(c for _, c in cols) + " | state + year"
    m = pf.feols(fml, data=d, vcov={"CRV1": cluster})
    tidy = m.tidy()                                  # index = coef name, has Estimate, CIs
    out = []
    for kk, col in cols:
        r = tidy.loc[col]
        out.append({"k": kk, "coef": r["Estimate"], "lo": r["2.5%"], "hi": r["97.5%"]})
    # add the pinned-to-zero baseline so the plot shows the normalization explicitly
    out.append({"k": drop_k, "coef": 0.0, "lo": 0.0, "hi": 0.0})
    return pd.DataFrame(out).sort_values("k").reset_index(drop=True)

es_clean = event_study(panel)
print("Event-study coefficients (clean panel, no injected pre-trend):")
print(es_clean.to_string(index=False, float_format=lambda v: f"{v:7.2f}"))

4.1 Reading the canonical event-study plot

Plot \hat\beta_k against event time, with 95% confidence bands and the k=-1 baseline pinned at zero. This is the DiD figure. We expect, on this clean panel:

  • Leads (k<0): small — an order of magnitude below the $70 effect — and showing no systematic slope, consistent with the treated and control states moving together before the rule. (They are not all exactly zero: a single treated state has its own year-to-year wiggle, and — as Section 6 stresses — with one treated cluster the confidence bands are not fully trustworthy. We read the pattern, not each point.)
  • Lags (k\ge0): jumping up to roughly the true $70 at k=0 and staying there — an immediate, permanent effect, which is what our DGP put in.
def plot_event_study(es, title, true_att=ATT_TRUE):
    """Canonical event-study figure: coefficients +/- 95% CI vs event time."""
    fig, ax = plt.subplots()
    pre  = es[es["k"] < 0]
    post = es[es["k"] >= 0]
    yerr_pre  = [pre["coef"] - pre["lo"],  pre["hi"] - pre["coef"]]
    yerr_post = [post["coef"] - post["lo"], post["hi"] - post["coef"]]
    ax.errorbar(pre["k"],  pre["coef"],  yerr=yerr_pre,  fmt="o", color="steelblue",
                capsize=3, label="leads (pre-treatment)")
    ax.errorbar(post["k"], post["coef"], yerr=yerr_post, fmt="s", color="crimson",
                capsize=3, label="lags (post-treatment)")
    ax.axhline(0.0, color="black", lw=1)
    ax.axvline(-0.5, color="gray", ls="--", lw=1, label="treatment starts (k=0)")
    ax.axhline(true_att, color="green", ls=":", label=f"true ATT = ${true_att:.0f}")
    ax.set_xlabel(r"event time $k = t - t^*$  (k = -1 normalized to 0)")
    ax.set_ylabel(r"$\hat\beta_k$  (treated-vs-control gap, dollars)")
    ax.set_title(title)
    ax.legend(loc="upper left", fontsize=9)
    fig.tight_layout()
    return fig

_ = plot_event_study(es_clean, "Clean panel: flat leads, lags at the true $70")

leads_clean = es_clean[es_clean["k"] < -1]   # genuine pre-periods (exclude pinned k=-1)
lags_clean  = es_clean[es_clean["k"] >= 0]
mean_abs_lead = leads_clean["coef"].abs().mean()
mean_lag      = lags_clean["coef"].mean()
print(f"Mean |lead coefficient| (k < -1): {mean_abs_lead:.2f}  (an order of magnitude below $70)")
print(f"Mean lag coefficient   (k >= 0): {mean_lag:.2f}  (near true $70)")
# Pattern test: leads small relative to the effect, and no strong systematic slope.
assert mean_abs_lead < 0.2 * ATT_TRUE, "clean-panel leads should be small relative to the $70 effect"
assert abs(mean_lag - ATT_TRUE) < 10, "clean-panel lags should sit near the true $70"
print("\nCHECK PASSED: clean leads are small vs the effect; lags recover ~$70.")

6. Inference: cluster by state, and the few-treated-clusters caveat

The point estimate is settled; the chapter’s hard question is the standard error. Bertrand, Duflo & Mullainathan (2004) showed that classical and even heteroskedasticity-robust (HC1) standard errors are far too small for DiD, because a state’s residuals are serially correlated across years and the treatment dummy is highly persistent within a unit. The fix is to cluster by the unit of treatment (here, state), letting each state’s residuals correlate freely across all its years.

We estimate the same \beta three ways and watch the SE — not the point estimate — move. But here we hit Priya’s bind head-on: she has one treated state, and with so few treated clusters the cluster-robust formula is itself unreliable. In a well-powered DiD (30–50 clusters, several treated) the clustered SE comes out larger than classical, correcting BDM’s over-confidence. With one treated cluster it can come out smaller and erratic — there is simply not enough independent treated variation to estimate the cluster correction. Watch which way it goes, and treat the number with suspicion.

# Same TWFE beta, three vcov flavors. pyfixest takes the vcov spec directly.
m_classical = pf.feols("premium ~ D | state + year", data=panel, vcov="iid")     # classical OLS
m_hc1       = pf.feols("premium ~ D | state + year", data=panel, vcov="HC1")     # heteroskedasticity-robust
m_cluster   = pf.feols("premium ~ D | state + year", data=panel, vcov={"CRV1": "state"})  # clustered by state

def se_of(m): return float(m.tidy().loc["D", "Std. Error"])
def b_of(m):  return float(m.coef().iloc[0])

rows = []
for name, m, note in [("Classical (iid)", m_classical, "No - ignores within-state persistence"),
                      ("Robust (HC1)",    m_hc1,       "No - fixes variance, not serial correlation"),
                      ("Clustered by state", m_cluster, "Right structure, but UNRELIABLE: 1 treated cluster")]:
    b, se = b_of(m), se_of(m)
    rows.append({"SE flavor": name, "beta": b, "Std. error": se, "t-stat": b / se, "Honest?": note})
se_table = pd.DataFrame(rows)
print(se_table.to_string(index=False, float_format=lambda v: f"{v:.3f}"))

# The point estimate must NOT move across vcov flavors -- that is the BDM/Ch.2.4 lesson.
assert np.allclose(se_table["beta"], se_table["beta"].iloc[0]), "beta must be identical across vcov flavors"
print("\nCHECK PASSED: beta is identical across all three; only the standard error changes.")
print(f"\nWith ONE treated cluster the clustered SE here ({se_of(m_cluster):.2f}) is actually SMALLER and")
print("erratic, not larger -- the few-treated-clusters pathology. Do NOT read its t-stat as honest;")
print("with many treated clusters it would instead correct the too-small classical SE upward.")
print("The placebo experiment below is the right tool when treated clusters are few.")

6.1 The BDM placebo experiment, and why one treated state breaks inference

To reproduce the spirit of Bertrand–Duflo–Mullainathan, we need two ingredients their experiment had: serially correlated residuals (so classical SEs are genuinely too small) and many states (so a rejection rate is meaningful). We build a fresh panel with 40 control states and an AR(1) error (\rho=0.6, so a state’s premium shock persists year to year) — the exact disease BDM identified.

Now the placebo experiment: assign a fake treatment to a control state that never got the rule (true effect = zero), refit, and check significance. With honest inference a placebo should be flagged significant about 5% of the time. We run this placebo on every control state in turn and tally the rejection rate under classical versus clustered SEs. Watch two things:

  1. Classical SEs over-reject — well above 5% — because they ignore the serial correlation. This is the BDM disease.
  2. Clustering does not save a single-treated-cluster design. Each placebo has exactly one fake-treated state, so the cluster correction has only one treated cluster to work with and is unreliable — it can over-reject even more. This is Priya’s bind made quantitative: you cannot cluster your way out of one treated unit.

So what can Priya do? Permutation / placebo inference. The spread of the fake (true-zero) effects is an empirical null distribution; ask how unusual her real $70 looks against it. That works even with one treated cluster, because it never needs a trustworthy analytic SE.

N_BDM = 40                              # many control states: a rejection RATE becomes meaningful
RHO   = 0.6                             # AR(1) serial correlation: the BDM disease

# Fresh panel with serial correlation. State 0 is the real treated state; 1..N_BDM are controls.
panel_bdm = build_panel(pretrend_slope=0.0, rho=RHO, n_ctrl=N_BDM)

# Priya's REAL estimate on this panel (one treated state, the true $70 effect).
m_real  = pf.feols("premium ~ D | state + year", data=panel_bdm, vcov="iid")
beta_real = float(m_real.coef().iloc[0])

def placebo_run(df, fake_treated_state, vcov):
    """Assign a FAKE post-2019 treatment to one control state (true effect = 0), refit TWFE.
    Returns (beta_placebo, pvalue) under the requested vcov. One fake-treated state = one treated cluster."""
    d = df[df["state"] != 0].copy()                 # drop the real treated state
    d["Dfake"] = ((d["state"] == fake_treated_state) & (d["k"] >= 0)).astype(int)
    m = pf.feols("premium ~ Dfake | state + year", data=d, vcov=vcov)
    t = m.tidy().loc["Dfake"]
    return float(t["Estimate"]), float(t["Pr(>|t|)"])

control_states = list(range(1, 1 + N_BDM))
placebo_betas, sig_classical, sig_cluster = [], 0, 0
for s in control_states:
    b_c, p_c = placebo_run(panel_bdm, s, vcov="iid")             # classical SE
    _,   p_k = placebo_run(panel_bdm, s, vcov={"CRV1": "state"}) # clustered SE (1 treated cluster)
    placebo_betas.append(b_c)                                    # beta is identical across vcov
    sig_classical += int(p_c < 0.05)
    sig_cluster   += int(p_k < 0.05)

placebo_betas = np.array(placebo_betas)
n = len(control_states)
print(f"Placebo (true effect = 0) on {n} control states, with AR(1) rho={RHO} serial correlation:")
print(f"  classical SE flagged significant : {sig_classical}/{n}  "
      f"({100*sig_classical/n:.0f}%)  <- should be ~5%; over-rejects = the BDM disease")
print(f"  clustered SE flagged significant : {sig_cluster}/{n}  "
      f"({100*sig_cluster/n:.0f}%)  <- 1 treated cluster per placebo: clustering does NOT rescue it")
print(f"\nPlacebo effect spread (the empirical null): mean {placebo_betas.mean():+.2f}, "
      f"sd {placebo_betas.std():.2f}, range [{placebo_betas.min():+.2f}, {placebo_betas.max():+.2f}]")

# Permutation p-value: how often a FAKE effect is at least as big as the REAL one. Needs no analytic SE.
perm_p = float((np.abs(placebo_betas) >= abs(beta_real)).mean())
print(f"Priya's REAL estimate = ${beta_real:.1f}; permutation p-value vs the placebo null = {perm_p:.3f}")

# The BDM signatures we MUST reproduce: classical over-rejects; clustering with 1 treated cluster does not fix it.
assert sig_classical / n > 0.05, "classical SE should over-reject placebos under serial correlation (BDM)"
assert perm_p < 0.05, "the real $70 effect should be unusual against the placebo null"
print("\nCHECK PASSED: classical SE over-rejects (BDM); permutation inference still flags the real effect.")

# Visualize the placebo distribution with the real estimate marked.
fig, ax = plt.subplots()
ax.hist(placebo_betas, bins=12, color="lightsteelblue", edgecolor="gray", label="placebo betas (true = 0)")
ax.axvline(0.0, color="black", lw=1)
ax.axvline(beta_real, color="crimson", lw=2, label=f"real estimate = ${beta_real:.0f}")
ax.set_xlabel("estimated effect (dollars)")
ax.set_ylabel("count")
ax.set_title("Placebo (permutation) inference: real $70 vs. the spread of fake effects")
ax.legend()
fig.tight_layout()
print("With one treated cluster, lean on this placebo/permutation distribution, not an analytic t-test.")

Your Turn — tune the pre-trend, then add more treated states

You now have a God’s-eye DiD lab. Two extensions, each a one-line change to the functions above.

1. Sweep the pre-trend slope \delta and watch the bias grow. Loop pretrend_slope over, say, [0, 4, 8, 12, 16]. For each, refit the static TWFE and record \hat\beta. Predict before you run: the bias should be roughly linear in \delta, because a differential trend of \delta/year over the post window adds a near-mechanical amount to the treated-control gap. Plot \hat\beta - 70 against \delta — that line is “how much of a trend gets mislabeled as a treatment effect.”

slopes = [0, 4, 8, 12, 16]
betas = []
for d in slopes:
    p = build_panel(pretrend_slope=d)
    m = pf.feols("premium ~ D | state + year", data=p, vcov={"CRV1": "state"})
    betas.append(float(m.coef().iloc[0]))
# plot betas vs slopes; the line should rise steadily above the true 70

2. Vary the number of treated states and re-examine the leads and the placebo. Add a n_treat argument to build_panel so the first n_treat states are treated (all at the same t^* — still a clean, single-treatment-time design, not staggered adoption, which is Chapter 4.2). With more treated units, the clean-panel leads should get flatter (less idiosyncratic single-state noise), the clustered SE should tighten, and the few-treated-clusters problem eases. This is exactly why researchers reach for many treated units — and a preview of why the staggered-adoption setting of the next chapter is both more powerful and more dangerous.

Check your understanding. (a) In Section 5, why did injecting a pre-trend make the lags overshoot $70 and not just tilt the leads? (b) In Section 5.2, you saw the same real pre-trend become undetectable as noise rose — explain in one sentence why a failure to reject a flat pre-trend is weak evidence. (c) In Section 6.1, why does clustering by state not rescue Priya’s inference when she has only one treated state, and what does the placebo distribution give her instead?