---
title: Some Demos of Bootstrapping
output: slidy_presentation
date: 12 February 2019
author: 36-402, Spring 2019, Section A
---

```{r, include=FALSE}
# General set up

# Load packages we may use later
library(knitr)

# Set knitr options for knitting code into the report:
# - Save results so that code blocks aren't re-run unless code changes (cache),
# _or_ a relevant earlier code block changed (autodep), but don't re-run if the
# only thing that changed was the comments (cache.comments)
# - Don't clutter R output with messages or warnings (message, warning)
  # This _will_ leave error messages showing up in the knitted report
opts_chunk$set(cache=TRUE, autodep=TRUE, cache.comments=FALSE,
               message=FALSE, warning=FALSE)
```

## An example: Let's predict the stock market! {.smaller}

```{r}
# Load some data
  # Loads as a time-series structure which we'll undo in a moment
require(pdfetch)
sp <- pdfetch_YAHOO("SPY", fields="adjclose",
  from=as.Date("1993-02-09"), to=as.Date("2019-02-09"))
# Compute log returns for each day
  # as.numeric() gets rid of the time-series structure
sp <- as.numeric(diff(log(sp)))
# need to drop the initial NA which makes difficulties later
sp <- sp[-1]
# Make into a two-column data frame, one column for day $t$, other for day $t+1$
  # So columns should have all but the last entry in sp vs. all but the first
sp.df <- data.frame(Today=head(sp,-1), Tomorrow=tail(sp,-1))
# Fit a linear model, see the coefficients
sp.lm <- lm(Tomorrow ~ Today, data=sp.df)
coefficients(sp.lm)
```

```{r}
# Plot tomorrow versus today
plot(Tomorrow ~ Today, data=sp.df, pch=19, col="darkgrey")
# Add a nice horizontal line
abline(h=0,lty=3)
# And the estimated model
abline(sp.lm,lwd=2)
```

Let's try putting some error bars on this, with bootstrapping



## What's going on in bootstrapping?

We want to get the distribution of $\tau(X)$ when $X \sim P$ by simulating $\tilde{X} \sim \hat{P}$, and calculating $\tau(\tilde{X})$ repeatedly

#### Three parts to every bootstrapping exercise

1. A way to simulate $\hat{P}$ to get $\tilde{X}$
2. A way to calculate $\tau$ on $\tilde{X}$
3. A way to summarize the distribution of $\tau(\tilde{X})$



## Some functions for bootstrapping {.smaller}

```{r, include=FALSE}
source("http://www.stat.cmu.edu/~cshalizi/uADA/19/lectures/bootstrap.R")
```

```{r}
bootstrap
```

```{r}
rboot
```

```{r}
bootstrap.se
```

## Start with `rboot`

```{r}
rboot
```

- What this calls for:
    + A `simulator` function to make pseudo-data (generates $\tilde{X})$
	+ A `statistic` function to calculate on each simulation (calculates $\tau)$
- We need to write these

## Simulate the linear model of stock returns

- We'll resample residuals

```{r}
resample
```

```{r}
# After resample.residuals.penn in the textbook
resample.residuals.sp <- function() {
    new.frame <- sp.df
    new.Tomorrows <- fitted(sp.lm) + resample(residuals(sp.lm))
    new.frame$Tomorrow <- new.Tomorrows
    return(new.frame)
}
```

## Example with that simulator

```{r}
plot(Tomorrow ~ Today, data=sp.df, pch=19, col="darkgrey")
```



## Example with that simulator

```{r}
plot(Tomorrow ~ Today, data=sp.df, pch=19, col="darkgrey")
sim.1 <- resample.residuals.sp()
points(x=sim.1$Today, y=sim.1$Tomorrow, pch=19, col="orange", cex=0.5)
```



## Example with that simulator

```{r}
plot(Tomorrow ~ Today, data=sp.df, pch=19, col="darkgrey")
sim.1 <- resample.residuals.sp()
points(x=sim.1$Today, y=sim.1$Tomorrow, pch=19, col="orange", cex=0.5)
sim.2 <- resample.residuals.sp()
points(x=sim.2$Today, y=sim.2$Tomorrow, pch=19, col="darkorange", cex=0.5)
```



## What do we want to calculate? {.smaller}

