Diff-in-Diff and Event Studies

Two-way fixed effects and dynamic treatment effect estimation

What You'll Learn

Setting up panel data for diff-in-diff estimation
Creating treatment timing and post-treatment variables
Two-way fixed effects (TWFE) regression
Event study specifications with leads and lags
Visualizing treatment effects over time

Background: No-Fault Divorce and Female Suicide

In this tutorial, we replicate a classic diff-in-diff analysis examining how no-fault divorce reform affected female suicide rates across US states.

The Research Question

Did no-fault divorce laws reduce female suicide rates? The hypothesis is that women trapped in unhappy or abusive marriages might experience better mental health outcomes when divorce becomes easier to obtain.

What is No-Fault Divorce?

Before no-fault divorce, ending a marriage required proving that one spouse was at "fault" (adultery, abandonment, abuse, etc.). No-fault divorce allows couples to divorce simply by declaring the marriage "irretrievably broken."

States adopted no-fault divorce at different times (1969-1985)
Some states never adopted it during our sample period
This staggered adoption creates our diff-in-diff design

The Data

We use a balanced panel from 1964 to 1996 covering 49 US states:

asmrs - Age-standardized female suicide mortality rate (outcome Y)
pcinc - Per capita income (control)
asmrh - Age-standardized homicide rate (control)
cases - Number of cases (control)
stfips - State FIPS code
year - Year

Setting Up the Analysis

Load and Explore the Data

* Load the data
use bacon_example_diff_in_diff_review.dta, clear

* Explore the structure
tab year
tab stfips

* Examine the outcome and controls
desc asmrs        // Female suicide rate (Y)
desc pcinc asmrh cases  // Controls

# Load packages
library(haven)
library(tidyverse)
library(fixest)

# Load the data
df <- read_dta("bacon_example_diff_in_diff_review.dta")

# Explore the structure
table(df$year)
table(df$stfips)

# Examine the outcome and controls
summary(df[c("asmrs", "pcinc", "asmrh", "cases")])

import pandas as pd
import pyfixest as pf

# Load the data
df = pd.read_stata("bacon_example_diff_in_diff_review.dta")

# Explore the structure
print(df['year'].value_counts().sort_index())
print(df['stfips'].value_counts().sort_index())

# Examine the outcome and controls
df[['asmrs', 'pcinc', 'asmrh', 'cases']].describe()

Create the Treatment Timing Variable

First, we need to identify when each state adopted no-fault divorce. States that never adopted act as our control group.

* Create variable for year state adopted no-fault divorce
* Missing = never adopted (control state)
gen _nfd = .

* Fill in adoption years for each state
replace _nfd = 1971 if stfips == 1   // Alabama
replace _nfd = 1973 if stfips == 4   // Arizona
replace _nfd = 1970 if stfips == 6   // California
replace _nfd = 1971 if stfips == 8   // Colorado
replace _nfd = 1973 if stfips == 9   // Connecticut
replace _nfd = 1977 if stfips == 11  // DC
replace _nfd = 1971 if stfips == 12  // Florida
replace _nfd = 1973 if stfips == 13  // Georgia
replace _nfd = 1971 if stfips == 16  // Idaho
replace _nfd = 1984 if stfips == 17  // Illinois
replace _nfd = 1973 if stfips == 18  // Indiana
replace _nfd = 1970 if stfips == 19  // Iowa
replace _nfd = 1969 if stfips == 20  // Kansas
replace _nfd = 1972 if stfips == 21  // Kentucky
replace _nfd = 1973 if stfips == 23  // Maine
replace _nfd = 1975 if stfips == 25  // Massachusetts
replace _nfd = 1972 if stfips == 26  // Michigan
replace _nfd = 1974 if stfips == 27  // Minnesota
replace _nfd = 1973 if stfips == 29  // Missouri
replace _nfd = 1975 if stfips == 30  // Montana
replace _nfd = 1972 if stfips == 31  // Nebraska
replace _nfd = 1973 if stfips == 32  // Nevada
replace _nfd = 1971 if stfips == 33  // New Hampshire
replace _nfd = 1971 if stfips == 34  // New Jersey
replace _nfd = 1973 if stfips == 35  // New Mexico
replace _nfd = 1971 if stfips == 38  // North Dakota
replace _nfd = 1974 if stfips == 39  // Ohio
replace _nfd = 1973 if stfips == 41  // Oregon
replace _nfd = 1980 if stfips == 42  // Pennsylvania
replace _nfd = 1976 if stfips == 44  // Rhode Island
replace _nfd = 1969 if stfips == 45  // South Carolina
replace _nfd = 1985 if stfips == 46  // South Dakota
replace _nfd = 1974 if stfips == 48  // Texas
replace _nfd = 1973 if stfips == 53  // Washington
replace _nfd = 1977 if stfips == 55  // Wisconsin
replace _nfd = 1977 if stfips == 56  // Wyoming

