Event-study DiD with staggered treatment design

This is another post in my series of attempts to learn the data.table package for R and to get more comfortable with base R graphics. Today, I reimplement the data generating process and one of the models in this excellent blog post by Andrew Baker: How to create relative time indicators. I also strongly recommend watching Andrew’s talk on Difference-in-Differences staggered treatment designs.

For this exercise, we will only use two libraries:

library(data.table)
library(fixest)

As in the original post, the data generating process includes unit and time fixed effects, as well as unit-specific treatment effects.

𝑦_{𝑖 𝑡} = 𝛼_{𝑖} + 𝛼_{𝑡} + 𝜏_{𝑖 𝑡} + 𝜀_{𝑖 𝑡}

𝛼_{𝑖}, 𝛼_{𝑡} \sim 𝑁 (0, 1)

𝜀_{𝑖 𝑡} \sim 𝑁 (0, 0.5)

𝜇_{𝑖 𝑡} \sim 𝑁 (0.3, {0.2}^{2})

and $𝜏_{𝑖 𝑡}$ is equal to the sum of each firm-specific $𝜇_{𝑖 𝑡}$ since the start of the treatment period.

We have 1000 firms split into 4 cohorts which receive treatment in 1986, 1992, 1998, and 2004.

make_data = function() {

  # unit fixed effects
  unit = data.table(
    unit    = 1:1000,
    unit_fe = rnorm(1000, 0, 1),                # unit fe
    state   = sample(1:40, 1000, replace=TRUE), # state id
    mu      = rnorm(1000, .3, .2)               # treatment effect
  )                

  # year fixed effects
  year = data.table(
    year    = 1980:2010,
    year_fe = rnorm(31, 0, 1)
  )

  # treatment groups by state
  treat_taus = data.table(
    state       = sample(1:40, 40),
    cohort_year = sort(rep(c(1986, 1992, 1998, 2004), 10))
  )

  # full interaction of unit X year
  out = CJ(
    unit = 1:1000,
    year = 1980:2010
  )
  out = merge(out, unit, by="unit")
  out = merge(out, year, by="year")
  out = merge(out, treat_taus, by="state")

  # error term, treatment indicator, and treatment effects
  out[, error := rnorm(31000, 0, .5)][
      , treat := fifelse(year >= cohort_year, 1, 0)][
      , tau   := fifelse(treat == 1, mu, 0)]

  # cumulative treatment effects
  setkey(out, unit, year)
  out[, tau_cum := cumsum(tau), by="unit"]

  # dependent variable
  out[, dep_var := unit_fe + year_fe + tau_cum + error]

  return(out)

}

Plot data from a single simulation

dat = make_data()

# empty plot window
plot(0,
     type = "n",
     ylim = c(-8, 20),
     xlim = c(1980, 2010),
     xlab = "",
     ylab = "Outcome variable")

# one line per firm
for (i in 1:1000) {
  lines(x   = dat[unit==i]$year,
        y   = dat[unit==i]$dep_var,
        col = adjustcolor("grey", alpha=.1))
}

# means by treatment cohort
cohorts = dat[, .(dep_var = mean(dep_var)), by=.(cohort_year, year)]

years <- c(
  "#CC6677" = 1986,
  "#332288" = 1992,
  "#DDCC77" = 1998,
  "#117733" = 2004
) 

for (y in years) {
  abline(v   = y,
         col = names(years)[years==y],
         lwd = 3)
  lines(x   = cohorts[cohort_year==y]$year,
        y   = cohorts[cohort_year==y]$dep_var,
        col = names(years)[years==y],
        lwd = 3)
}

# simulate data
dat = make_data()

# drop 2004 treatment
dat = dat[year <= 2003]

# years since/to treatment
dat[, rel_year := year - cohort_year]
min_year = dat[cohort_year==2004, min(rel_year)]
max_year = dat[cohort_year==2004, max(rel_year)]

# identification requires dropping two dummies
# dat[, rel_year := fifelse(!rel_year %in% c(min_year, -1), as.character(rel_year), "omitted")][
    # , rel_year := factor(rel_year)][
    # , rel_year := relevel(rel_year, ref=c("omitted")]
dat[, rel_year := factor(rel_year, sort(unique(rel_year)))]

# regression model
f = dep_var ~ rel_year | unit + year
mod = feols(f, data=dat)

# clean results
out = data.table("term" = names(coef(mod)), "estimate" = coef(mod))
out = out[ , term := as.numeric(gsub(".*year", "", term))]
out = out[term %in% -5:5]

