mat `xxinv' = e(V)

tempname betavec v_naive

di in gr "Number of unique values of `varlist' is  " in ye _result(18)

foreach var in pairs_bs rad_res rad_unres webb_res {

foreach var in pairs_bs rad_res webb_res {

foreach var in `varlist' {

		confirm new var `generate'

gen `var'_deviation = `var' - `var'_avg_y

gen `var'_deviation = ds_`var' - yavg_ds_`var'

		di "    If some variable in your options list contains the string 'robust', you will"

		di "Calculating cov part for variables: `grouplist' (`plusminus')"

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

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

reg lnwage policy age age2 yrseduc , vce(bootstrap , cluster(statefip) reps(`bsreps') saving(`bootout') )

/* note that vce(robust) will do cluster robust on the i-level (statefip) cluster */
xtreg `lhs' `rhs'   , pa corr(ar1) vce(robust) noconstant

/* Row 1, default standard erorrs */
reg `reglhs' `regrhs' , vce(ols) noconstant

/* Row 3, Cluster on state */
reg `reglhs' `regrhs' , vce(cluster statefip) noconstant

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

	local vcetype =         e(vcetype)

end commented for version 2.3.1 fix of missing cluster values issue

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