--- 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) ```