sdid and event-study plots

[Lukas Lang]

Hi,

First of all a big 'thank you' to the authors of the sdid command, it really is easy to use.

I have a query about the event-study plot that is described in the Stata paper introducing this new command by Clarke et al. (2023).

My problem is that in my application the estimated 'overall' ATT (the one in the printed output table) is 3 times higher than the values that I obtain when plotting the ATT by year using the event-study plot mentioned above.

What I have noticed is that while using controls makes a big difference in terms of the estimated 'overall' ATT, it does not make that much of a difference when producing the study event plot.

To make this point I will use the same application used in the Stata paper.

Let's start with the scenario that does not control for covariates:

webuse set www.damianclarke.net/stata/
webuse quota_example.dta, clear
egen m=min(year) if quota==1, by(country) //indicator for the year of adoption
egen mm=mean(m), by(country)
keep if mm==2002 | mm==. //keep only one time of adoption
drop if lngdp==.

sdid womparl country year quota, vce(noinference) graph g2_opt(ylab(-5(5)20) ///
        ytitle("Women in Parliament") scheme(sj))
matrix lambda = e(lambda)[1..12,1] //save lambda weight
matrix yco = e(series)[1..12,2] //control baseline
matrix ytr = e(series)[1..12,3] //treated baseline
matrix aux = lambda'*(ytr - yco) //calculate the pre-treatment mean
scalar meanpre_o = aux[1,1]
matrix difference = e(difference)[1..26,1..2] // Store Ytr-Yco
svmat difference
ren (difference1 difference2) (time d)
replace d = d - meanpre_o // Calculate vector in (8)

gen y=time>=2002 & time!=.
bys y: egen d_mean=mean(d)
sort time d
preserve
keep if time==2002
di d_mean
restore

In this case the estimated 'overall' ATT (6.85377) is essentially identical to the average of the ATTs by year (6.8537698)

However, in the scenario with controls:

webuse set www.damianclarke.net/stata/
webuse quota_example.dta, clear
egen m=min(year) if quota==1, by(country) //indicator for the year of adoption
egen mm=mean(m), by(country)
keep if mm==2002 | mm==. //keep only one time of adoption
drop if lngdp==.

sdid womparl country year quota, vce(noinference) graph g2_opt(ylab(-5(5)20) ///
        ytitle("Women in Parliament") scheme(sj)) ///
        covariates(lngdp lnmmrt, projected)
matrix lambda = e(lambda)[1..12,1] //save lambda weight
matrix yco = e(series)[1..12,2] //control baseline
matrix ytr = e(series)[1..12,3] //treated baseline
matrix aux = lambda'*(ytr - yco) //calculate the pre-treatment mean
scalar meanpre_o = aux[1,1]
matrix difference = e(difference)[1..26,1..2] // Store Ytr-Yco
svmat difference
ren (difference1 difference2) (time d)
replace d = d - meanpre_o // Calculate vector in (8)

gen y=time>=2002 & time!=.
bys y: egen d_mean=mean(d)
sort time d
preserve
keep if time==2002
di d_mean
restore

The estimated 'overall' ATT is bigger (7.11653) than the average of the ATTs by year (6.8687816) which is almost as the same level as the scenario without covariates (6.8537698).

In my application this issue is exacerbated, i.e. the average of the ATTs by year is a lot lower than the estimated 'overall' ATT.

Am I misunderstanding anything?

Should we actually expect the estimated 'overall' ATT to look similar to the average of the ATTs by year?

Perhaps Damian Clarke can help with this?

Many thanks,

Lukas

[Jared Greathouse]

In their paper, they note that ATT is a weighted average of ATTS in staggered adoption settings. An average of averages, as far as I understand

[Lukas Lang]

Hi Jared, sorry I do not think I understand your answer. The example above is not on staggered DiD.

[Lukas Lang]

Any thoughts anyone?

[Jared Greathouse]

I guess I'm not understanding the issue. Are you saying that the ATT as output by SDID differs from the ATT suggested by the event study plot?

Am I understanding the issue well? Or am I missing something? If so, please provide your data using dataex (even of the results) so we can actually see what's what.


EDIT: beyond this, all I can see is

What I have noticed is that while using controls makes a big difference in terms of the estimated 'overall' ATT, it does not make that much of a difference when producing the study event plot.

Well... yeah. Adding in controls will sometimes change your results by a lot. Is that wrong? Or unexpected? I guess that's my thing, what results are you getting that you're not expecting to get? And are they because of the included covariates?

[Lukas Lang]

No, I think you've got it right now.

Let me use your own words.

When you do not use controls, ATT from sdid and from event study plot are the same.

However, when you add controls, ATT from sdid is different from ATT from event study plot.

Why is that the case?

I would expect ATT from sdid and from event study plot to remain similar also when adding controls.

[Jared Greathouse]

Are you sure you're calculating the ATT correct? Once; Daniel gave me this code to estimate the ATT. Beyond ATT calculation errors, I'm not sure why they'd be different

[Lukas Lang]

