* Excel file to record influences of individual variables -- this part just calculates difference between US and OECD mean; later add regression coefficients
use USA oecd year $notimputable $imputable if year==2005 & oecd using paneldata.dta, clear
collapse $notimputable $imputable, by(USA)
foreach var of varlist $notimputable $imputable {
	replace `var'=-`var' if !USA
	}
collapse (sum) $notimputable $imputable
export excel using results/contributions.xlsx, replace firstrow(var) cell(B1)
putexcel set results/contributions.xlsx, modify
putexcel A2=("USminusOECD") A3=("hom") A4=("pris") A5=("d_hom") A6=("d_pris") // one actually needs to add the calculations of row 5 and 6 by hand


* MI estimates (incl saving contributions)
use MIdataset.dta, clear

_eststo NRviol3: mi estimate, post imputations(1/`=$imputations'): regress ln_IHMEviolentdeaths USA $notimputable $imputable
_eststo viol3: mi estimate, post imputations(1/`=$imputations'): regress ln_IHMEviolentdeaths USA $notimputable $imputable, robust
matrix results=e(b)
matrix results=results[1,2..`=1+`: word count $imputable' + `: word count $notimputable'']
putexcel B3=matrix(results)
test $notimputable $imputable
estadd scalar pF=r(p)

_eststo NRpris3: mi estimate, post imputations(`=1+$imputations'/`=2*$imputations'): regress i_ln_pris USA $notimputable $imputable
_eststo pris3: mi estimate, post imputations(`=1+$imputations'/`=2*$imputations'): regress i_ln_pris USA $notimputable $imputable, robust
matrix results=e(b)
matrix results=results[1,2..`=1+`: word count $imputable' + `: word count $notimputable'']
putexcel B4=matrix(results)
test $notimputable $imputable
estadd scalar pF=r(p)

* "naive" OLS (casewise deletion)
mi extract 0, clear

_eststo NRviol1: reg ln_IHME USA $notimputable $imputable
_eststo viol1: reg ln_IHME USA $notimputable $imputable, robust
test $notimputable $imputable
estadd scalar pF=r(p)
_eststo NRpris1: reg i_ln_pris USA $notimputable $imputable
_eststo pris1: reg i_ln_pris USA $notimputable $imputable, robust
test $notimputable $imputable
estadd scalar pF=r(p)

* FIML

_eststo NRviol4: sem (ln_IHMEviolentdeaths <- USA $notimputable $imputable) if mi(ln_IHME)==0, method(mlmv) nocapslatent
_eststo viol4: sem (ln_IHMEviolentdeaths <- USA $notimputable $imputable) if mi(ln_IHME)==0, method(mlmv) nocapslatent vce(robust)
test $notimputable $imputable
estadd scalar pF=r(p)
_eststo NRpris4: sem (i_ln_pris <- USA $notimputable $imputable) if mi(i_ln_pris)==0, method(mlmv) nocapslatent
_eststo pris4: sem (i_ln_pris <- USA $notimputable $imputable) if mi(i_ln_pris)==0, method(mlmv) nocapslatent vce(robust)
test $notimputable $imputable
estadd scalar pF=r(p)

* OLS+ 

local imputed ""
local complete ""
foreach var of varlist $imputable {
	quietly misstable sum `var'
	if r(N_eq_dot)!=. local imputed "`imputed' `var'"
	else local complete "`complete' `var'"
	}
foreach indepvar in `imputed' {
	gen md_`indepvar'=mi(`indepvar')
	recode `indepvar' (.=0)
	}


_eststo NRviol2: reg ln_IHME USA $notimputable $imputable md*
_eststo viol2: reg ln_IHME USA $notimputable $imputable md*, robust
test $notimputable $imputable
estadd scalar pF=r(p)
_eststo NRpris2: reg i_ln_pris USA $notimputable $imputable md*
_eststo pris2: reg i_ln_pris USA $notimputable $imputable md*, robust
test $notimputable $imputable
estadd scalar pF=r(p)


*** Scatter plots

tempvar residpris
tempvar residhom
tempvar residhomUS

quietly: reg ln_IHME $notimputable $imputable md*, robust
predict `residhom', residuals
gen `residhomUS'=`residhom' if USA
replace `residhom'=. if USA
quietly: reg i_ln_pris $notimputable $imputable md*, robust
predict `residpris', residuals

scatter `residhom' `residhomUS' `residpris', msymbol(i i) mlabel(ISO ISO) mlabgap(*-.5 *-.5) mlabcolor(gray red) mlabsize(. *1.5) xtitle("log-residual incarceration rate") legend(off) ///
	ytitle("log-residual homicide rate") title("Log-residual homicide and incarceration rates around the world", size(automatic)) subtitle("2005") ///
	note("Residuals from models equivalent to models 2 and 6 of table 3 except that they do not" "contain a separate US dummy.") scheme(s1color) ///
	saving(results/figure3, replace)

drop `residhom' `residpris'

sum i_ln_pris if mi(ln_IHMEviolentdeath)==0
scalar lowerp=r(min)
scalar upperp=r(max)
sum ln_IHMEviolentdeath if mi(i_ln_pris)==0
scalar lowerh=r(min)
scalar upperh=r(max)

tempvar hom OECDhom UShom
gen `hom'= ln_IHMEviolentdeath if !oecd
gen `OECDhom'=ln_IHMEviolentdeath if oecd & !USA
gen `UShom'=ln_IHMEviolentdeath if USA

scatter `hom' `OECDhom' `UShom' i_ln_pris if i_ln_pris>=lowerp & i_ln_pris<=upperp & ln_IHMEviolentdeath>=lowerh & ln_IHMEviolentdeath<=upperh, ///
	msymbol(i i i) mlabel(ISO ISO ISO) mlabgap(*-.5 *-.5 *-.5) mlabcolor(gray black red) mlabsize(. . *1.5) ///
	xtitle("Log incarceration rate")  ylabel(-1/5) ytitle("Log homicide rate") legend(off) ///
	title("Log homicide and incarceration rates around the world", size(automatic)) subtitle("log-rates per 100,000, 2005") ///
	note("Countries included: union of estimation samples of models 2 and 6 of table 3." "OECD countries in solid black." "Source: IHME-GBD (homicides), ICPS (incarceration).") ///
	scheme(s1color) saving(results/figure1, replace)

replace `hom'=exp(ln_IHMEviol) if !oecd
replace `OECDhom'=exp(ln_IHMEviol) if oecd & !USA
replace `UShom'=exp(ln_IHMEviol) if USA
replace ln_IHMEviolentdeath = exp(ln_IHMEviol)
replace i_ln_pris=exp(i_ln_pris) // I replace the old variables so I can just copy the graph code above; will clear the data set anyway after this

sum i_ln_pris if mi(ln_IHMEviolentdeath)==0
scalar lowerp=r(min)
scalar upperp=r(max)
sum ln_IHMEviolentdeath if mi(i_ln_pris)==0
scalar lowerh=r(min)
scalar upperh=r(max)

scatter `hom' `OECDhom' `UShom' i_ln_pris if i_ln_pris>=lowerp & i_ln_pris<=upperp & ln_IHMEviolentdeath>=lowerh & ln_IHMEviolentdeath<=upperh, ///
	msymbol(i i i) mlabel(ISO ISO ISO) mlabgap(*-.5 *-.5 *-.5) mlabcolor(gray black red) mlabsize(. . *1.5) ///
	xtitle("Prisoners (ICPS)")  ytitle("Homicides (IHME-GBD)") legend(off) ///
	title("Homicide and incarceration rates around the world") subtitle("levels per 100k inhabitants, 2005") ///
	note("Countries included: union of estimation samples of models 2 and 6 of table 2." "OECD countries in solid black") ///
	scheme(s1color) saving(results/figure1_alt_levels, replace)	

clear

*** regression output (problem: STATA's esttab chokes on the SEM output when included in the big table; couldn't fix that)

esttab viol1 viol2 viol3 pris1 pris2 pris3 using results/table3_OLS_OLSplus_MI, ///
	replace csv label b(2) se(2) r2(2) scalars(pF) sfmt(3) star(* 0.10 ** 0.05 *** 0.01) nogaps nonumbers drop(md* _cons) order(USA $notimputable $imputable) ///
	title(Table: Homicides and Incarceration) ///
	mtitles(OLS OLS+ MI OLS OLS+ MI) mgroups("Homicide rate (2005)" "Incarceration rate (2005)", pattern(1 0 0  1 0 0 )) ///
	coeflabels($variabledescriptions) nonotes ///
	addnotes("Robust standard errors in parentheses." "OLS+ sets missing values to zero and adds a dummy variable equal to one for those observations." ///
		"MI uses multiple imputation to fill in missing values. FIML derives the joint likelihood.")

esttab NRviol1 NRviol2 NRviol3 NRpris1 NRpris2 NRpris3 using results/table3_nonrobust, replace csv label b(2) se(2) r2(2) star(* 0.10 ** 0.05 *** 0.01) ///
	mtitles(OLS OLS+ MI OLS OLS+ MI) mgroups("Homicide rate (2005)" "Incarceration rate (2005)", pattern(1 0 0  1 0 0 )) ///
	nogaps nonumbers drop(md* _cons) order(USA $notimputable $imputable) coeflabels($variabledescriptions) nonotes ///
	addnotes("non-robust standard errors")

		
foreach depvar in viol pris {
esttab `depvar'4 NR`depvar'4 using results/table3_FIML`depvar', ///
	replace csv label b(2) se(2) scalars(pF) sfmt(3) star(* 0.10 ** 0.05 *** 0.01) nogaps nonumbers drop( _cons) ///
	mtitles(FIML) coeflabels($variabledescriptions)
	}