* States with missing _nfd never adopted (control group)

# Create adoption year variable
# NA = never adopted (control state)
adoption_years <- c(
  `1` = 1971, `4` = 1973, `6` = 1970, `8` = 1971, `9` = 1973,
  `11` = 1977, `12` = 1971, `13` = 1973, `16` = 1971, `17` = 1984,
  `18` = 1973, `19` = 1970, `20` = 1969, `21` = 1972, `23` = 1973,
  `25` = 1975, `26` = 1972, `27` = 1974, `29` = 1973, `30` = 1975,
  `31` = 1972, `32` = 1973, `33` = 1971, `34` = 1971, `35` = 1973,
  `38` = 1971, `39` = 1974, `41` = 1973, `42` = 1980, `44` = 1976,
  `45` = 1969, `46` = 1985, `48` = 1974, `53` = 1973, `55` = 1977,
  `56` = 1977
)

df <- df %>%
  mutate(nfd = adoption_years[as.character(stfips)])

# States with NA never adopted (control group)

# Create adoption year variable
# NaN = never adopted (control state)
adoption_years = {
    1: 1971, 4: 1973, 6: 1970, 8: 1971, 9: 1973, 11: 1977,
    12: 1971, 13: 1973, 16: 1971, 17: 1984, 18: 1973, 19: 1970,
    20: 1969, 21: 1972, 23: 1973, 25: 1975, 26: 1972, 27: 1974,
    29: 1973, 30: 1975, 31: 1972, 32: 1973, 33: 1971, 34: 1971,
    35: 1973, 38: 1971, 39: 1974, 41: 1973, 42: 1980, 44: 1976,
    45: 1969, 46: 1985, 48: 1974, 53: 1973, 55: 1977, 56: 1977
}

df['nfd'] = df['stfips'].map(adoption_years)

# States with NaN never adopted (control group)

Key Insight: Staggered Adoption

Notice that treatment is staggered - different states adopt at different times. Kansas and South Carolina adopted first (1969), while South Dakota adopted last (1985). This variation is what allows identification.

Two-Way Fixed Effects (TWFE)

Create the Post-Treatment Indicator

For the basic diff-in-diff, we need a binary indicator for whether a state-year is in the post-treatment period:

* Create post-treatment indicator
* = 1 if year >= year state adopted no-fault divorce
* = 0 otherwise (including never-adopters)
gen treat_post = (year >= _nfd)

# Create post-treatment indicator
# = 1 if year >= year state adopted no-fault divorce
# = 0 otherwise (including never-adopters)
df <- df %>%
  mutate(treat_post = ifelse(!is.na(nfd) & year >= nfd, 1, 0))

# Create post-treatment indicator
# = 1 if year >= year state adopted no-fault divorce
# = 0 otherwise (including never-adopters)
df['treat_post'] = ((df['nfd'].notna()) & (df['year'] >= df['nfd'])).astype(int)

Run the TWFE Regression

* Two-way fixed effects regression
* Y = female suicide rate
* X = post-treatment indicator
* FEs = state and year
* Controls = income, homicide rate, cases
reg asmrs treat_post i.stfips i.year pcinc asmrh cases, vce(cluster stfips)

# Two-way fixed effects regression using fixest
# Y = female suicide rate
# X = post-treatment indicator
# FEs = state and year
# Controls = income, homicide rate, cases
model_twfe <- feols(asmrs ~ treat_post + pcinc + asmrh + cases | stfips + year,
                    cluster = ~stfips, data = df)
summary(model_twfe)

# Two-way fixed effects regression using pyfixest
# Y = female suicide rate
# X = post-treatment indicator
# FEs = state and year
# Controls = income, homicide rate, cases
model_twfe = pf.feols('asmrs ~ treat_post + pcinc + asmrh + cases | stfips + year',
                      vcov={'CRV1': 'stfips'}, data=df)
print(model_twfe.summary())

Quick Check: TWFE

Question: In the TWFE regression, what does the coefficient on treat_post capture?

Interpreting the Coefficient

