Agenda

The sampling distribution of a statistic tells us about statistical uncertainty (standard errors, biases, confidence sets)
The bootstrap principle: approximate the sampling distribution by simulating from a good model, and treating the simulation run just like the data
Sometimes we simulate from the model we’re estimating (model-based or “parametric” bootstrap)
Sometimes we simulate by re-sampling the original data (resampling or “nonparametric” bootstrap)
Stronger assumptions \(\Rightarrow\) less uncertainty if we’re right

Statistical Uncertainty

Re-run the experiment (survey, census, …) and get different data
\(\therefore\) everything we calculate from data (estimates, test statistics, policies, …) will change from trial to trial as well
This variability is (the source of) statistical uncertainty
Quantifying this is a way to be honest about what we actually know

Measures of Uncertainty

Standard error = standard deviation of an estimator
- (could equally well use median absolute deviation, etc.)
\(p\)-value = Probability we’d see a signal this big if there was just noise
Confidence region = All the parameter values we can’t reject at low error rates:

Either the true parameter is in the confidence region
or we are very unlucky
or our model is wrong

etc., etc.

The Sampling Distribution Is the Source of All Knowledge

Data \(X \sim P_{X,\theta_0}\), for some true \(\theta_0\)
We calculate a statistic \(T = \tau(X)\) so it has distribution \(P_{T,\theta_0}\)
If we knew \(P_{T,\theta_0}\), we could calculate
- \(\Var{T}\) (and so standard error)
- \(\Expect{T}\) (and so bias)
- quantiles (and so confidence intervals or \(p\)-values), etc.

The Problems

Most of the time, \(P_{X,\theta_0}\) is very complicated
Most of the time, \(\tau\) is a very complicated function
We certainly don’t know \(\theta_0\)

Upshot: We don’t know \(P_{T,\theta_0}\) and can’t use it to calculate anything

The Solutions

Classically (\(\approx 1900\)–\(\approx 1975\)): Restrict the model and the statistic until you can calculate the sampling distribution, at least for very large \(n\)
- This is the asymptotic theory we did in Lectures 14 and 15
Modern (\(\approx 1975\)–): Use complex models and statistics, but simulate calculating the statistic on the model

The Bootstrap Principle

Find a good estimate \(\hat{P}\) for \(P_{X,\theta_0}\)
Generate a simulation \(\tilde{X}\) from \(\hat{P}\), set \(\tilde{T} = \tau(\tilde{X})\)
Use the simulated distribution of the \(\tilde{T}\) to approximate \(P_{T,\theta_0}\)

Refinements:

improving the initial estimate \(\hat{P}\)
reducing the number of simulations or speeding them up
transforming \(\tau\) so the final approximation is more stable

First step: find a good estimate \(\hat{P}\) for \(P_{X,\theta_0}\)

Model-based Bootstrap

If we are using a model, our best guess at \(P_{X,\theta_0}\) is \(P_{X,\hat{\theta}}\), with our best estimate \(\hat{\theta}\) of the parameters

The Model-based Bootstrap

Get data \(X\), estimate \(\hat{\theta}\) from \(X\)
Repeat \(m\) times:
- Simulate \(\tilde{X}\) from \(P_{X,\hat{\theta}}\) (simulate data of same size/“shape” as real data)
- Calculate \(\tilde{T} = \tau(\tilde{X}\)) (treat simulated data the same as real data)
Use empirical distribution of \(\tilde{T}\) as \(P_{T,\theta_0}\)

Example: lynxes in Canada

Annual time series of number of lynxes trapped for the fur company in the Mackenzie River area of Canada:

data(lynx)
plot(lynx, xlab="Year", ylab="Lynxes")

Fitting an AR model to the example data

This oscillates, so an AR(1) isn’t appropriate — maybe an AR(2)

lynx.ar2 <- ar.ols(lynx, order.max=2, aic=FALSE, demean=FALSE, intercept=TRUE)
lynx.ar2$ar

## , , 1
## 
##           [,1]
## [1,]  1.152423
## [2,] -0.606229

