cd "${output_path}\tmp\"

use "Linked_data.dta", clear
*Restrict sample to children aged 25-40
keep if age>=25 & age<=40

xtset pid syear

*This is our outcome
gen fbirth=f.birth
gen age2=age*age
gen age3=age2*age

gen vage2=vage*vage
gen vage3=vage2*vage
gen mage2=mage*mage

*These are the discontinuity indicators for fathers and mothers
gen vr60=cond(vage>=60 & vage!=.,1,cond(vage!=.,0,.))
gen mr60=cond(mage>=60 & age!=.,1,cond(mage!=.,0,.))

*Reform indicators for fathers and mothers based on year of birth
gen vtreat=.
replace vtreat=0 if vgebjahr<=1936
replace vtreat=1 if vgebjahr>=1937 & vgebjahr<=1941
replace vtreat=2 if vgebjahr>=1942 & vgebjahr<=1945
replace vtreat=3 if vgebjahr>=1946 & vgebjahr<=1948
replace vtreat=4 if vgebjahr>=1949 & vgebjahr<=1951
replace vtreat=5 if vgebjahr>=1952


gen mtreat=.
replace mtreat=0 if mgebjahr<=1939
replace mtreat=1 if mgebjahr>=1940 & mgebjahr<1945
replace mtreat=2 if mgebjahr>=1945 & mgebjahr<1952
replace mtreat=3 if mgebjahr>1952


*Standardizing age
foreach var of varlist vage vage2 mage mage2{
sum `var'
gen `var'_sd=(`var'-r(mean))/r(sd)
}


*These are the additional isntruments created by interacting the discontinuity with the reform indicators
gen mr60_treat1=cond(mtreat==1,mr60,0)
gen mr60_treat2=cond(mtreat==2,mr60,0)
gen mr60_treat3=cond(mtreat==3,mr60,0)

gen vr60_treat1=cond(vtreat==1,vr60,0)
gen vr60_treat2=cond(vtreat==2,vr60,0)
gen vr60_treat3=cond(vtreat==3,vr60,0)
gen vr60_treat4=cond(vtreat==4,vr60,0)
gen vr60_treat5=cond(vtreat==5,vr60,0)

*Cohort-specific age trend for fathers and mothers
gen vage_treat1=cond(vtreat==1,vage,0)
gen vage_treat2=cond(vtreat==2,vage,0)
gen vage_treat3=cond(vtreat==3,vage,0)
gen vage_treat4=cond(vtreat==4,vage,0)
gen vage_treat5=cond(vtreat==5,vage,0)

gen vage2_treat1=cond(vtreat==1,vage2,0)
gen vage2_treat2=cond(vtreat==2,vage2,0)
gen vage2_treat3=cond(vtreat==3,vage2,0)
gen vage2_treat4=cond(vtreat==4,vage2,0)
gen vage2_treat5=cond(vtreat==5,vage2,0)

gen mage_treat1=cond(mtreat==1,mage,0)
gen mage_treat2=cond(mtreat==2,mage,0)
gen mage_treat3=cond(mtreat==3,mage,0)

gen mage2_treat1=cond(mtreat==1,mage2,0)
gen mage2_treat2=cond(mtreat==2,mage2,0)
gen mage2_treat3=cond(mtreat==3,mage2,0)


cd "${output_path}\Results\"

set logtype text
cap log close
log using "IV_Validity.txt", replace 


preserve
*We first estimate the model to restrict the sample to our working sample. The remainder of the code follows the example code provided by Mourifié and Wan. We use 500 grid points for the search and include father's age (quadratic trend) and children's age (quadratic trend) as control variables.
eststo: xi: xtivreg2 fbirth age age2 vage vage2 (vretired= vr60) if vage>=50 & vage<=70,fe cluster(pid) first savefirst savefprefix(_all_)
keep if e(sample)==1


egen c1=mean(vretired)
gen c0=1-c1

scalar lb=0
scalar ub=1

gen T=lb+0.002*_n
replace T=. if T>ub

gen T2=T
gen T3=T
gen T4=T
gen T5=T

gen l1=c1*vretired*(1-vr60)-c0*vretired*vr60
gen l0=c0*(1-vretired)*vr60-c1*(1-vretired)*(1-vr60)


clrtest (l1 fbirth vage vage2 age age2 T T2 T3 T4 T5) (l0 fbirth vage vage2 age age2 T T2 T3 T4 T5), level(0.9,0.95,0.99) met("parametric")
restore

eststo: xi: xtivreg2 fbirth age age2 mage mage2 (mretired= mr60) if mage>=50 & mage<=70,fe cluster(pid) first savefirst savefprefix(_all_)
keep if e(sample)==1


egen c1=mean(mretired)
gen c0=1-c1

scalar lb=0
scalar ub=1

gen T=lb+0.002*_n
replace T=. if T>ub

gen T2=T
gen T3=T
gen T4=T
gen T5=T

gen l1=c1*mretired*(1-mr60)-c0*mretired*mr60
gen l0=c0*(1-mretired)*mr60-c1*(1-mretired)*(1-mr60)


clrtest (l1 fbirth mage mage2 age age2 T T2 T3 T4 T5) (l0 fbirth mage mage2 age age2 T T2 T3 T4 T5), level(0.9,0.95,0.99) met("parametric")


log close