I am using exactly the same example and code in the Stata journal paper. You could cross-check the code in the first post if you want.
Perhaps we could also involve Daniel PV in this discussion?

[Jared Greathouse]

You should actually email Damian. He usually gets back with you pretty quickly, in my experience with talking with him. And... once you have the answer, post it back here. That way, everyone could learn as to why the results are this way😂

[Noah Spencer]

I have also encountered this issue and my solution was to manually calculate residualized outcomes and then run SDID using those residualized outcomes.

My solution is below:

webuse set www.damianclarke.net/stata/
webuse quota_example.dta, clear
egen m=min(year) if quota==1, by(country) //indicator for the year of adoption
egen mm=mean(m), by(country)
keep if mm==2002 | mm==. //keep only one time of adoption
drop if lngdp==.


*** My changes start:

encode country, gen(country_code)

* Estimate beta^
reg womparl lngdp lnmmrt i.country_code i.year if quota == 0

* Generate Y^{res}
gen y_res = womparl - (lngdp*e(b)[1,1] + lnmmrt*e(b)[1,2])

* Run SDID on Y^res
sdid y_res country year quota, vce(noinference) graph

*** my changes end

matrix lambda = e(lambda)[1..12,1] //save lambda weight
matrix yco = e(series)[1..12,2] //control baseline
matrix ytr = e(series)[1..12,3] //treated baseline
matrix aux = lambda'*(ytr - yco) //calculate the pre-treatment mean
scalar meanpre_o = aux[1,1]
matrix difference = e(difference)[1..26,1..2] // Store Ytr-Yco
svmat difference
ren (difference1 difference2) (time d)
replace d = d - meanpre_o // Calculate vector in (8)
gen y=time>=2002 & time!=.
bys y: egen d_mean=mean(d)
sort time d
preserve
keep if time==2002
di d_mean
restore

/// ATT = 7.1180406 (much closer to ATT estimate from standard SDID with controls)

[Noah Spencer]

Note that my solution above only works for the "projected" covariates approach. I am not sure what regression to run to estimate beta^ with non-projected covariates.

[Damian Clarke]

Thanks very much Jared for suggesting Lukas write me. I unfortunately often miss relevant things on Statalist, so thanks for the heads up! And thanks so much Noah for your response. As Noah suggests, the point that Lukas correctly raises is that with controls, you can get some inconsistencies with these two methods. The reason is that what we are doing in the event study laid out in the Stata Journal paper is to take the baseline difference as Y_co-Y_tr using the original outcome variable, rather than using the residualised version where the effect of controls have been concentrated out. This is laid out as equations (9) and (10) in the paper linked to above. Internally, when covariates are included, sdid correctly always works with the residualised version of outcomes when calculating ATTs. Presumably, what you are after Lukas is the more logical version of the event study where the baseline difference is calculated as Y^{resid}_co-Y^{resid}_tr. This is an excellent point, and one solution, noted by Noah, is to just residualise the variable first, and then proceed as standard.

While this is a nice solution, Noah is right that this will only work for the "projected" controls. The optimized controls method suggested in the original synthetic diff-in-diff paper is slightly more involved, and we did not previously return the residualised quantities from sdid. We have just issued an update on github (soon to be released on the SSC) where these values are just directly returned, meaning that you can essentially do what Noah lays out with an automatic procedure with the optimized control method. This requires doing everything you do above, but rather than using the e(series) matrix, work with the new e(series_resid) matrix. Specifically, based on the example you lay out this would be:

clear all
set more off
webuse set www.damianclarke.net/stata/
webuse quota_example.dta, clear
set more off
egen m=min(year) if quota==1, by(country) //indicator for the year of adoption
egen mm=mean(m), by(country)
keep if mm==2002 | mm==. //keep only one time of adoption
drop if lngdp==.


sdid womparl country year quota, vce(noinference) covariates(lngdp lnmmrt)
matrix lambda = e(lambda)[1..12,1] //save lambda weight
matrix yco = e(series_resid)[1..12,2] //control baseline -- NOTE USE OF series_resid
matrix ytr = e(series_resid)[1..12,3] //treated baseline -- NOTE USE OF series_resid
matrix aux = lambda'*(ytr - yco) //calculate the pre-treatment mean
scalar meanpre_o = aux[1,1]
matrix difference = e(difference)[1..26,1..2] // Store Ytr-Yco
svmat difference
ren (difference1 difference2) (time d)
replace d = d - meanpre_o // Calculate vector in (8)
gen y=time>=2002 & time!=.
bys y: egen d_mean=mean(d)
sort time d
preserve
keep if time==2002
di d_mean
restore

/// ATT = 6.9159718 (identical to SDID output)

This should all work as standard, and replicates the results as expected by Lukas in his question.

As a final point, I would just like to note that it may be worth exploring the very recently released sdid_event (developed by Diego Ciccia in interaction with sdid) which makes the process of conducting event studies based on SDID much more easy! This code which is available at the same github link where sdid is developed (a preliminary version is also available on SSC) also allows for the direct generation of event study estimates in a quite straightforward way. Diego has put together a really nice example at the link above, and in particular, this can be used to replicate the process which Noah describes in his post:

