
/***********************/
/*** 1-line ES graph ***/
/***********************/

capture program drop ES_graph

program define ES_graph

set more off

syntax, ci(real) [ xti(string) yti(string) lab1(string) outname(string) rawti(string) ]

preserve

gen fig_upper=b+`ci'*se
gen fig_lower=b-`ci'*se
gen fig_t = t-101
rename b fig_b

twoway (scatter fig_b fig_t if  fig_t<6 & fig_t>-5, c(l) xline(-0.0, lc(gs11) lp(shortdash)) yline(0, lc(gs11) lp(shortdash))) ///
(rcap fig_upper fig_lower fig_t  if fig_t<6 & fig_t>-5 , lpattern(solid) lcolor(dknavy) c(l) m(i) ), ///
plotregion(fcolor(white) lcolor(white)) ///
graphregion(fcolor(white) lcolor(white)) ///
xtitle("`xti'") ///
ytitle("`yti'") ///
xlab(-4(1)5) ///
title("`rawti'", color(black)) ///
legend(label(1 "`lab1'") label(2 "95% CI"))
	
graph save "$figdir/raw/`outname'.gph", replace

restore

end


/***********************/
/*** 2-line ES graph ***/
/***********************/

capture program drop ES_graph2

program define ES_graph2

set more off
	
syntax, ci(real) [ xti(string) yti(string) outname(string) lab1(string) lab2(string) rawti(string) ]

preserve

gen fig_upper1=b1+`ci'*se1
gen fig_lower1=b1-`ci'*se1
gen fig_upper2=b2+`ci'*se2
gen fig_lower2=b2-`ci'*se2
gen fig_t = t-101
rename b1 fig_b1
rename b2 fig_b2

	
*Figure in gph format with title (only if specified with rawti option)
twoway (scatter fig_b1 fig_t if  fig_t<6 & fig_t>-5, c(l) lc(dknavy) mcolor(dknavy) msymbol(smtriangle) xline(-0.0, lc(gs11) lp(shortdash)) yline(0, lc(gs11) lp(shortdash))) ///
(rcap fig_upper1 fig_lower1 fig_t  if fig_t<6 & fig_t>-5 , lpattern(solid) lcolor(dknavy) c(l) m(i) ) ///
(scatter fig_b2 fig_t if  fig_t<6 & fig_t>-5, c(l) lpattern(dash) lc(maroon) mcolor(maroon)  msymbol(smdiamond) ) ///
(rcap fig_upper2 fig_lower2 fig_t  if fig_t<6 & fig_t>-5 , lpattern(dash) lcolor(maroon) c(l) m(i) ), ///
plotregion(fcolor(white) lcolor(white)) ///
graphregion(fcolor(white) lcolor(white)) ///
xtitle("`xti'") ///
ytitle("`yti'") ///
xlab(-4(1)5) ///
title("`rawti'", color(black)) ///
legend(label(1 "`lab1'") label(3 "`lab2'") label(2 "95% CI") label(4 "95% CI"))
	
graph save "$figdir/raw/`outname'.gph", replace

restore

end


/***********************/
/*** Program specifically for processing the output from ES_interactions ***/
/***********************/

capture program drop ES_interactions_process
program define ES_interactions_process
syntax, fn(string)


import delim "$dumpdir/`fn'.csv", clear
drop if inlist(mod(_n, 11), 0, 1) /* Don't want binned time periods in graph, so remove these obs; NOTE THIS DEPENDS ON WHAT WE INCLUDE IN THE RAW FILE SO CHECK */
local new_obs = _N + 1
set obs `new_obs' /* Add one more obs so we have an additional time period for t = -1, which is normalized to 0 */

g t = _n - 5 if _n < 4
replace t = _n - 4 if missing(t)
replace t = -1 in `new_obs'
sort t
replace coef = 0 if missing(coef) /* Normalized period (t = 1 -> 0) */
replace stderr = 0 if missing(stderr)

g coef2 = coef[_n+10] if t < -1 /* Now moving the other group into its own column(s) */
replace coef2 = coef[_n+9] if t > -1
replace coef2 = 0 if t == -1

