***************************************************
*Date created: 11 Oct, 2016
*Purpose: to provide the analysis program to DR 
*note: Stata version 13.1 
****************************************************

*set the file location
local mypath="C:\DR\Reg"  

global ddo="`mypath'\do"
global dlog="`mypath'\log"
global dgph="`mypath'\gph"

global ddta="L:\LAD"   


*****************************************************

set more off
set linesize 105
capture log close

log using "${dlog}\analysis.log", replace 

use "${ddta}\cohort1992.dta", clear
qui append using "${ddta}\cohort2002.dta"


/* variable label====================
lin__i		"LAD individual ID"
cohort		"1992 and 2002"
time		"time from separation"
marstat 	"1=married  2=cohabiting"
Quebec		"1=Quebec 0=other provinces"
xtircf		"before-tax total family income" 
aftaxf		"after-tax total family income" 
eqincbt		"before-tax family-equivalent total income"
eqincat		"after-tax family-equivalent total income"
age__i		"age"
agecat		"age category: 1=25-29 2=30-34 3=35-39 4=40-44" 
work		"flag_working"
tkid0_17	"total number of children under 18"
gnumkid		"tkid0_17 category: 0=no child 1=one child 2=two children 3=three children 4=four and more children" 
recup  		"flag_re-partnered"
almi_f		"alimony income received, family"
wgt2_i		"LAD 20% sample weights"
======================================== */

tab time cohort [iw=wgt2_i], mis

tab marstat Quebec if cohort==1992 & time==0 [iw=wgt2_i], mis
tab marstat Quebec if cohort==2002 & time==0 [iw=wgt2_i], mis


*=================================================
* Figure 1  
*=================================================
//=== Replacement Rate (RR=current_inc/baseyr_inc) ===* 
*1) total income before tax
by lin__i cohort (time), sort: gen bftotinc0=eqincbt[1] 
gen bftotinc_rr=eqincbt/bftotinc0  if bftotinc0>0

*2) total income after tax 
by lin__i cohort (time), sort: gen aftotinc0=eqincat[1] 
gen aftotinc_rr=eqincat/aftotinc0  if aftotinc0>0

foreach T in bf af { 
			
			preserve			
			keep if `T'totinc0>0    // all years
			 			
			collapse `T'totinc_rr (median) med`T'totinc_rr=`T'totinc_rr [iw=wgt2_i], by(cohort marstat time Quebec) 
			
			foreach C of num 1992  2002 {								
						tabdisp time marstat if cohort==`C' & Quebec, c(med`T'totinc_rr)  format(%9.2f)
						tabdisp time marstat if cohort==`C' & !Quebec, c(med`T'totinc_rr)  format(%9.2f)	
								
						for new m`C'b1-m`C'b2 \ num 1 2 : qui gen X=med`T'totinc_rr if marstat==Y  & Quebec  & cohort==`C' 
						for new m`C'b3-m`C'b4 \ num 1 2 : qui gen X=med`T'totinc_rr if marstat==Y  & !Quebec & cohort==`C' 
						
						scatter m`C'b1-m`C'b4 time, ///  
								c(l l l l) m(o t oh th) sort xlab(0(1)5) ylab() ///
								legend(lab(1 "married_Que") lab(2 "comlaw_Que") lab(3 "married_oth")  lab(4 "comlaw_oth") c(4)) /// lab(1 "married") lab(2 "comlaw")
								yti("") xti("Time") /// 
								ti("Cohort`C' (median)", box bexpand bc(ltblue)) ///
								name(coht`C'_med)
								*saving("${dgph}\coht`C'_`T'RR_med.gph", replace)
			}
									
			gr combine coht1992_med  coht2002_med,  ///
				iscale(.5) cap("`T' tax family income replacement rate(Median)", position(11))  
				saving("${dgph}\bytcoht_`T'RR_med_w.gph", replace)
			
			graph drop _all
					
			restore 
}		

