program define example4
version 2.1
*   Monte Carlo simulation of bootstrap confidence intervals for
*   a misspecified (heteroscedastic, nonnormal errors) regression
*   model.  Generates 100 Monte Carlo samples, and resamples each
*   of them 2,000 times.
*
drop _all
set more 1
set maxobs 2100
set seed 9999
capture erase example4.log
macro define _mcit=1
while %_mcit<101 {
      quietly drop _all
      quietly set obs 80
      quietly generate X=(invnorm(uniform()))^2
      quietly generate Y=3*X+X*((invnorm(uniform())^2)-1)
      *   Previous two lines define the true model.
      quietly regress Y X
      macro define _orb=_b[X]
      macro define _orSE=%_orb/sqrt(_result(6))
      *   Perform the original-sample regression, storing slope
      *   as _orb and standard error _orSE.
      quietly generate XX=.
      quietly generate YY=.
      quietly generate randnum=.
      macro define _bsample=1
      capture erase bstemp.log
      log using bstemp.log
      log off
      while %_bsample<2001 {
            *   Begin bootstrap iterations, indexed by _bsample.
            quietly replace randnum=int(_N*uniform())+1
            quietly replace XX=X[randnum]
            quietly replace YY=Y[randnum]
            quietly regress YY XX
            *   Data resampling, not assuming i.i.d. errors.
            log on
            display %_orb
            display %_orSE
            display _b[XX]
            display (_b[XX]-%_orb)/(_b[XX]/sqrt(_result(6)))
            *   Previous four lines display the original-sample
            *   regression slope (_orb), original-sample standard
            *   error (_orSE), bootstrap regression slope, and
            *   bootstrap studentized slope.
            log off
            macro define _bsample=%_bsample+1
            }
      log close
      drop _all
      infile orb orSE bootb stub using bstemp.log
      save mcit%_mcit, replace
      *   Each Monte Carlo iteration saves a dataset containing
      *   results from 2,000 resamplings.
      macro define _mcit=%_mcit+1
      }
macro define _mcit=1
capture erase mctemp.log
log using mctemp.log
log off
while %_mcit<101 {
      quietly use mcit%_mcit, clear
      log on
      display %_mcit
      quietly summ orb
      display _result(3)
      quietly summ orSE
      display _result(3)
      quietly summ bootb, detail
      display _result(3)
      display sqrt(_result(4))
      display _result(7)
      display _result(13)
      quietly summ stub, detail
      display _result(7)
      display _result(13)
      display
      log off
      *   Preceding lines display the Monte Carlo sample number,
      *   original-sample slope and standard error, bootstrap
      *   mean and standard deviation of regression slope; and
      *   bootstrap percentiles and studentized percentiles.
      macro define _mcit=%_mcit+1
      }
log close
drop _all
#delimit ;
infile iterate orb orSE bootb bootSE p05 p95 t05 t95
       using mctemp.log;
#delimit cr
label variable orb "original sample b coefficient"
label variable orSE "original sample SE of b"
label variable bootb "mean bootstrap b coefficient"
label variable bootSE "SD of bootstrap b coefficient"
label define correct 0 "wrong" 1 "correct"
*
generate stanlo=orb-1.665*orSE
generate stanhi=orb+1.665*orSE
generate stanwide=stanhi-stanlo
label variable stanwide "width standard t interval"
generate stancor=0
replace stancor=1 if stanlo<3 & 3<stanhi
label variable stancor "standard 90% c.i. correct?"
label values stancor correct
*
generate perwide=p95-p05
label variable perwide "width 5th-95th percentile"
generate percor=0
replace percor=1 if p05<3 & 3<p95
label variable percor "percentile 90% c.i. correct?"
label values percor correct
*
generate hyblo=orb-(p95-orb)
generate hybhi=orb-(p05-orb)
generate hybwide=hybhi-hyblo
label variable hybwide "width hybrid percentile interval"
generate hybcor=0
replace hybcor=1 if hyblo<3 & 3<hybhi
label variable hybcor "hybrid 90% c.i. correct?"
label values hybcor correct
*
generate stulo=orb-t95*bootSE
generate stuhi=orb-t05*bootSE
generate stuwide=stuhi-stulo
label variable stuwide "width percentile-t interval"
generate stucor=0
replace stucor=1 if stulo<3 & 3<stuhi
label variable stucor "student 90% c.i. correct?"
label values stucor correct
*
label data "n=80 bootstrap/Monte Carlo"
label variable orb "Y=3X + Xe,e~chi2-1,X~chi2"
save boot4.dta, replace
end