*===============================================================================
* analysis.do
* Staggered Synthetic Difference-in-Differences (SDID): Gender quotas & women in
* parliament.  Companion script for the Stata tutorial post `stata_sdid_staggered`.
*
* Run (batch):
*   "/Applications/Stata 18.0/StataMP.app/Contents/MacOS/stata-mp" -b do analysis.do
*
* Requires (all from SSC): sdid, sdid_event, panelview, reghdfe (+ ftools, unique)
*   . ssc install sdid
*   . ssc install sdid_event          // Ciccia, Clarke & Pailanir (2024)
*   . ssc install panelview
*   . ssc install reghdfe
*
* Data: quota_example.dta (Bhalotra, Clarke, Gomes & Venkataramani 2023, shipped
*       with the sdid package).  119 countries, 1990-2015, balanced panel.
*         outcome   womparl   = % women in the national parliament
*         treatment quota     = 1 once a country adopts a gender quota (absorbing)
*         control   lngdp     = log GDP per capita
*       Treatment is STAGGERED: 9 countries adopt across 7 cohorts
*       (2000, 2002, 2003, 2005, 2010, 2012, 2013); 110 never-treated controls.
*       Estimand: ATT -- the effect of adopting a quota on the adopting countries'
*       women-in-parliament share, aggregated across cohorts (Clarke et al. 2024).
*===============================================================================

clear all
set more off
set scheme s2color
set linesize 100

cd "/Users/carlosmendez/Documents/GitHub/starter-academic-v501/content/post/stata_sdid_staggered"   // students: set this to your local copy of this post folder
capture mkdir web_app
capture mkdir web_app/data

* Site colour palette (RGB strings for Stata graphs)
global TREAT "217 119 87"   // #d97757 warm orange  -> treated / adopting cohorts
global CTRL  "106 155 204"  // #6a9bcc steel blue    -> control / synthetic / never-treated
global TEAL  "0 212 200"    // #00d4c8 teal          -> SDID estimate / aggregate
global INK   "20 20 19"     // #141413 ink

* Packages (capture so re-runs do not error). sdid_event is not always on SSC.
capture ssc install sdid, replace
capture ssc install sdid_event, replace
capture which sdid_event
if _rc net install sdid_event, from("https://raw.githubusercontent.com/DiegoCiccia/sdid/main/sdid_event") replace
capture ssc install panelview, replace
capture ssc install reghdfe, replace
capture ssc install ftools, replace

capture log close
log using "analysis.log", replace text

*-------------------------------------------------------------------------------
* 1. DATA -- load, persist, document the staggered structure
*-------------------------------------------------------------------------------
capture confirm file quota_example.dta
if _rc {
    webuse set www.damianclarke.net/stata/
    webuse quota_example, clear
    save quota_example.dta, replace
}
use quota_example.dta, clear
label variable quota "Parliamentary gender quota"

describe womparl quota lngdp country year quotaYear
summarize womparl quota lngdp
encode country, gen(id)
xtset id year
xtdescribe

* Adoption cohort = first year a country is treated (quotaYear ships with the data).
* Verify it equals the first treated year, then tabulate cohort sizes.
bysort country (year): egen firsttreat = min(cond(quota==1, year, .))
gen byte evertreat = !missing(firsttreat)
capture assert firsttreat == quotaYear if evertreat   // sanity check (quotaYear ships with the data)
di as result "firsttreat == quotaYear check rc = " _rc
di as result _n "=== Adoption cohorts (one row per country) ==="
preserve
    keep country firsttreat
    duplicates drop
    tab firsttreat, missing
restore
count if quota==1
di as result "treated country-year observations = " r(N)

