clear all
set mem 32m
set more on

cap log close
log using bsw-example, replace

* eliminating simulation bias
sysuse auto, clear
bstrap , reps(150) : mean mpg
estat bootstrap, all
more

bsweights bsw, reps(148) n(50) balanced nosvy
bs4rw , rw(bsw*) : mean mpg
estat bootstrap, all
more

cap noi use http://www.stata-press.com/data/r10/nhanes2.dta, clear
* use locally if no Internet connection
if _rc use nhanes2, clear

* is bs4rw available?
cap noi which bs4rw
if _rc==111 {
   net from http://www.stata.com/users/jpitblado/
   net install bs4rw
}
help bs4rw
more

cap noi which bsweights
if _rc==111 {
   net from http://web.missouri.edu/~kolenikovs/stata/
   net install bsweights
}

help bsweights
more

* define some analysis variables
gen bmi = 10000*weight/(height*height)

* the analysis I want to run
* local myestim reg bmi sex age
local myestim logit highbp female age bmi

* basic correct estimation
svy : `myestim'
estimates note : "Original design"
est store correct
more

svy , subpop(black) : `myestim'
est store correct_subpop
more

* compute the totals to use in calibration example
gen byte ones = 1
svy : total ones, over(sex)
mat Totals = e(b)

* create a quasi-design
egen fakepsu = group(strata psu)

recode strata (1/2 = 1) (3/5 = 2) (6/8 = 3) (9/12=4) (13/15=5) ///
  (16/19=6) (20/21=7) (22/25=8) (26/29=9) (30/32=10), gen(fakestrata)

svyset fakepsu [pw=finalwgt], strata(fakestrata)
svydes
more

* basic linearized fake estimation
svy : `myestim'
estimates note : "Faked design, linearization"
est store fake_lin

* simple bootstrap
bsweights bsw , n(3) reps(100) replace dots
bs4rw , rw(bsw*) : `myestim' [pw=final]
estimates note : "Faked design, simple bootstrap, 100 replications"
est store fake_bs1

est tab correct fake*, b(%6.4f) se(%8.4f)
more

* balanced bootstrap: note the particular number of replications!
bsweights bsw , n(-3) reps(100) replace dots balanced
drop bsw*
set seed 56098
bsweights bsw , n(-3) reps(96) replace dots balanced
bs4rw , rw(bsw*) : `myestim' [pw=final]
estimates note : "Faked design, balanced bootstrap, 96 replications"
est store fake_bs2

est tab correct fake*, b(%6.4f) se(%8.4f)
more

* QMC bootstrap: note the particular number of replications!
drop bsw*
bsweights bsw , n(3) reps(96) qmcstrat shuffle replace dots float
bs4rw , rw(bsw*) : `myestim' [pw=final]
estimates note : "Faked design, QMC stratified bootstrap, 320 replications"
est store fake_bs3

est tab correct fake*, b(%6.4f) se(%8.4f)
more

* another flavor of QMC
drop bsw*
bsweights bsw , n(3) reps(120) qmcmat shuffle balance replace dots float
bs4rw , rw(bsw*) : `myestim' [pw=final]
estimates note : "Faked design, QMC 2D, 120 replications"
est store fake_bs4

est tab correct fake*, b(%6.4f) se(%8.4f)
more

* calibrated bootstrap
cap program drop CalGender
program define CalGender
   args wvar

   total ones [pw=`wvar'], over(sex)
   replace `wvar' = `wvar'*Totals[1,1]/_b[ones:Male]
   replace `wvar' = `wvar'*Totals[1,2]/_b[ones:Female]

end

set seed 56098
drop bsw*
bsweights bsw , n(-1) reps(96) replace dots balanced calibrate(CalGender @)
bs4rw , rw(bsw*) : `myestim' [pw=final]
estimates note : "Faked design, balanced bootstrap, 96 replications, calibrated for gender"
est store fake_bs5

est tab correct fake*, b(%6.4f) se(%8.4f)
more

* subpopulation
cap program drop CalSubpop
program define CalSubpop
   args wvar
   replace `wvar' = `wvar' * black
end

drop bsw*
bsweights bsw , n(-1) reps(96) replace dots balanced calibrate(CalSubpop @)
bs4rw , rw(bsw*) : `myestim' [pw=final*black ]
estimates note : "Faked design, balanced bootstrap, 96 replications, black subpopulation"
est store fake_subpop1

svy , subpop(black) : `myestim'
est store fake_subpop2

est tab *subpop*, se
more

log close

exit