g stderr2 = stderr[_n+10] if t < -1 /* Now moving the other group into its own column(s) */
replace stderr2 = stderr[_n+9] if t > -1
replace stderr2 = 0 if t == -1

drop if t > 5 /* Now drop that group */
replace t = t + 101 /* For figure program compatibility */

*Now order as we want in figure
rename (coef coef2 stderr stderr2) ///
	   (b2 b1 se2 se1) /* Reordering so that we have variable that is part of the indicated group, then other */

end



/****************************************************************************************/
/*** Average impact of a state corporate tax cut on different outcomes (in time = t): ***/
/****************************************************************************************/
capture program drop Event_outcome
program define Event_outcome
syntax [if] [in], yvar(varname) event(varname)  [ controls(string) outname(string) ]

set more off

preserve

*Keep if specified
if "`if'" != "" {
	keep `if'
}

*Setup event
g e = year if `event' == 1
qui sum year
local years = `r(max)' - `r(min)' + 1

*Fill up:
local i = 0
while `i' < `years' {
	qui replace e = e[_n+1] if fips_state == fips_state[_n+1] & abs(year - e[_n+1]) < abs(year - e)
	local ++i
}
*Fill down: 
local i = 0
while `i' < `years' {
	qui replace e = e[_n-1] if fips_state == fips_state[_n-1] & abs(year - e[_n-1]) < abs(year - e)
	local ++i
}

*DEFINE EVENT
*binned
gen Dn5=year<=e-5 
gen Dn4=year==e-4
gen Dn3=year==e-3
gen Dn2=year==e-2
gen D0=year==e
gen D1=year==e+1
gen D2=year==e+2
gen D3=year==e+3
gen D4=year==e+4
gen D5=year==e+5
*binned
gen D6=year>=e+6


*Proper GDP weight (4 years prior to event; in cases with no event, take GDP in median year as weight)
sort fips_state year
bysort fips_state: gen GDP_wgt = GDP[_n-4] if e == year/* If event 4 years in the future, capture gdp */
*Conditional if year too soon (i.e. if event occurs in 1981, we don't have GDP prior to 1980, so capture 1980 GDP)
bysort fips_state: replace GDP_wgt = GDP[1] if e == year & missing(GDP_wgt)
*Now within that event (i.e. every year for which e == that), replace GDP_wgt as such; if no event, then take median year's GDP (assumed 1995, can generalize later)
sort fips_state year e
g e_id = e
replace e_id = 0 if missing(e_id)
egen tempgroup = group(fips_state e_id)
levelsof tempgroup, local(groups)
foreach group in `groups' {
	qui sum GDP_wgt if tempgroup == `group'
	if `r(N)' != 0 {
		qui replace GDP_wgt = `r(mean)' if tempgroup == `group'
	}
	else {
		qui sum GDP if tempgroup == `group' & year == 1995
		qui replace GDP_wgt = `r(mean)' if tempgroup == `group'
	}
}


*Base Controls (if specified)
if "`controls'" != "" {
	local x = 1
	foreach control in `controls' {
		g temp_e = year if ch_`control' != 0
		*Fill up:
		local i = 0
		while `i' < `years' {
			qui replace temp_e = temp_e[_n+1] if fips_state == fips_state[_n+1] & abs(year - temp_e[_n+1]) < abs(year - temp_e)
			local ++i
		}
		*Fill down: 
		local i = 0
		while `i' < `years' {
			qui replace temp_e = temp_e[_n-1] if fips_state == fips_state[_n-1] & abs(year - temp_e[_n-1]) < abs(year - temp_e)
			local ++i
		}
		*gen variable and increment
		g control`x' = year <= temp_e-5
		local ++x
		g control`x' = year == temp_e-4
		local ++x
		g control`x' = year == temp_e-3
		local ++x
		g control`x' = year == temp_e-2
		local ++x
		g control`x' = year == temp_e
		local ++x
		g control`x' = year == temp_e+1
		local ++x
		g control`x' = year == temp_e+2
		local ++x
		g control`x' = year == temp_e+3
		local ++x
		g control`x' = year == temp_e+4
		local ++x
		g control`x' = year == temp_e+5
		local ++x
		g control`x' = year >= temp_e+6
		local ++x
		
		drop temp_e
	}
	local --x
}


*MAIN REG
areg `yvar' Dn5-D6 i.fips_state, absorb(year) cluster(fips_state)