* cohort sizes (countries per adoption year) for later merges + the web app
preserve
    keep if evertreat
    bysort country: keep if _n==1
    contract firsttreat, freq(n_treated)
    rename firsttreat cohort
    list, noobs
    tempfile csize
    save `csize', replace
restore

*-------------------------------------------------------------------------------
* 2. EDA -- panelview (staggered structure) + a site-coloured trend figure
*-------------------------------------------------------------------------------
* (2a) Treatment-timing heatmap: the staggered "staircase"
panelview womparl quota, i(country) t(year) type(treat) bytiming               ///
    xtitle("Year") ytitle("Country (sorted by adoption timing)")               ///
    title("Staggered adoption of parliamentary gender quotas", size(medium))   ///
    ylabdist(10) xlabdist(5)
graph export "stata_sdid_staggered_panelview_treat.png", replace width(2400)

* (2b) Outcome trajectories, treated (orange) vs control (blue)
panelview womparl quota, i(country) t(year) type(outcome)                      ///
    xtitle("Year") ytitle("% women in parliament")                             ///
    title("Women in parliament: treated vs. control trajectories", size(medium))
graph export "stata_sdid_staggered_panelview_outcome.png", replace width(2400)

* (2c) Mean outcome: ever-adopting vs never-adopting (site colours)
preserve
    collapse (mean) womparl, by(evertreat year)
    reshape wide womparl, i(year) j(evertreat)
    label var womparl1 "Ever-adopting countries (mean)"
    label var womparl0 "Never-adopting countries (mean)"
    twoway (line womparl1 year, lcolor("$TREAT") lwidth(thick))                 ///
           (line womparl0 year, lcolor("$CTRL")  lwidth(medthick) lpattern(dash)), ///
           ytitle("% women in parliament") xtitle("") xlabel(1990(5)2015)       ///
           legend(order(1 "Ever-adopting (mean)" 2 "Never-adopting (mean)")      ///
                  pos(11) ring(0) cols(1) size(small))                          ///
           title("Raw outcome trends by treatment group", size(medium))         ///
           note("Adoption is staggered (2000-2013); a single group mean blurs the timing.")
    graph export "stata_sdid_staggered_raw_trends.png", replace width(2400)
restore

*-------------------------------------------------------------------------------
* 3. BASELINE -- static two-way fixed-effects DiD (the biased foil)
*    Under staggered timing with heterogeneous effects, this TWFE coefficient is
*    a contaminated weighted average that uses already-treated units as controls
*    (Goodman-Bacon 2021; de Chaisemartin & D'Haultfoeuille 2020).  Reported only
*    as a benchmark, NOT as a credible ATT.
*-------------------------------------------------------------------------------
capture reghdfe womparl quota, absorb(id year) vce(cluster id)
if _rc {
    di as error "reghdfe unavailable (rc=" _rc "); falling back to xtreg."
    xtreg womparl quota i.year, fe vce(cluster id)
}
scalar twfe_att = _b[quota]
scalar twfe_se  = _se[quota]
di as result "Static TWFE 'ATT' (biased foil) = " twfe_att "  (cluster SE " twfe_se ")"

*-------------------------------------------------------------------------------
* 4. MAIN STAGGERED SDID -- bootstrap inference + cohort-specific effects
*-------------------------------------------------------------------------------
use quota_example.dta, clear
label variable quota "Parliamentary gender quota"

sdid womparl country year quota, vce(bootstrap) seed(1213)
scalar sdid_att = e(ATT)
scalar sdid_se  = e(se)
scalar sdid_cil = e(ATT_l)
scalar sdid_cir = e(ATT_r)
di as result "SDID staggered ATT = " sdid_att "  SE = " sdid_se ///
             "  95% CI = [" sdid_cil ", " sdid_cir "]"

* ---- e(tau): cohort-specific ATTs (Tau, Std.Err., Time) ----
matrix Tau = e(tau)
matrix list Tau
preserve
    clear
    svmat Tau                                   // Tau1=tau, Tau2=se, Tau3=cohort year
    rename (Tau1 Tau2 Tau3) (tau se cohort)
    gen lci = tau - 1.96*se
    gen uci = tau + 1.96*se
    gen t_post = 2015 - cohort + 1
    merge 1:1 cohort using `csize', nogen
    * aggregation weight = treated (unit x post-period) share -> reproduces overall ATT
    gen w_raw = n_treated * t_post
    egen w_tot = total(w_raw)
    gen agg_weight = w_raw / w_tot
    gsort cohort
    order cohort tau se lci uci n_treated t_post agg_weight
    list cohort tau se n_treated t_post agg_weight, noobs
    * verification: weighted average of cohort taus = overall ATT
    gen wtau = agg_weight*tau
    egen check_att = total(wtau)
    di as result "Sum of weighted cohort taus = " check_att "  (should match overall ATT " sdid_att ")"
    drop w_raw w_tot wtau check_att
    export delimited cohort tau se lci uci n_treated t_post agg_weight ///
        using "web_app/data/cohorts.csv", replace