gr combine "${dgph}\bytcoht_bfRR_med_w.gph" "${dgph}\bytcoht_afRR_med_w.gph", cols(1) ///
	saving("${dgph}\bytcoht_RR_med.gph", replace)


*====================================================================
*  Table 1 and 2a  : sample means and T-tests
*====================================================================

capture drop g_*
for var agecat gnumkid recup: qui tab X, gen(g_X) mis   
			
global tvar0 "g_agecat1 g_agecat2 g_agecat3 g_agecat4 age__i work g_gnumkid1 g_gnumkid2 g_gnumkid3 g_gnumkid4 g_gnumkid5 tkid0_17 eqincbt eqincat"
global tvar5 "work g_gnumkid1 g_gnumkid2 g_gnumkid3 g_gnumkid4 g_gnumkid5 tkid0_17 g_recup2"


// ==Mean and T test - first combination=====================================
capture mat drop _all
capture scalar drop _all
capture itab clear

mark married if marstat==1
tab marstat married, mis

foreach c of num 1992 2002 {
	foreach y of num 1 0 {

		qui tabstat $tvar0 if time==0 & cohort==`c' & Quebec==`y' & marstat==1 [aw=wgt2_i] , stats(mean) col(stat) save f(%12.4g) 
		mat Mean1_`y'_`c' =r(StatTotal)'
		qui tabstat $tvar0 if time==0 & cohort==`c' & Quebec==`y' & marstat==2 [aw=wgt2_i] , stats(mean) col(stat) save f(%12.4g) 
		mat Mean2_`y'_`c' =r(StatTotal)'
		mat Mean_`y'_`c'=(., .\ Mean1_`y'_`c', Mean2_`y'_`c')
		
		qui tabstat $tvar5 if time==5 & cohort==`c' & Quebec==`y' & marstat==1 [aw=wgt2_i] , stats(mean) col(stat) save f(%12.4g) 
		mat Mean1_`y'_`c'_5 =r(StatTotal)'
		qui tabstat $tvar5 if time==5 & cohort==`c' & Quebec==`y' & marstat==2 [aw=wgt2_i] , stats(mean) col(stat) save f(%12.4g) 
		mat Mean2_`y'_`c'_5 =r(StatTotal)'
		mat Mean_`y'_`c'=(Mean_`y'_`c'\ ., .\ Mean1_`y'_`c'_5, Mean2_`y'_`c'_5)
			
		mat Ttest_`y'_`c'=(., ., ., .)
		foreach x of global tvar0 {
			qui reg `x' i.married if time==0 & cohort==`c' & Quebec==`y' [iw=wgt2_i]
			mat T=r(table)'
			mat ST=T[2, 1..4]
			mat rownames ST=`x'
			mat Ttest_`y'_`c'=(Ttest_`y'_`c'\ST) 
		}	
		
		mat Ttest_`y'_`c'=(Ttest_`y'_`c'\ ., ., ., .)
		foreach x of global tvar5 {
			qui reg `x' i.married if time==5 & cohort==`c' & Quebec==`y' [iw=wgt2_i]
			mat T=r(table)'
			mat ST=T[2, 1..4]
			mat rownames ST=`x'
			mat Ttest_`y'_`c'=(Ttest_`y'_`c'\ST) 
		}	
		
		mat Meantest_`y'_`c'=(Mean_`y'_`c', Ttest_`y'_`c') 		
		mat colnames Meantest_`y'_`c'=married_`y'_`c' cohabiting_`y'_`c' diff std_err t_stat p_val
		mat rownames Meantest_`y'_`c'=Var_Time0 age25-29 age30-34 age35-39 age40-44 meanAGE work numkid0 numkid1 numkid2 numkid3 numkid4+ meanKID bt_faminc at_faminc Var_Time5 work numkid0 numkid1 numkid2 numkid3 numkid4+ meanKID recup1
		mat list Meantest_`y'_`c', format(%12.3f)
	}
}

mat FstMeantest=(Meantest_1_1992, Meantest_0_1992, Meantest_1_2002, Meantest_0_2002)  
mat list FstMeantest, format(%12.3f)


// ==Mean and T test - second combination=====================================
capture mat drop _all
capture scalar drop _all

foreach c of num 1992 2002 {
	foreach y of num 1 2 {

		qui tabstat $tvar0 if time==0 & cohort==`c' & marstat==`y' & Quebec==1 [aw=wgt2_i] , stats(mean) col(stat) save f(%12.4g) 
		mat Mean1_`y'_`c' =r(StatTotal)'
		qui tabstat $tvar0 if time==0 & cohort==`c' & marstat==`y' & Quebec==0 [aw=wgt2_i] , stats(mean) col(stat) save f(%12.4g) 
		mat Mean2_`y'_`c' =r(StatTotal)'
		mat Mean_`y'_`c'=(., .\ Mean1_`y'_`c', Mean2_`y'_`c')
		
		qui tabstat $tvar5 if time==5 & cohort==`c' & marstat==`y' & Quebec==1 [aw=wgt2_i] , stats(mean) col(stat) save f(%12.4g) 
		mat Mean1_`y'_`c'_5 =r(StatTotal)'
		qui tabstat $tvar5 if time==5 & cohort==`c' & marstat==`y' & Quebec==0 [aw=wgt2_i] , stats(mean) col(stat) save f(%12.4g) 
		mat Mean2_`y'_`c'_5 =r(StatTotal)'
		mat Mean_`y'_`c'=(Mean_`y'_`c'\ ., .\ Mean1_`y'_`c'_5, Mean2_`y'_`c'_5)

		
		mat Ttest_`y'_`c'=(., ., ., .)
		foreach x of global tvar0 {
			qui reg `x' i.Quebec if time==0 & cohort==`c' & marstat==`y' [iw=wgt2_i]
			mat T=r(table)'
			mat ST=T[2, 1..4]
			mat rownames ST=`x'
			mat Ttest_`y'_`c'=(Ttest_`y'_`c'\ST) 
		}	
		
		mat Ttest_`y'_`c'=(Ttest_`y'_`c'\ ., ., ., .)
		foreach x of global tvar5 {
			qui reg `x' i.Quebec if time==5 & cohort==`c' & marstat==`y' [iw=wgt2_i]
			mat T=r(table)'
			mat ST=T[2, 1..4]
			mat rownames ST=`x'
			mat Ttest_`y'_`c'=(Ttest_`y'_`c'\ST) 
		}	
		
		mat Meantest_`y'_`c'=(Mean_`y'_`c', Ttest_`y'_`c') 		
		mat colnames Meantest_`y'_`c'=Quebec_`y'_`c' oth_`y'_`c' diff std_err t_stat p_val
		mat rownames Meantest_`y'_`c'=Var_Time0 age25-29 age30-34 age35-39 age40-44 meanAGE work numkid0 numkid1 numkid2 numkid3 numkid4+ meanKID bt_faminc at_faminc Var_Time5 work numkid0 numkid1 numkid2 numkid3 numkid4+ meanKID recup1 
		mat list Meantest_`y'_`c', format(%12.3f)
	}
}

mat SndMeantest=(Meantest_1_1992, Meantest_2_1992, Meantest_1_2002, Meantest_2_2002)  
mat list SndMeantest, format(%12.3f)


// ==Mean and T test - third combination=====================================
capture mat drop _all
capture scalar drop _all

mark c1992 if cohort==1992
tab cohort c1992, mis

foreach c of num 1 0 {
	foreach y of num 1 2 {

		qui tabstat $tvar0 if time==0 & Quebec==`c' & marstat==`y' & cohort==1992 [aw=wgt2_i] , stats(mean) col(stat) save f(%12.4g) 
		mat Mean1_`y'_`c' =r(StatTotal)'
		qui tabstat $tvar0 if time==0 & Quebec==`c' & marstat==`y' & cohort==2002 [aw=wgt2_i] , stats(mean) col(stat) save f(%12.4g) 
		mat Mean2_`y'_`c' =r(StatTotal)'
		mat Mean_`y'_`c'=(., .\ Mean1_`y'_`c', Mean2_`y'_`c')
		
		qui tabstat $tvar5 if time==5 & Quebec==`c' & marstat==`y' & cohort==1992 [aw=wgt2_i] , stats(mean) col(stat) save f(%12.4g) 
		mat Mean1_`y'_`c'_5 =r(StatTotal)'
		qui tabstat $tvar5 if time==5 & Quebec==`c' & marstat==`y' & cohort==2002 [aw=wgt2_i] , stats(mean) col(stat) save f(%12.4g) 
		mat Mean2_`y'_`c'_5 =r(StatTotal)'
		mat Mean_`y'_`c'=(Mean_`y'_`c'\ ., .\ Mean1_`y'_`c'_5, Mean2_`y'_`c'_5)

		mat Ttest_`y'_`c'=(., ., ., .)
		foreach x of global tvar0 {
			qui reg `x' i.c1992 if time==0 & Quebec==`c' & marstat==`y' [iw=wgt2_i]
			mat T=r(table)'
			mat ST=T[2, 1..4]
			mat rownames ST=`x'
			mat Ttest_`y'_`c'=(Ttest_`y'_`c'\ST) 
		}	
		
		mat Ttest_`y'_`c'=(Ttest_`y'_`c'\ ., ., ., .)
		foreach x of global tvar5 {
			qui reg `x' i.c1992 if time==5 & Quebec==`c' & marstat==`y' [iw=wgt2_i]
			mat T=r(table)'
			mat ST=T[2, 1..4]
			mat rownames ST=`x'
			mat Ttest_`y'_`c'=(Ttest_`y'_`c'\ST) 
		}	

		mat Meantest_`y'_`c'=(Mean_`y'_`c', Ttest_`y'_`c') 		
		mat colnames Meantest_`y'_`c'=c1992_`y'_`c' c2002_`y'_`c' diff std_err t_stat p_val
		mat rownames Meantest_`y'_`c'=Var_Time0 age25-29 age30-34 age35-39 age40-44 meanAGE work numkid0 numkid1 numkid2 numkid3 numkid4+ meanKID bt_faminc at_faminc Var_Time5 work numkid0 numkid1 numkid2 numkid3 numkid4+ meanKID recup1 
		mat list Meantest_`y'_`c', format(%12.3f)
	}
}

mat TrdMeantest=(Meantest_1_1, Meantest_2_1, Meantest_1_0, Meantest_2_0)  
mat list TrdMeantest, format(%12.3f)


*===============================================================================
* Table 2b: mean and SD
*===============================================================================
for num 1992 2002: sum eqincbt if time==5 & marstat==1 & cohort==X [iw=wgt2_i] \ scalar mean1_X = r(mean) \ scalar sd1_X = r(sd)
for num 1992 2002: sum eqincbt if time==5 & marstat==2 & cohort==X [iw=wgt2_i] \ scalar mean2_X = r(mean) \ scalar sd2_X = r(sd)
capture {
	mat define tab=(mean1_1992, mean2_1992, mean1_2002, mean2_2002\ sd1_1992, sd2_1992, sd1_2002, sd2_2002)
	mata: tab=st_matrix("tab")
	mata: tab=round(tab, 100)
	mata: st_matrix("tab", tab)
	mat rownames tab=eqincbtMean eqincbtSD
	mat canada=canada\tab
}
mat list tab

foreach Y of num 0/1 {
	for num 1992 2002: sum eqincbt if time==5 & marstat==1 & cohort==X & Quebec==`Y' [iw=wgt2_i] \ scalar mean1_X = r(mean) \ scalar sd1_X = r(sd)
	for num 1992 2002: sum eqincbt if time==5 & marstat==2 & cohort==X & Quebec==`Y' [iw=wgt2_i] \ scalar mean2_X = r(mean) \ scalar sd2_X = r(sd)
	mat define tab`Y'=(mean1_1992, mean2_1992, mean1_2002, mean2_2002\ sd1_1992, sd2_1992, sd1_2002, sd2_2002)
	mata: tab`Y'=st_matrix("tab`Y'")
	mata: tab`Y'=round(tab`Y', 100)
	mata: st_matrix("tab`Y'", tab`Y')
	mat rownames tab`Y'=eqincbtMean eqincbtSD
	mat Quebec`Y'=Quebec`Y'\tab`Y'
	mat list tab`Y'
}


for num 1992 2002: sum eqincat if time==5 & marstat==1 & cohort==X [iw=wgt2_i] \ scalar mean1_X = r(mean) \ scalar sd1_X = r(sd)
for num 1992 2002: sum eqincat if time==5 & marstat==2 & cohort==X [iw=wgt2_i] \ scalar mean2_X = r(mean) \ scalar sd2_X = r(sd)
capture {
	mat define tab=(mean1_1992, mean2_1992, mean1_2002, mean2_2002\ sd1_1992, sd2_1992, sd1_2002, sd2_2002)
	mata: tab=st_matrix("tab")
	mata: tab=round(tab, 100)
	mata: st_matrix("tab", tab)
	mat rownames tab=eqincatMean eqincatSD
	mat canada=canada\tab
}
mat list tab

foreach Y of num 0/1 {
	for num 1992 2002: sum eqincat if time==5 & marstat==1 & cohort==X & Quebec==`Y' [iw=wgt2_i] \ scalar mean1_X = r(mean) \ scalar sd1_X = r(sd)
	for num 1992 2002: sum eqincat if time==5 & marstat==2 & cohort==X & Quebec==`Y' [iw=wgt2_i] \ scalar mean2_X = r(mean) \ scalar sd2_X = r(sd)
	mat define tab`Y'=(mean1_1992, mean2_1992, mean1_2002, mean2_2002\ sd1_1992, sd2_1992, sd1_2002, sd2_2002)
	mata: tab`Y'=st_matrix("tab`Y'")
	mata: tab`Y'=round(tab`Y', 100)
	mata: st_matrix("tab`Y'", tab`Y')
	mat rownames tab`Y'=eqincatMean eqincatSD
	mat Quebec`Y'=Quebec`Y'\tab`Y'
	mat list tab`Y'
}

matrix list canada, format(%9.3g)
matrix list Quebec1, format(%9.3g)
matrix list Quebec0, format(%9.3g)


*===============================================================================
* Table 3: alimony income 
*===============================================================================
*proportion of women receiving alimony : mean(alimonyf)
*share of alimony income in before tax family income: mean(ali_share)

mark alimonyf if almi_f>0
gen double ali_share=almi_f/xtircf if xtircf>0  // only can be calculated for positive income

global tvar2 "alimonyf ali_share"

// ==Mean and T test - Fst combination=====================================
capture mat drop _all
capture scalar drop _all

foreach c of num 1992 2002 {
	foreach y of num 1 2 {

		qui tabstat $tvar2 if time==2 & cohort==`c' & marstat==`y' & Quebec==1 [aw=wgt2_i] , stats(mean) col(stat) save f(%12.4g) 
		mat Mean1_`y'_`c' =r(StatTotal)'
		qui tabstat $tvar2 if time==2 & cohort==`c' & marstat==`y' & Quebec==0 [aw=wgt2_i] , stats(mean) col(stat) save f(%12.4g) 
		mat Mean2_`y'_`c' =r(StatTotal)'
		mat Mean_`y'_`c'=(., .\ Mean1_`y'_`c', Mean2_`y'_`c')
				
		mat Ttest_`y'_`c'=(., ., ., .)
		foreach x of global tvar2 {
			qui reg `x' i.Quebec if time==2 & cohort==`c' & marstat==`y' [iw=wgt2_i]
			mat T=r(table)'
			mat ST=T[2, 1..4]
			mat rownames ST=`x'
			mat Ttest_`y'_`c'=(Ttest_`y'_`c'\ST) 
		}	
				
		mat Meantest_`y'_`c'=(Mean_`y'_`c', Ttest_`y'_`c') 		
		mat colnames Meantest_`y'_`c'=Quebec_`y' oth_`y' diff std_err t_stat p_val
		mat rownames Meantest_`y'_`c'=Time2_`c' alimonyf ali_share  
		mat list Meantest_`y'_`c', format(%12.3f)
	}
}
 
mat T3FstMeantest=(Meantest_1_1992, Meantest_2_1992 \ Meantest_1_2002, Meantest_2_2002)  
mat list T3FstMeantest, format(%12.3f)


// ==Mean and T test - Snd combination=====================================
capture mat drop _all
capture scalar drop _all

foreach c of num 1992 2002 {
	foreach y of num 1 0 {

		qui tabstat $tvar2 if time==2 & cohort==`c' & Quebec==`y' & marstat==1 [aw=wgt2_i] , stats(mean) col(stat) save f(%12.4g) 
		mat Mean1_`y'_`c' =r(StatTotal)'
		qui tabstat $tvar2 if time==2 & cohort==`c' & Quebec==`y' & marstat==2 [aw=wgt2_i] , stats(mean) col(stat) save f(%12.4g) 
		mat Mean2_`y'_`c' =r(StatTotal)'
		mat Mean_`y'_`c'=(., .\ Mean1_`y'_`c', Mean2_`y'_`c')
		
			
		mat Ttest_`y'_`c'=(., ., ., .)
		foreach x of global tvar2 {
			qui reg `x' i.married if time==2 & cohort==`c' & Quebec==`y' [iw=wgt2_i]
			mat T=r(table)'
			mat ST=T[2, 1..4]
			mat rownames ST=`x'
			mat Ttest_`y'_`c'=(Ttest_`y'_`c'\ST) 
		}	
		
			
		mat Meantest_`y'_`c'=(Mean_`y'_`c', Ttest_`y'_`c') 		
		mat colnames Meantest_`y'_`c'=married_`y' cohabiting_`y' diff std_err t_stat p_val
		mat rownames Meantest_`y'_`c'=Time2_`c' alimonyf ali_share  
		mat list Meantest_`y'_`c', format(%12.3f)

	}
}
 
mat T3SndMeantest=(Meantest_1_1992, Meantest_0_1992\ Meantest_1_2002, Meantest_0_2002)  
mat list T3SndMeantest, format(%12.3f)

capture mat drop _all
capture scalar drop _all
*====================================================================================
* set up for a fixed effect model
*====================================================================================
capture estimates clear 

mark sepcom if marstat==2 & time>0			
lab var sepcom "a flag for sep from comlaw"	
tab time sepcom if cohort==1992, mis

mark sepcom2 if marstat==2

mark everptn if recup==1 & time>0
lab var everptn "a flag for ever re-partnered"	
tab sepcom everptn if cohort==1992, mis


*===================================================================
* Table 4 
*===================================================================

for num 1992  2002: xtreg eqincat ib0.time##ib1.marstat if cohort==X  & Quebec, fe \ est sto M1_at_trX_Que
for num 1992  2002: xtreg eqincat ib0.time##ib1.marstat if cohort==X  & !Quebec, fe \ est sto M1_at_trX_oth


for num 1992  2002: xtreg eqincat ib0.time##ib1.marstat i.gnumkid i.work i.everptn if cohort==X  & Quebec, fe \ est sto T4_at_trX_Que
for num 1992  2002: xtreg eqincat ib0.time##ib1.marstat i.gnumkid i.work i.everptn if cohort==X  & !Quebec, fe \ est sto T4_at_trX_oth


*==========================
* with interactions
*==========================
*count persons repeating in cohort 2002 and 1992, lin__i is not unique when time==0
duplicates report lin__i if time==0 
duplicates tag lin__i if time==0, gen(tag)
by lin__i:  egen duptag=max(tag)
tab duptag time , mis
tab duptag marstat if time==0 & cohort==1992 [iw=wgt2_i], mis col
tab duptag marstat if time==0 & cohort==2002 [iw=wgt2_i], mis col

tab cohort marstat if time==0 & duptag==1 [iw=wgt2_i], mis col

egen pid=group(lin__i cohort)
duplicates report pid cohort
xtset pid time

*1) time##cohabiting##Quebec 
for num 1992 2002: xtreg eqincat ib0.time##i.marstat##i.Quebec i.gnumkid i.work i.everptn if cohort==X, fe  \ est sto IntQue_at_trX

*2) time##cohabiting##cohort
for num 0 1: xtreg eqincat ib0.time##i.marstat##i.cohort i.gnumkid i.work i.everptn if Quebec==X, fe cluster(lin__i)  \ est sto Intcoh_at_trX