*Alternative reg with controls if specified:
if "`controls'" != "" {
	qui ds control*
	areg `yvar' Dn5-D6 control1-control`x' i.fips_state, absorb(year) cluster(fips_state)
}

*Raw regression output
regsave  
drop if substr(var,1,1) != "D"
drop N r2
export delimited $dumpdir/`outname'.csv, replace

restore

end


/*******************************/
/*** Kernal Density by decade **/
/*******************************/

capture program drop kdensity_decade

program kdensity_decade

set more off
	
syntax, var(varname) xmax(real) [ xti(string) yti(string) INCrement(real 1) outname(string) ]

preserve

keep if inlist(year, 1980, 1990, 2000, 2010)
drop if fips_state == 0 | fips_state == 11 | fips_state > 56 /* US, DC, Foreign */

*construct x axis
sort year fips_state /* This is necessary for x axis construction */
gen xaxis_max=`xmax'
gen id= mod(_n - 1 , 50) + 1
replace id=id/50
g xaxis=id* xaxis_max

tempfile data2
save `data2', replace

foreach year in 1980 1990 2000 2010 {
	use `data2', clear
	keep if year == `year'
	kdensity `var', gen(x_`var' kd_`var'_`year') n(50) at(xaxis) nograph
	keep fips_state year kd*
	tempfile data_`year'
	save `data_`year'', replace
}
use `data2'
merge 1:1 fips_state year using `data_1980', nogen
merge 1:1 fips_state year using `data_1990', nogen
merge 1:1 fips_state year using `data_2000', nogen
merge 1:1 fips_state year using `data_2010', nogen
	
