***************************************************************
// Filename: test_sun_shapiro.do
// Author: eddie yu
// Date: 2023-10-16
// Task: Sun and Shapiro Test
//
***************************************************************

clear all
snapshot erase _all
est clear

cd "$output/Work2024"

use "$data/Replication Data/figure6_data_A.dta", clear

summ pr_link_i if base_emp > 0, d 
global iqr = r(p75)-r(p25) // iqr = 1 because it is already standardized in figure6_clean_pr_link_raw.do
di $iqr

* Create Decile bins of PR Exposure
xtile bin_pr_link = pr_link_i, nq(10)

* create bin1, bin2, ..., bin10 where bin1 is totally unaffected group
* Check if bin1's PR exposure = 0
tab bin_pr_link, summ(pr_link_i) // 25% of sample has PR exposure = 0

// * create sample of regression to run DiD w/ TWFE
// forval v = 3/10 {
// 	gen reg_sample`v' = 0
// 	replace reg_sample`v' = 1 if bin_pr_link == `v' | bin_pr_link == 1
// 	tab bin_pr_link if reg_sample`v' == 1 // check if we have totally unaffected & each group
// }



// Regressions
keep if year == 1990 | year == 1995 | year == 2012
xi i.year|pr_link_i , noomit
drop _IyeaXp*1995

xi i.year*i.bin_pr_link, noomit
drop _IyeaXbin_*_1
drop _IyeaXbin_1995_*

local loopvalues = "3/10"

forval i = `loopvalues' {
	label var _IyeaXbin_2012_`i' "Bin `i'"
}


forval i = `loopvalues' {
	eststo bin_`i':  reghdfe emp_growth _IyeaX* [aw=wgt] if bin_pr_link == 1 | bin_pr_link == `i', a(i.pid i.year#i.industr ) cl(fips_state indust ) nocons
	lincom _IyeaXbin_2012_`i'
		estadd local beta_pr = string(r(estimate), "%8.3fc")
		estadd local se_pr   = "(" + string(r(se), "%8.3fc") + ")"
		estadd local p_pr    = "[" + string(r(p), "%8.3fc") + "]"
		estadd local ci_pr   = "[" + string(r(lb), "%8.3fc") + "," + string(r(ub), "%8.3fc") + "]"
		estadd local lb_pr   = string(r(lb), "%8.3fc")
		estadd local ub_pr   = string(r(ub), "%8.3fc")
	summ pr_link_i if bin_pr_link == `i', d
		estadd local avg_pr_link = string(r(mean), "%8.3fc")
		estadd local beta_pr_scaled = string( _b[_IyeaXbin_2012_`i'] / r(mean), "%8.3fc")
		estadd local iqr_adj_pr_link = string(r(mean)/$iqr, "%8.3fc")
		estadd local beta_pr_scaled_iqr = string( _b[_IyeaXbin_2012_`i'] / (r(mean)/$iqr ), "%8.3fc")
}


esttab bin_* using sun_shapiro_pr_bin_longdiff.tex, drop(*) stats(beta_pr se_pr p_pr ci_pr avg_pr_link iqr_adj_pr_link beta_pr_scaled  beta_pr_scaled_iqr N, /// 
	labels("Long Diff Estimate" "s.e." "p-val" "95\% CI" ///
			  "Mean PR Exposure" "Mean PR Exposure (IQR adjusted)" "Estimate / Mean PR Exposure" "Estimate / IQR adjusted PR Exposure")) /// 
	mlabels(, nodepvars) ///
	nostar nonotes nonum nogap replace label nocons collabels(none) varwidth(40) 
	

	
	
	

********************************************************************************
* All-in-one regression [START]
********************************************************************************
eststo emp_1: reghdfe emp_growth _IyeaX* [aw=wgt], a(i.pid i.year#i.industr ) cl(fips_state indust ) nocons
	forval i = `loopvalues' {
		summ pr_link_i if pr_link_i > 0 & _Ibin_pr_li_`i'==1
		estadd local avg_pr_exp_`i' = string(r(mean), "%8.3fc")
		estadd local beta_pr_scaled_`i' = string( _b[_IyeaXbin_2012_`i'] / r(mean), "%8.3fc")
		estadd local beta_pr_scaled_iqr_`i' = string( _b[_IyeaXbin_2012_`i'] / $iqr, "%8.3fc")
	}
	

