global madedat = "../data"

graph twoway (line Rademacher t_stat if t_stat >= `lower' & t_stat <= `upper') (line Webb t_stat if t_stat >= `lower' & t_stat <= `upper') , xline(`main_t') ti("CDFs of Bootstrapped t-distributions") 
	note("Note: 6 Clusters.  999 bootstrap replications.  Vertical line at main t-statistic.  Zoomed in near main t-staistic.")

/* do nonparametric bootstraps: for standard errors */
use `main_data' , replace

use statefip year `lhs'* `keyrhs'* using `main_data' , replace

merge 1:m statefip using `main_data' , assert(match) keep(match) nogenerate

merge 1:m pseudo_cluster using `main_data' , assert(match) keep(match) nogenerate

/* Rows 4 and 5 */
CR23_IK_CSS , betavec(`betavec') lhs(`reglhs') rhs(`regrhs') key_rhs(`regrhs') 
	noconstant(1) main_data("`main_data'")

qui save `main_data' , replace

syntax  , [lhs(varlist) rhs(varlist) betavec(namelist) key_rhs(varlist) 
	noconstant(integer 0) main_data(string) ]

tempfile main_data

/* get ready to run the monte carlo */
tempfile dummies to_be_sampled tobsample main_data bsout pairs_bs_out

local main_se = _se[policy]

/* Do percentile-T bootstraps */
local main_t = (`b' - `beta_hypothesis') / `se_clu'

/* Do percentile-T bootstraps, for rows 7 and 8  */
local main_t = (`b' - `beta_hypothesis') / `se_clu'

/* next we do a wild bootstrap, Rademacher weights (+-1), and show that the resulting distribution of t-statistics takes only a few values
	(and hence the p-value of the main t-staistic is only interval-identified We will do the same with the Webb 6-point distribution, to see 
	that the t-stats has a much more "continuous" distribution.  */
use `maindata' , clear

qui save `maindata'

cap prog drop make_deviations

do makedata

mata: quickreg_1("`lhs'_deviation","`rhs'_deviation","`clusterid'")

cap mata: mata drop quickreg_1()

clear mata

local matabeta = mybeta

local matase = myse

merge 1:m statefip using "`sourcedata'" , assert(match using) keep(match) nogenerate

merge 1:m pseudo_cluster using `main_data' , assert(match) keep(match) nogenerate

merge 1:m statefip using "`sourcedata'" , assert(match using) keep(match) nogenerate

merge 1:m pseudo_cluster using `main_data' , assert(match) keep(match) nogenerate

			    	di " -> NOTE: Raw estimated variance matrix was non positive semi-definite."

		    di " -> NOTE: Raw estimated variance matrix was non-positive semi-definite."

/* another way to to random effects model 
(compare with xtreg output)
(not in table)							*/
set matsize 6000

set matsize 10000

local toprank = r(max)

local clusize_max = r(max)

local maxrank = r(N)     

di "Number MC replications," _column(25) %6.0f `mc_reps'

di "Number MC replications," _column(35) %6.0f `mc_reps'

/* the main loop */
qui forvalues mc = 1/`mcreps' {

local mc_reps = r(N)

local mc_reps_bigG = 1000

local mc_reps_smallG = 4000

/* prepare the state list, and ensure that all files are sorted on statefip  */
tempfile statelist mcout_all

use `mcout_`ns'' , replace

qui save "`mcout_`ns''" , replace

tempfile mcout_`ns'

qui postfile mcoutput numobs numstates K numbsreps `postlist' using "t4_post_G`numstates'_`mydate'.dta"  , replace

post mcoutput (`numobs') (`numstates') (`K') (`bsreps') (`b') (`se_rob') (`se_clu') (`se_CR2') (`se_CR3') (`dof_CR2_IK') (`dof_CSS') (`se_clu_bs') (`bs_misreps')
	(`p_pairs_bs') (`p_rad_res') (`p_rad_res_high') (`p_rad_res_low') (`p_rad_unres') (`p_webb_res')

cap postclose mcoutput

postclose mcoutput

do_the_monte_carlo , statelist("`statelist'") sourcedata("`mysourcedata'") savedata("`mcout_`ns''")
		mcreps(`num_mc_reps') bsreps(`num_bs_reps') numstates(`ns')

syntax  , [mcreps(integer 3) bsreps(integer 4) numstates(integer 6) statelist(string) sourcedata(string) savedata(string)]

collapse (mean) lnwage lnwage_sy (rawsum) wtsupp [pw=wtsupp] , by(statefip year)

egen `keylhs'_avg_y = mean(`keylhs') , by(year)

egen policy_avg_y = mean(policy) , by(year)

egen `lhs'_avg_y = mean(`lhs') , by(year)

egen `keyrhs'_avg_y = mean(`keyrhs') , by(year)

egen y_wild_rad_restricted_avg_y = mean(y_wild_rad_restricted) , by(year)

egen y_wild_webb_restricted_avg_y = mean(y_wild_webb_restricted) , by(year)

local mean_p_`mthd' = r(mean)

local IK_dof_mean = r(mean)

local nobs_mean = r(mean)

local mean_bsreps_miss = r(mean)

di "CSS dof (mean)," _column(25) %6.1f `CSS_dof_mean'

local yrmean = r(mean)

local rej_rate_`mthd' = r(mean)

local bs_missed_b = r(mean)

local yrmean = r(mean)

local mean_bsreps_miss = r(mean)

summ `touse', meanonly

local meanrank = r(mean)

	di "Note: +/- means the corresponding matrix is added/subtracted"

		    		di "          See Cameron, Gelbach & Miller, "

local varlist "req ex min(1)"

local myp = 2 * min(`pctile' , (1-`pctile'))

local myp2 = 2 * min(`pctile_bot' , (1-`pctile_bot'))

local myp = 2 * min(`pctile' , (1-`pctile'))