*Local assignment
local fignum = 9

clear
set more off

/********************************************************/
/***Analysis 3: Figures of Heterogeneous Elasticities ***/
/********************************************************/


global baserules "rec_val sales_wgt Losscarryback Losscarryforward FranchiseTax FedIncomeTaxDeductible FederalIncomeasStateTaxBase AllowFedAccDep ACRSDepreciation FederalBonusDepreciation throwback combined investment_credit incr_ma incr_fixed"

use "$dtadir/state_taxes_analysis.dta", clear
drop if State == ""
set more off 


bys fips: egen mean_GDP = mean(GDP)

local i = 1
foreach xvar of varlist $baserules  {
	 sum `xvar'  [aw=mean_GDP], d
	 replace `xvar'=(`xvar'-r(mean))/r(sd)
	 
	 gen s`i' = `xvar'
	 local ++i
}

label var rec_val "R\&D Tax Credit"

** Gen main varaibles 
gen log_rev = log( rev_corp)
gen log_corp = log( 1-corpo/100 ) 
gen log_rev_gdp = log(rev_corp/GDP)
*replace r_g = r_g*100000
gen r_gos = rev_corptax/GOS_s
gen log_rev_gos = log(rev_corp/GOS_s)
gen log_gos = log(GOS_s)
gen log_gdp = log(GDP)
gen gos = GOS_s/100000

gen r_g = rev_corptax/GDP*10^4

