Hi, is there a way to simulate data and plot all the regression coefficients from the different DiD commands and plot them on an event plot? I tried to store the results as a matrix and in some cases, I just use the eststo command. However, when I try to use the event_plot command, I keep getting the error estimates not found. It is the coefficients that are stored as matrix that are causing the issue. I am attaching a simulated data set and the codes that I have tried.
Code:
clear all
set seed 10
global T = 15
global I = 400
set obs `=$I*$T'
g i = int((_n-1)/$T )+1
g t = mod((_n-1),$T )+1
xtset i t
// Randomly generate treatment rollout periods uniformly across Ei=10..16
// (note that periods t>=16 would not be useful since all units are treated by then)
*Period when unit is first treated
g Ei = ceil(runiform()*7)+$T -6 if t==1
bys i (t): replace Ei = Ei[1]
*Relative time, i.e. number of periods since treated (could be missing if never treated)
g K = t-Ei
*Treatment indicator
g D = K>=0 & Ei!=.
*Heterogeneous treatment effects
g tau = cond(D==1, (t-12.5), 0)
*Error term
g eps = rnormal()
*Outcome
g Y = i + 3*t + tau*D + eps
********************************************************************************
*DiD_imputation of Borusyak et al. (2021) *
********************************************************************************
did_imputation Y i t Ei, allhorizons pretrends(5)
eststo bjs
********************************************************************************
*did_multiplegt of de Chaisemartin and D'Haultfoeuille (2020) *
********************************************************************************
did_multiplegt Y i t D, robust_dynamic dynamic(5) placebo(5) longdiff_placebo breps(100) cluster(i)
matrix dcdh_b = e(estimates)
matrix dcdh_v = e(variances)
********************************************************************************
*csdid of Callaway and Sant'Anna (2020)
********************************************************************************
g gvar = cond(Ei>15, 0, Ei) // group variable as required for the csdid command
csdid Y, ivar(i) time(t) gvar(gvar) agg(event)
estat event, window(-10 10) estore(csdd)
********************************************************************************
*Eventstudyinteract of Sun and Abraham (2020)*
********************************************************************************
*Preparation
sum Ei
g lastcohort = Ei==r(max) // dummy for the latest- or never-treated cohort
forvalues l = 0/5 {
g L`l'event = K==`l'
}
forvalues l = 1/14 {
g F`l'event = K==-`l'
}
drop F1event // normalize K=-1 (and also K=-15) to zero
eventstudyinteract Y L*event F*event, vce(cluster i) absorb(i t) cohort(Ei) control_cohort(lastcohort)
matrix sa_b = e(b_iw)
matrix sa_v = e(V_iw)
********************************************************************************
// did2s of Gardner (2021)
********************************************************************************
did2s Y, first_stage(i.i i.t) second_stage(F*event L*event) treatment(D) cluster(i)
matrix did2s_b = e(b)
matrix did2s_v = e(V)
********************************************************************************
// stackedev of Cengiz et al. (2019)
********************************************************************************
*Create a variable equal to the time that the unit first received treatment. It must be missing for never treated units.
gen treat_year=.
replace treat_year=Ei if Ei!=16
*Create never treated indicator that equals one if a unit never received treatment and zero if it did.
gen no_treat= (Ei==16)
*Preparation
cap drop F*event L*event
sum Ei
forvalues l = 0/5 {
g L`l'event = K==`l'
replace L`l'event = 0 if no_treat==1
}
forvalues l = 1/14 {
g F`l'event = K==-`l'
replace F`l'event = 0 if no_treat==1
}
drop F1event // normalize K=-1 (and also K=-15) to zero
*Run stackedev
preserve
stackedev Y F*event L*event, cohort(treat_year) time(t) never_treat(no_treat) unit_fe(i) clust_unit(i)
restore
matrix stackedev_b = e(b)
matrix stackedev_v = e(V)
restore
********************************************************************************
// TWFE Estimation
********************************************************************************
eststo TWFE: reghdfe Y F*event L*event, absorb(i t) cl(i)
// Constructing a vector of the true average treatment effects by number of periods since treatment
matrix btrue = J(1,6,.)
matrix colnames btrue = tau0 tau1 tau2 tau3 tau4 tau5
qui forvalues h = 0/5 {
sum tau if K==`h'
matrix btrue[1,`h'+1]=r(mean)
}
// Combining all the plots using the stored estimates (5 leads and lags around event)
event_plot bjs dcdh_b#dcdh_v csdd sa_b#sa_v did2s_b#did2s_v stackedev_b#stackedev_v TWFE, ///
stub_lag(tau# tau# Effect_# T+# L#event L#event L#event L#event) stub_lead(pre# pre# Placebo_# T-# F#event F#event F#event F#event) ///
plottype(scatter) ciplottype(rcap) ///
together perturb(-0.325(0.1)0.325) trimlead(5) noautolegend ///
graph_opt(title("Event study estimators in a simulated panel (400 units, 15 periods)", size(med)) ///
xtitle("Periods since the event", size(small)) ytitle("Average causal effect", size(small)) xlabel(-5(1)5, nogrid) ///
legend(order(1 "True value" 2 "Borusyak et al." 4 "de Chaisemartin-D'Haultfoeuille" ///
6 "Callaway-Sant'Anna" 8 "Sun-Abraham" 10 "Gardner" 12 "Cengiz et al." 14 "TWFE OLS") rows(2) position(6) region(style(none))) ///
xline(-0.5, lcolor(gs8) lpattern(dash)) yline(0, lcolor(gs8)) graphregion(color(white)) bgcolor(white) ylabel(, angle(horizontal))) ///
lag_opt1(msymbol(+) color(black)) lag_ci_opt1(color(black)) ///
lag_opt2(msymbol(O) color(cranberry)) lag_ci_opt2(color(cranberry)) ///
lag_opt3(msymbol(Dh) color(navy)) lag_ci_opt3(color(navy)) ///
lag_opt4(msymbol(Th) color(forest_green)) lag_ci_opt4(color(forest_green)) ///
lag_opt5(msymbol(Sh) color(dkorange)) lag_ci_opt5(color(dkorange)) ///
lag_opt6(msymbol(Th) color(blue)) lag_ci_opt6(color(blue)) ///
lag_opt7(msymbol(Dh) color(red)) lag_ci_opt7(color(red)) ///
lag_opt8(msymbol(Oh) color(purple)) lag_ci_opt8(color(purple))
