clear all
use "/Users/matthijskalmijn/Dropbox/OG trends/OG_full.dta", clear

gen div=divfam
label variable div "Parents divorced"
label variable agediv "Age at separation"
label variable sex "Daughter vs son"
label define z 0 "intact families" 1 "separated families"
label variable retrobiasdiv "Length of time since divorce"
label values div z
keep if div==1
tabl div
* Create same N
drop if div==1 & (relfatherdiv==. | relmotherdiv==. | relparentsdiv==. |  agediv==.) // 330 cases dropped

gen divyear=birthyear+agediv
recode divyear (1938/1969=1)(1970/1979=2)(1980/1989=3)(1990/2009=4), gen(divcoh)
label define divcoh 1 "< 1970" 2 "1970s" 3 "1980s" 4 "1990s"
label values divcoh divcoh
tab divcoh

recode relmotherdiv (1=2) // 5 cases added to bad

tab divcoh relfatherdiv, row 
tab divcoh relmotherdiv, row 
tab divcoh relparentsdiv, row 
tab divcoh relfatherdiv, row nofreq
tab divcoh relmotherdiv, row nofreq
tab divcoh relparentsdiv, row nofreq

bys divcoh: corr rel*div 
egen mby=mean(birthyear)
egen mdy=mean(divyear)
gen by=(birthyear-mby)/10
gen dy=(divyear-mdy)/10
label variable by "Birth year"
label variable dy "Divorce year"

recode eduparents (1=1)(2 3=0)(else=.), gen(ed1)
recode eduparents (2=1)(1 3=0)(else=.), gen(ed2)
recode eduparents (3=1)(1 2=0)(else=.), gen(ed3)
label variable ed1 "Parents lower education"
label variable ed2 "Parents middle education"
label variable ed3 "Parents higher education"

gen div_by=div*by
gen div_dy=div*dy

gen ed1_by=ed1*by
gen ed2_by=ed2*by
gen ed3_by=ed3*by

gen ed1_dy=ed1*dy
gen ed2_dy=ed2*dy
gen ed3_dy=ed3*dy

gen ed1_div=ed1*div
gen ed2_div=ed2*div
gen ed3_div=ed3*div

gen ed1_div_by=ed1*div*by
gen ed2_div_by=ed2*div*by
gen ed3_div_by=ed3*div*by

recode divyear (min/1969=1)(1970/1979=2)(1980/1989=3)(1990/max=4), gen(coh4)
tab coh4 relfatherdiv, row

* Logistic regression analyses divorced parents
recode relfatherdiv (1=1)(2 3 4=0), gen(fnone)
recode relfatherdiv (1=.)(2=1)(3=2)(4=3), gen(fqual)
recode relmotherdiv (1=.)(2=1)(3=2)(4=3), gen(mqual)
recode relparentsdiv (1=1)(2 3 4=0), gen(pnone)
recode relparentsdiv (1=.)(2=1)(3=2)(4=3), gen(pqual)
label variable fnone "No contact father"
label variable fqual "Quality relation father"
label variable mqual "Quality relation mother"
label variable pnone "No contact parents"
label variable pqual "Quality relation parents"
drop if div~=1

eststo: estpost summarize sex agediv divyear retrobiasdiv ed1 ed2 ed3 fnone fqual 
esttab using "/Users/matthijskalmijn/Dropbox/og trends/tables-revised2-final", ///
 cells("mean(fmt(2)) sd(fmt(2)) min(fmt(2)) max(fmt(2)) count(fmt(0))") rtf nodepvar label ///
 title("Table 1: Descriptive statistics of variables in sample") replace compress nomtitle nonumber ///
 substitute("\fonttbl{\f0\fnil Times New Roman" "\fonttbl{\f0\fnil Calibri")

eststo a1: logit fnone dy sex agediv retrobiasdiv ed2 ed3 
eststo a2: logit fnone dy sex agediv retrobiasdiv ed2 ed3 ed2_dy ed3_dy
eststo a3: ologit fqual dy sex agediv retrobiasdiv ed2 ed3 
eststo a4: ologit fqual dy sex agediv retrobiasdiv ed2 ed3 ed2_dy ed3_dy

esttab a* using "/Users/matthijskalmijn/Dropbox/og trends/tables-revised2-final", b(3) ar2 nogap parentheses rtf append star( ~ 0.10 * 0.05) ///
  nolz not p title("Table 2: Logit regression of post-divorce relations fathers") label ///
  note(Note: P-values in parentheses. ///
             Year of separation centered and divided by 10.) ///
  stats(N chi2 ll,labels("N" "Chi-2" "Log-likelihood") fmt(%8.0f %8.1f %8.1f)) 

save test, replace
use test, clear
sum sex agediv retrobiasdiv ed2 ed3 ed2_dy ed3_dy
mi set wide
mi register imputed  ed2 ed3 ed2_dy ed3_dy
mi impute chained (regress) ed2 ed3 ed2_dy ed3_dy  ///
 = fnone sex agediv retrobiasdiv, add(30) rseed(03062019)
eststo a1: mi estimate, post: logit fnone dy sex agediv retrobiasdiv ed2 ed3 ed2_dy ed3_dy

use test, clear
sum sex agediv retrobiasdiv ed2 ed3 ed2_dy ed3_dy
drop if fqual==.
mi set wide
mi register imputed  ed2 ed3 ed2_dy ed3_dy
mi impute chained (regress) ed2 ed3 ed2_dy ed3_dy  ///
 = fqual sex agediv retrobiasdiv, add(30) rseed(03062019)
eststo a1: mi estimate, post: ologit fqual dy sex agediv retrobiasdiv ed2 ed3 ed2_dy ed3_dy

use test, clear
labe define edunew 1 "Low education" 2 "Middle education" 3 "High education"
label values eduparents edunew

logit fnone c.divyear i.sex c.agediv c.retrobiasdiv i.eduparents i.eduparents#c.divyear
margins, at(divyear=(1950 1960 1970 1980 1990 2000) agediv retrobiasdiv sex) noestimcheck post over(eduparents)
set scheme mrc
keep if divyear>=1950
marginsplot, xdimension(divyear) plotdimension(eduparents) name(test,replace) ///
 title(" ") ylab(0(.10).80, labsize(small)) xlab(1950(10)2000, labsize(small)) plotregion(color(white)) legend(rows(1) size(small)) ///
 xtitle(Divorce cohort,size(small)) ytitle(Proportion no contact, size(small))  plot1opts(msymbol(S)) plot2opts(msymbol(O)) plot3opts(msymbol(B)) 
graph export "/users/matthijskalmijn/dropbox/og trends/figure2-final.pdf", replace as(pdf) 

use test, clear
ologit fqual c.divyear i.sex c.agediv c.retrobiasdiv i.eduparents i.eduparents#c.divyear
margins,  predict(outcome(3)) at(divyear=(1950 1960 1970 1980 1990 2000) agediv retrobiasdiv sex) noestimcheck post over(eduparents)
set scheme mrc
keep if divyear>=1950
marginsplot, xdimension(divyear) plotdimension(eduparents) name(test,replace) ///
 title(" ") ylab(0(.10).70, labsize(small)) xlab(1950(10)2000, labsize(small)) plotregion(color(white)) legend(rows(1) size(small)) ///
 xtitle(Divorce cohort,size(small)) ytitle(Proportion good contact,size(small))  plot1opts(msymbol(S)) plot2opts(msymbol(O)) plot3opts(msymbol(B)) 
graph export "/users/matthijskalmijn/dropbox/og trends/figure3-final.pdf", replace as(pdf)