Key: <unit, year>
       state  year  unit   unit_fe         mu    year_fe cohort_year
       <int> <int> <int>     <num>      <num>      <num>       <num>
    1:    20  1980     1 0.2470796 0.02171754 -1.2023995        2004
    2:    20  1981     1 0.2470796 0.02171754 -0.7647634        2004
    3:    20  1982     1 0.2470796 0.02171754 -0.2958321        2004
    4:    20  1983     1 0.2470796 0.02171754  0.2016919        2004
    5:    20  1984     1 0.2470796 0.02171754 -0.1343072        2004
   ---                                                              
23996:    36  1999  1000 0.5774815 0.26184108 -0.6427118        1998
23997:    36  2000  1000 0.5774815 0.26184108  0.6224544        1998
23998:    36  2001  1000 0.5774815 0.26184108 -1.3553682        1998
23999:    36  2002  1000 0.5774815 0.26184108 -0.7718564        1998
24000:    36  2003  1000 0.5774815 0.26184108  2.3014549        1998
             error treat       tau   tau_cum     dep_var rel_year
             <num> <num>     <num>     <num>       <num>    <num>
    1: -0.25827408     0 0.0000000 0.0000000 -1.21359402      -24
    2:  0.01859074     0 0.0000000 0.0000000 -0.49909312      -23
    3: -0.82524439     0 0.0000000 0.0000000 -0.87399691      -22
    4: -0.35391091     0 0.0000000 0.0000000  0.09486055      -21
    5: -0.15934624     0 0.0000000 0.0000000 -0.04657391      -20
   ---                                                           
23996: -0.71740800     1 0.2618411 0.5236822 -0.25895615        1
23997: -0.48937463     1 0.2618411 0.7855233  1.49608451        2
23998:  0.36482017     1 0.2618411 1.0473643  0.63429776        3
23999:  0.52880334     1 0.2618411 1.3092054  1.64363379        4
24000: -0.22777875     1 0.2618411 1.5710465  4.22220417        5

Key: <unit, year>
Index: <cohort_year>
       state  year  unit   unit_fe         mu    year_fe cohort_year
       <int> <int> <int>     <num>      <num>      <num>       <num>
    1:    20  1980     1 0.2470796 0.02171754 -1.2023995        2004
    2:    20  1981     1 0.2470796 0.02171754 -0.7647634        2004
    3:    20  1982     1 0.2470796 0.02171754 -0.2958321        2004
    4:    20  1983     1 0.2470796 0.02171754  0.2016919        2004
    5:    20  1984     1 0.2470796 0.02171754 -0.1343072        2004
   ---                                                              
23996:    36  1999  1000 0.5774815 0.26184108 -0.6427118        1998
23997:    36  2000  1000 0.5774815 0.26184108  0.6224544        1998
23998:    36  2001  1000 0.5774815 0.26184108 -1.3553682        1998
23999:    36  2002  1000 0.5774815 0.26184108 -0.7718564        1998
24000:    36  2003  1000 0.5774815 0.26184108  2.3014549        1998
             error treat       tau   tau_cum     dep_var rel_year
             <num> <num>     <num>     <num>       <num>   <fctr>
    1: -0.25827408     0 0.0000000 0.0000000 -1.21359402      -24
    2:  0.01859074     0 0.0000000 0.0000000 -0.49909312      -23
    3: -0.82524439     0 0.0000000 0.0000000 -0.87399691      -22
    4: -0.35391091     0 0.0000000 0.0000000  0.09486055      -21
    5: -0.15934624     0 0.0000000 0.0000000 -0.04657391      -20
   ---                                                           
23996: -0.71740800     1 0.2618411 0.5236822 -0.25895615        1
23997: -0.48937463     1 0.2618411 0.7855233  1.49608451        2
23998:  0.36482017     1 0.2618411 1.0473643  0.63429776        3
23999:  0.52880334     1 0.2618411 1.3092054  1.64363379        4
24000: -0.22777875     1 0.2618411 1.5710465  4.22220417        5

The variables 'rel_year12', 'rel_year13', 'rel_year14', 'rel_year15',
'rel_year16' and 'rel_year17' have been removed because of collinearity (see
$collin.var).

Plot results

# event study estimates
plot(x    = out$term,
     y    = out$estimate, 
     pch=20,
     xlim = c(-5, 5),
     ylim = c(-.5, 1.5),
     xlab = "Relative time",
     ylab = "Estimate")

# truth
truth = data.table(term=-5:5)
truth[, estimate := fifelse(term >= 0, (term + 1) * .3, 0)]
lines(x   = truth$term,
      y   = truth$estimate,
      col = "#CC6677",
      lty = "dashed")

     term estimate
    <int>    <num>
 1:    -5      0.0
 2:    -4      0.0
 3:    -3      0.0
 4:    -2      0.0
 5:    -1      0.0
 6:     0      0.3
 7:     1      0.6
 8:     2      0.9
 9:     3      1.2
10:     4      1.5
11:     5      1.8