---
title: "The Bootstrap"
author: "Eric Zivot"
date: "Monday, May 11, 2015"
output: slidy_presentation
---

## Set options and load packages

```{r}
options(digits=3, width=70)
library(boot)  
library(IntroCompFinR)
library(PerformanceAnalytics)

```

## R function <code>sample()</code>

Consider random sampling from the vector <code>r=c(0.1,0.05,-0.02,0.03,-0.04)</code>. 

* Method 1: Use <code>sample()</code> to create a random sample from the integers 1 to 5 and use it to subset the vector <code>r</code>

```{r}
r = c(0.1,0.05,-0.02,0.03,-0.04)
set.seed(123)
idx = sample(5, replace=TRUE)
idx
r[idx]
```

## R function <code>sample()</code>

* Method 2: Use <code>sample()</code> directly to automate the previous two step process:

```{r}
set.seed(123)
sample(r, replace=TRUE)

```


```{r, echo=FALSE}
# get data silently
data(msftDailyPrices, sp500DailyPrices, sbuxDailyPrices)
msftPrices = to.monthly(msftDailyPrices, OHLC=FALSE)
sp500Prices = to.monthly(sp500DailyPrices, OHLC=FALSE)
sbuxPrices = to.monthly(sbuxDailyPrices, OHLC=FALSE)
# set sample to match book chapter
smpl = "1998-01::2012-05"
msftPrices = msftPrices[smpl]
sp500Prices = sp500Prices[smpl]
sbuxPrices = sbuxPrices[smpl]
# calculate returns
msftRetS = Return.calculate(msftPrices, method="simple")
sp500RetS = Return.calculate(sp500Prices, method="simple")
sbuxRetS = Return.calculate(sbuxPrices, method="simple")
msftRetS = msftRetS[-1]
sp500RetS = sp500RetS[-1]
sbuxRetS = sbuxRetS[-1]
msftRetC = log(1 + msftRetS)
sp500RetC = log(1 + sp500RetS)
sbuxRetC = log(1 + sbuxRetS)
# merged data set
cerRetS = merge(msftRetS, sbuxRetS, sp500RetS)
cerRetC = merge(msftRetC, sbuxRetC, sp500RetC)
colnames(cerRetS) = colnames(cerRetC) = c("MSFT", "SBUX", "SP500")
```

## Example Data: monthly cc returns on MSFT

```{r}
my.panel <- function(...) {
  lines(...)
  abline(h=0)
}
plot.zoo(msftRetC, col="blue", lwd=2, main="", panel=my.panel, type="l")

```

## Put <code>xts</code> data in matrix for use with R package boot

```{r}
returns.mat = as.matrix(cerRetC)
MSFT = returns.mat[,"MSFT", drop=FALSE]
```

## Estimate CER model parameters for MSFT

Compute CER model estimates for MSFT

```{r}
n.obs = nrow(MSFT)
muhat.MSFT = mean(MSFT)
sigmahat.MSFT = sd(MSFT)
estimate = c(muhat.MSFT, sigmahat.MSFT)
```

Calculate analytic standard errors

```{r}
se.muhat.MSFT = sigmahat.MSFT/sqrt(n.obs)
se.sigmahat.MSFT = sigmahat.MSFT/sqrt(2*n.obs)
se = c(se.muhat.MSFT, se.sigmahat.MSFT)
```

Show estimates with analytic standard errors

```{r}
ans = rbind(estimate, se)
colnames(ans) = c("Mu", "Sigma")
ans
```

## Brute force bootstrap

Same idea as Monte Carlo Simulation but instead of
generating random data from an assumed distribution or model,
you generate random data by sampling with replacement
from the observed data using <code>sample()</code>.

Example: Bootstrap $\hat{\mu}_{MSFT}$ 

```{r}
# Create bootstrap simulation for muhat
B = 999
muhat.boot = rep(0, B)
n.obs = nrow(MSFT)
set.seed(123)
for (i in 1:B) {
  boot.data = sample(MSFT, n.obs, replace=TRUE)
  muhat.boot[i] = mean(boot.data)
}


```

## Bootstrap bias estimate

Calculate bootstrap bias estimate

```{r}
mean(muhat.boot) - muhat.MSFT
```

Bootstrap bias is close to zero 

## Bootstrap SE

Calculate bootstrap SE and compare with analytic SE

```{r}
sd(muhat.boot)
# analytic SE
se.muhat.MSFT
```

bootstrap SE is very close to analytic SE.

