Generate the "Baker" simulated dataset for staggered treatment effects

baker.R

set.seed(1234)

# Create the base dataset as a cross-join of 1,000 firms over 30 periods
baker = expand.grid(n = 1:30, id = 1:1000)

# Add additional columns
baker = 
  baker |>
  within({
    year       = n + 1980 - 1
    state      = 1 + (id-1) %/% 25
    firms      = runif(id*year, 0, 5)
    group      = 1 + (state-1) %/% 10
    treat_date = 1980 + group*6
    time_til   = year - treat_date
    treat      = time_til>=0
    e          = rnorm(id*year, 0, 0.5^2)
    te         = rnorm(id*year, 10-2*(group-1), 0.2^2)
    y          = firms + n + treat*te*(year - treat_date + 1) + e 
    y2         = firms + n + te*treat + e
    })

Estimation and plot step using @kylebutts 's excellent did2s:

setDT(baker)
# subset to before last cohort treated; set treatment time for never treated to 0
baker = baker[year <= 2003][, treat_date := ifelse(treat_date == 2004, NA, treat_date)]
# estimation 
library(did2s)
es_estimates = event_study(baker, yname = 'y', idname = 'id', tname = 'year', gname = 'treat_date')
library(ggplot2)
p1 = ggplot(es_estimates, aes(term, estimate)) +
  geom_pointrange(
    aes(ymin = estimate - 1.96 * `std.error`, ymax = estimate + 1.96 * `std.error`, colour = estimator), size = 0.2, alpha = 0.6) +
  geom_hline(yintercept = 0,    linetype = 'dashed') + geom_vline(xintercept = -0.5, linetype = 'dashed') +
  facet_wrap(~ estimator) + theme(legend.pos = "None")

Are the BJS + SA estimators doing something unexpected here? Non-negligible biases (unless i messed something up in the DGP / subsetting).

So, I don't have too much time to dig into this all right now, buttttt

1. I think my Sun and Abraham code is wrong and it's on my todo list to fix for event_study, so I think that might be user-error
2. BJS and did2s will produce identical estimates in the post-period. I'm not sure exactly the reason, but BJS estimates the pre-period in a separate regression:

pre_est <- fixest::feols(pre_formula, data[data$zz000treat == 0, ], weights = weights_vector[data$zz000treat == 0], warn=FALSE, notes=FALSE) 

source: https://github.com/kylebutts/didimputation/blob/0330e81be3267d8c759297b0dd1ebc3480aac7db/R/did_imputation.R#L270

I think that is causing the weird level shift, but can't really place a finger on what the problem is there