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.