cd "${output_path}\tmp\"


*Figure 5 and 6 are estimated based on parent level data, rather than children-parent level data. Therefore, we use individual-level observation (dataset child.dta) and restrict observations to fathers and mothers for whom we also observe the child in our main estimation sample.
use "Linked_data.dta", clear
keep if age>=25 & age<=40
keep vnr
duplicates drop
rename vnr pid
save "IDs_fathers.dta", replace

use "Linked_data.dta", clear
keep if age>=25 & age<=40
keep mnr
duplicates drop
rename mnr pid
save "IDs_mothers.dta", replace

append using "IDs_fathers.dta"

save "Complete_IDs_parents.dta", replace

graph drop _all




use "child.dta", replace
merge m:1 pid using "Complete_IDs_parents.dta", keep(match) gen(merge_childsample)

*Restrict to men aged 50-70
keep if sex==1 & age>=50 & age<=70

xtset pid syear

gen age2=age*age


gen vr60=cond(age>=60 & age!=.,1,cond(age!=.,0,.))

*Generate treatment indicator based on birth cohort. 

gen vtreat=.
replace vtreat=0 if gebjahr<=1936
replace vtreat=1 if gebjahr>=1937 & gebjahr<=1941
replace vtreat=2 if gebjahr>=1942 & gebjahr<=1945
replace vtreat=3 if gebjahr>=1946 & gebjahr<=1948
replace vtreat=4 if gebjahr>=1949 & gebjahr<=1951
replace vtreat=5 if gebjahr>=1952


gen age_int=floor(age)
bys vtreat age_int: egen mretired=mean(retired)

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



twoway (scatter mretired age_int if vtreat==0, msymbol("circle")) ///
(lpoly retired age_int if vtreat==0 & age_int<60, lpattern("solid") lcolor("black")) (lpoly retired age_int if vtreat==0 & age_int>=60, lpattern("solid") lcolor("black")), xline(60) scheme(s1mono) ///
legend(label(1 "Means") label(2 "Local polynomial fit") order(1 2)) title("A. Pre-reform cohorts - ≤1936") xtitle("Age") ytitle("% retired") name(part1) yscale(range(0 1)) ylabel(0(0.2)1)

twoway (scatter mretired age_int if vtreat==1, msymbol("circle")) ///
(lpoly retired age_int if vtreat==1 & age_int<60, lpattern("solid") lcolor("black")) (lpoly retired age_int if vtreat==1 & age_int>=60, lpattern("solid") lcolor("black")), xline(60) scheme(s1mono) ///
legend(label(1 "Means") label(2 "Local polynomial fit") order(1 2)) title("B. Reform I: Extended scope - 1937-1941") xtitle("Age") ytitle("% retired") name(part2) yscale(range(0 1)) ylabel(0(0.2)1)

twoway (scatter mretired age_int if vtreat==2, msymbol("circle")) ///
(lpoly retired age_int if vtreat==2 & age_int<60, lpattern("solid") lcolor("black")) (lpoly retired age_int if vtreat==2 & age_int>=60, lpattern("solid") lcolor("black")), xline(60) scheme(s1mono) ///
legend(label(1 "Means") label(2 "Local polynomial fit") order(1 2)) title("C. Reform II: Deductions of 18% - 1942-1945") xtitle("Age") ytitle("% retired") name(part3) yscale(range(0 1)) ylabel(0(0.2)1)

twoway (scatter mretired age_int if vtreat==3, msymbol("circle")) ///
(lpoly retired age_int if vtreat==3 & age_int<60, lpattern("solid") lcolor("black")) (lpoly retired age_int if vtreat==3 & age_int>=60, lpattern("solid") lcolor("black")), xline(60) scheme(s1mono) ///
legend(label(1 "Means") label(2 "Local polynomial fit") order(1 2)) title("D. Reform III: Threshold increases - 1946-1948") xtitle("Age") ytitle("% retired") name(part4) yscale(range(0 1)) ylabel(0(0.2)1)

twoway (scatter mretired age_int if vtreat==4, msymbol("circle")) ///
(lpoly retired age_int if vtreat==4 & age_int<60, lpattern("solid") lcolor("black")) (lpoly retired age_int if vtreat==4 & age_int>=60, lpattern("solid") lcolor("black")), xline(60) scheme(s1mono) ///
legend(label(1 "Means") label(2 "Local polynomial fit") order(1 2)) title("E. Reform IV: Early retirement at 63 - 1949-1951") xtitle("Age") ytitle("% retired") name(part5) yscale(range(0 1)) ylabel(0(0.2)1)

twoway (scatter mretired age_int if vtreat==5, msymbol("circle")) ///
(lpoly retired age_int if vtreat==5 & age_int<60, lpattern("solid") lcolor("black")) (lpoly retired age_int if vtreat==5 & age_int>=60, lpattern("solid") lcolor("black")), xline(60) scheme(s1mono) ///
legend(label(1 "Means") label(2 "Local polynomial fit") order(1 2)) title("F. Reform V: Scheme abolished - ≥1952") xtitle("Age") ytitle("% retired") name(part6) yscale(range(0 1)) ylabel(0(0.2)1)

grc1leg part1 part2 part3 part4 part5 part6, scheme(s1mono) altshrink title("Retirement by age and birth cohort")
graph export "Fig5.tif", width(2400) replace 
graph export "Fig5.eps", replace 


