clear
clear matrix

graph drop _all
cd "${output_path}\tmp\"

use "Linked_data.dta", clear
set more off

gen age2=age*age
gen age3=age2*age

gen vage2=vage*vage
gen vage3=vage2*vage


gen vdisc60=cond((vage+0.001)>=(60+1/12),1,0)
replace vdisc60=cond((vage+0.001)>=(60+2/12),1,0) if vgebjahr==1952 & vgebmonat==1
replace vdisc60=cond((vage+0.001)>=(60+3/12),1,0) if vgebjahr==1952 & vgebmonat==2
replace vdisc60=cond((vage+0.001)>=(60+4/12),1,0) if vgebjahr==1952 & vgebmonat==3
replace vdisc60=cond((vage+0.001)>=(60+5/12),1,0) if vgebjahr==1952 & vgebmonat==4
replace vdisc60=cond((vage+0.001)>=(60+6/12),1,0) if vgebjahr==1952 & vgebmonat==5
replace vdisc60=cond((vage+0.001)>=(60+7/12),1,0) if vgebjahr==1952 & vgebmonat>=6 & vgebmonat<=12
replace vdisc60=cond((vage+0.001)>=(60+8/12),1,0) if vgebjahr==1953 
replace vdisc60=cond((vage+0.001)>=(60+9/12),1,0) if vgebjahr==1954
replace vdisc60=cond((vage+0.001)>=(60+10/12),1,0) if vgebjahr==1955
replace vdisc60=cond((vage+0.001)>=(60+11/12),1,0) if vgebjahr==1956
replace vdisc60=cond((vage+0.001)>=61,1,0) if vgebjahr==1957
 

sort pid syear


preserve
clear
gen beta=.
gen lower=.
gen upper=.
gen depvar=""

save "Timing_birth_fathers_5y.dta", replace
restore



foreach var of varlist  l7.birth l6.birth l5.birth l4.birth l3.birth l2.birth l.birth birth f.birth f2.birth f3.birth f4.birth f5.birth f6.birth f7.birth{

eststo: xi: xtivreg2 `var' age age2 vage vage2 i.month i.syear (vretired = vdisc60) if vage>=55 & vage<=65 & age>=25 & age<=40 , fe cluster(pid) partial(i.month i.syear)

preserve
    use "Timing_birth_fathers_5y.dta", clear
    local nplus = _N +1   
    set obs `nplus' 
	
	replace beta=_b[vretired] in l
	replace  lower= _b[vretired]-invnormal(0.975)*_se[vretired] in l
	replace  upper= _b[vretired]+invnormal(0.975)*_se[vretired] in l
	replace depvar=e(depvar) in l
    save "Timing_birth_fathers_5y.dta", replace
restore

}


use "Timing_birth_fathers_5y.dta", clear

gen sortnr=_n
gen t=sortnr-8





sort t

twoway (rspike upper lower t , lpattern(solid) lcolor(black) lwidth(medium) yline(0))  (scatter beta t , msize(small)) , ytitle("Effect size") title("a. Paternal retirement") ///
xtitle("Time to retirement") legend(label(1 "95% CI")  label(2 "Point estimate")) scheme(s1mono) xscale(range(-7 7)) yscale(range(-0.6 0.6)) xlabel(-7(1)7) xline(-0.1) ylabel(-0.6(0.1)0.6) name(fathers)






clear
clear matrix
cap log close
set logtype text

use "Linked_data.dta", clear
set more off


gen age2=age*age
gen age3=age2*age

gen mage2=mage*mage
gen mage3=mage2*mage



gen mdisc60=cond((mage+0.001)>=(60+1/12),1,0)
replace mdisc60=cond((mage+0.001)>=(60+2/12),1,0) if mgebjahr==1952 & mgebmonat==1
replace mdisc60=cond((mage+0.001)>=(60+3/12),1,0) if mgebjahr==1952 & mgebmonat==2
replace mdisc60=cond((mage+0.001)>=(60+4/12),1,0) if mgebjahr==1952 & mgebmonat==3
replace mdisc60=cond((mage+0.001)>=(60+5/12),1,0) if mgebjahr==1952 & mgebmonat==4
replace mdisc60=cond((mage+0.001)>=(60+6/12),1,0) if mgebjahr==1952 & mgebmonat==5
replace mdisc60=cond((mage+0.001)>=(60+7/12),1,0) if mgebjahr==1952 & mgebmonat>=6 & mgebmonat<=12
replace mdisc60=cond((mage+0.001)>=(60+8/12),1,0) if mgebjahr==1953 
replace mdisc60=cond((mage+0.001)>=(60+9/12),1,0) if mgebjahr==1954
replace mdisc60=cond((mage+0.001)>=(60+10/12),1,0) if mgebjahr==1955
replace mdisc60=cond((mage+0.001)>=(60+11/12),1,0) if mgebjahr==1956
replace mdisc60=cond((mage+0.001)>=61,1,0) if mgebjahr==1957
 


sort pid syear


preserve
clear
gen beta=.
gen lower=.
gen upper=.
gen depvar=""

save "Timing_birth_mothers_5y.dta", replace
restore



foreach var of varlist  l7.birth l6.birth l5.birth l4.birth l3.birth l2.birth l.birth birth f.birth f2.birth f3.birth f4.birth f5.birth f6.birth f7.birth{

eststo: xi: xtivreg2 `var' age age2 mage mage2 i.month i.syear (mretired = mdisc60) if mage>=55 & mage<=65 & age>=25 & age<=40, fe cluster(pid) partial(i.month i.syear)

preserve
    use "Timing_birth_mothers_5y.dta", clear
    local nplus = _N +1   
    set obs `nplus' 
	
	replace beta=_b[mretired] in l
	replace  lower= _b[mretired]-invnormal(0.975)*_se[mretired] in l
	replace  upper= _b[mretired]+invnormal(0.975)*_se[mretired] in l

	replace depvar=e(depvar) in l
    save "Timing_birth_mothers_5y.dta", replace
restore

}


use "Timing_birth_mothers_5y.dta", clear

gen sortnr=_n
gen t=sortnr-8





sort t

twoway (rspike upper lower t, lpattern(solid) lcolor(black) lwidth(medium) yline(0) xscale(range(-7 7)))  ///
(scatter beta t , msize(small) xscale(range(-7 7))) , ytitle("Effect size") title("b. Maternal retirement") ///
xtitle("Time to retirement") legend(label(1 "95% CI") label(2 "Point estimate")) scheme(s1mono) xscale(range(-7 7)) xlabel(-7(1)7) xline(-0.1) yscale(range(-0.6 0.6)) ylabel(-0.6(0.1)0.6) name(mothers)




graph combine fathers mothers, altshrink scheme(s1mono)
graph export "${output_path}\Results\FigA3_B.tif", replace width(3200) height(1200)

