/** compute standardized mortality rate for FEM mortality using CDC life tables as reference population
*/

global scen : env SCEN

log using logs/${scen}_cdc_smr.log, replace

* prepare CDC life tables for comparison
clear
forvalues yr=2005/2016 {
          append using ../../base_data/cdc_nvsr_lifetable`yr'.dta
}

desc

* only doing gender categories right now
drop if !missing(black) | !missing(hispanic)

drop if age < 55 | age > 99
rename age agegrp

* collapse to 2-year mortality in even years
gen birthcohort = year - age

sort male birthcohort age

gen pdied2yr_cdc = 1 - (1 - qx[_n-1])*(1 - qx) if mod(year, 2)==0 & age > 55

# 2-year mortality rate = (num. of deaths in years 1 and 2) / (approx. exposure time in years 1 and 2)
gen mx2yr =  (dx[_n-1] + (lx[_n-1]*(1-qx[_n-1])*qx)) / (Lx[_n-1] + lx[_n-1]*(1-qx[_n-1])*(1-qx/2)) if mod(year, 2)==0 & age > 55

* drop odd years for merging with FEM results
drop if mod(year,2)==1

tempfile cdc_rates
save `cdc_rates'

* merge FEM and CDC rates for comparison
use data/${scen}_mortality_rates.dta
drop if year < 2006
sort year male agegrp
desc
li

merge 1:1 year male agegrp using `cdc_rates'
keep if _merge==3

li

gen pdied2yr = n_died / startpop

gen excess_died  = n_died - (startpop * pdied2yr_cdc)

gen excess_rate = pdied2yr / pdied2yr_cdc

* save age-specific 2-year mortality probabilities for plotting
save data/${scen}_cdc_mortbyage.dta, replace

gen expdied = exposure * mx2yr

egen agebin = cut(agegrp), at(40,55,65,75,85,125)

collapse (sum) n_died expdied, by(male year agebin)

gen smr = n_died / expdied

sort male agebin year

li

save data/${scen}_cdc_smr.dta, replace