> - Coefficients?
> - Predictions for each observed value?
> - Predictions at _arbitrary_ values?
> - Whatever it is, write a function to calculate it
>    + Make it read in the output from the simulator

## An example of a `statistic`

```{r}
selected.today.returns <- data.frame(Today=seq(from=-0.1,
                                               to=0.1,
                                               length.out=5))

sp.selected.preds <- function(data) {
    mdl <- lm(Tomorrow ~ Today, data=data)
    preds <- predict(mdl, newdata=selected.today.returns)
}
```

```{r}
signif(sp.selected.preds(sp.df), 3) # How'd you check this?
signif(sp.selected.preds(sim.1), 3)
```


## Put these together

```{r}
# You don't need to run this inside system.time()
system.time(my.first.bootstrap <- rboot(statistic=sp.selected.preds,
                                    simulator=resample.residuals.sp,
                                    B=800))
dim(my.first.bootstrap)
# Rows are different dimensions of the statistic
# Columns are different bootstrap runs
```

## What does this look like?

```{r}
plot(Tomorrow ~ Today, data=sp.df, pch=19, col="darkgrey",
     ylim=range(my.first.bootstrap))
abline(sp.lm, lwd=2)
for (repnum in 1:ncol(my.first.bootstrap)) {
    lines(x=selected.today.returns$Today,
           y=my.first.bootstrap[,repnum], lwd=0.1,
           col="blue")
}    
```

> - It's a huge pile of numbers
> - We need summaries
> - The summaries could be standard deviations, confidence intervals, etc.
>    + We need a function to calculate summaries

##

```{r}
sp.selected.preds.ses <- bootstrap(tboots = my.first.bootstrap,
                                summarizer=sd)
signif(sp.selected.preds.ses[1,],3)
# Why is this symmetric around 0?
```

##

```{r}
plot(Tomorrow ~ Today, data=sp.df, pch=19, col="darkgrey")
abline(sp.lm, lwd=2)
segments(x0=selected.today.returns$Today,
         y0=sp.selected.preds(sp.df)-2*sp.selected.preds.ses,
         x1=selected.today.returns$Today,
         y1=sp.selected.preds(sp.df)+2*sp.selected.preds.ses,
         lwd=2)
```


## We can use a different simulator

```{r}
resample.data.frame
```

```{r}
resample.sp <- function() { resample.data.frame(sp.df) }
```


## Example of simulating by resampling

```{r}
plot(Tomorrow ~ Today, data=sp.df, pch=19, col="darkgrey")
```

## Example of simulating by resampling

```{r}
plot(Tomorrow ~ Today, data=sp.df, pch=19, col="darkgrey")
sim.3 <- resample.sp()
# Jitter just a little to see multiple points
points(x=jitter(sim.3$Today), y=jitter(sim.3$Tomorrow), pch=19, col="lightgreen", cex=0.5)
```

## Example of simulating by resampling

```{r}
plot(Tomorrow ~ Today, data=sp.df, pch=19, col="darkgrey")
sim.3 <- resample.sp()
# Jitter just a little to see multiple points
points(x=jitter(sim.3$Today), y=jitter(sim.3$Tomorrow), pch=19, col="lightgreen", cex=0.5)
sim.4 <- resample.sp()
points(x=jitter(sim.4$Today), y=jitter(sim.4$Tomorrow), pch=19, col="darkgreen",
       cex=0.3)
```

## DISCUSSION

- How would we get standard errors for those five predicted values using the
resampling simulator?

## SOLUTION

```{r}
system.time(my.second.bootstrap <- rboot(statistic=sp.selected.preds,
                                    simulator=resample.sp,
                                    B=800))
sp.selected.preds.ses.resampled <- bootstrap(tboots = my.second.bootstrap,
                                summarizer=sd)

signif(sp.selected.preds.ses.resampled, 3)
# Why are those not symmetric around 0?
```

##

```{r}
plot(Tomorrow ~ Today, data=sp.df, pch=19, col="darkgrey")
abline(sp.lm)
segments(x0=selected.today.returns$Today,
         y0=sp.selected.preds(sp.df)-2*sp.selected.preds.ses.resampled,
         x1=selected.today.returns$Today,
         y1=sp.selected.preds(sp.df)+2*sp.selected.preds.ses.resampled,
         lwd=2)
```