lynx.ar2$x.intercept

## [1] 710.1056

so the estimated model is \[ X(t) = 710.1055888 + 1.1524226 X(t-1) - 0.606229 X(t-2) + \epsilon(t) \]

Simulating from the fitted model

What do some time series from this look like?
- (We wrote ar2.sim last time)

lynx.sims <- replicate(5, ar2.sim(n=length(lynx),
                                   a=lynx.ar2$x.intercept,
                                   b=lynx.ar2$ar,
                                   sd.innovation=sd(lynx.ar2$resid, na.rm=TRUE),
                                   x.start=lynx[1:2]))
plot(lynx, xlab="Year", ylab="Lynxes caught")
for (i in 1:ncol(lynx.sims)) {
    lines(1821:1934, lynx.sims[,i], col="grey")
}

Uncertainty in the coefficients

Write a function which estimates the coefficients

# Estimate the coefficients of an AR(2) model
# Inputs: a time series (x)
# Outputs: vector containing intercept and 2 AR slopes
ar2.estimator <- function(x) {
    # Treat the input exactly the same way we treated the real data
    fit <- ar.ols(x, order.max=2, aic=FALSE, demean=FALSE, intercept=TRUE)
    return(c(a=fit$x.intercept, b1=fit$ar[1], b2=fit$ar[2]))
}

Check that this gives the right answer on the data, and a different answer on a simulation:

ar2.estimator(lynx)

##          a         b1         b2 
## 710.105589   1.152423  -0.606229

ar2.estimator(lynx.sims[,1])

##           a          b1          b2 
## 762.6063220   1.1102328  -0.5743644

Uncertainty in the coefficients

Now do lots of simulations, and apply this function to each simulation:

lynx.coef.sims <- replicate(1000, ar2.estimator( ar2.sim(n=length(lynx),
                                   a=lynx.ar2$x.intercept,
                                   b=lynx.ar2$ar,
                                   sd.innovation=sd(lynx.ar2$resid, na.rm=TRUE),
                                   x.start=lynx[1:2])))
dim(lynx.coef.sims)

## [1]    3 1000

lynx.coef.sims[1:3, 1:5]

##           [,1]        [,2]        [,3]       [,4]        [,5]
## a  830.8435812 724.7291944 598.4342755 761.794216 885.9718640
## b1   1.1771914   1.1892204   1.2535834   1.109654   0.9910510
## b2  -0.7153194  -0.6591029  -0.6459415  -0.578512  -0.5180509

Uncertainty in the coefficients

Standard errors \(\approx\) standard deviation across simulations:

apply(lynx.coef.sims, 1, sd)

##            a           b1           b2 
## 118.50209000   0.07622429   0.07498753

95% confidence intervals \(\approx\) quantiles across simulations:

apply(lynx.coef.sims, 1, quantile, prob=c(0.025, 0.975))

##              a        b1         b2
## 2.5%  508.0321 0.9839003 -0.7473846
## 97.5% 963.5467 1.2820643 -0.4422757

For Gaussian AR\((p)\) models…

… the bootstrapping isn’t needed, because there’s an asymptotic theory which tells us about standard errors and confidence intervals
- Like what you did for the AR(1) in homework 8, only with even more book-keeping
But we can use this approach for any generative time-series model
- E.g., make \(X(t)\) have a Poisson distribution rather than a Gaussian
- All we have to do is swap out the simulation

Prediction

You can apply these ideas when the “statistic” is \(X(t)\) for a fixed \(t\)
Or when it’s \(X(t) | X(t-1) = x\) for a fixed \(x\)
- Since the model tells us how to calculate that!
Expected values across simulations = best guess for the point prediction
Variances and intervals incorporate modeling uncertainty

Resampling

Problem: Suppose we don’t have a trust-worthy model
Resource: We do have data, which tells us a lot about the distribution
Solution: Resampling, treat the sample like a whole population

The Resampling (“Nonparameteric”) Bootstrap

