Difference-in-Differences Basics

41.1 Two bad comparisons and one good one (40 min)

NorthMart ran a campaign in every Ontario and Quebec store from week 40. Two comparisons come to mind first, and both mislead:

Difference-in-differences (DiD) subtracts one from the other: the treated stores' change, minus the other stores' change. Permanent store differences cancel in each change; common shocks cancel in the difference. What remains is the effect, if the treated stores would have changed like the others without the campaign. That's parallel trends, the assumption the whole design rests on. It's about trends, not levels: Ontario stores can be bigger, as long as they'd have grown the same.

41.2 DiD from four averages (40 min)

Create causalkit/stores.py from the provided file above, then day41/two_by_two.py:

# day41/two_by_two.py   (python -m day41.two_by_two)
"""
NorthMart ran a campaign in every Ontario and Quebec store from week 40.
Three ways to estimate its effect from four averages; only one removes both store and time differences.
"""
import polars as pl

from causalkit.stores import campaign, public, truth

full = campaign()
df = public(full)
print(df.head(3))

means = (df.group_by("treated_group", "post").agg(pl.col("sales").mean().round(2))
         .pivot(on="post", index="treated_group", values="sales").sort("treated_group")
         .select("treated_group", pl.col("0").alias("before (weeks 1-39)"), pl.col("1").alias("after (weeks 40-60)")))
print(means)

m = {(g, p): v for g, p, v in df.group_by("treated_group", "post").agg(pl.col("sales").mean()).iter_rows()}
after_only = m[1, 1] - m[0, 1]
before_after = m[1, 1] - m[1, 0]
did = (m[1, 1] - m[1, 0]) - (m[0, 1] - m[0, 0])
print(f"\ntreated vs control, after only:   {after_only:+.2f}   (mixes in store differences)")
print(f"treated, after vs before:        {before_after:+.2f}   (mixes in what changed for everyone)")
print(f"difference-in-differences:       {did:+.2f}")
print(f"Truth: {truth(full)['ATT']:+.2f} ($ thousands per store-week)")

python -m day41.two_by_two

│ treated_group ┆ before (weeks 1-39) ┆ after (weeks 40-60) │
│ 0             ┆ 100.85              ┆ 101.17              │
│ 1             ┆ 103.54              ┆ 106.6               │

treated vs control, after only:   +5.43   (mixes in store differences)
treated, after vs before:        +3.06   (mixes in what changed for everyone)
difference-in-differences:       +2.74
Truth: +3.00 ($ thousands per store-week)

Sales are in $ thousands per store per week. The after-only comparison nearly doubles the effect; DiD lands within 0.3 of the truth.

41.3 DiD as a fixed-effects regression (45 min)

Concept. The four-average DiD is a regression: sales ~ treated_group * post, where the interaction coefficient is the DiD. The general form, which works for any number of groups and periods, is the two-way fixed effects (TWFE) regression:

sales_it = store effect_i + week effect_t + β · d_it + error

d_it is 1 when store i is treated in week t. Store effects absorb permanent differences; week effects absorb common shocks. We use pyfixest, a fast fixed-effects package with R's fixest syntax: variables after | are absorbed as fixed effects.

41.4 Standard errors: the part most people get wrong (50 min)

Create day41/did_regression.py:

# day41/did_regression.py   (python -m day41.did_regression)
"""
The same DiD as a regression: store fixed effects + week fixed effects + the treatment dummy (TWFE).
Same estimate every time; four different standard errors.
"""
import pyfixest as pf

from causalkit.stores import REGIONS, campaign, public, truth

full = campaign()
df = public(full).to_pandas()
df["region_id"] = df.region.map({r: i for i, r in enumerate(REGIONS)}).astype("int64")

simple = pf.feols("sales ~ treated_group * post", df, vcov={"CRV1": "store_id"})
print(f"2x2 regression (group x post): {simple.coef()['treated_group:post']:.2f}")
for name, vcov in [("i.i.d.", "iid"), ("clustered by store", {"CRV1": "store_id"}),
                   ("clustered by region (8)", {"CRV1": "region_id"})]:
    fit = pf.feols("sales ~ d | store_id + week", df, vcov=vcov)
    print(f"TWFE, {name:<24} {fit.coef()['d']:.2f}  SE {fit.se()['d']:.2f}  p = {fit.pvalue()['d']:.4f}")
wild = fit.wildboottest(param="d", reps=9999, seed=1)          # uses the region clusters of the last fit
print(f"TWFE, wild cluster bootstrap by region:     p = {wild['Pr(>|t|)']:.4f}")
print(f"Truth: {truth(full)['ATT']:.2f}")

python -m day41.did_regression

2x2 regression (group x post): 2.74
TWFE, i.i.d.                   2.74  SE 0.13  p = 0.0000
TWFE, clustered by store       2.74  SE 0.30  p = 0.0000
TWFE, clustered by region (8)  2.74  SE 0.18  p = 0.0000
TWFE, wild cluster bootstrap by region:     p = 0.0156
Truth: 3.00

Same estimate, very different uncertainty. Why?

