/* calls program "table4_np_boot", which is defined in t34_programs.do.  That program re-does the "deviations from means" and then 
re-runs the regression, to get the betas and (clustered) standard errors. */
qui bootstrap b=r(mybeta) se=r(myse) , reps(`bsreps') nodots cluster(pseudo_cluster) idcluster(pseudo_cluster_2) 
	seed(10101) saving(`pairs_bs_out' , double)
	: table4_np_boot , lhs(`lhs') rhs(`rhs') keyrhs(`keyrhs') clusterid(pseudo_cluster_2)

local rej_rate_`mthd' = r(mean)

local IK_dof_p5 = r(p5)

local CSS_dof_mean = r(mean)

local mc_reps = r(N)

local mean_bsreps_miss = r(mean)

local maxrank = r(N)     

local toprank = r(max)

local G_ar1 = r(r)

local yrmean = r(mean)

local se_fe_clu_bs = r(sd)

local meanrank = r(mean)

local mc_reps = r(N)

local bs_missed_b = r(mean)

local yrmean = r(mean)

local bs_misreps_t = r(sum)

local meanrank = r(mean)

local clusize_max = r(max)

	tempname N_clust N df_m F r2 mss rss r2_a ll ll_0 

*	ereturn scalar df_r =              `df_r'
	ereturn scalar r2 =                `r2'

	scalar `r2_a' =         e(r2_a)

di "Avg Rad Res p-avg, " _column(25) %6.4f `mean_p_rad_res'

di "Avg Rad Res p-low, " _column(25) %6.4f `mean_p_rad_res_low'

di "percentile T Rad2, ,"											_column(35) %6.4f `p_rad_res'

merge t_stat using `rad'

tempfile rad

foreach mthd in rob_Tnk clu_Tnk clu_Tg1 clu_Tg2 CR2_Tg1 CR3_Tg1 CR2_Tik clu_Tcss bsse_Tg1 pairs_bs rad_res_low rad_res rad_res_high webb_res rad_unres {

foreach var in pairs_bs rad_res webb_res {

foreach mthd in def_Tnk clu_Tnk clu_Tg1 clu_Tg2 bsse_Tg1 pairs_bs rad_res webb_res {

foreach mthd in rob_Tnk clu_Tnk clu_Tg1 clu_Tg2 CR2_Tg1 CR3_Tg1 CR2_Tik clu_Tcss bsse_Tg1 pairs_bs rad_res rad_res_high rad_res_low webb_res rad_unres {

foreach mthd in rob_Tnk clu_Tnk clu_Tg1 clu_Tg2 CR2_Tg1 CR3_Tg1 CR2_Tik clu_Tcss bsse_Tg1 pairs_bs rad_res rad_res_high rad_res_low webb_res rad_unres {

foreach mthd in rob_Tnk clu_Tnk clu_Tg1 clu_Tg2 CR2_Tg1 CR3_Tg1 CR2_Tik clu_Tcss bsse_Tg1 pairs_bs rad_res rad_res_high rad_res_low webb_res rad_unres {

foreach var in pairs_bs rad_res rad_unres webb_res {

keep Rademacher t_rad_res

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.")

di "Random effects"

gen rank = .

gen rank = .

gen rank = .

summ rank if abs(t_`var' - `main_t') < 0.0001  

summ rank if abs(t_`var' - `main_t') < 0.0001  

local K = e(rank)

local K = e(rank)

di "Rej. rate `mthd',"	_column(35) %6.3f `rej_rate_`mthd''

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

	*now reconstruct variance matrix using spectral decomposition formula (e.g., Def A.16 in Greene, 6th)
	tempname raw_running_sum

		    di "            in e(rawcovmat)."

global rawdat = "../data"

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

prog def table3_np_boot , rclass

xtreg lnwage policy age age2 yrseduc , re

/* Random Effects model with cluster bootstrap for standard errors */
xtreg lnwage policy age age2 yrseduc , re vce(bootstrap , reps(399) seed(10101) )

	label var `count' "Records per `varlist'"

local command = "reg `lhs' `rhs'  , cluster(pseudo_cluster) "

local command = "reg `reglhs' `regrhs' , noconstant cluster(pseudo_cluster) "

		local options = regexr("`options'", "robust", "")

local reglhs = "`lhs'_deviation"

	ereturn local estat_cmd =         "regress_estat"

local regrhs = "`rhs'_deviation"

di "Rej. rate `mthd',"	_column(25) %6.4f `rej_rate_`mthd''

qui gen rej = .

qui gen rej = .

qui summ rej

local rej_rate_`mthd' = r(mean)

		di " -> Removing string 'robust' from your options line: it's unnecessary as an option,"

use "${rawdat}\cps_00003" , replace

/* generate state policy variable */
use statefip using "${madedat}/CPS_2012_micro" , replace

use statefip using `to_be_sampled' , replace

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

use pseudo_cluster using `to_be_sampled' , replace

use "`sourcedata'" , replace

use `bsout' , replace

use "`sourcedata'" , replace