global scen : env SCEN

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

global nreps : env NREPS

include ../../fem_env.do

/*** load FEM and average over reps (so we count reps with zero population) ***/

use ../../output/$scen/${scen}_by_rep.dta, clear

keep year rep endpop*

* drop data years beyond Census projections
drop if year > 2050

gen stockpopsize_m = .
gen stockpopsize_f = .
forvalues yy=2004(2)2050 {
	  local stockagemin = 51 + `yy' - 2004
	  foreach sex in m f {
	            replace stockpopsize_`sex' = endpop_`sex'_`stockagemin'p if year==`yy'
          }
}

drop endpop*
reshape long stockpopsize, i(year rep) j(sexlab) string
gen sex = strupper(substr(sexlab,2,1))
drop sexlab

collapse (sum) stockpopsize, by(year sex)
replace stockpopsize = stockpopsize / $nreps

rename stockpopsize stockpopsize_fem
preserve

/*** append Census projections for merging ***/
clear
foreach yr in 2004 /* 2012 2017 */ {
	local projdata `projdata' data/census_proj`yr'.dta
}
append using `projdata', generate(projectionyear)
recode projectionyear (1=2004) (2=2012) (3=2017)
tempfile censusproj
save `censusproj'

restore
preserve
merge 1:m year sex using `censusproj'
tab _merge
replace stockpopsize_fem = 0.0 if _merge==2
drop _merge

sort year projectionyear sex

* rescale to 100,000s (previously scaled to 1,000s)
replace stockpopsize_census = stockpopsize_census / 100.0
replace stockpopsize_fem = stockpopsize_fem / 100.0

bys sex projectionyear: egen init_stockpopsize = max(stockpopsize_census)
sort year sex projectionyear
li year sex projectionyear stockpopsize_census init_stockpopsize if year <= 2018

* difference in FEM population and Census projection
gen excess_pop = 100.0 * (stockpopsize_fem - stockpopsize_census) / init_stockpopsize
* gen excess_pop = stockpopsize_fem - stockpopsize_census
gen pop_ratio = stockpopsize_fem / stockpopsize_census

gen series = "FEM vs. Census projection"

li if sex=="M"

li if sex=="F"

keep year sex excess_pop pop_ratio series

tempfile fem_v_proj
save `fem_v_proj'

/*** merge Census estimates ***/

restore
drop if year > 2016
merge 1:1 year sex using data/census_popsizeest.dta
tab _merge
replace stockpopsize_fem = 0.0 if _merge==2
drop _merge

* rescale to 100,000s (FEM previously scaled to 1,000s)
replace stockpopsizeest = stockpopsizeest / 100000.0
replace stockpopsize_fem = stockpopsize_fem / 100.0

sort sex year
by sex: egen init_stockpopsize = max(stockpopsizeest)
sort year sex
li year sex stockpopsizeest init_stockpopsize if year <= 2018

* difference in FEM population and Census projection
gen excess_pop = 100.0 * (stockpopsize_fem - stockpopsizeest) / init_stockpopsize
* gen excess_pop = stockpopsize_fem - stockpopsizeest
gen pop_ratio = stockpopsize_fem / stockpopsizeest

gen series = "FEM vs. Census estimate"

li if sex=="M"

li if sex=="F"

keep year sex excess_pop pop_ratio series

append using `fem_v_proj'

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