clear all
set more off
* This program illustrates the use of the Residual Bootstrap
* to improve the size of a heteroskedasticity robust t-test
* under homoskedasticity; bootstrap errors are generated with
* the bsample command

* Define main parameters 
local n = 50
local rep = 1000
local b = 99

* Initialize series
local smpl=max(`n',`rep',`b')
set obs `smpl'
qui g t=0 if _n<=`rep'
qui g pb=0 if _n<=`rep'
qui g byte tr=0 if _n<=`rep'
qui g byte trb=0 if _n<=`rep'
qui g x=0 if _n<=`n'
qui g e=0 if _n<=`n'
qui g y=0 if _n<=`n'
qui g yb=0 if _n<=`n'

* Start outer loop over different values of the parameters
foreach m of numlist 0 0.25 {
set seed 31052011

* Start main loop
forvalues i = 1/`rep' {
* Generate data
qui replace x = rchi2(2) if _n<=`n'
qui replace e = rnormal(0,1) if _n<=`n'
qui replace y=1+`m'*x+e   if _n<=`n'
* Run regression under the null and save results
qui reg y if _n<=`n'
qui predict fit if _n<=`n', xb
qui predict res if _n<=`n', residuals
* Run regression and compute t statistic
qui reg y x, robust  
qui replace t=_b[x]/_se[x] if _n==`i'
local at=(_b[x]/_se[x]) 
local count=0
qui save z.dta, replace

* Start bootstrap loop
      forvalues j = 1/`b' {
	    qui keep res
            qui bsample if _n<=`n'
            qui rename res eb
            qui merge 1:1 _n using z.dta, nogenerate noreport
	qui replace yb=fit+0*x+eb  if _n<=`n'
	qui reg yb x, robust 
	qui local count=`count'+ ((_b[x]/_se[x])>=(`at')) 
	drop eb 
        * End bootstrap loop
        }
* Compute bootstrap results
drop res fit	  
qui replace pb=2*min(`count',(`b'-`count'))/`b' if _n==`i'

* End main loop
}
qui replace tr=abs(t)>1.96 if _n<=`rep'
qui replace trb=pb<0.05 if _n<=`rep'
di
di as result "Results for m = " `m'
su tr trb if _n<=`rep'

* End outer loop 
}