## Bootstrap distribution

```{r}
# plot bootstrap distribution and show normal qq-plot
par(mfrow=c(1,2))
  hist(muhat.boot, col="cornflowerblue")
  abline(v=muhat.MSFT, col="white", lwd=2)
  qqnorm(muhat.boot, col="cornflowerblue")
  qqline(muhat.boot)
par(mfrow=c(1,1))
```

Bootstrap distribution looks like normal distribution.

## Bootstrap 95% confidence interval

Calculate bootstrap 95% confidence interval based on normal distribution

```{r}
se.boot = sd(muhat.boot)
lower = muhat.MSFT - 2*se.boot
upper = muhat.MSFT + 2*se.boot
width = upper - lower
c(lower, upper, width)
```


## R package boot

* Implements a variety of bootstrapping functions
* Background material is book by Davidson and Hinkley, *Bootstrap Methods and Their
Application*, Cambridge University Press, 1997.
* Main functions are:
    + <code>boot()</code> bootstrap user supplied function
    + <code>boot.ci()</code> compute bootstrap confidence interval
    
## Using the <code>boot()</code> function

The function <code>boot()</code> requires user-supplied functions that take
two arguments: data and an index. The index is created
by the <code>boot()</code> function and represents random resampling with
replacement

```{r}
# function for bootstrapping sample mean
mean.boot = function(x, idx) {
# arguments:
# x      data to be resampled
# idx    vector of scrambled indices created by boot() function
# value:
# ans		 mean value computed using resampled data
     ans = mean(x[idx])
     ans
}
```

## Using the <code>boot()</code> function

```{r}
# bootstrapping muhat
MSFT.muhat.boot = boot(MSFT, statistic = mean.boot, R=999)
class(MSFT.muhat.boot)
names(MSFT.muhat.boot)
MSFT.muhat.boot
```

Here, "original" is $\hat{\mu}$; "bias" is bootstrap bias; "std. error" is the bootstrap SE.
Recall, the analytic SE is `r se.muhat.MSFT`. 

## Plot method for "<code>boot</code>" objects

Plot method shows histogram of bootstrap values and normal qq-plot.

```{r}
plot(MSFT.muhat.boot)

```

## Bootstrap confidence intervals

Use the function <code>boot.ci()</code> to compute bootstrap confidence intervals

```{r}
boot.ci(MSFT.muhat.boot, conf = 0.95, type = c("norm","perc"))

```

## Bootstrapping $\hat{\sigma}$ for MSFT

```{r}
# function for bootstrapping sample standard deviation
sd.boot = function(x, idx) {
# arguments:
# x   	data to be resampled
# idx		vector of scrambled indices created by boot() function
# value:
# ans		sd value computed using resampled data
     ans = sd(x[idx])
     ans
}
```

```{r}
# bootstrap sigmahat
MSFT.sigmahat.boot = boot(MSFT, statistic = sd.boot, R=999)
MSFT.sigmahat.boot
```

Recall, the analytic SE for sigmahat is `r se.sigmahat.MSFT`. Here the bootstrap SE should be more accurate because the analytic SE is only an approximation based on the central limit theorem.

## Bootstrap distribution

```{r}
# plot bootstrap distribution
plot(MSFT.sigmahat.boot)
```

Bootstrap distribution looks pretty normal

## Bootstrap 95% confidence intervals

```{r}
# compute 95% confidence intervals
boot.ci(MSFT.sigmahat.boot, conf=0.95, type=c("norm", "perc"))

```

## Bootstrapping normal VaR for MSFT

Create function to compute normal Value-at-Risk to be passed to <code>boot()</code>

```{r}
ValueAtRisk.boot = function(x, idx, p=0.05, w=100000) {
# x     data to be resampled
# idx		vector of scrambled indices created by boot() function
# p		  probability value for VaR calculation
# w		  value of initial investment
# value:
# ans		Value-at-Risk computed using resampled data

	q = mean(x[idx]) + sd(x[idx])*qnorm(p)
	VaR = (exp(q) - 1)*w
	VaR
}

```

## Bootstrapping normal VaR

```{r}
MSFT.VaR.boot = boot(MSFT, statistic = ValueAtRisk.boot, R=999)
MSFT.VaR.boot
```

## Bootstrapping normal VaR

```{r}
plot(MSFT.VaR.boot)
```

```{r}
boot.ci(MSFT.VaR.boot, conf=0.95, type=c("norm", "perc"))
```

95% confidence interval for 5% VaR is fairly wide.
