Unit and time fixed effects
With many units and periods, the 2×2 regression generalises to two-way fixed effects:
Y_it = α_i + λ_t + δ·D_it + ε_it
- α_i (unit fixed effects) absorb everything constant about a municipality: geography, industry mix.
- λ_t (time fixed effects) absorb every shock common to all units in a period: recessions, national policy.
- δ is identified from units whose treatment status changes, relative to those whose does not.
pyfixest (a Python port of R's fixest) estimates such models fast, with formula syntax
y ~ x | fe1 + fe2.
!pip install -q pyfixest numpy pandas matplotlib
import numpy as np, pandas as pd
import pyfixest as pf
import matplotlib.pyplot as plt
rng = np.random.default_rng(42)
TWFE with pyfixest
A minimum wage increase in one province from month 12, observed for 24 months in 60 municipalities. The effect builds up: −0.02 in the first months, reaching −0.05 after six months.
U, T, T0 = 60, 24, 12
muni = np.repeat(np.arange(U), T); month = np.tile(np.arange(T), U)
treated = (muni >= 30).astype(int)
rel = np.where(treated == 1, month - T0, -99) # months since the policy (treated only)
effect = np.where((treated == 1) & (rel >= 0), -0.02 - 0.03 * np.minimum(rel, 6) / 6, 0.0)
y = (rng.normal(0.6, 0.05, U)[muni] + 0.03 * np.sin(month / 3) + effect + rng.normal(0, 0.015, len(muni)))
df = pd.DataFrame({"muni": muni, "month": month, "treated": treated, "rel": rel,
"D": ((treated == 1) & (month >= T0)).astype(int), "employment": y})
twfe = pf.feols("employment ~ D | muni + month", data=df, vcov={"CRV1": "muni"})
twfe.summary()
print(f"True average effect in the post period: {effect[df.D == 1].mean():.4f}")
Event studies: effects by time since treatment
A single δ averages over all post-treatment periods. An event study replaces D with a set of dummies for each period relative to treatment (leads before, lags after), omitting one reference period (usually −1):
Y_it = α_i + λ_t + Σ_k β_k·1[t − T0 = k] + ε_it
- β_k for k < 0 (leads) should be near zero if pre-trends are parallel.
- β_k for k ≥ 0 trace out how the effect evolves.
df["rel_binned"] = df.rel.clip(-8, 10) # bin distant periods
df.loc[df.treated == 0, "rel_binned"] = -1 # never-treated sit in the reference period
es = pf.feols("employment ~ i(rel_binned, ref=-1) | muni + month", data=df, vcov={"CRV1": "muni"})
coefs = es.tidy()
print(coefs[["Estimate", "Std. Error"]].round(4).head(10))
Reading and plotting the event study
tab = es.tidy().reset_index()
tab["k"] = tab["Coefficient"].str.extract(r"\[(?:T\.)?(-?\d+(?:\.\d+)?)\]")[0].astype(float)
tab = pd.concat([tab, pd.DataFrame({"k": [-1.0], "Estimate": [0.0], "Std. Error": [0.0]})]).sort_values("k")
true_path = {k: (-0.02 - 0.03 * min(k, 6) / 6 if k >= 0 else 0.0) for k in tab.k}
plt.figure(figsize=(10, 5))
plt.errorbar(tab.k, tab.Estimate, yerr=1.96 * tab["Std. Error"], fmt="o", capsize=3, label="Estimate ± 95% CI")
plt.plot(list(true_path), list(true_path.values()), "r--", label="True effect")
plt.axhline(0, color="grey"); plt.axvline(-0.5, color="black", ls=":")
plt.xlabel("Months relative to the policy"); plt.ylabel("Effect on employment rate")
plt.title("Event study: leads are flat, lags show the effect building up"); plt.legend(); plt.show()
leads = [c for c in es.coef().index if "-" in c and "-1]" not in c]
print("Joint test that all leads are zero:")
print(es.wald_test(R=np.eye(len(es.coef()))[[list(es.coef().index).index(c) for c in leads]]))
A flat pre-period supports parallel trends; a joint test of the leads makes the check formal. Two caveats (Roth, 2022): pre-trend tests have low power, so "not significant" is weak reassurance; and conditioning your analysis on passing the test can itself bias the results. Report the plot and think about the size of plausible violations (the HonestDiD approach of Rambachan and Roth).
Add an anticipation effect: employment in the treated province starts falling two months before the policy (firms adjust in advance). Re-run the event study. What do the leads show, and how should you redefine the reference period?
Show a solution
y_anti = df.employment + np.where((df.treated == 1) & df.rel.between(-2, -1), -0.02, 0)
es_a = pf.feols("employment ~ i(rel_binned, ref=-3) | muni + month", data=df.assign(employment=y_anti),
vcov={"CRV1": "muni"})
print(es_a.tidy()[["Estimate", "Std. Error"]].round(4).head(10))
# With ref=-1 the anticipation contaminates the baseline. Move the reference to before anticipation
# starts (here -3) and treat periods -2 and -1 as part of the effect.
Checklist and recap
- I can estimate TWFE with
pyfixestand cluster by unit. - I can build and plot an event study with a reference period.
- I read leads as evidence on pre-trends, with their limits in mind.