clear
clear matrix
cd "${output_path}\tmp\"

gen sim=.
gen beta_fathers=.
gen se_fathers=.


gen beta_mothers=.
gen se_mothers=.

save "placebo_simulation.dta", replace


use "Linked_data.dta", replace
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



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 fbirth=f.birth




egen y5=cut(syear), at(1984, 1989, 1994, 1999, 2004, 2009, 2017)

egen age5=cut(age), at(25,30,35,41)
egen mage5=cut(mage), at(47.5,52.5,57.5,62.5,67.5,72.5)
egen vage5=cut(vage), at(47.5,52.5,57.5,62.5,67.5,72.5)




bys age5 sex y5 vage5: egen vbirthrate=mean(fbirth)
bys age5 sex y5 mage5: egen mbirthrate=mean(fbirth)




foreach sim of numlist 1/1000{

preserve

gen vplacebo=rbinomial(1,vbirthrate)


gen mplacebo=rbinomial(1,mbirthrate)


qui xi: xtivreg2 vplacebo age age2 vage vage2 i.month i.syear (vretired = vdisc60) if vage>=50 & vage<=70 & age>=25 & age<=40, fe cluster(pid) partial(i.month i.syear)
estimates store fathers


qui xi: xtivreg2 mplacebo age age2 mage mage2 i.month i.syear (mretired = mdisc60) 	if mage>=50 & mage<=70 & age>=25 & age<=40, fe cluster(pid) partial(i.month i.syear)
estimates store mothers


use "placebo_simulation.dta", clear

			local nplus = _N +1   
			set obs `nplus' 
	
			replace sim=`sim' in l
			
			estimates restore fathers
			replace beta_fathers=_b[vretired] in l
			replace se_fathers=_se[vretired] in l
			
		
			
			estimates restore mothers
			replace beta_mothers=_b[mretired] in l
			replace se_mothers=_se[mretired] in l
			
	
			
save "placebo_simulation.dta", replace
restore
}
