comparing estimators for event studies

event_study_comparison.R

rm(list = ls())
library(LalRUtils)
libreq(data.table, fixest, did2s, ggplot2)
theme_set(lal_plot_theme())

# %%
simulate = function() {
  # dimensions
  N = 250
  T = 35
  # factor loadings
  λ_1i = rnorm(N); λ_2i = rnorm(N)
  # factors with stochastic drift
  a_t = rnorm(T)
  for (k in 2:length(a_t)) {
    a_t[k] = a_t[k - 1] + a_t[k]
  }
  f_1t = a_t + 0.1 * 1:T + 3
  f_2t = rnorm(T)
  # time fixed effects with stochastic drift
  a_t = rnorm(T)
  for (k in 2:length(a_t)) {
    a_t[k] = a_t[k - 1] + a_t[k]
  }
  ξ_t = a_t + 0.1 * 1:T + 3

  # pre-treatment periods
  α_i = rnorm(N); ω_i = rnorm(N)
  tr_i = λ_1i + λ_2i + α_i + ω_i
  tr_i = rank(tr_i) / length(tr_i)
  T_0i = fcase(tr_i <= .5, 35,
               .5 < tr_i & tr_i <= .6, 32,
               .6 < tr_i & tr_i <= .7, 29,
               .7 < tr_i & tr_i <= .8, 26,
               .8 < tr_i & tr_i <= .9, 23,
               .9 < tr_i, 20)
  # unit-specific data
  unit = data.table(unit = 1:N, T_0i, α_i, ω_i, λ_1i, λ_2i)

  # time-specific data
  time = data.table( time = 1:T, ξ_t, f_1t, f_2t)
  # panel data
  dat = unit[CJ(unit, time = 1:T), on = .(unit)]
  dat = time[dat, on = .(time)]
  dat[, ε  := rnorm(N * T)][
      , X1 := rnorm(N * T)][
      , X2 := rnorm(N * T)][
      , D  := fifelse(T_0i < time, 1, 0)][
      , δ  := 0.2 * (time - T_0i) + rnorm(N * T)][
      , Y  := 5 +
              δ * D +
              1 * X1 +
              3 * X2 +
              λ_1i * f_1t +
              λ_2i * f_2t +
              α_i +
              ξ_t +
              ε]
  dat[, `Time after treatment` := time - T_0i]
  return(dat)
}
# %%
set.seed(1000)
dat = simulate()
att = dat[, .(Truth = mean(δ)), by = `Time after treatment`]
att[`Time after treatment` < 0, Truth := 0]

# %% ife
library(fect)
mod_ife = fect(
  Y ~ D + X1 + X2,
  data = dat,
  index = c("unit", "time"),
  method = "ife")

tmp = data.table(
  `Time after treatment` = mod_ife$time,
  `LWX: Interacted Fixed Effects` = mod_ife$att)
att = tmp[att, on = .(`Time after treatment`)]

# %% MC
mod_mc = fect(
  Y ~ D + X1 + X2,
  data = dat,
  index = c("unit", "time"),
  method = "mc")

tmp = data.table(
  `Time after treatment` = mod_mc$time,
  `LWX: Matrix completion` = mod_mc$att)
att = tmp[att, on = .(`Time after treatment`)]

# %%
tmp = melt(att, id.vars = "Time after treatment",
           variable.name = "Model", value.name = "ATT")
# %%
dat[, treat := D][
    , first.treat := fifelse(T_0i == 35, 0, T_0i)]
out = event_study(data = dat, yname = "Y", idname = "unit",
    xformla = ~ X1 + X2,
    tname = "time", gname = "first.treat")
# %%
out2 = as.data.table(out)
colnames(out2) = c("Time after treatment", "ATT", "SE", "Model")
outp = rbind(tmp, out2, fill = T)
# %%
ggplot(outp, aes(`Time after treatment`, y = ATT, colour = as.factor(Model))) +
  geom_vline(xintercept = -0.5, alpha = 0.5) + geom_hline(yintercept = 0, alpha = 0.5) +
  geom_line() +
  scale_x_continuous(breaks = seq(-30, 15, by = 5)) +
  facet_wrap(~ Model, nrow = 4)

# %%