*3) time##cohabiting##Quebec##cohort
xtreg eqincat ib0.time##i.marstat##i.Quebec##i.cohort i.gnumkid i.work i.everptn, fe cluster(lin__i)  
est sto Intall_at_tr


// estimation out=======================================================

*1) Table 4 without other controls
estout M1_at_tr1992_Que M1_at_tr1992_oth M1_at_tr2002_Que M1_at_tr2002_oth , ///
       cells("b(star fmt(%9.2f) drop(*0b.*,relax)) se(par fmt(%-9.2fc) drop(*0b.*,relax))") stats(sigma_u sigma_e rho N, fmt(2 2 2 0) labels(sigma_u sigma_e rho N)) ///
       starlevels(^ 0.1 * 0.05 ** 0.01 *** 0.001)  varl( _cons constant) ///
       margin legend collabels(coef. Std.Err.) ml() varwidth(20) 


*2) Table 4 with other controls
estout T4_at_tr1992_Que T4_at_tr1992_oth T4_at_tr2002_Que T4_at_tr2002_oth , ///
       cells("b(star fmt(%9.2f) drop(*0b.*,relax)) se(par fmt(%-9.2fc) drop(*0b.*,relax))") stats(sigma_u sigma_e rho N, fmt(2 2 2 0) labels(sigma_u sigma_e rho N)) ///
       starlevels(^ 0.1 * 0.05 ** 0.01 *** 0.001)  varl( _cons constant) ///
       margin legend collabels(coef. Std.Err.) ml() varwidth(20) 