Get data \(X\), containing \(n\) samples, find \(T=\tau(X)\)
Fix a block length \(k\) and divide \(X\) into (overlapping) blocks of length \(k\)
- So first block is \((X(1), \ldots X(k))\), second block is \((X(2), \ldots X(k+1))\), etc.
Repeat \(m\) times:
- Generate \(\tilde{X}\) by drawing \(n/k\) blocks from \(X\) with replacement (resample the data)
- Concatenate the blocks to make \(\tilde{X}\)
- Calculate \(\tilde{T} = \tau(\tilde{X})\) (treat simulated data the same as real data)
Use empirical distribution of \(\tilde{T}\) as \(P_{T,\theta}\)

Example: Back to the lynxes

# Convert a time series into a data frame of lagged values
# Input: a time series, maximum lag to use, whether older values go on the right
  # or the left
# Output: a data frame with (order+1) columns, named lag0, lag1, ... , and
  # length(ts)-order rows
design.matrix.from.ts <- function(ts,order,right.older=TRUE) {
  n <- length(ts)
  x <- ts[(order+1):n]
  for (lag in 1:order) {
    if (right.older) {
      x <- cbind(x,ts[(order+1-lag):(n-lag)])
    } else {
      x <- cbind(ts[(order+1-lag):(n-lag)],x)
    }
  }
  lag.names <- c("lag0",paste("lag",1:order,sep=""))
  if (right.older) {
    colnames(x) <- lag.names
  } else {
    colnames(x) <- rev(lag.names)
  }
  return(as.data.frame(x))
}

head(lynx)

## [1]  269  321  585  871 1475 2821

head(design.matrix.from.ts(lynx,2))

##   lag0 lag1 lag2
## 1  585  321  269
## 2  871  585  321
## 3 1475  871  585
## 4 2821 1475  871
## 5 3928 2821 1475
## 6 5943 3928 2821

Example: Back to the lynxes

# Simple block bootstrap
# Inputs: time series (ts), block length, length of output
# Output: one resampled time series
# Presumes: ts is a univariate time series
rblockboot <- function(ts,block.length,len.out=length(ts)) {
  # chop up ts into blocks
  the.blocks <- as.matrix(design.matrix.from.ts(ts,block.length-1,
    right.older=FALSE))
    # look carefully at design.matrix.from.ts to see why we need the -1
  # How many blocks is that?
  blocks.in.ts <- nrow(the.blocks)
  # Sanity-check
  stopifnot(blocks.in.ts == length(ts) - block.length+1)
  # How many blocks will we need (round up)?
  blocks.needed <- ceiling(len.out/block.length)
  # Sample blocks with replacement
  picked.blocks <- sample(1:blocks.in.ts,size=blocks.needed,replace=TRUE)
  # put the blocks in the randomly-selected order
  x <- the.blocks[picked.blocks,]
  # convert from a matrix to a vector and return
    # need transpose because R goes from vectors to matrices and back column by
    # column, not row by row
  x.vec <- as.vector(t(x))
    # Discard uneeded extra observations at the end silently
  return(x.vec[1:len.out])
}

head(lynx)

## [1]  269  321  585  871 1475 2821

head(rblockboot(lynx, 3))

## [1] 1676 2251 1426  473  358  784

Resampling the lynxes

We can do everything the same as before, just changing how we simulate

lynx.coef.sims.resamp <- replicate(1000, ar2.estimator(rblockboot(lynx, 3)))

Standard errors:

apply(lynx.coef.sims.resamp, 1, sd)

##            a           b1           b2 
## 198.36971077   0.09812245   0.07406839

95% confidence intervals:

apply(lynx.coef.sims.resamp, 1, quantile, prob=c(0.025, 0.975))

##               a        b1          b2
## 2.5%   598.7163 0.3864529 -0.35398540
## 97.5% 1365.3888 0.7632487 -0.06106865

(Anything funny here?)

Why should this work?