• A store's sales this week are strongly correlated with its sales last week (serial correlation). The i.i.d. formula treats 14,400 store-weeks as 14,400 independent observations; really there are 240 stores, each giving one noisy before-after change. Its SE is less than half the right size.
• Clustering by store lets each store's errors be correlated over time. It's the minimum for any panel.
• The campaign was assigned by region. If regions have their own shocks, the right cluster is the region, but there are only 8 regions, and clustered SEs with few clusters are too small (notice: 0.18 < 0.30).
• The wild cluster bootstrap (Cameron, Gelbach & Miller, 2008) is the standard fix for few clusters: resample the region-level residuals with random sign flips.

Rule: cluster at the level at which treatment was assigned, and use the wild cluster bootstrap when there are fewer than about 40 clusters.

Try it. Estimate the DiD using only weeks 39 and 40 (one week before, one week after). Create day41/one_week.py. What happens to the estimate and its SE?

# day41/one_week.py   (python -m day41.one_week)
import pyfixest as pf

from causalkit.stores import campaign, public

df = public(campaign()).to_pandas()
for name, d in [("weeks 39 and 40 only", df[df.week.isin([39, 40])]), ("all 60 weeks", df)]:
    fit = pf.feols("sales ~ d | store_id + week", d, vcov={"CRV1": "store_id"})
    print(f"{name:<22} {fit.coef()['d']:.2f}  (SE {fit.se()['d']:.2f})")

weeks 39 and 40 only   3.86  (SE 0.43)
all 60 weeks           2.74  (SE 0.30)

41.5 Day 41 mini-project: fake campaigns (75 min)

Bertrand, Duflo & Mullainathan (2004) showed that most published DiD papers of the time over-rejected badly. Repeat their experiment: give two random regions a fake campaign in data where nothing happened, 200 times, and count how often each method calls it significant. Two worlds: store noise only, and stores that also share regional shocks (region_sd=0.1).

Create day41/clustered_se.py (it takes about 4 minutes):

# day41/clustered_se.py   (python -m day41.clustered_se)
"""
Mini-project (Bertrand, Duflo & Mullainathan, 2004): give two random regions a FAKE campaign in data where
nothing happened, 200 times. A correct 5% test rejects about 5% of the time. Which standard error does?
Two worlds: store noise only, and stores that also share regional shocks.
"""
import numpy as np
import polars as pl
import pyfixest as pf

from causalkit.stores import REGIONS, campaign, public

methods = ["i.i.d.", "clustered by store", "clustered by region", "wild bootstrap by region"]
rows = []
for world, region_sd in [("store noise only", 0.0), ("plus regional shocks", 0.1)]:
    rng = np.random.default_rng(0)
    reject = dict.fromkeys(methods, 0)
    n = 200
    for s in range(n):
        fake = tuple(rng.choice(REGIONS, 2, replace=False))
        df = public(campaign(seed=7000 + s, effect=0.0, treated_regions=fake, region_sd=region_sd)).to_pandas()
        df["region_id"] = df.region.map({r: i for i, r in enumerate(REGIONS)}).astype("int64")
        for name, vcov in zip(methods[:3], ["iid", {"CRV1": "store_id"}, {"CRV1": "region_id"}]):
            fit = pf.feols("sales ~ d | store_id + week", df, vcov=vcov)
            reject[name] += fit.pvalue()["d"] < 0.05
        reject[methods[3]] += fit.wildboottest(param="d", reps=999, seed=s)["Pr(>|t|)"] < 0.05
    rows.append({"world": world} | {m: f"{r / n:.1%}" for m, r in reject.items()})
print("False positive rate of a 5% test (should be 5%):")
with pl.Config(tbl_width_chars=140, fmt_str_lengths=30):
    print(pl.DataFrame(rows))

python -m day41.clustered_se

False positive rate of a 5% test (should be 5%):
│ world                ┆ i.i.d. ┆ clustered by store ┆ clustered by region ┆ wild bootstrap by region │
│ store noise only     ┆ 38.0%  ┆ 8.0%               ┆ 15.0%               ┆ 2.5%                     │
│ plus regional shocks ┆ 57.0%  ┆ 24.0%              ┆ 19.0%               ┆ 2.5%                     │

• i.i.d. standard errors find a "significant" effect of a campaign that never happened 38–57% of the time.
• Store clustering is nearly right when there are no regional shocks (8%), and badly wrong when there are (24%).
• Region clustering with 8 clusters is wrong in both worlds.
• The wild cluster bootstrap is conservative (2.5%) but never over-rejects.

In real data you don't know whether regional shocks exist, so assume they do.

41.6 Save your work with Git (20 min)

git add .
git commit -m "Day 41: store simulator, 2x2 DiD, TWFE, clustered SEs, wild bootstrap"
git push

Day 41 checklist

• I can explain the two bad comparisons DiD fixes, and state parallel trends
• I can compute a 2×2 DiD from four averages and as a TWFE regression with pyfixest
• I cluster standard errors at the level of treatment assignment
• I use the wild cluster bootstrap when there are few clusters
• I can show by simulation which standard errors give a 5% false positive rate