The R package boot implements a variety of bootstrappingtechniques including the basic non-parametric bootstrap describedabove. The boot package was written to accompany the textbook Bootstrap Methods and Their Application by (Davison and Hinkley 1997).The two main functions in boot are boot() andboot.ci(), respectively. The boot() function implementsthe bootstrap for a statistic computed from a user-supplied function.The boot.ci() function computes bootstrap confidence intervalsgiven the output from the boot() function.
The arguments to boot() are:
library(boot)args(boot)## function (data, statistic, R, sim = "ordinary", stype = c("i", ## "f", "w"), strata = rep(1, n), L = NULL, m = 0, weights = NULL, ## ran.gen = function(d, p) d, mle = NULL, simple = FALSE, ..., ## parallel = c("no", "multicore", "snow"), ncpus = getOption("boot.ncpus", ## 1L), cl = NULL) ## NULLwhere data is a data object (typically a vector or matrixbut can be a time series object too), statistic is a use-specifiedfunction to compute the statistic of interest, and R is thenumber of bootstrap replications. The remaining arguments are notimportant for the basic non-parametric bootstrap. The function assignedto the argument statistic has to be written in a specificform, which is illustrated in the next example. The boot()function returns an object of class boot for whichthere are print and plot methods.
The arguments to boot.ci() are
args(boot.ci)## function (boot.out, conf = 0.95, type = "all", index = 1L:min(2L, ## length(boot.out$t0)), var.t0 = NULL, var.t = NULL, t0 = NULL, ## t = NULL, L = NULL, h = function(t) t, hdot = function(t) rep(1, ## length(t)), hinv = function(t) t, ...) ## NULLwhere boot.out is an object of class boot,conf specifies the confidence level, and type isa subset from c("norm", "basic", "stud", "perc", "bca") indicating the type of confidence interval to compute.The choices norm, perc, and “bca” compute thenormal confidence interval (8.16), the percentileconfidence interval (8.17), and the BCA confidence interval, respectively. The remaining arguments are not important for the computation of thesebootstrap confidence intervals.
Example 2.39 (Nonparametric bootstrap of GWN model using the R package boot)To use the boot() function to implement the bootstrap for\(\hat{\mu},\) a function must be specified to compute \(\hat{\mu}\)for each bootstrap sample. The function must have two arguments: xand idx. Here, x represents the original data andidx represents the random integer index (created internallyby boot()) to subset x for each bootstrap sample.For example, a function to be passed to boot() for \(\hat{\mu}\) is
mean.boot = function(x, idx) { # arguments: # x data to be resampled # idxvector of scrambled indices created by boot() function # value: # ansmean value computed using resampled data ans = mean(x[idx]) ans }To implement the nonparametric bootstrap for \(\hat{\mu}\) with 999samples use
set.seed(123)muhat.boot = boot(msftRetS, statistic = mean.boot, R=999) class(muhat.boot) ## [1] "boot"names(muhat.boot)## [1] "t0""t" "R" "data" "seed" "statistic"## [7] "sim""call" "stype" "strata""weights"The returned object muhat.boot is of class boot.The component t0 is the sample estimate \(\hat{\mu}\), andthe component t is a \(999\times1\) matrix containing the bootstrapestimates \(\{\hat{\mu}_{1}^{*},\ldots,\hat{\mu}_{999}^{*}\}\). Theprint method shows the sample estimate, the bootstrap bias and thebootstrap standard error:
muhat.boot## ## ORDINARY NONPARAMETRIC BOOTSTRAP## ## ## Call:## boot(data = msftRetS, statistic = mean.boot, R = 999)## ## ## Bootstrap Statistics :## originalbiasstd. error## t1* 0.00915 0.000482 0.00758These statistics can be computed directly from the components of muhat.boot
c(muhat.boot$t0,mean(muhat.boot$t) - muhat.boot$t0, sd(muhat.boot$t))## [1] 0.009153 0.000482 0.007584A histogram and nornal QQ-plot of the bootstrap values, shown in Figure8.3, can be created using the plot method
plot(muhat.boot)Figure 8.3: plot method for objects of class boot
Normal, percentile, and bias-adjusted (bca) percentile 95% confidence intervals can be computed using
boot.ci(muhat.boot, conf=0.95, type=c("norm", "perc", "bca"))## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS## Based on 999 bootstrap replicates## ## CALL : ## boot.ci(boot.out = muhat.boot, conf = 0.95, type = c("norm", ## "perc", "bca"))## ## Intervals : ## Level Normal PercentileBCa ## 95%(-0.0062, 0.0235 )(-0.0057, 0.0245 )(-0.0062, 0.0239 ) ## Calculations and Intervals on Original ScaleBecause the bootstrap distribution looks normal, the normal, percentile, and bca confidence intervals are very similar.
The GWN model estimate \(\hat{\sigma}\) can be “bootstrapped” ina similar fashion. First, we write a function to compute \(\hat{\sigma}\)for each bootstrap sample
sd.boot = function(x, idx) {ans = sd(x[idx]) ans }Then we call boot() with statistic=sd.boot
set.seed(123)sigmahat.boot = boot(msftRetS, statistic = sd.boot, R=999) sigmahat.boot## ## ORDINARY NONPARAMETRIC BOOTSTRAP## ## ## Call:## boot(data = msftRetS, statistic = sd.boot, R = 999)## ## ## Bootstrap Statistics :## originalbiasstd. error## t1*0.101 -0.000783 0.00764Figure 8.4: Bootstrap distribution for \(\hat{\sigma}\)
The bootstrap distribution, shown in Figure 8.4,looks a bit non-normal. The 95% confidences intervals are computed using
boot.ci(sigmahat.boot, conf=0.95, type=c("norm", "perc", "bca"))## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS## Based on 999 bootstrap replicates## ## CALL : ## boot.ci(boot.out = sigmahat.boot, conf = 0.95, type = c("norm", ## "perc", "bca"))## ## Intervals : ## Level Normal PercentileBCa ## 95%( 0.0873, 0.1173 )( 0.0864, 0.1162 )( 0.0892, 0.1200 ) ## Calculations and Intervals on Original Scale## Some BCa intervals may be unstableWe see that all of these intervals are quite similar.
\(\blacksquare\)
Example 2.44 (Nonparametric bootstrap for example functions)The real power of the bootstrap procedure comes when we apply it toplug-in estimates of functions like (8.1)-(8.4) for which thereare no easy analytic formulas for estimating bias and standard errors. In this situation, the bootstrap procedure easily computes numerical estimates of the bias,standard error, and 95% confidence interval. All we need is an R function for computingthe plug-in function estimate in a form suitable for boot(). The R functions for the example functions are:
f1.boot = function(x, idx, alpha=0.05) { q = mean(x[idx]) + sd(x[idx])*qnorm(alpha) q}f2.boot = function(x, idx, alpha=0.05, w0=100000) { q = mean(x[idx]) + sd(x[idx])*qnorm(alpha) VaR = -w0*q VaR}f3.boot = function(x, idx, alpha=0.05, w0=100000) { q = mean(x[idx]) + sd(x[idx])*qnorm(alpha) VaR = -(exp(q) - 1)*w0 VaR }f4.boot = function(x, idx, r.f=0.0025) { SR = (mean(x[idx]) - r.f)/sd(x[idx]) SR}The functions \(f_1\), \(f_2\), and \(f_3\) have the additional argument alphawhich specifies the tail probability for the quantile, and the functions \(f_2\), and \(f_3\) have the additional argument w0 for the initial wealth invested. The function \(f_4\) has the additional argument r.f for the risk-free rate.
To compute the nonparametric bootstrap for \(f_1\) with \(\alpha=0.05\) use:
set.seed(123)f1.bs = boot(msftRetS, statistic = f1.boot, R=999) f1.bs## ## ORDINARY NONPARAMETRIC BOOTSTRAP## ## ## Call:## boot(data = msftRetS, statistic = f1.boot, R = 999)## ## ## Bootstrap Statistics :## original biasstd. error## t1*-0.158 0.00177 0.0124Here, the bootstrap estimate of bias is small and the bootstrap estimated standard error is similar to the delta method and jackknife standard errors.
plot(f1.bs)The bootstrap distribution looks normal, and the three methods for computing 95% confidence intervals are essentially the same:
boot.ci(f1.bs, type=c("norm", "perc", "bca"))## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS## Based on 999 bootstrap replicates## ## CALL : ## boot.ci(boot.out = f1.bs, type = c("norm", "perc", "bca"))## ## Intervals : ## Level Normal PercentileBCa ## 95%(-0.184, -0.135 )(-0.181, -0.132 )(-0.189, -0.138 ) ## Calculations and Intervals on Original Scale## Some BCa intervals may be unstableTo compute the nonparametric bootstrap for \(f_2\) with \(\alpha=0.05\) and \(W_0=\$100,000\) use:
set.seed(123)f2.bs = boot(msftRetS, statistic = f2.boot, R=999) f2.bs## ## ORDINARY NONPARAMETRIC BOOTSTRAP## ## ## Call:## boot(data = msftRetS, statistic = f2.boot, R = 999)## ## ## Bootstrap Statistics :## original biasstd. error## t1*15780-1771244The bootstrap estimate of bias is small and the bootstrap estimated standard error is similar to the delta method and jackknife standard errors.
plot(f2.bs)The bootstrap distribution looks normal, and the three bootstrap 95% confidence intervals are all very similar:
boot.ci(f2.bs, type=c("norm", "perc", "bca"))## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS## Based on 999 bootstrap replicates## ## CALL : ## boot.ci(boot.out = f2.bs, type = c("norm", "perc", "bca"))## ## Intervals : ## Level Normal PercentileBCa ## 95%(13518, 18395 )(13161, 18089 )(13801, 18940 ) ## Calculations and Intervals on Original Scale## Some BCa intervals may be unstableTo compute the nonparametric bootstrap for \(f_3\) with \(\alpha=0.05\) and \(W_0=\$100,000\) use:
set.seed(123)f3.bs = boot(msftRetC, statistic = f3.boot, R=999) f3.bs## ## ORDINARY NONPARAMETRIC BOOTSTRAP## ## ## Call:## boot(data = msftRetC, statistic = f3.boot, R = 999)## ## ## Bootstrap Statistics :## original biasstd. error## t1*14846-1821223The bias is small and the bootstrap estimated standard error is larger than the delta method standard error and is closer to the jackknife standard error.
plot(f3.bs)The bootstrap distribution looks a little asymmetric and the BCA confidence interval is a bit different from the normal and percentile intervals:
boot.ci(f3.bs, type=c("norm", "perc", "bca"))## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS## Based on 999 bootstrap replicates## ## CALL : ## boot.ci(boot.out = f3.bs, type = c("norm", "perc", "bca"))## ## Intervals : ## Level Normal PercentileBCa ## 95%(12631, 17426 )(12295, 17155 )(12994, 18061 ) ## Calculations and Intervals on Original Scale## Some BCa intervals may be unstableFinally, to compute the bootstrap for \(f_4\) with \(r_f = 0.0025\) use:
set.seed(123)f4.bs = boot(msftRetS, statistic = f4.boot, R=999) f4.bs## ## ORDINARY NONPARAMETRIC BOOTSTRAP## ## ## Call:## boot(data = msftRetS, statistic = f4.boot, R = 999)## ## ## Bootstrap Statistics :## original biasstd. error## t1*0.0655 0.00391 0.0739The bootstrap estimated bias is small and the standard error is similar to the delta method and the jackknife.
plot(f4.bs)Here, the bootstrap distribution is slightly asymmetric.
boot.ci(f4.bs, conf=0.95, type=c("norm", "perc", "bca"))## BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS## Based on 999 bootstrap replicates## ## CALL : ## boot.ci(boot.out = f4.bs, conf = 0.95, type = c("norm", "perc", ## "bca"))## ## Intervals : ## Level Normal PercentileBCa ## 95%(-0.0832, 0.2065 )(-0.0821, 0.2132 )(-0.0897, 0.2030 ) ## Calculations and Intervals on Original ScaleNotice that the Percentile and BCA confidence intervals are a bit different from the normal approximation interval.
\(\blacksquare\)