SDID event study graph, Loop over all permutation issue

[Pierangela Peruzzo]

Hi,

I'm trying to replicate this event study graph in Clark et al., 2023 using Synthetic difference in difference. However, I am encountering some issues in the part where I loop over all permutation replications.
I have no idea how to solve it and why it is not working.
Can someone help me?
Thanks a lot!

/*=======================================================================
    Store Equation (8) data series ========================================================================*/  
* *Prepare data
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==. //keep if covariates observed

* Run SDiD
sdid womparl country year quota, vce(noinference) graph

* Calculate second term in Equation 8
matrix lambda = e(lambda)[1..12,1] // time weights
matrix yco = e(series)[1..12,2] // pre-treatment values for treated unit
matrix ytr = e(series)[1..12,3] // average pre-treatment values for control units
matrix aux = lambda'*(ytr - yco) // calculate the pre-treatment baseline difference
scalar meanpre_o = aux[1,1]        // store this value
di meanpre_o

* Store pre-made output on first term of Equation 8
matrix difference = e(difference)[1..26,1..2]  // store Ytr - Yco
svmat difference
keep difference1 difference2
ren (difference1 difference2) (year trt_ctrl_diff)

* Calculate the vector in equation 8
gen eqn8 = trt_ctrl_diff - meanpre_o
drop trt_ctrl_diff
drop if year == .

* Save
save "eqn8_values", replace




/*=======================================================================
    Do permutation procedure to generate confidence intervals
========================================================================*/

* Prepare data
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==. //keep if covariates observed


* Drop treated units
drop if country == "Djibouti" | country == "Morocco"

* Reset state codes so they're all in order
encode country, gen(country_code)


* Set permutation counter and number of reps
local p = 1
local P = 100


* Loop over all permutation replications
while (`p' <= `P') {
    drop quota
    // Need to choose 2 placebo states.  Let's just take max and min of rand number
    // This could probably be done in less lines, but the below should work fine
    gen c = rnormal()
    bys country: egen chooser = min(c)
    qui sum chooser
    gen treated = chooser==r(min) | chooser==r(max)
    gen quota = (year >= 2002 & treated==1)
    * Run SDiD
    qui: sdid womparl country year quota, vce(noinference) graph
    * Collect the quantities you stored above to calculate the second term in Equation 8
    matrix lambda_p = e(lambda)[1..12,1]
    matrix yco_p = e(series)[1..12,2]
    matrix ytr_p = e(series)[1..12,3]
    matrix aux_p = lambda_p'*(ytr_p - yco_p)
    matrix meanpre_p = J(26,1,aux_p[1,1])
    * Collect Equation 8 values (and store in a matrix)
    matrix d`p'=e(difference)[1..26,2] - meanpre_p
    * Update loop counter
    local ++p
    drop c chooser treated
}

* Save all placebo estimates of Equation 8 values
sort country year
forval p=1/`P' {
    svmat d`p'
}

* Calculate standard deviation of estimates across each time period
egen rsd = rowsd(d11 - d`P'1)
keep year rsd
drop if rsd == .

* Merge in true Equation 8 values
merge 1:1 year using "eqn8_values", nogen

* Generate confidence intervals
gen LCI = eqn8 + invnormal(0.025)*rsd
gen UCI = eqn8 + invnormal(0.975)*rsd


/*=======================================================================
    Generate plot
========================================================================*/


*generate plot
tw rarea UCI LCI year, color(gray%40) || scatter eqn8 year, color(blue) m(d) xtitle("") ytitle("") legend(order(2 "Point Estimate" 1 "95% CI") pos(12) col(2)) xline(2002, lc(black) lp(solid)) yline(0, lc(red) lp(shortdash)) scheme(sj)
graph export sdidplacebo.pdf, replace

[Daniel Schaefer]

When you say it's "not working" what do you mean? Does it keep going seemingly forever? Does it produce an error message? Does it produce incorrect results?

[Pierangela Peruzzo]

It produces an error "<= invalid name" when running this command:

while (`p' <= `P') {
    drop quota
    // Need to choose 2 placebo states.  Let's just take max and min of rand number
    // This could probably be done in less lines, but the below should work fine
    gen c = rnormal()
    bys country: egen chooser = min(c)
    qui sum chooser
    gen treated = chooser==r(min) | chooser==r(max)
    gen quota = (year >= 2002 & treated==1)
    * Run SDiD
    qui: sdid womparl country year quota, vce(noinference) graph
    * Collect the quantities you stored above to calculate the second term in Equation 8
    matrix lambda_p = e(lambda)[1..12,1]
    matrix yco_p = e(series)[1..12,2]
    matrix ytr_p = e(series)[1..12,3]
    matrix aux_p = lambda_p'*(ytr_p - yco_p)
    matrix meanpre_p = J(26,1,aux_p[1,1])
    * Collect Equation 8 values (and store in a matrix)
    matrix d`p'=e(difference)[1..26,2] - meanpre_p
    * Update loop counter
    local ++p
    drop c chooser treated
}

[Daniel Schaefer]

Okay. It sounds to me like you're trying to run the while loop on its own without also instantiating the two locals, `p' and `P'. Locals are created and destroyed at the start and end of execution respectively. If you just highlight the while loop and run that without also highlighting these lines:

local p = 1
local P = 100

You will run into an error. So you should make sure you are running the relevant locals at the same time.

Additionally, from a style perspective, this should probably be a forvalues loop, not a while loop. Something like this might be a bit better, and you would avoid having to create the two locals p and P in the two lines above:

forvalues p = 1/100 {
    drop quota
    // Need to choose 2 placebo states. Let's just take max and min of rand number
    // This could probably be done in less lines, but the below should work fine
    gen c = rnormal()
    bys country: egen chooser = min(c)
    qui sum chooser
    gen treated = chooser==r(min) | chooser==r(max)
    gen quota = (year >= 2002 & treated==1)
    * Run SDiD
    qui: sdid womparl country year quota, vce(noinference) graph
    * Collect the quantities you stored above to calculate the second term in Equation 8
    matrix lambda_p = e(lambda)[1..12,1]
    matrix yco_p = e(series)[1..12,2]
    matrix ytr_p = e(series)[1..12,3]
    matrix aux_p = lambda_p'*(ytr_p - yco_p)
    matrix meanpre_p = J(26,1,aux_p[1,1])
    * Collect Equation 8 values (and store in a matrix)
    matrix d`p'=e(difference)[1..26,2] - meanpre_p
    * Don't need to update on the next line with a forvalues loop. I've commented it out for demonstration purposes, but you can delete.
    * local ++p
    drop c chooser treated
}

I would normally test all of this on my end first, but I have my own technical issue I need to deal with first. Definitely let us know if this doesn't work.

[Pierangela Peruzzo]

Thank you very much! It works!