replace log_corp = corporate_rate
*** Define program for bootstrap 
capture: program drop cent_ME
program define cent_ME, rclass
	args y
	reg `y' log_corp i.year i.fips c.s1 c.s2 c.s3 c.s4 c.s5 c.s6 c.s7 c.s8 c.s9 c.s10 c.s11 c.s12 c.s13 c.s14 c.s15 ///	
	c.log_corp#(c.s1 c.s2 c.s3 c.s4 c.s5 c.s6 c.s7 c.s8 c.s9 c.s10 c.s11 c.s12 c.s13 c.s14 c.s15)[aw=mean_GDP]	
	gen temp_prog = _b[log_corp]+_b[c.log_corp#c.s1]*s1 ///
				 +_b[c.log_corp#c.s2]*s2+_b[c.log_corp#c.s3]*s3 ///   
				 +_b[c.log_corp#c.s4]*s4+_b[c.log_corp#c.s5]*s5 ///
				 +_b[c.log_corp#c.s6]*s6+_b[c.log_corp#c.s7]*s7 ///
				 +_b[c.log_corp#c.s8]*s8+_b[c.log_corp#c.s9]*s9 ///
				 +_b[c.log_corp#c.s10]*s10+_b[c.log_corp#c.s11]*s11 ///
				 +_b[c.log_corp#c.s12]*s12+_b[c.log_corp#c.s13]*s13 ///
				 +_b[c.log_corp#c.s14]*s14+_b[c.log_corp#c.s15]*s15 if e(sample) 
	centile temp_prog if e(sample) , centile(10(10)90)
	
	mat c = (r(c_1),r(c_2),r(c_3),r(c_4),r(c_5),r(c_6),r(c_7),r(c_8),r(c_9))
	return matrix c = c
	return scalar c1 = r(c_1)
	return scalar c2 = r(c_2)
	return scalar c3 = r(c_3)
	return scalar c4 = r(c_4)
	return scalar c5 = r(c_5)
	return scalar c6 = r(c_6)
	return scalar c7 = r(c_7)
	return scalar c8 = r(c_8)
	return scalar c9 = r(c_9)
	
	centile temp_prog if e(sample) & year < 1990, centile(10(10)90)
	mat c = (r(c_1),r(c_2),r(c_3),r(c_4),r(c_5),r(c_6),r(c_7),r(c_8),r(c_9))
	return matrix c = c
	return scalar c1_1980 = r(c_1)
	return scalar c2_1980 = r(c_2)
	return scalar c3_1980 = r(c_3)
	return scalar c4_1980 = r(c_4)
	return scalar c5_1980 = r(c_5)
	return scalar c6_1980 = r(c_6)
	return scalar c7_1980 = r(c_7)
	return scalar c8_1980 = r(c_8)
	return scalar c9_1980 = r(c_9)

	centile temp_prog if e(sample) & inrange(year,1990,1999) , centile(10(10)90)
	mat c = (r(c_1),r(c_2),r(c_3),r(c_4),r(c_5),r(c_6),r(c_7),r(c_8),r(c_9))
	return matrix c = c
	return scalar c1_1990 = r(c_1)
	return scalar c2_1990 = r(c_2)
	return scalar c3_1990 = r(c_3)
	return scalar c4_1990 = r(c_4)
	return scalar c5_1990 = r(c_5)
	return scalar c6_1990 = r(c_6)
	return scalar c7_1990 = r(c_7)
	return scalar c8_1990 = r(c_8)
	return scalar c9_1990 = r(c_9)	

	centile temp_prog if e(sample) & inrange(year,2000,2010) , centile(10(10)90)
	mat c = (r(c_1),r(c_2),r(c_3),r(c_4),r(c_5),r(c_6),r(c_7),r(c_8),r(c_9))
	return matrix c = c
	return scalar c1_2000 = r(c_1)
	return scalar c2_2000 = r(c_2)
	return scalar c3_2000 = r(c_3)
	return scalar c4_2000 = r(c_4)
	return scalar c5_2000 = r(c_5)
	return scalar c6_2000 = r(c_6)
	return scalar c7_2000 = r(c_7)
	return scalar c8_2000 = r(c_8)
	return scalar c9_2000 = r(c_9)		
	
	drop temp_prog
end

bootstrap  c1=r(c1) c2=r(c2) c3=r(c3) c4=r(c4) c5=r(c5) c6=r(c6) c7=r(c7) c8=r(c8) c9=r(c9) ///
		   c1_1980=r(c1_1980) c2_1980=r(c2_1980) c3_1980=r(c3_1980) c4_1980=r(c4_1980) c5_1980=r(c5_1980) c6_1980=r(c6_1980) c7_1980=r(c7_1980) c8_1980=r(c8_1980) c9_1980=r(c9_1980)  ///
		   c1_1990=r(c1_1990) c2_1990=r(c2_1990) c3_1990=r(c3_1990) c4_1990=r(c4_1990) c5_1990=r(c5_1990) c6_1990=r(c6_1990) c7_1990=r(c7_1990) c8_1990=r(c8_1990) c9_1990=r(c9_1990)  ///
		   c1_2000=r(c1_2000) c2_2000=r(c2_2000) c3_2000=r(c3_2000) c4_2000=r(c4_2000) c5_2000=r(c5_2000) c6_2000=r(c6_2000) c7_2000=r(c7_2000) c8_2000=r(c8_2000) c9_2000=r(c9_2000)  ///
		   , reps(50) : cent_ME r_g
		   		   
		   
test (c1_1980 = c1_2000) /// 		   
	(c2_1980 = c2_2000) /// 
	(c8_1980 = c8_2000) /// 
	(c9_1980 = c9_2000) 
scalar pval = r(p) 

*** Ouput for graphs 
mat q = (10,20,30,40,50,60,70,80,90)'
mat ci = e(ci_normal) 	
mat b = e(b) 	

mat c = b[1,1..9]'
mat c_1980 = b[1,10..18]'
mat c_1990 = b[1,19..27]'
mat c_2000 = b[1,28...]'

mat lb = ci[1,1..9]'
mat ub = ci[2,1..9]'
mat lb_1980 = ci[1,10..18]'
mat ub_1980 = ci[2,10..18]'
mat lb_1990 = ci[1,19..27]'
mat ub_1990 = ci[2,19..27]'
mat lb_2000 = ci[1,28...]'
mat ub_2000 = ci[2,28...]'

foreach mat in q c c_1980 c_1990 c_2000 lb ub lb_1980 lb_1990 lb_2000 ub_1980 ub_1990 ub_2000 {
	capture: drop `mat'1
	svmat `mat'
}

reg r_g log_corp i.year i.fips c.s1 c.s2 c.s3 c.s4 c.s5 c.s6 c.s7 c.s8 c.s9 c.s10 c.s11 c.s12 c.s13 c.s14 c.s15 ///	
	c.log_corp#(c.s1 c.s2 c.s3 c.s4 c.s5 c.s6 c.s7 c.s8 c.s9 c.s10 c.s11 c.s12 c.s13 c.s14 c.s15)[aw=mean_GDP]	
	gen temp_prog = _b[log_corp]+_b[c.log_corp#c.s1]*s1 ///
				 +_b[c.log_corp#c.s2]*s2+_b[c.log_corp#c.s3]*s3 ///   
				 +_b[c.log_corp#c.s4]*s4+_b[c.log_corp#c.s5]*s5 ///
				 +_b[c.log_corp#c.s6]*s6+_b[c.log_corp#c.s7]*s7 ///
				 +_b[c.log_corp#c.s8]*s8+_b[c.log_corp#c.s9]*s9 ///
				 +_b[c.log_corp#c.s10]*s10+_b[c.log_corp#c.s11]*s11+_b[c.log_corp#c.s12]*s12+_b[c.log_corp#c.s13]*s13 ///
				 +_b[c.log_corp#c.s14]*s14+_b[c.log_corp#c.s15]*s15 if e(sample) 
			 
* create time series with different percentiles			 
preserve 
collapse (mean) temp_prog (p75) p75_temp=temp_prog (p25) p25_temp=temp_prog ///
(p10) p10_temp=temp_prog (p90) p90_temp=temp_prog , by(year) 
keep if temp_prog !=. 
twoway (scatter temp_prog year , connect(line)) ///
, graphregion(fcolor(white)) graphregion(color(white))  xtitle("Year")  ytitle("Estimated Average Effect by Year") ///
 yline(0, lcolor(black)) 


graph export "$figdir/figure`fignum'.pdf", replace 

restore