*3) interactions 

estout IntQue_at_tr1992 IntQue_at_tr2002, ///
       cells("b(star fmt(%9.2f) drop(*0b.*,relax)) se(par fmt(%-9.2fc) drop(*0b.*,relax))") stats(sigma_u sigma_e rho N, fmt(2 2 2 0) labels(sigma_u sigma_e rho N)) ///
       starlevels(^ 0.1 * 0.05 ** 0.01 *** 0.001)  varl( _cons constant) ///
       margin legend collabels(coef. Std.Err.) ml() varwidth(20)  


estout Intcoh_at_tr1 Intcoh_at_tr0 , ///
       cells("b(star fmt(%9.2f) drop(*0b.*,relax)) se(par fmt(%-9.2fc) drop(*0b.*,relax))") stats(sigma_u sigma_e rho N, fmt(2 2 2 0) labels(sigma_u sigma_e rho N)) ///
       starlevels(^ 0.1 * 0.05 ** 0.01 *** 0.001)  varl( _cons constant) ///
       margin legend collabels(coef. Std.Err.) ml() varwidth(20)  


estout Intall_at_tr , ///
       cells("b(star fmt(%9.2f) drop(*0b.*,relax)) se(par fmt(%-9.2fc) drop(*0b.*,relax))") stats(sigma_u sigma_e rho N, fmt(2 2 2 0) labels(sigma_u sigma_e rho N)) ///
       starlevels(^ 0.1 * 0.05 ** 0.01 *** 0.001)  varl( _cons constant) ///
       margin legend collabels(coef. Std.Err.) ml() varwidth(20)  
	   
estimates clear
log close