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
* 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 = 12>
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
* 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 = 12>
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 = 12>
else return scalar CIp2 = 0
mysim 100 200
* 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
* 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.