## Wednesday, October 10, 2012

### Bootstrapping Percentile Count Vs. Standard Error Approximation

* A typical method for calculating confidence intervals is by using the estimates from some number of redraws and taking the standard error from that.

* Then using the standard error to calculate the confidence interval drawing from the t-distribution.

* An alternative method is to use the draws themselves as the sample and to calculate the confidence interval as the range of values that 1-alpha part of draws fall within.

* This might seem fishy but it is really the theoretically superior method (I believe).

* This simple simulation will see how well both methods do.

* First let's program a data generating command.

cap program drop mysim
program define mysim, rclass

clear
* The first argument that the user puts into the mysim function will be the population size.
set obs 1'

* Now let's generate some simple data.
gen x = rnormal()
gen u = rnormal()

* The coefficient of interest is the coefficient on x.
* Specifically, a good 95% confidence interval will have the simulated value of the coefficient (2) fall within it 95% of the time.
gen y = 5+2*x+u

* Now let's bootstrap the standard errors.
bs, rep(2') saving(bstemp.dta, replace): reg y x
* Simple OLS need not use bootstrapping to derive a useful confidence interval.
* However, for this exercise because OLS is fast and speed is useful in Monte Carlo aproaches.

matrix CIp = e(ci_normal)
scalar CInlower = CIp[1,1]
scalar CInupper = CIp[2,1]

if (CInlower<2 amp="amp" nupper="nupper">2) return scalar CIn = 1
else return scalar CIn = 0

* Now let's see if we can't use the empirical draws to calculate a CI.
* First let's load the BS data data that was generated from the commands.
use bstemp.dta, clear

* Sort the beta coefficients
sort _b_x

* Now let's calculate the values that we will use for our CI.
* We want to start looking at the sorted estimates at the 5% observation or greater to the 95% observation or less and see how frequently the true value (2) falls within that.

* The number of observations in 2'.
scalar lower1 = floor(2'*.025)
* This makes it so that it rounds down at the bottom and up at the top.  This is less conservative.
scalar upper1 = ceil(2'*.975)
di "CIp1: " _b_x[lower1] " & "_b_x[upper1]
if (_b_x[lower1]<2 amp="amp" b_x="b_x" upper1="upper1">2) return scalar CIp1 = 1
else return scalar CIp1 = 0

scalar lower2 = ceil(2'*.025)
* This makes it so that it rounds up at the bottom and down at the top.  This is more conservative.
scalar upper2 = floor(`2'*.975)
di "CIp2: " _b_x[lower2] " & "_b_x[upper2]
if (_b_x[lower2]<2 amp="amp" b_x="b_x" upper2="upper2">2) return scalar CIp2 = 1
else return scalar CIp2 = 0

end

mysim 100 200
return list
* One run produces results are effectively identical because the difference in methods is small.

* Let's instead try 1000 simulations.
simulate CIn=r(CIn) CIp1=r(CIp1) CIp2=r(CIp2), reps(10000):  mysim 100 200

sum

* Normally I leave results to be run by blog readers.  However, this took a long time and many repetitions is generally important.
/*   Variable |       Obs        Mean    Std. Dev.       Min        Max
-------------+--------------------------------------------------------
CIn |      9992    .9377502    .2416208          0          1
CIp1 |      9992    .9334468    .2492592          0          1
CIp2 |      9992    .9334468    .2492592          0          1
*/

* It seems that the percentile method is not working as well as the standard error approximation method.

* However, the difference is slight though surprising.  I suspect that I might not be choosing the correct lower and upper range combination.  Of course the difference is slight.  Only 4 draws out of 10,000.

* The percentile method has other values as well.  It can also be used to generate a joint rejection rate. This can be very useful in cracking difficult problems.