The coefficient on treat_post represents the average effect of no-fault divorce on female suicide rates, controlling for state fixed effects, year fixed effects, and time-varying controls. A negative coefficient suggests the reform reduced suicide rates.

Event Study

The TWFE gives us one average effect, but we often want to see how the effect evolves over time. An event study lets us:

Test the parallel trends assumption (pre-treatment coefficients should be near zero)
See how the effect builds over time after treatment
Identify any anticipation effects before treatment

Create Time-to-Treatment Variable

* Create time relative to treatment
* Negative = years before treatment
* 0 = year of treatment
* Positive = years after treatment
gen time_to_treat = year - _nfd

* For never-treated states, set to 0 (they'll be absorbed by FEs)
replace time_to_treat = 0 if missing(_nfd)

# Create time relative to treatment
# Negative = years before treatment
# 0 = year of treatment
# Positive = years after treatment
df <- df %>%
  mutate(time_to_treat = year - nfd,
         # For never-treated states, set to 0
         time_to_treat = ifelse(is.na(time_to_treat), 0, time_to_treat))

# Create time relative to treatment
# Negative = years before treatment
# 0 = year of treatment
# Positive = years after treatment
df['time_to_treat'] = df['year'] - df['nfd']

# For never-treated states, set to 0
df['time_to_treat'] = df['time_to_treat'].fillna(0).astype(int)

Run the Event Study Regression

* Handle Stata's factor variable limitation (no negatives)
summ time_to_treat
local min_val = r(min)

* Shift so minimum is 0
gen shifted_ttt = time_to_treat - `min_val'

* Find the value that corresponds to t=-1 (our reference period)
summ shifted_ttt if time_to_treat == -1
local true_neg1 = r(mean)

* Event study regression
* ib`true_neg1' sets t=-1 as the reference category
reg asmrs ib`true_neg1'.shifted_ttt pcinc asmrh cases i.stfips i.year, vce(cluster stfips)

# Event study using fixest's i() syntax
# ref = -1 sets t=-1 as the reference category
model_es <- feols(asmrs ~ i(time_to_treat, ref = -1) + pcinc + asmrh + cases | stfips + year,
                  cluster = ~stfips, data = df)
summary(model_es)

# Create event study plot
iplot(model_es,
      xlab = "Time to Treatment",
      main = "Event Study: No-Fault Divorce and Female Suicide")

# Event study using pyfixest's i() syntax
# Create factor variable with -1 as reference
df['time_to_treat_factor'] = pd.Categorical(df['time_to_treat'])

# Event study regression
model_es = pf.feols('asmrs ~ i(time_to_treat, ref=-1) + pcinc + asmrh + cases | stfips + year',
                    vcov={'CRV1': 'stfips'}, data=df)
print(model_es.summary())

# Create event study plot
model_es.iplot()

Quick Check: Event Studies

Question: In an event study plot, what should pre-treatment coefficients (t < 0) look like if parallel trends holds?

Why t = -1 as Reference?

We use t = -1 (one year before treatment) as our reference period because it's the last pre-treatment period. All coefficients are then interpreted relative to this period. This is standard practice in event studies.

Reading an Event Study Plot

Pre-treatment coefficients (left of 0): Should be close to zero and statistically insignificant if parallel trends holds
Post-treatment coefficients (right of 0): Show how the treatment effect evolves over time
Reference period (t = -1): Always zero by construction

Interpreting Results

What to Look For

Parallel Trends

In the event study plot, pre-treatment coefficients should be statistically indistinguishable from zero. If you see a clear trend before treatment, the parallel trends assumption may be violated.

Common Issues

Pre-trends: If coefficients are trending before treatment, your estimate may be biased
Anticipation effects: If coefficients become significant before t=0, states may have responded before official adoption
Heterogeneous treatment effects: With staggered adoption, TWFE can be biased if treatment effects vary over time

Modern DiD Methods

Recent econometrics research (Goodman-Bacon 2021, Callaway & Sant'Anna 2021, Sun & Abraham 2021) has shown that TWFE with staggered treatment can produce biased estimates. For more advanced analysis, consider using:

Stata: csdid, did_imputation, eventstudyinteract
R: did, did2s, fixest::sunab()
Python: pyfixest, differences

Tutorial Complete!

You've learned the essentials of diff-in-diff analysis and event studies. Key takeaways:

TWFE provides an average treatment effect but hides dynamics
Event studies reveal how effects evolve over time
Pre-treatment coefficients test the parallel trends assumption
With staggered treatment, consider modern DiD methods

← Back to Tutorials