/* Process bootstrap simulation results and do the averaging for "nested" simulations */

local reps 500
local nested 20
local simulations psid_bootstrap_pop psid_bootstrap_coh
local measures p_diabe p_cancre p_hearte p_hibpe p_lunge p_stroke p_smoken m_endpop a_bmi
local subpop all m f


foreach sim in `simulations' {
	use ../output/`sim'/`sim'_by_rep.dta, replace

	* ID for each bs sample
	gen bs = ceil((rep+1)/`nested')

	foreach mes in `measures' {
		foreach sub in `subpop' {
			preserve
			local var `mes'_`sub'
			* collapse to the BS level
			collapse `var', by(bs year)
			* collapse to the summary level
			collapse `var' (p2) `var'_p2=`var' (p98) `var'_p98=`var' (sem) `var'_sem=`var' (sd) `var'_sd=`var', by(year)
			
			gen `var'_upper_se = `var' + `var'_sem
			gen `var'_lower_se = `var' - `var'_sem
			
			gen `var'_upper_sd = `var' + 1.96*`var'_sd
			gen `var'_lower_sd = `var' - 1.96*`var'_sd
			
			* Show CI as 2nd and 98th percentiles
			#d ;
			twoway line `var' `var'_p2  `var'_p98 year, 
			lpattern(solid dash dash)
			saving(bs_results/`sim'_`var'_02_98.gph, replace)
			;
			#d cr
			
			* Show CI as standard error about the mean (this is very tight ...)
			#d ;
			twoway line `var' `var'_lower_se `var'_upper_se year, 
			lpattern(solid dash dash)
			saving(bs_results/`sim'_`var'_SEM.gph, replace)
			;
			#d cr
			
			* Show CI as 1.96 standard deviations
			#d ;
			twoway line `var' `var'_lower_sd `var'_upper_sd year, 
			lpattern(solid dash dash)
			saving(bs_results/`sim'_`var'_SD.gph, replace)
			;
			#d cr
			
			restore
		}
	}
}













capture log close