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)

keep if sex==2 & age>=50 & age<=70

xtset pid syear

gen age2=age*age


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

*Generate cohort indicators
gen mtreat=.
replace mtreat=0 if gebjahr<=1939
replace mtreat=1 if gebjahr>=1940 & gebjahr<1945
replace mtreat=2 if gebjahr>=1945 & gebjahr<1952
replace mtreat=3 if gebjahr>1952



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


cd "${output_path}\Results\"




twoway (scatter mretired age_int if mtreat==0, msymbol("circle")) (lpoly retired age_int if mtreat==0 & age_int<60, lpattern("solid") lcolor("black")) (lpoly retired age_int if mtreat==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 - ≤1939") xtitle("Age") ytitle("% retired") yscale(range(0 1)) ylabel(0(0.2)1) name(pre)

twoway (scatter mretired age_int if mtreat==1, msymbol("circle")) (lpoly retired age_int if mtreat==1 & age_int<60, lpattern("solid") lcolor("black")) (lpoly retired age_int if mtreat==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: Step-wise increase in deductions - 1940 to 1944") xtitle("Age") ytitle("% retired") yscale(range(0 1)) ylabel(0(0.2)1) name(reform1)

twoway (scatter mretired age_int if mtreat==2, msymbol("circle")) (lpoly retired age_int if mtreat==2 & age_int<60, lpattern("solid") lcolor("black")) (lpoly retired age_int if mtreat==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: Stable deductions of 18% - 1945 to 1951") xtitle("Age") ytitle("% retired") yscale(range(0 1)) ylabel(0(0.2)1) name(reform2)

twoway (scatter mretired age_int if mtreat==3, msymbol("circle")) (lpoly retired age_int if mtreat==3 & age_int<60, lpattern("solid") lcolor("black")) (lpoly retired age_int if mtreat==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. Post-reform: Scheme abolished - ≥1952") xtitle("Age") ytitle("% retired") yscale(range(0 1)) ylabel(0(0.2)1) name(post)

grc1leg pre reform1 reform2 post, scheme(s1mono) altshrink title("Retirement by age and birth cohort")
graph export "Fig6.tif", width(2400) replace 
graph export "Fig6.eps", replace

