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
replace father_near=l.father_near if syear==1992 | syear==2002 | syear==2007 | syear==2012 | syear==2017
replace father_near=l2.father_near if syear==1993 | syear==2003 | syear==2008 | syear==2013
replace father_near=f.father_near if syear==1990 | syear==2000 | syear==2005 | syear==2010 | syear==2015
replace father_near=f2.father_near if syear==1989 | syear==1999 | syear==2004 | syear==2009 | syear==2014

 
 



preserve
clear
gen beta=.
gen lower=.
gen upper=.
gen depvar=""
gen dist=.
save "Timing_birth_fathers_hetero1_distance.dta", replace
restore


foreach num of numlist 0/1{
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 & father_near==`num', fe cluster(pid) partial(i.month i.syear)

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

}
}

use "Timing_birth_fathers_hetero1_distance.dta", clear


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



replace t=t+0.1 if dist==1


sort t

twoway (rspike upper lower t if dist==0, lpattern(solid) lcolor(black) lwidth(medium) yline(0)) (rspike upper lower t if dist==1, lpattern(dash) lcolor(black) lwidth(medium)) ///
(scatter beta t if dist==0, msize(small)) (scatter beta t if dist==1, msize(small)), ytitle("Effect size") title("a. Paternal retirement") name(fathers) ///
xtitle("Time to retirement") legend(label(1 "95% CI - >1 hour") label(2 "95% CI - ≤1 hour") label(3 ">1 hour") label(4 "≤1 hour")) 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
replace mother_near=l.mother_near if syear==1992 | syear==2002 | syear==2007 | syear==2012 | syear==2017
replace mother_near=l2.mother_near if syear==1993 | syear==2003 | syear==2008 | syear==2013
replace mother_near=f.mother_near if syear==1990 | syear==2000 | syear==2005 | syear==2010 | syear==2015
replace mother_near=f2.mother_near if syear==1989 | syear==1999 | syear==2004 | syear==2009 | syear==2014



preserve
clear
gen beta=.
gen lower=.
gen upper=.
gen depvar=""
gen dist=.
save "Timing_birth_mothers_hetero1_distance.dta", replace
restore


foreach num of numlist 0/1{
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 & mother_near==`num', fe cluster(pid) partial(i.month i.syear)
preserve
    use "Timing_birth_mothers_hetero1_distance.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 dist=`num' in l
	replace depvar=e(depvar) in l
    save "Timing_birth_mothers_hetero1_distance.dta", replace
restore

}
}

use "Timing_birth_mothers_hetero1_distance.dta", clear


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



replace t=t+0.1 if dist==1


sort t

twoway (rspike upper lower t if dist==0, lpattern(solid) lcolor(black) lwidth(medium) yline(0) xscale(range(-7 7))) (rspike upper lower t if dist==1, lpattern(dash) lcolor(black) lwidth(medium) xscale(range(-7 7))) ///
(scatter beta t if dist==0, msize(small) xscale(range(-7 7))) (scatter beta t if dist==1, msize(small) xscale(range(-7 7))), ytitle("Effect size") title("b. Maternal retirement") name(mothers) ///
xtitle("Time to retirement") legend(label(1 "95% CI - >1 hour") label(2 "95% CI - ≤1 hour") label(3 ">1 hour") label(4 "≤1 hour")) 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\FigA4.tif", replace 