// 	bys reg_sample*: summ pr_link_i if pr_link_i > 0 & reg_sample`i', d
// 	estadd summ, iqr
//eststo emp_1:  reghdfe emp_growth _IyeaX* reg_sample* [aw=wgt], a(i.pid  ///
//  i.year#i.industr#reg_sample4 i.year#i.industr#reg_sample5 i.year#i.industr#reg_sample6 ///
//  i.year#i.industr#reg_sample7 i.year#i.industr#reg_sample8 i.year#i.industr#reg_sample9 ///
//  i.year#i.industr#reg_sample10) cl(fips_state indust ) nocons

esttab emp_1, drop(*) stats(avg_pr_exp_3 beta_pr_scaled_3 beta_pr_scaled_iqr_3 /// 
							avg_pr_exp_4 beta_pr_scaled_4 beta_pr_scaled_iqr_4 )

esttab emp_1, keep(_IyeaXbin_2012_*) nostar noobs ///
	cells(b(fmt(3)) se(fmt(3) par) p(fmt(3) par([ ]))) ///  // ci(fmt(3) par([ , ])))
	mlabel("2012 Estimate") label replace
	
esttab emp_1, keep(_IyeaXbin_2012_4) stats() nostar noobs ///
	cells(b(fmt(3)) se(fmt(3) par) p(fmt(3) par([ ])) ci(fmt(3) par([ , ]))) ///
	mlabel("2012 Estimate") append
	
********************************************************************************
* All-in-one regression [END]
********************************************************************************	



forval i = 3/10 {
	preserve
// 	local i = 3
	keep if reg_sample`i'==1 
	tab bin_pr_link
	di `avg_pr_link'
	eststo emp_`i':  reghdfe emp_growth _IyeaX* [aw=wgt], a(i.pid i.year#i.industr ) cl(fips_state indust ) nocons
	lincom _IyeaXpr__2012
		estadd local beta_pr = string(r(estimate), "%8.3fc")
		estadd local se_pr   = "(" + string(r(se), "%8.3fc") + ")"
		estadd local p_pr    = "[" + string(r(p), "%8.3fc") + "]"
		estadd local ci_pr   = "[" + string(r(lb), "%8.3fc") + "," + string(r(ub), "%8.3fc") + "]"
		estadd local lb_pr   = string(r(lb), "%8.3fc")
		estadd local ub_pr   = string(r(ub), "%8.3fc")
	summ pr_link_i if pr_link_i > 0, d
		estadd local avg_pr_link = string(r(mean), "%8.3fc")
		estadd local beta_pr_scaled = string( _b[_IyeaXpr__2012] / r(mean), "%8.3fc")
		estadd local iqr_adj_pr_link = string(r(mean)/(r(p75)-r(p25)), "%8.3fc")
		estadd local beta_pr_scaled_iqr = string( _b[_IyeaXpr__2012] / (r(mean)/(r(p75)-r(p25))), "%8.3fc")
// 	est save ster/emp_reg_`i', replace
	restore
}



esttab emp_* using sun_shapiro_pr_bin_longdiff.tex, drop(*) stats(beta_pr se_pr p_pr ci_pr avg_pr_link beta_pr_scaled iqr_adj_pr_link beta_pr_scaled_iqr N, /// 
	labels("Long Diff Estimate" "s.e." "p-val" "95\% CI" ///
			  "Mean PR Exposure" "Estimate / Mean PR Exposure" "Mean PR Exposure (IQR adjusted)" "Estimate / IQR adjusted PR Exposure")) /// 
	mtitles("Bin3" "Bin4" "Bin5" "Bin6" "Bin7" "Bin8" "Bin9" "Bin10") ///
	nostar nonotes nonum nogap replace label nocons collabels(none) varwidth(40)