Difference-in-Differences & Event Studies

Two-way fixed effects with leanfe

Overview

This tutorial shows how to use leanfe for:

Two-way fixed effects (TWFE) DiD estimation
Event study regressions with lead/lag indicators
Clustered standard errors

Example Data

Panel of 100 states over 20 years, with 40 states treated starting in 2015:

import numpy as np
import polars as pl
from leanfe import leanfe

np.random.seed(42)

# Setup
n_states = 100
years = list(range(2005, 2025))
n_years = len(years)
treatment_year = 2015
treated_states = set(range(1, 41))
true_effect = 3.0

# Generate fixed effects
state_fes = {s: np.random.normal(0, 5) for s in range(1, n_states + 1)}
year_fes = {y: 0.2 * (y - 2005) for y in years}

# Generate panel
data = []
for state in range(1, n_states + 1):
    is_treated = state in treated_states
    for year in years:
        event_time = year - treatment_year if is_treated else None
        post = 1 if year >= treatment_year else 0
        treated_post = float(is_treated and post)
        y = (state_fes[state] + year_fes[year] + 
             true_effect * treated_post + np.random.normal(0, 1))
        data.append({
            "state_id": state, "year": year, "event_time": event_time,
            "treated": float(is_treated), "treated_post": treated_post, "y": y
        })

df = pl.DataFrame(data)
print(f"Panel: {df.shape[0]:,} obs ({n_states} states × {n_years} years)")
print(f"Treated: {len(treated_states)} states starting {treatment_year}")

Panel: 2,000 obs (100 states × 20 years)
Treated: 40 states starting 2015

TWFE Estimation

result = leanfe(
    formula="y ~ treated_post | state_id + year",
    data=df,
    vcov="cluster",
    cluster_cols=["state_id"]
)

print(f"True effect: {true_effect}")
print(f"Estimated:   {result['coefficients']['treated_post']:.4f}")
print(f"SE:          {result['std_errors']['treated_post']:.4f}")
print(f"Clusters:    {result['n_clusters']}")

True effect: 3.0
Estimated:   3.1737
SE:          0.1239
Clusters:    100

Event Study

Create indicators for each period relative to treatment:

# Event time dummies (-5 to +9, omitting t=-1)
event_window = list(range(-5, 10))
reference_period = -1

for et in event_window:
    if et == reference_period:
        continue
    col_name = f"et_m{abs(et)}" if et < 0 else f"et_p{et}"
    df = df.with_columns(
        pl.when(pl.col("event_time") == et).then(1.0).otherwise(0.0).alias(col_name)
    )

et_vars = [f"et_m{abs(et)}" if et < 0 else f"et_p{et}" 
           for et in event_window if et != reference_period]
formula = f"y ~ {' + '.join(et_vars)} | state_id + year"

result_es = leanfe(formula=formula, data=df, vcov="cluster", cluster_cols=["state_id"])

# Extract coefficients
coefs = []
for et in event_window:
    if et == reference_period:
        coefs.append({"event_time": et, "coef": 0.0, "se": 0.0})
    else:
        col_name = f"et_m{abs(et)}" if et < 0 else f"et_p{et}"
        coefs.append({
            "event_time": et,
            "coef": result_es['coefficients'][col_name],
            "se": result_es['std_errors'][col_name]
        })

coef_df = pl.DataFrame(coefs).sort("event_time")

Event Study Plot

import matplotlib.pyplot as plt

fig, ax = plt.subplots(figsize=(10, 6))

et = coef_df["event_time"].to_numpy()
coef = coef_df["coef"].to_numpy()
se = coef_df["se"].to_numpy()

# Plot with confidence bands
ax.fill_between(et, coef - 1.96*se, coef + 1.96*se, alpha=0.2, color='steelblue')
ax.plot(et, coef, 'o-', color='steelblue', markersize=8, linewidth=2, label='Estimate ± 95% CI')

# Reference lines
ax.axhline(y=0, color='gray', linestyle='-', linewidth=0.8)
ax.axhline(y=true_effect, color='darkgreen', linestyle='--', linewidth=1.5, 
           alpha=0.7, label=f'True effect = {true_effect}')
ax.axvline(x=-0.5, color='red', linestyle='--', linewidth=1.5, alpha=0.7)

ax.annotate('Treatment', xy=(-0.5, true_effect + 0.3), fontsize=10, ha='center', color='red')

ax.set_xlabel('Years relative to treatment', fontsize=12)
ax.set_ylabel('Coefficient (relative to t = -1)', fontsize=12)
ax.set_title('Event Study', fontsize=14)
ax.set_xticks(event_window)
ax.set_xlim(-6, 10)
ax.legend(loc='upper left')
ax.grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

Comparing Standard Errors

result_iid = leanfe(data=df, formula="y ~ treated_post | state_id + year", vcov="iid")
result_cluster = leanfe(data=df, formula="y ~ treated_post | state_id + year",
                        vcov="cluster", cluster_cols=["state_id"])

print("Standard Error Comparison:")
print(f"  IID:       {result_iid['std_errors']['treated_post']:.4f}")
print(f"  Clustered: {result_cluster['std_errors']['treated_post']:.4f}")

Standard Error Comparison:
  IID:       0.0926
  Clustered: 0.1239

Quick Reference

Task	Formula
Basic DiD	`y ~ treated_post \\| unit + time`
Event study	`y ~ et_m5 + et_m4 + ... + et_p9 \\| unit + time`
With covariates	`y ~ treated_post + x1 + x2 \\| unit + time`