clear all
set more off
webuse set www.damianclarke.net/stata/
webuse quota_example.dta, clear
set more off
egen m=min(year) if quota==1, by(country) //indicator for the year of adoption
egen mm=mean(m), by(country)
keep if mm==2002 | mm==. //keep only one time of adoption
drop if lngdp==.
sdid_event womparl country year quota, vce(off) placebo(all) covariates(lngdp lnmmrt)

///ATT = 7.118041 as in Noah's projected method above

We will shortly update both sdid and sdid_event on the SSC so that e(series_resid) is available in sdid, and the covariates() option in sdid_event is also available there.

[Jared Greathouse]

Damian Clarke you're quite welcome. may I email you about my FDID estimator? We're (me and my two mentors who I cajoled into working with me) working on a SJ paper that describes FDID. I wanted your thoughts on extending FDID to the event study setting (I've also spoken a little with Diego about this too). I think it's the missing jewel from the crown, if you will.

[Damian Clarke]

Hi Jared, of course, I'm very interested to hear about how all this work is going!! I'm very happy to follow up on email if easiest?

[Lukas Lang]

Thanks very much Jared for suggesting Lukas write me. I unfortunately often miss relevant things on Statalist, so thanks for the heads up! And thanks so much Noah for your response. As Noah suggests, the point that Lukas correctly raises is that with controls, you can get some inconsistencies with these two methods. The reason is that what we are doing in the event study laid out in the Stata Journal paper is to take the baseline difference as Y_co-Y_tr using the original outcome variable, rather than using the residualised version where the effect of controls have been concentrated out. This is laid out as equations (9) and (10) in the paper linked to above. Internally, when covariates are included, sdid correctly always works with the residualised version of outcomes when calculating ATTs. Presumably, what you are after Lukas is the more logical version of the event study where the baseline difference is calculated as Y^{resid}_co-Y^{resid}_tr. This is an excellent point, and one solution, noted by Noah, is to just residualise the variable first, and then proceed as standard.

While this is a nice solution, Noah is right that this will only work for the "projected" controls. The optimized controls method suggested in the original synthetic diff-in-diff paper is slightly more involved, and we did not previously return the residualised quantities from sdid. We have just issued an update on github (soon to be released on the SSC) where these values are just directly returned, meaning that you can essentially do what Noah lays out with an automatic procedure with the optimized control method. This requires doing everything you do above, but rather than using the e(series) matrix, work with the new e(series_resid) matrix. Specifically, based on the example you lay out this would be:

clear all
set more off
webuse set www.damianclarke.net/stata/
webuse quota_example.dta, clear
set more off
egen m=min(year) if quota==1, by(country) //indicator for the year of adoption
egen mm=mean(m), by(country)
keep if mm==2002 | mm==. //keep only one time of adoption
drop if lngdp==.


sdid womparl country year quota, vce(noinference) covariates(lngdp lnmmrt)
matrix lambda = e(lambda)[1..12,1] //save lambda weight
matrix yco = e(series_resid)[1..12,2] //control baseline -- NOTE USE OF series_resid
matrix ytr = e(series_resid)[1..12,3] //treated baseline -- NOTE USE OF series_resid
matrix aux = lambda'*(ytr - yco) //calculate the pre-treatment mean
scalar meanpre_o = aux[1,1]
matrix difference = e(difference)[1..26,1..2] // Store Ytr-Yco
svmat difference
ren (difference1 difference2) (time d)
replace d = d - meanpre_o // Calculate vector in (8)
gen y=time>=2002 & time!=.
bys y: egen d_mean=mean(d)
sort time d
preserve
keep if time==2002
di d_mean
restore

/// ATT = 6.9159718 (identical to SDID output)

This should all work as standard, and replicates the results as expected by Lukas in his question.

As a final point, I would just like to note that it may be worth exploring the very recently released sdid_event (developed by Diego Ciccia in interaction with sdid) which makes the process of conducting event studies based on SDID much more easy! This code which is available at the same github link where sdid is developed (a preliminary version is also available on SSC) also allows for the direct generation of event study estimates in a quite straightforward way. Diego has put together a really nice example at the link above, and in particular, this can be used to replicate the process which Noah describes in his post:

clear all
set more off
webuse set www.damianclarke.net/stata/
webuse quota_example.dta, clear
set more off
egen m=min(year) if quota==1, by(country) //indicator for the year of adoption
egen mm=mean(m), by(country)
keep if mm==2002 | mm==. //keep only one time of adoption
drop if lngdp==.
sdid_event womparl country year quota, vce(off) placebo(all) covariates(lngdp lnmmrt)

///ATT = 7.118041 as in Noah's projected method above

We will shortly update both sdid and sdid_event on the SSC so that e(series_resid) is available in sdid, and the covariates() option in sdid_event is also available there.

A very big thank you to both Damian and Noah for their help with this.