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
gen dretired=d.vretired
bys pid: egen total=total(abs(dretired))


gen pre_kids=children if f.dretired==1 & total==1
gen pre_inc=hghinc if f.dretired==1 & total==1


bys pid: egen kids_ind=mean(pre_kids) if total==1
bys pid: egen hinc_ind=mean(pre_inc) if total==1

gen vkids_cat=cond(kids_ind==0,0,cond(kids_ind==1,1,cond(kids_ind>=2 & kids_ind!=.,2,.)))
sum hinc_ind, detail
gen vinc_cat=cond(hinc_ind>r(p75) & hinc_ind!=.,1,cond(hinc_ind<=r(p25) & hinc_ind!=.,0,.))
 
 
preserve
clear
gen beta=.
gen lower=.
gen upper=.
gen depvar=""
gen kidscat=.
save "Timing_birth_fathers_hetero4_kids.dta", replace
restore


foreach num of numlist 0/2{
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>=50 & vage<=70 & age>=25 & age<=40 & vkids_cat==`num', fe cluster(pid) partial(i.month i.syear)


preserve
    use "Timing_birth_fathers_hetero4_kids.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 kidscat=`num' in l
	replace depvar=e(depvar) in l
    save "Timing_birth_fathers_hetero4_kids.dta", replace
restore
}
}


use "Timing_birth_fathers_hetero4_kids.dta", clear


bys kidscat: gen sortnr=_n
gen t=sortnr-8



replace t=t+0.1 if kidscat==1
replace t=t+0.2 if kidscat==2


sort t

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





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
gen dretired=d.mretired
bys pid: egen total=total(abs(dretired))


gen pre_kids=children if f.dretired==1 & total==1
gen pre_inc=hghinc if f.dretired==1 & total==1


bys pid: egen kids_ind=mean(pre_kids) if total==1
bys pid: egen hinc_ind=mean(pre_inc) if total==1

gen mkids_cat=cond(kids_ind==0,0,cond(kids_ind==1,1,cond(kids_ind>=2 & kids_ind!=.,2,.)))
sum hinc_ind, detail
gen minc_cat=cond(hinc_ind>r(p75) & hinc_ind!=.,1,cond(hinc_ind<=r(p25) & hinc_ind!=.,0,.))


preserve
clear
gen beta=.
gen lower=.
gen upper=.
gen depvar=""
gen kidscat=.
save "Timing_birth_mothers_hetero4_kids.dta", replace
restore




foreach num of numlist 0/2{
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>=50 & mage<=70 & age>=25 & age<=40 & mkids_cat==`num', fe cluster(pid) partial(i.month i.syear)
preserve
    use "Timing_birth_mothers_hetero4_kids.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 kidscat=`num' in l
	replace depvar=e(depvar) in l
    save "Timing_birth_mothers_hetero4_kids.dta", replace
restore

}
}


use "Timing_birth_mothers_hetero4_kids.dta", clear


bys kidscat: gen sortnr=_n
gen t=sortnr-8



replace t=t+0.1 if kidscat==1
replace t=t+0.2 if kidscat==2



sort t

twoway (rspike upper lower t if kidscat==0, lpattern(solid) lcolor(black) lwidth(medium) yline(0)) (rspike upper lower t if kidscat==1, lpattern(dash) lcolor(black) lwidth(medium)) (rspike upper lower t if kidscat==2, lpattern(dot) lcolor(black) lwidth(medium)) ///
(scatter beta t if kidscat==0, msize(small)) (scatter beta t if kidscat==1, msize(small)) (scatter beta t if kidscat==2, msize(small)), ytitle("Effect size") title("b. Maternal retirement") name(mothers) ///
xtitle("Time to retirement") legend(label(1 "95% CI - 0") label(2 "95% CI - 1") label(3 "95% CI - 2+") label(4 "0") label(5 "1") label(6 "2+")) scheme(s1mono) xscale(range(-7 7)) yscale(range(-0.6 0.6)) ylabel(-0.6(0.1)0.6) xlabel(-7(1)7) xline(-0.1)
grc1leg fathers mothers, altshrink scheme(s1mono) col(1)
graph export "${output_path}\Results\FigA6.tif", replace width(2400) 


