*=============================================================================== * analysis.do * Synthetic Difference-in-Differences (SDID): California's Proposition 99 * Companion script for the Stata tutorial post `stata_sdid`. * * Run (batch): * "/Applications/Stata/StataSE.app/Contents/MacOS/stata-se" -b do analysis.do * * Requires (all from SSC): sdid, synth, synth2 * . ssc install sdid * . ssc install synth2 * . ssc install synth * * Data: prop99_example.dta (Arkhangelsky et al. 2021 / Abadie et al. 2010). * 39 US states, 1970-2000, outcome = cigarette packs per capita. * California is the single treated unit; Proposition 99 takes effect 1989. * The panel is OUTCOME-ONLY (no covariates), so synthetic control and SDID * see exactly the same information set -- an apples-to-apples comparison. * Estimand: ATT -- the effect of Proposition 99 on California, 1989-2000. *=============================================================================== clear all set more off set scheme s2color set linesize 100 cd "/Users/carlos/GitHub/starter-academic-v501/content/tutorials/stata_sdid" // 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 -> California / observed global CTRL "106 155 204" // #6a9bcc steel blue -> control / synthetic global TEAL "0 212 200" // #00d4c8 teal -> SDID global INK "20 20 19" // #141413 *------------------------------------------------------------------------------- * 1. DATA *------------------------------------------------------------------------------- capture confirm file prop99_example.dta if _rc { webuse set www.damianclarke.net/stata/ webuse prop99_example.dta, clear save prop99_example.dta, replace } use prop99_example.dta, clear describe tab treated summarize packspercapita encode state, gen(id) xtset id year summ id if state=="California", meanonly local ca = r(mean) scalar ca_id = `ca' di as result "California id = `ca'" *------------------------------------------------------------------------------- * 2. EDA: California vs. the simple average of the 38 control states *------------------------------------------------------------------------------- preserve gen byte iscal = state=="California" collapse (mean) packs=packspercapita, by(iscal year) reshape wide packs, i(year) j(iscal) label var packs1 "California" label var packs0 "Average of 38 controls" twoway (line packs1 year, lcolor("$TREAT") lwidth(thick)) /// (line packs0 year, lcolor("$CTRL") lwidth(medthick) lpattern(dash)), /// xline(1989, lcolor(gs10) lpattern(solid)) /// ytitle("Cigarette packs per capita") xtitle("") /// xlabel(1970(5)2000) /// legend(order(1 "California" 2 "Average of 38 controls") pos(1) ring(0) cols(1) size(small)) /// title("California vs. the raw control average", size(medium)) /// note("Vertical line: Proposition 99 takes effect (1989).") graph export "stata_sdid_raw_trends.png", replace width(2000) restore *------------------------------------------------------------------------------- * 3. The original difference-in-differences: a raw 2x2 *------------------------------------------------------------------------------- gen byte cal = state=="California" gen byte post = year>=1989 quietly summ packspercapita if cal==1 & post==0 scalar m_ca_pre = r(mean) quietly summ packspercapita if cal==1 & post==1 scalar m_ca_post = r(mean) quietly summ packspercapita if cal==0 & post==0 scalar m_co_pre = r(mean) quietly summ packspercapita if cal==0 & post==1 scalar m_co_post = r(mean) scalar did2x2 = (m_ca_post - m_ca_pre) - (m_co_post - m_co_pre) di as result "2x2 DiD = " did2x2 * identical to the interaction in a saturated regression: reg packspercapita i.cal##i.post drop cal post *------------------------------------------------------------------------------- * 4. The original synthetic control (Abadie et al. 2010) via synth2 * Match on the full pre-period outcome PATH (each pre-year a separate * predictor) -- the fair analogue to SDID's unit weights. *------------------------------------------------------------------------------- local preds "" forvalues y = 1970/1988 { local preds "`preds' packspercapita(`y')" } synth2 packspercapita `preds', trunit(`ca') trperiod(1989) frame(sc2) symbol(2) nofigure scalar sc_att = e(att) scalar sc_rmse = e(rmse) matrix scU = e(U_wt) matrix list scU * export synth2 donor weights (rownames carry the donor names) local rn : rownames scU file open fh using "web_app/data/sc_omega.csv", write replace file write fh "state,weight" _n local i = 0 foreach s of local rn { local ++i local w = scU[`i',1] file write fh "`s',`w'" _n } file close fh * pull California's synthetic path + gap out of the synth2 frame (by var label) tempfile scser frame change sc2 keep if id==`ca' foreach v of varlist _all { local lbl : variable label `v' if strpos("`lbl'","prediction") rename `v' sc_synth if strpos("`lbl'","treatment effect") rename `v' sc_effect } keep year sc_synth sc_effect save `scser', replace frame change default * SC inference: in-space placebo (RMSPE-ratio) test -> p-value (optional) capture noisily synth2 packspercapita `preds', trunit(`ca') trperiod(1989) placebo(unit) nofigure if _rc==0 { matrix scPV = e(pval) di as result "--- synth2 in-space placebo p-values e(pval) ---" matrix list scPV } *------------------------------------------------------------------------------- * 5. Synthetic difference-in-differences (Arkhangelsky et al. 2021) via sdid * Run A: point estimate + canonical figure + returned unit/time weights *------------------------------------------------------------------------------- use prop99_example.dta, clear encode state, gen(id) sdid packspercapita state year treated, method(sdid) vce(noinference) graph g1on /// returnweights mattitles /// g1_opt(ylabel(-110(20)50) xtitle("")) /// g2_opt(ylabel(0(25)150) ytitle("Packs per capita")) graph export "stata_sdid_sdid_main.png", replace width(2000) scalar sdid_att = e(ATT) di as result "SDID ATT = " sdid_att * capture the treated + SDID synthetic outcome trajectories matrix S = e(series) tempfile sdidser preserve clear svmat S, names(col) rename (Yco1989 Ytr1989) (sdid_synth ca_actual) keep year sdid_synth ca_actual save `sdidser', replace restore * export SDID unit weights (omega) and time weights (lambda) preserve keep state omega1989 duplicates drop rename omega1989 omega gsort -omega export delimited using "web_app/data/sdid_omega.csv", replace restore tempfile lamf preserve keep year lambda1989 duplicates drop rename lambda1989 lambda export delimited using "web_app/data/sdid_lambda.csv", replace save `lamf', replace restore capture drop omega1989 lambda1989 * Run B: placebo inference (the valid choice with ONE treated unit) sdid packspercapita state year treated, vce(placebo) seed(1213) scalar sdid_se = e(se) scalar sdid_cil = e(ATT_l) scalar sdid_cir = e(ATT_r) di as result "SDID ATT = " e(ATT) " SE = " sdid_se " 95% CI = [" sdid_cil "," sdid_cir "]" *------------------------------------------------------------------------------- * 6. The three estimators in one framework: method(did|sc|sdid) *------------------------------------------------------------------------------- sdid packspercapita state year treated, method(did) vce(noinference) graph g1on /// g1_opt(ylabel(-110(20)50) xtitle("")) /// g2_opt(ylabel(0(25)150) ytitle("Packs per capita")) graph export "stata_sdid_did_panel.png", replace width(2000) scalar did_att = e(ATT) sdid packspercapita state year treated, method(sc) vce(noinference) graph g1on /// g1_opt(ylabel(-110(20)50) xtitle("")) /// g2_opt(ylabel(0(25)150) ytitle("Packs per capita")) graph export "stata_sdid_sc_panel.png", replace width(2000) scalar sc_sdidframe_att = e(ATT) di as result "DiD(sdid)=" did_att " SC(sdid)=" sc_sdidframe_att " SDID=" sdid_att *------------------------------------------------------------------------------- * 7. Putting them side by side: one chart of California vs. every counterfactual *------------------------------------------------------------------------------- preserve keep if state!="California" collapse (mean) ctrl_mean=packspercapita, by(year) tempfile cm save `cm', replace restore use `sdidser', clear merge 1:1 year using `scser', nogen merge 1:1 year using `cm', nogen merge 1:1 year using `lamf', nogen * DiD counterfactual: California's pre-level shifted by the controls' change summ ca_actual if year<=1988, meanonly local capre = r(mean) summ ctrl_mean if year<=1988, meanonly local copre = r(mean) gen did_cf = `capre' + (ctrl_mean - `copre') * SDID counterfactual: e(series) Yco matches TRENDS not levels (unit FE absorbs * the offset), so anchor it by the lambda-WEIGHTED pre-period gap -- exactly the * baseline SDID differences against. Post-period (ca_actual - sdid_cf) then * averages to the SDID ATT. gen double pregap = lambda*(sdid_synth - ca_actual) if year<=1988 & lambda<. egen double offset_sdid = total(pregap) gen sdid_cf = sdid_synth - offset_sdid summ offset_sdid, meanonly di as result "SDID lambda-weighted pre-period gap (anchor) = " r(mean) drop pregap offset_sdid order year ca_actual did_cf sc_synth sdid_cf sdid_synth ctrl_mean sc_effect lambda export delimited using "web_app/data/series.csv", replace list year ca_actual did_cf sc_synth sdid_cf, sepby(year) noobs twoway (line ca_actual year, lcolor("$TREAT") lwidth(thick)) /// (line did_cf year, lcolor(gs7) lpattern(dash)) /// (line sc_synth year, lcolor("$CTRL") lpattern(shortdash) lwidth(medthick)) /// (line sdid_cf year, lcolor("$TEAL") lpattern(solid) lwidth(medthick)), /// xline(1989, lcolor(gs10)) /// ytitle("Cigarette packs per capita") xtitle("") xlabel(1970(5)2000) /// legend(order(1 "California (observed)" 2 "DiD counterfactual" /// 3 "Synthetic control (synth2)" 4 "SDID counterfactual") /// rows(2) pos(6) size(small)) /// title("Four counterfactuals for California", size(medium)) /// note("SDID counterfactual anchored to California by its {&lambda}-weighted pre-period gap.") graph export "stata_sdid_compare_paths.png", replace width(2000) * synth2 path + gap figures (built from the merged series) twoway (line ca_actual year, lcolor("$TREAT") lwidth(thick)) /// (line sc_synth year, lcolor("$CTRL") lpattern(dash) lwidth(medthick)), /// xline(1989, lcolor(gs10)) ytitle("Packs per capita") xtitle("") /// xlabel(1970(5)2000) /// legend(order(1 "California" 2 "Synthetic California") pos(1) ring(0) cols(1) size(small)) /// title("Synthetic control fit (synth2)", size(medium)) graph export "stata_sdid_sc_path.png", replace width(2000) twoway (line sc_effect year, lcolor("$INK") lwidth(medthick)), /// yline(0, lcolor(gs10)) xline(1989, lcolor(gs10) lpattern(dash)) /// ytitle("Gap: California - synthetic") xtitle("") xlabel(1970(5)2000) /// title("Estimated gap (synth2)", size(medium)) graph export "stata_sdid_sc_gap.png", replace width(2000) * time-weight bar chart (SDID puts all pre-weight on 1986-1988) preserve import delimited "web_app/data/sdid_lambda.csv", clear twoway (bar lambda year if year<=1988, color("$CTRL") barwidth(0.8)), /// ytitle("SDID time weight ({&lambda})") xtitle("") /// xlabel(1970(2)1988) ylabel(0(.1).5) legend(off) /// title("Where SDID looks: pre-period time weights", size(medium)) /// note("Pre-period weight (1970-1988): zero until 1986, then 0.37, 0.21, 0.43 on 1986-1988." /// "Post-1989 years are omitted; SDID weights them uniformly at 1/12.") graph export "stata_sdid_lambda.png", replace width(2000) restore *------------------------------------------------------------------------------- * 8. Inference figure: SDID in-space placebo distribution * Drop California; assign each control as the placebo-treated unit at 1989; * re-estimate SDID; collect placebo ATTs. p = share with |placebo| >= |obs|. *------------------------------------------------------------------------------- use prop99_example.dta, clear drop if state=="California" levelsof state, local(ctrls) tempname pf postfile `pf' str20 pstate double ptau using "_placebo.dta", replace foreach s of local ctrls { preserve gen byte ptreat = (state=="`s'") & (year>=1989) capture sdid packspercapita state year ptreat, vce(noinference) if _rc==0 post `pf' ("`s'") (e(ATT)) restore } postclose `pf' use "_placebo.dta", clear local aatt = sdid_att count if abs(ptau) >= abs(`aatt') local pc = r(N) count local pn = r(N) scalar sdid_pperm = `pc'/`pn' di as result "SDID placebo permutation p-value = " sdid_pperm " (n=`pn')" export delimited using "web_app/data/placebo.csv", replace local aatts : di %4.1f `aatt' local pvals : di %5.3f sdid_pperm twoway (histogram ptau, width(2) color("${CTRL}%70") lcolor(white)), /// xline(`aatt', lcolor("$TREAT") lwidth(thick)) /// xtitle("Placebo ATT among control states") ytitle("Density") /// title("SDID placebo distribution", size(medium)) /// note("Orange line: California's estimated ATT (`aatts'). Permutation p = `pvals'.") graph export "stata_sdid_placebo_hist.png", replace width(2000) *------------------------------------------------------------------------------- * 9. Summary table of every ATT estimate -> web_app/data/atts.csv *------------------------------------------------------------------------------- clear set obs 5 gen str30 method = "" gen double att = . gen double se = . gen double ci_l = . gen double ci_r = . gen double note_pval = . replace method = "Raw 2x2 DiD" in 1 replace att = did2x2 in 1 replace method = "DiD (TWFE, sdid)" in 2 replace att = did_att in 2 replace method = "Synthetic control (synth2)" in 3 replace att = sc_att in 3 replace method = "SC (sdid framework)" in 4 replace att = sc_sdidframe_att in 4 replace method = "SDID" in 5 replace att = sdid_att in 5 replace se = sdid_se in 5 replace ci_l = sdid_cil in 5 replace ci_r = sdid_cir in 5 replace note_pval = sdid_pperm in 5 export delimited using "web_app/data/atts.csv", replace list, noobs *------------------------------------------------------------------------------- * Clean up scratch file *------------------------------------------------------------------------------- capture erase "_placebo.dta" di as result _n "==================== KEY NUMBERS ====================" di as result "Raw 2x2 DiD = " %7.2f did2x2 di as result "DiD (sdid framework) = " %7.2f did_att di as result "Synthetic control (synth2) = " %7.2f sc_att " (RMSE " %4.2f sc_rmse ")" di as result "SC (sdid framework) = " %7.2f sc_sdidframe_att di as result "SDID = " %7.2f sdid_att di as result "SDID placebo SE = " %7.2f sdid_se di as result "SDID 95% CI = [" %6.2f sdid_cil ", " %6.2f sdid_cir "]" di as result "SDID permutation p-value = " %5.3f sdid_pperm di as result "====================================================="