twoway (scatter kd_`var'_1980 xaxis , c(l) msymbol(none) lcolor(black) lpattern(solid) lwidth(thick)) ///
	(scatter kd_`var'_1990 xaxis , c(l) msymbol(none) lcolor(black) lpattern(dash) lwidth(thick)) ///
	(scatter kd_`var'_2000 xaxis , c(l) msymbol(none) lcolor(gs9) lpattern(solid)) ///
	(scatter  kd_`var'_2010 xaxis, c(l) msymbol(none) lcolor(gs9) lpattern(dash)), ///
	graphregion(fcolor(white)) graphregion(color(white)) ///
	xtitle("`xti'", size(medlarge)) ///
	ytitle("`yti'", size(medlarge)) ///
	legend(order(1 2 3 4) size(medlarge) label(1 "1980") label(2 "1990") ///
	label(3 "2000") label(4 "2010")) ///
	ylab(, nogrid) ///
	xscale(range(0 `xmax') ) ///
	xlabel(0(`increment')`xmax')
	
graph export "$figdir/`outname'.pdf", replace

restore

end



/*******************************/
/*** Kernal Density by region **/
/*******************************/

capture program drop kdensity_region

program kdensity_region

set more off
	
syntax, var(varname) xmax(real) [ xti(string) yti(string) outname(string) ]

preserve

tempfile data
save `data', replace

*construct x axis
sort year fips_state /* This is necessary for x axis construction */
gen xaxis_max=`xmax'
tempfile data_temp
save `data_temp', replace
foreach region in 1 2 3 4 {
	use `data_temp', clear
	keep if region_broad == `region'
	sort `var'
	qui count
	gen id_`region'= mod(_n - 1 , `r(N)') + 1 
	replace id_`region'=id_`region'/`r(N)'
	g xaxis_`region'=id_`region'*xaxis_max
	tempfile data_temp`region'
	save `data_temp`region'', replace
}
use `data_temp', clear
merge 1:1 fips_state year using `data_temp1', nogen
merge 1:1 fips_state year using `data_temp2', nogen
merge 1:1 fips_state year using `data_temp3', nogen
merge 1:1 fips_state year using `data_temp4', nogen

tempfile data_temp
save `data_temp', replace

foreach region in 1 2 3 4 {
	use `data_temp', clear
	keep if region_broad == `region'
	kdensity `var', gen(x_`var' kd_`var'_`region') at(xaxis_`region') nograph
	keep fips_state region_broad year kd*
	tempfile data_`region'
	save `data_`region'', replace
}
use `data_temp', clear
merge 1:1 fips_state region_broad year using `data_1', nogen
merge 1:1 fips_state region_broad year using `data_2', nogen
merge 1:1 fips_state region_broad year using `data_3', nogen
merge 1:1 fips_state region_broad year using `data_4', nogen

gen xaxis = xaxis_1
replace xaxis = xaxis_2 if missing(xaxis)
replace xaxis = xaxis_3 if missing(xaxis)
replace xaxis = xaxis_4 if missing(xaxis)

sort `var' xaxis*
	
twoway (scatter kd_`var'_1 xaxis_1, c(l) msymbol(none) lcolor(black) lpattern(solid) lwidth(thick) ) ///
	(scatter kd_`var'_2 xaxis_2, c(l) msymbol(none) lcolor(black) lpattern(dash) lwidth(thick) ) ///
	(scatter kd_`var'_3 xaxis_3, c(l) msymbol(none) lcolor(gs9) lpattern(solid) ) ///
	(scatter  kd_`var'_4 xaxis_4, c(l) c(l) msymbol(none) lcolor(gs9) lpattern(dash) ), ///
	graphregion(fcolor(white)) graphregion(color(white)) ///
	xtitle(`xti', size(medlarge)) ///
	ytitle(`yti', size(medlarge)) ///
	legend(order(1 2 3 4) size(medlarge) label(1 "North East") label(2 "Midwest") ///
	label(3 "South") label(4 "West")) ///
	ylab(, nogrid) 
	
graph export "$figdir/`outname'.pdf", replace

restore

end


/*******************************/
/*** Kernal Density by base **/
/*******************************/

capture program drop kdensity_base

program kdensity_base

set more off
	
syntax, var(varname) xmax(real) [ xti(string) yti(string) outname(string) ]

preserve

tempfile data
save `data', replace

*construct x axis
sort year fips_state /* This is necessary for x axis construction */
gen xaxis_max=`xmax'
tempfile data_temp
save `data_temp', replace
foreach split in broad med narrow {
	use `data_temp', clear
	keep if base_`split' == 1
	sort `var'
	qui count
	gen id_`split'= mod(_n - 1 , `r(N)') + 1 
	replace id_`split'=id_`split'/`r(N)'
	g xaxis_`split'=id_`split'*xaxis_max
	tempfile data_temp`split'
	save `data_temp`split'', replace
}
use `data_temp', clear
merge 1:1 fips_state year using `data_tempbroad', nogen
merge 1:1 fips_state year using `data_tempmed', nogen
merge 1:1 fips_state year using `data_tempnarrow', nogen

tempfile data_temp
save `data_temp', replace

foreach split in broad med narrow {
	use `data_temp', clear
	keep if base_`split' == 1
	kdensity `var', gen(x_`var' kd_`var'_`split') at(xaxis_`split') nograph
	keep fips_state base_`split' year kd*
	tempfile data_`split'
	save `data_`split'', replace
}
use `data_temp', clear
merge 1:1 fips_state year using `data_broad', nogen
merge 1:1 fips_state year using `data_med', nogen
merge 1:1 fips_state year using `data_narrow', nogen

sort `var' xaxis*
	
twoway (scatter kd_`var'_broad xaxis_broad, c(l) msymbol(none) lcolor(black) lpattern(solid) lwidth(thick) ) ///
	(scatter kd_`var'_med xaxis_med, c(l) msymbol(none) lcolor(black) lpattern(dash) lwidth(thick) ) ///
	(scatter kd_`var'_narrow xaxis_narrow, c(l) msymbol(none) lcolor(gs9) lpattern(solid) ), ///
	graphregion(fcolor(white)) graphregion(color(white)) ///
	xtitle(`xti', size(medlarge)) ///
	ytitle(`yti', size(medlarge)) ///
	legend(order(1 2 3) size(medlarge) label(1 "Broad Base") label(2 "Medium Base") ///
	label(3 "Narrow Base")) ///
	ylab(, nogrid) 
	
graph export "$figdir/`outname'.pdf", replace

restore

end