Preserves the distribution within blocks
Independence across blocks
So \(\tilde{X} \approx\) stationary, with the same short-range distribution as \(X\) but no longer-range dependence
Should work well if \(X\) is stationary without long-range dependence
Want to let \(k\) grow with \(n\), but still have a growing number of blocks, so \(n/k_n \rightarrow \infty\), or \(k=o(n)\)
Some theory suggests \(k \propto n^{1/3}\)
- Politis and White (2004) provides an automatic procedure for picking \(k\)
- R implementation at [http://www.math.ucsd.edu/~politis/SOFT/PPW/ppw.R]
- …which suggests a block-length of 3 for the lynxes

What about non-stationary time series?

Model-based: OK if the model is non-stationary
Otherwise:
- Estimate the trend
- Bootstrap the deviations from the trend (by model or resampling)
- Add these back to the trend
- \(\tilde{X} =\) trend + bootstrapped deviations

What about spatial data?

Model-based: simulate the model
Resampling: use spatial blocks (e.g., rectangles)

Sources of Error in Bootstrapping

Simulation Using only \(m\) bootstrap replicates
- Make this small by letting \(m\rightarrow\infty\)
- Costs computing time
- Diminishing returns: error is generally \(O(1/\sqrt{m})\)
Approximation Using \(\hat{P}\) instead of \(P_{X,\theta_0}\)
- Make this small by careful statistical modeling
Estimation Only a finite number of samples
- Make this small by being careful about what we simulate (example in backup slides)

Generally: for fixed \(n\), resampling boostrap shows more uncertainty than model-based bootstraps, but is less at risk to modeling mistakes (yet another bias-variance tradeoff)

Summary

Bootstrapping is applying the Monte Carlo idea to the sampling distribution
- Find a good approximation to the data-generating distribution
- Simulate it many times
- Treat each simulation run just like the original data
- Distribution over simulations \(\approx\) sampling distribution
Choice of how to simulate:
- Based on a model
- By resampling blocks
Many, many refinements, but this is the core of it

Backup: The “surrogate data” method

You think that \(X\) comes from some interesting class of processes \(\mathcal{H}_1\)
Your skeptical friend thinks maybe it’s just some boring process in \(\mathcal{H}_0\)
Fit a model to \(X\) using only the boring processes in \(\mathcal{H}_0\), say \(\hat{Q}\)
Now simulate \(\hat{Q}\) to get many “surrogate data” sets \(\tilde{X}\)
Look at the distribution of \(\tau(\tilde{X})\) over your simulations
If \(\tau(X)\) is really unlikely under that simulation, it’s evidence against \(\mathcal{H}_0\)
- When would it be evidence for \(\mathcal{H}_1\)?
The term “surrogate data” comes out of physics (Kantz and Schreiber 2004)
- testing the null hypothesis \(\mathcal{H}_0\) against the alternative \(\mathcal{H}_1\)

Backup: Incorporating uncertainty in the initial conditions

Simulating a dynamical model needs an initial value \(X(0)\) (or \(X(1)\), or maybe \(X(1), X(2)\), etc.)
What if this isn’t known precisely?
Easy case: it follows a probability distribution
- Then draw \(X(0)\) from that distribution independently on each run
- Aggregate across runs as usual
Less easy: \(X(0)\) within known limits but no distribution otherwise
- My advice: draw \(X(0)\) uniformly from inside the region
- Rationale 1: Maximum entropy distribution compatible with constraints
- Rationale 2: If the process is rapidly mixing, memory of the distribution of \(X(0)\) decays quickly anyway

Backup: Improving on the Crude Confidence Interval

Crude CI use distribution of \(\tilde{\theta}\) under \(\hat{\theta}\)

But really we want the distribution of \(\hat{\theta}\) under \(\theta\)

Mathematical observation: Generally speaking, (distributions of) \(\tilde{\theta} - \hat{\theta}\) is closer to \(\hat{\theta}-\theta_0\) than \(\tilde{\theta}\) is to \(\hat{\theta}\)

(“centering” or “pivoting”)

Backup: The Basic, Pivotal CI

quantiles of \(\tilde{\theta}\) = \(q_{\alpha/2}, q_{1-\alpha/2}\)

\[ \begin{eqnarray*} 1-\alpha & = & \Probwrt{q_{\alpha/2} \leq \tilde{\theta} \leq q_{1-\alpha/2}}{\hat{\theta}} \\ & = & \Probwrt{q_{\alpha/2} - \hat{\theta} \leq \tilde{\theta} - \hat{\theta} \leq q_{1-\alpha/2} - \hat{\theta}}{\hat{\theta}} \\ & \approx & \Probwrt{q_{\alpha/2} - \hat{\theta} \leq \hat{\theta} - \theta_0 \leq q_{1-\alpha/2} - \hat{\theta}}{\theta_0}\\ & = & \Probwrt{q_{\alpha/2} - 2\hat{\theta} \leq -\theta_0 \leq q_{1-\alpha/2} - 2\hat{\theta}}{\theta_0}\\ & = & \Probwrt{2\hat{\theta} - q_{1-\alpha/2} \leq \theta_0 \leq 2\hat{\theta}-q_{\alpha/2}}{\theta_0} \end{eqnarray*} \]

Basically: re-center the simulations around the empirical data

2*ar2.estimator(lynx)-apply(lynx.coef.sims.resamp, 1, quantile, prob=0.975)

##         a        b1        b2 
## 54.822397  1.541596 -1.151389

2*ar2.estimator(lynx)-apply(lynx.coef.sims.resamp, 1, quantile, prob=0.025)

##           a          b1          b2 
## 821.4949207   1.9183923  -0.8584727

Backup: Further reading

Bootstrap begins with Efron (1979) (drawing on an older idea from the 1940s called the “jackknife”)
Block bootstrap for time series comes from Künsch (1989)
Bühlmann (2002) is a good review on time series bootstraps
Davison and Hinkley (1997) is the best over-all textbook on the bootstrap
Lahiri (2003) is very theoretical

References

Bühlmann, Peter. 2002. “Bootstraps for Time Series.” Statistical Science 17:52–72. https://doi.org/10.1214/ss/1023798998.

Davison, A. C., and D. V. Hinkley. 1997. Bootstrap Methods and Their Applications. Cambridge, England: Cambridge University Press.

Efron, Bradley. 1979. “Bootstrap Methods: Another Look at the Jackknife.” Annals of Statistics 7:1–26. https://doi.org/10.1214/aos/1176344552.

Kantz, Holger, and Thomas Schreiber. 2004. Nonlinear Time Series Analysis. Second. Cambridge, England: Cambridge University Press.

Künsch, Hans R. 1989. “The Jackknife and the Bootstrap for General Stationary Observations.” Annals of Statistics 17:1217–41. https://doi.org/10.1214/aos/1176347265.

Lahiri, S. N. 2003. Resampling Methods for Dependent Data. New York: Springer-Verlag.

Politis, Dimitris N., and Halbert White. 2004. “Automatic Block-Length Selection for the Dependent Bootstrap.” Econometric Reviews 23:53–70. https://doi.org/10.1081/ETC-120028836.

Simulation for Inference I — The Bootstrap

In our last episode

Agenda

Statistical Uncertainty

Measures of Uncertainty

The Sampling Distribution Is the Source of All Knowledge

The Problems

The Solutions

The Bootstrap Principle

Model-based Bootstrap

The Model-based Bootstrap

Example: lynxes in Canada

Example: lynxes in Canada

Fitting an AR model to the example data

Simulating from the fitted model

Uncertainty in the coefficients

Uncertainty in the coefficients

Uncertainty in the coefficients

For Gaussian AR\((p)\) models…

Prediction

Resampling

The Resampling (“Nonparameteric”) Bootstrap

Example: Back to the lynxes

Example: Back to the lynxes

Resampling the lynxes

Why should this work?

What about non-stationary time series?

What about spatial data?

Sources of Error in Bootstrapping

Summary

Backup: The “surrogate data” method

Backup: Incorporating uncertainty in the initial conditions

Backup: Improving on the Crude Confidence Interval

Backup: The Basic, Pivotal CI

Backup: Further reading

References