* Set directory 
****************

local user Guangyi
// local user Justin

if "`user'"=="Guangyi" {
	global main "/Users/guangyi/Library/CloudStorage/Box-Box/School segregation/NHIS Segregation/Manuscripts/Journal Submission/DiD methods paper/2 Epidemiology/R&R/simulation"
	cd "$main"
	}
	

if "`user'"=="Justin" {
	global main "~/Dropbox (BOSTON UNIVERSITY)/Research/School segregation/Analysis"
	global code "$main/code"
	cd "$main"
	}

	

* Run simulation for all scenarios (HTE + PTA violation: steeper slope for early treated)
******************************************************************************************
local scenario "1a 1b 1c 2a 2b 2c 3a 3b 3c 3d 3e 3f"
foreach x of local scenario {
	global scenario "`x'"
	run "code/DiD_simulation_for Epidemiology submission.do"  
}


// local scenario "1a 1b 1c 2a 2b 2c 3a 3b 3c 3d 3e 3f"
// foreach x of local scenario {
// 	global scenario "`x'"
// 	run "$code/DiD_simulation_500_30_60_02-24-2024.do"
// }


* Create figures and tables 
******************************************************
	
* calculate bias & RMSE *
local scenario "1a 1b 1c 2a 2b 2c 3a 3b 3c 3d 3e 3f"
foreach x of local scenario {
	use "data/sim_estimates_`x'.dta", clear

	* Calculate bias and RMSE *
	foreach var in beta_twfe beta_cs beta_bjs beta_sa beta_jw {
		* % bias 
		gen bias_`var' = `var' - beta_true
		gen bias_percent_`var' = bias_`var' * 100 / beta_true // % bias
		* RMSE
		gen sqr_bias_`var' = bias_percent_`var'^2
		egen mean_sqr_bias_`var' = mean(sqr_bias_`var')
		gen rmse_`var' = sqrt(mean_sqr_bias_`var')
	}

	
* plot *
	gr hbox bias_percent*, noout ///
					legend(order(1 "TWFE" 2 "CS" 3 "BJS" 4 "SA" 5 "JW") row(1)) ///
					box(1, color(gs12)) box(2, color(gs9)) box(3, color(gs6)) box(4, color(gs3)) box(5, color(gs1)) ///
					yline(0, lcolor(gs0)) ///
					graphregion(color(white)) bgcolor(white) ///
					ytitle("Bias (%)") ///
					title("Scenario `x'") ///
					ylabel(-160(40)60)
					
	gr save "figures/scenario_`x'.gph", replace
}


* combine figures *
// figure 1: HTE: no PTA violation //
grc1leg "figures/scenario_1a.gph" ///
		"figures/scenario_1b.gph" ///
		"figures/scenario_1c.gph" ///
		"figures/scenario_2a.gph" ///
		"figures/scenario_2b.gph" ///
		"figures/scenario_2c.gph"		
gr save "figures/scenario_HTE_PTA_1a-2c.gph", replace 


// figure 2: HTE: PTA violation //
grc1leg "figures/scenario_3a.gph" ///
		"figures/scenario_3b.gph" ///
		"figures/scenario_3c.gph" ///
		"figures/scenario_3d.gph" ///
		"figures/scenario_3e.gph" ///
		"figures/scenario_3f.gph"	
gr save "figures/scenario_HTE_violated_PTA_3a-3f.gph", replace


* tables *
eststo clear
local scenario "1a 1b 1c 2a 2b 2c 3a 3b 3c 3d 3e 3f"
foreach x of local scenario {
	use "data/sim_estimates_`x'.dta", clear

	* Calculate bias and RMSE *
	foreach var in beta_twfe beta_cs beta_bjs beta_sa beta_jw {
		* % bias 
		gen bias_`var' = `var' - beta_true
		gen bias_percent_`var' = bias_`var' * 100 / beta_true // % bias
// 		drop if bias_percent_`var' < -1000 | bias_percent_`var' > 1000 // drop outliers
		summarize bias_percent_`var', detail
		keep if inrange(bias_percent_`var', r(p1), r(p99)) // drop outliers
		* RMSE
		gen sqr_bias_`var' = bias_`var'^2
		egen mean_sqr_bias_`var' = mean(sqr_bias_`var')
		gen rmse_`var' = sqrt(mean_sqr_bias_`var')
	}

	* sum bias and RMSE *
	eststo: estpost summarize bias_percent* rmse_*
}
esttab using "tables/table3_simulation_1a-3f_v2.csv", replace ///
			  cells("mean(fmt(2))") noobs ///
			  title(Table 3. Simulation results: bias(%) and RMSE for estiation of the ATT) ///
			  mtitle("1a" "1b" "1c" "2a" "2b" "2c" "3a" "3b" "3c" "3d" "3e" "3f")