clear all
set more off
* This program illustrates the use of the Wild Bootstrap
* to improve the size of a heteroskedasticity robust t-test

* 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 double t=0 if _n<=`rep'
qui g double pb=0 if _n<=`rep'
qui g byte tr=0 if _n<=`rep'
qui g byte trb=0 if _n<=`rep'
qui g double tb=0 if _n<=`b'
qui g double x=0 if _n<=`n'
qui g double e=0 if _n<=`n'
qui g double y=0 if _n<=`n'
qui g double u=0 if _n<=`n'
qui g double 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+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 if _n<=`n', vce(hc3) 
qui replace t=_b[x]/_se[x] if _n==`i'
local at=(_b[x]/_se[x]) 

* Start bootstrap loop
      forvalues j = 1/`b' {
	  qui replace u = uniform()                    if _n<=`n'
      qui replace yb=fit+abs(res)*2*(0.5-(u<0.5))  if _n<=`n'
	  qui reg yb x if _n<=`n', vce(hc3)
	  qui replace tb=_b[x]/_se[x]  if _n==`j' 
	  * End bootstrap loop
	  }
* Compute bootstrap results	  
drop res fit	  
qui count if tb>=(`at') in 1/`b'
qui replace pb=2*min(r(N),(`b'-r(N)))/`b' if _n==`i'

* End main loop
}
* Compute and display results
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 
}