restore

* ---- cohort-ATT figure: tau_a with 95% CI, zero line, aggregate-ATT line ----
preserve
    clear
    svmat Tau                                   // Tau1=tau, Tau2=se, Tau3=cohort year
    rename (Tau1 Tau2 Tau3) (tau se cohort)
    gen lci = tau - 1.96*se
    gen uci = tau + 1.96*se
    twoway (rcap lci uci cohort, lcolor("$CTRL"))                               ///
           (scatter tau cohort, mcolor("$TREAT") msize(large) ms(d)),           ///
           yline(0, lcolor(gs10) lpattern(dash))                                ///
           yline(`=sdid_att', lcolor("$TEAL") lwidth(medthick))                 ///
           xlabel(2000 2002 2003 2005 2010 2012 2013, angle(45))                ///
           ytitle("Adoption-cohort ATT (pp)") xtitle("Adoption year (cohort)")   ///
           legend(off)                                                          ///
           title("Cohort-specific SDID effects", size(medium))                  ///
           note("Teal line: overall weighted ATT (8.0 pp). Cohorts range from -3.5 (2005) to +21.8 (2012).")
    graph export "stata_sdid_staggered_cohort_taus.png", replace width(2400)
restore

* ---- e(series): treated vs synthetic outcome path per cohort -> CSV ----
matrix S = e(series)
preserve
    clear
    svmat S, names(col)
    reshape long Yco Ytr, i(year) j(cohort)
    rename (Yco Ytr) (y_synth y_treated)
    drop if missing(y_treated) & missing(y_synth)
    order cohort year y_treated y_synth
    sort cohort year
    export delimited using "web_app/data/series_by_cohort.csv", replace
restore

* ---- e(lambda): pre-period time weights per cohort -> CSV ----
matrix L = e(lambda)
local cohlist 2000 2002 2003 2005 2010 2012 2013
preserve
    clear
    svmat L
    * last row holds adoption-year labels; last column (L8) is the calendar year
    local nr = rowsof(L)
    drop in `nr'
    rename L8 year
    local j = 0
    foreach c of local cohlist {
        local ++j
        rename L`j' lam`c'
    }
    reshape long lam, i(year) j(cohort)
    rename lam lambda
    drop if missing(lambda)
    order cohort year lambda
    sort cohort year
    export delimited using "web_app/data/lambda_by_cohort.csv", replace
restore

* ---- donor (unit) weights per cohort -> CSV ----
*      Use returnweights: the country name is native in the data, avoiding the
*      e(omega) matrix rownames (which mattitles fills with names that can carry
*      spaces).  A quick noinference re-fit returns the same weights.
sdid womparl country year quota, vce(noinference) returnweights
preserve
    keep country quotaYear omega2000 omega2002 omega2003 omega2005 omega2010 omega2012 omega2013
    duplicates drop
    keep if missing(quotaYear)                 // donors = never-treated countries
    drop quotaYear
    reshape long omega, i(country) j(cohort)
    drop if missing(omega) | omega==0          // keep nonzero donors
    order cohort country omega
    gsort cohort -omega
    export delimited using "web_app/data/omega_by_cohort.csv", replace
restore

* ---- treated-vs-synthetic path for the 2002 cohort (the worked example) ----
*      SDID matches the pre-period TREND, not the level (the unit fixed effect
*      absorbs the level gap).  To visualise the counterfactual we anchor the
*      synthetic to the treated cohort by its lambda-weighted pre-period gap
*      (exactly the baseline SDID differences against; see the sdid post).
preserve
    import delimited "web_app/data/lambda_by_cohort.csv", clear
    keep if cohort==2002
    keep year lambda
    tempfile l2002
    save `l2002', replace

    import delimited "web_app/data/series_by_cohort.csv", clear
    keep if cohort==2002
    keep year y_treated y_synth
    merge 1:1 year using `l2002', nogen
    replace lambda = 0 if missing(lambda)
    gen double pg = lambda*(y_synth - y_treated) if year<2002
    egen double offset = total(pg)
    gen y_synth_anch = y_synth - offset
    di as result "2002 cohort lambda-weighted pre-period gap (anchor) = " offset[1]
    twoway (line y_treated    year, lcolor("$TREAT") lwidth(thick))             ///
           (line y_synth_anch year, lcolor("$CTRL")  lwidth(medthick) lpattern(dash)), ///
           xline(2001.5, lcolor(gs10))                                          ///
           ytitle("% women in parliament") xtitle("") xlabel(1990(5)2015)       ///
           legend(order(1 "Treated cohort (2002)" 2 "Synthetic control (anchored)") ///
                  pos(11) ring(0) cols(1) size(small))                          ///
           title("SDID counterfactual for the 2002 cohort", size(medium))        ///
           note("Synthetic anchored to the treated cohort by its {&lambda}-weighted pre-2002 gap; the post-2002 divergence is the effect.")
    graph export "stata_sdid_staggered_cohort2002_path.png", replace width(2400)
restore

* ---- 2002-cohort pre-period time weights (lambda) bar chart ----
preserve
    import delimited "web_app/data/lambda_by_cohort.csv", clear
    keep if cohort==2002
    twoway (bar lambda year, color("$CTRL") barwidth(0.8)),                     ///
           ytitle("SDID time weight ({&lambda})") xtitle("")                    ///
           xlabel(1990(2)2001, angle(45)) legend(off)                           ///
           title("Where SDID looks: 2002-cohort pre-period time weights", size(medium)) ///
           note("Weight concentrates on the years just before 2002 -- the pre-period most like the post-period.")
    graph export "stata_sdid_staggered_lambda.png", replace width(2400)
restore

*-------------------------------------------------------------------------------
* 5. COVARIATES -- optimized (Arkhangelsky et al.) vs projected (Kranz 2022)
*    sdid needs a balanced panel, so drop the 104 obs with missing lngdp first.
*-------------------------------------------------------------------------------
use quota_example.dta, clear
label variable quota "Parliamentary gender quota"
drop if missing(lngdp)

sdid womparl country year quota, vce(bootstrap) seed(2022) covariates(lngdp, optimized)
scalar att_opt = e(ATT)
scalar se_opt  = e(se)
di as result "SDID + lngdp (optimized) ATT = " att_opt "  SE = " se_opt

sdid womparl country year quota, vce(bootstrap) seed(1213) covariates(lngdp, projected)
scalar att_prj = e(ATT)
scalar se_prj  = e(se)
di as result "SDID + lngdp (projected) ATT = " att_prj "  SE = " se_prj

*-------------------------------------------------------------------------------
* 6. EVENT STUDY -- sdid_event on the full staggered panel + the 2002 cohort
*-------------------------------------------------------------------------------
* (6a) Full staggered panel: aggregated ATT + cohort-aggregated dynamic effects
use quota_example.dta, clear
label variable quota "Parliamentary gender quota"
drop if missing(lngdp)
sdid_event womparl country year quota, vce(bootstrap) brep(100) effects(8) ///
    placebo(5) covariates(lngdp)
matrix Hfull = e(H)
di as result "sdid_event full-panel aggregated ATT = " Hfull[1,1] "  SE = " Hfull[1,2]

* (6b) Clean event study on the 2002 cohort (the package authors' worked example).
*      Effect_l = event time (l-1); Placebo_l = event time (-l).
use quota_example.dta, clear
label variable quota "Parliamentary gender quota"
keep if quotaYear==2002 | quotaYear==.
drop if missing(lngdp)
sdid_event womparl country year quota, vce(placebo) brep(100) placebo(all) covariates(lngdp)
matrix H = e(H)
matrix list H

local Lg  = 2015 - 2002 + 1     // 14 post-treatment (dynamic) effects
local Lpl = 2002 - 1990         // 12 pre-treatment placebos
preserve
    clear
    svmat H
    rename (H1 H2 H3 H4 H5) (coef se ci_l ci_u switchers)
    gen row = _n
    drop if row==1                                   // drop the aggregate ATT row
    gen event_time = .
    gen str4 period_type = ""
    replace event_time = row-2      if row>=2          & row<=1+`Lg'
    replace period_type = "post"    if row>=2          & row<=1+`Lg'
    replace event_time = -(row-1-`Lg') if row>=2+`Lg'  & row<=1+`Lg'+`Lpl'
    replace period_type = "pre"        if row>=2+`Lg'  & row<=1+`Lg'+`Lpl'
    keep event_time coef se ci_l ci_u period_type
    sort event_time
    list, noobs
    export delimited using "web_app/data/event_study.csv", replace

    * headline event-study figure
    twoway (rarea ci_l ci_u event_time, color("${CTRL}%35") lwidth(none))       ///
           (line coef event_time, lcolor("$TREAT") lwidth(medthick))            ///
           (scatter coef event_time, mcolor("$TREAT") msize(small) ms(O)),      ///
           yline(0, lcolor("$TEAL") lpattern(dash))                             ///
           xline(-0.5, lcolor(gs9) lpattern(solid))                            ///
           xlabel(-12(2)13) xtitle("Years relative to quota adoption (event time)") ///
           ytitle("Effect on women in parliament (pp)")                        ///
           legend(order(3 "Point estimate" 1 "95% CI") pos(11) ring(0) cols(1) size(small)) ///
           title("Event-study SDID for the 2002 cohort (sdid_event)", size(medium)) ///
           note("Pre-period placebos hug zero (parallel trends); post-period effects trace the dynamic ATT.")
    graph export "stata_sdid_staggered_event_study.png", replace width(2400)
restore

*-------------------------------------------------------------------------------
* 7. INFERENCE -- bootstrap vs placebo vs jackknife (paper's 2-cohort subsample)
*    Drop the five single-country cohorts so jackknife (needs >1 treated unit per
*    period) is defined; only the 2002 & 2003 cohorts remain.
*-------------------------------------------------------------------------------
use quota_example.dta, clear
label variable quota "Parliamentary gender quota"
drop if inlist(country,"Algeria","Kenya","Samoa","Swaziland","Tanzania")

sdid womparl country year quota, vce(bootstrap) seed(1213)
scalar b_att = e(ATT)
scalar b_se  = e(se)
scalar b_cil = e(ATT_l)
scalar b_cir = e(ATT_r)
sdid womparl country year quota, vce(placebo) seed(1213)
scalar p_att = e(ATT)
scalar p_se  = e(se)
scalar p_cil = e(ATT_l)
scalar p_cir = e(ATT_r)
sdid womparl country year quota, vce(jackknife)
scalar j_att = e(ATT)
scalar j_se  = e(se)
scalar j_cil = e(ATT_l)
scalar j_cir = e(ATT_r)

clear
set obs 3
gen str10 method = ""
gen double att = .
gen double se  = .
gen double ci_l = .
gen double ci_u = .
replace method="bootstrap" in 1
replace att=b_att in 1
replace se=b_se in 1
replace ci_l=b_cil in 1
replace ci_u=b_cir in 1
replace method="placebo" in 2
replace att=p_att in 2
replace se=p_se in 2
replace ci_l=p_cil in 2
replace ci_u=p_cir in 2
replace method="jackknife" in 3
replace att=j_att in 3
replace se=j_se in 3
replace ci_l=j_cil in 3
replace ci_u=j_cir in 3
gen tstat = att/se
gen pval  = 2*(1-normal(abs(tstat)))
list, noobs
export delimited using "web_app/data/inference.csv", replace

* forest plot of the three inference methods (same point estimate, different SEs)
gen order = _n
twoway (rcap ci_l ci_u order, horizontal lcolor("$CTRL"))                       ///
       (scatter order att, mcolor("$TREAT") msize(large) ms(d)),                ///
       xline(0, lcolor(gs10) lpattern(dash))                                    ///
       ylabel(1 "bootstrap" 2 "placebo" 3 "jackknife", angle(0))                ///
       ytitle("") xtitle("ATT on women in parliament (pp)")                     ///
       legend(off) ysc(reverse)                                                 ///
       title("Same ATT, three variance estimators (2002 & 2003 cohorts)", size(medium)) ///
       note("Point estimate is identical (10.3 pp); jackknife is most conservative, placebo tightest.")
graph export "stata_sdid_staggered_inference.png", replace width(2400)

*-------------------------------------------------------------------------------
* 8. SUMMARY TABLE -> web_app/data/atts.csv
*-------------------------------------------------------------------------------
clear
set obs 4
gen str28 spec = ""
gen double att = .
gen double se  = .
gen double ci_l = .
gen double ci_u = .
replace spec="Static TWFE (biased foil)" in 1
replace att=twfe_att in 1
replace se=twfe_se in 1
replace spec="SDID (no covariates)" in 2
replace att=sdid_att in 2
replace se=sdid_se in 2
replace ci_l=sdid_cil in 2
replace ci_u=sdid_cir in 2
replace spec="SDID + lngdp (optimized)" in 3
replace att=att_opt in 3
replace se=se_opt in 3
replace spec="SDID + lngdp (projected)" in 4
replace att=att_prj in 4
replace se=se_prj in 4
gen tstat = att/se
gen pval  = 2*(1-normal(abs(tstat)))
replace ci_l = att - 1.96*se if missing(ci_l)
replace ci_u = att + 1.96*se if missing(ci_u)
list, noobs
export delimited spec att se ci_l ci_u pval using "web_app/data/atts.csv", replace

*-------------------------------------------------------------------------------
* 9. KEY NUMBERS
*-------------------------------------------------------------------------------
di as result _n "==================== KEY NUMBERS ===================="
di as result "Static TWFE ATT (biased foil)  = " %7.2f twfe_att "  SE " %5.2f twfe_se
di as result "SDID staggered ATT             = " %7.2f sdid_att "  SE " %5.2f sdid_se
di as result "  95% CI                       = [" %5.2f sdid_cil ", " %5.2f sdid_cir "]"
di as result "SDID + lngdp (optimized)       = " %7.2f att_opt  "  SE " %5.2f se_opt
di as result "SDID + lngdp (projected)       = " %7.2f att_prj  "  SE " %5.2f se_prj
di as result "Inference subsample ATT        = " %7.2f b_att
di as result "  bootstrap SE                 = " %7.2f b_se
di as result "  placebo SE                   = " %7.2f p_se
di as result "  jackknife SE                 = " %7.2f j_se
di as result "====================================================="

log close