clear
estimates clear
use state_dataset.dta, clear
replace hhinc00=hhinc00/1000
local intllist "hhinc00 cs_race_bla gini cs_fam_wkidsinglemom"
local betterlist "hhinc00 cs00_seg_inc_pov25 gini99 scap_ski90pcm" // consider using cs_race_theil_2000 cs00_seg_inc
local combinedlist "hhinc00 cs_race_bla cs00_seg_inc_pov25 gini99 scap_ski90pcm"
pwcorr `intllist' `betterlist'

*************************************************************
****** state-level regressions for discussion part **********

foreach depvar in ln_inc ln_murderrate2000 crime {
	foreach controls in intllist betterlist combinedlist {
		_eststo: reg `depvar' ``controls'' if !mi(state_id) // the if statement makes sure that the US total is not in the regression
		test ``controls''
		estadd scalar pF=r(p)
		}
	}

esttab * using results/stateregs, replace csv label b(2) se(2) r2(2) ar2(2) scalar(pF) sfmt(2) star(* 0.10 ** 0.05 *** 0.01) nogaps nonumbers nonotes drop(_cons)

*************************************************************
***** States vs. countries: scatter plot *******************

* putting together data set
keep state stateabbrv ln_incarceration2003 ln_murderrate2005
tempfile state
save `state'

use ISO3 country ln_IHMEviol oecd year if year==2005 & oecd==1 using paneldata.dta, clear
drop oecd year
tempfile countries
save `countries'

use ISO3 country i_ln_pris oecd year if year==2003 & oecd==1 using paneldata.dta, clear
drop oecd year
merge 1:1 ISO3 country using `countries', nogenerate

append using `state'

label var i_ln_pris "Incarceration rate, 2003 (ICPS, logs)"
label var ln_IHMEviol "Violent death rate, 2005 (IHME, logs)"
label var ln_incarceration2003 "Incarceration rate, 2003 (Sourcebook, logs)"
label var ln_murderrate2005 "Murder rate, 2004-06 avg. (UCR, logs)"

* adjusting for reporting differences between IHME (doctors' reports!) and UCR (police reports)
sum ln_IHME if ISO=="USA"
local IHME=r(mean)
sum ln_murder if stateabbrv=="USA" // note: this is the log murderrate (UCR) for the entire US, i.e., "sum" here just references a single observation (the US total)
gen adj_ln_murderrate2005= ln_murderrate2005+`IHME'-`r(mean)'
label var adj_ln_murderrate2005 "ln_murderrate plus US difference IHME-UCR"


* scatter output
twoway (scatter ln_IHMEviol i_ln_pris if ISO!="USA", msymbol(i) mlabel(ISO) mlabgap(*-.5)) (scatter adj_ln_murder ln_inc, msymbol(i) mlabel(stateabbrv) mlabgap(*-.5)) ///
	(lfit ln_IHMEviol i_ln_pris if ISO!="USA", lcolor(green)) (lfit adj_ln_murder ln_inc, lcolor(red)), ///
	xtitle("Log incarceration rate")  ytitle("Log homicide rate") legend(off) ///
	title("US states vs OECD countries") subtitle("log-rates per 100,000; 2005/2003") ///
	note("Source: IHME-GBD and ICPS (OECD); UCR and BJS (US states). Linear trend lines.") ///
	scheme(s1color) saving(results/figure2.gph, replace)
