---
title: Model-Agnostic Prediction
bibliography: locusts.bib
output: slidy_presentation
---

\[
\newcommand{\Expect}[1]{\mathbb{E}\left[ #1 \right]}
\DeclareMathOperator*{\argmin}{argmin}
\newcommand{\Var}[1]{\mathrm{Var}\left[ #1 \right]}
\newcommand{\Prob}[1]{\mathbb{P}\left( #1 \right)}
\]


## Previously

- Linear autoregressive models
\[
X(t+1) = \epsilon(t+1) + \sum_{i=1}^{p}{\beta_i X(t+1-i)}
\]
$\Rightarrow$ obvious linear prediction $\sum_{i=1}^{p}{\beta_i X(t+1-i)}$
- We can fit an AR($p$) even when we don't think that's the right model
- Today: nonlinear predictions


## Agenda

- Reminder about regression functions
- Nonlinear autoregressions
    + Time-delay coordinates
    + "The method of analogs"
	+ Smoothing
	+ Additive autoregressions
- Conditional density estimation


## Nonlinear Prediction

- We see $X(1:t) \equiv (X(1), X(2), \ldots X(t-1), X(t))$
- We want a guess at $X(t+1)$
- This guess will be some function of $X(1:t)$
- We are now allowing it to be a _nonlinear_ function


## Regression and Conditional Expectations

> - Remember that if we use mean squared error, the optimal function is the conditional expectation
\[
\Expect{Y|Z} = \argmin_{m: ~ \mathcal{Z} \mapsto \mathbb{R}}{\Expect{(Y-m(Z))^2}}
\]
>    + Expected squared error would be $\Var{Y|Z}$
> - So the function we'd ideally like is
\[
\Expect{X(t+1)|X(1:t)}
\]
> - We could estimate this if we had many realizations of the process
>    + In general we don't get many realizations
> - If $X$ is a Markov process we'd only need $\Expect{X(t+1)|X(t)}$
>    + In general $X$ isn't Markov


## Limited Memory Approximations

- Even if $X$ isn't Markov, maybe the distant past doesn't matter very much any more
- So we _pick_ an **order** or **memory length** $p$, and try to predict $X(t+1)$ from $X(t-p+1), \ldots X(t-1), X(t)$
- Abbreviate $X(t-p+1, \ldots X(t)$ as $S(t)$, the **time-delay vector** (or **time-delay coordinates**)
    + Keeping $p$ implicit in the notation
	+ Come back to how to pick $p$
- We now want to regress $X(t+1)$ on $S(t)$
- This now looks like a pretty ordinary nonlinear regression problem


## Nearest Neighbors

- Nearest neighbor regression: we've got data $(z_1, y_1), \ldots (z_n, y_n)$, we see $Z=z$ and we want to predict $Y$
- Find the $z_i$ closest to $z$, say $z_{NN}$, then predict $y_{NN}$
- "What does this most remind me of?"
- In our context, $Z$ is $S(t)$ and $Y$ is $X(t+1)$
- For a stationary, ergodic process, [the nearest neighbor will get closer and closer to $S(t)$](http://www.stat.cmu.edu/~cshalizi/dm/20/lectures/03/lecture-03.html#31_convergence_of_the_nearest_neighbor), and the [expected squared error approaches $2\Var{X(t+1)|S(t)}$](http://www.stat.cmu.edu/~cshalizi/dm/20/lectures/03/lecture-03.html#4_putting_the_noise_back_in_for_1nn)
    + Remember that the true regression function would have mean squared error $\Var{X(t+1)|S(t)}$
- For time series, sometimes called the "method of analogs"


## $k$ Nearest Neighbors

- Why use just one historical case?  Why not use say $k$ of them, and take
the average?
- This is $k$ nearest neighbors (kNN)
- As $n\rightarrow\infty$, the expected squared error approaches $(1+\frac{1}{k})\Var{X(t+1)|S(t)}$
- But at finite $n$, larger $k$ always means more bias and less variance than smaller $k$
- We usually want to let $k\rightarrow\infty$ as $n\rightarrow\infty$
- Select $k$ at finite $n$ by cross-validation
    + Or some other way of approximating generalization error, like a suitable bootstrap [@Lunde-CRS-bootstrap-GEB]


## Kernel Smoothing

- Often nice to get smoother predictions than $kNN$ will give us
- Kernel regression/kernel smoothing/Nadaraya-Watson regression:
\[
\hat{y}(z) = \frac{\sum_{i=1}^{n}{y_i K(\|z-z_i\|/h)}}{\sum_{j=1}^{n}{K(\|z-z_j\|/h)}}
\]
where $K$ is the **kernel function**, e.g., a standard Gaussian pdf, and $h$ is the **bandwidth**
- This is a weighted average of the $y_i$s, where the weights are $\propto K(\|z-z_i\|/h)$
    + More weight on near-by points
	+ Weights change smoothly with distance, not dis-continuously as with kNN
- Big $h$: very smooth, little detail, high bias, low variance
- Small $h$: less smooth, lots of detail, low bias, high variance
- Again, let $h\rightarrow 0$ as $n\rightarrow\infty$, pick $h$ by CV


## Curse of Dimensionality

> - Fitting high-dimensional regressions is intrinsically hard
> - With $p$ dimensions, the MSE has to be at least $\Var{Y|Z} + O(n^{-4/(4+p)})$ _in general_
>     + The **minimax rate**, the one which minimizes the maximum excess error
> - You can do better at learning some regression functions only if you do worse at learning others
>    + E.g., if you use linear regression and the truth is linear, you get a rate of $O(p/n)$, but if the truth is non-linear you get a rate of $O(1)$
> - With $p=1$, the excess is $O(n^{-0.8})$ which isn't so bad
> - With $p=100$ the excess is $O(n^{-`r signif(4/(4+100), 2)`})$ which hardly moves at all
> - As $p$ grows, the amount of data you need grows exponentially
> - _unless_ we **either** have special structure **or** live with approximation error



## Additive Autoregressions

- Try fitting an **additive model**,
\[
\hat{X}(t+1) = \alpha + \sum_{i=1}^{p}{f_i(X(t+1-i))}
\]
- The functions $f_i$ are the **partial response** functions
- We can find each $f_i$ by nonparametric regression (kNN, kernel smoothing, splines...)
- Each $\hat{f}_i$ converges fairly quickly, $O(n^{-0.8})$
- If the real regression function isn't additive, there's still an $O(1)$ approximation error
- Additive autoregressions include linear autoregressions as a special case; they are almost as fast to fit; and they are just as easy to interpret

```{r, include=FALSE}
# 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))
}


# Fit an additive autoregressive model
  # additive model fitting is outsourced to mgcv::gam, with splines
# Inputs: time series (x), order of autoregression (order)
# Output: fitted GAM object
aar <- function(ts,order) {
  stopifnot(require(mgcv))
  # Automatically generate a suitable data frame from the time series
  # and a formula to go along with it
  fit <- gam(as.formula(auto.formula(order)),
    data=design.matrix.from.ts(ts,order))
  return(fit)
}

# Generate formula for an autoregressive GAM of a specified order
# Input: order (integer)
# Output: a formula which looks like
  # "lag0 ~ s(lag1) + s(lag2) + ... + s(lagorder)"
auto.formula <- function(order) {
  inputs <- paste("s(lag",1:order,")",sep="",collapse="+")
    form <- paste("lag0 ~ ",inputs)
    return(form)
}
```

## Fitting an AAR(2) to the Lynx Data

```{r, include=FALSE}
library(mgcv)
```

```{r, echo=FALSE}
lynx.aar2 <- aar(lynx,2)
plot(lynx.aar2,pages=1)
```

- Note: a high count two years ago predicts a _low_ count this year $\Rightarrow$ cycling


## Fitting an AAR(2) to the Lynx Data

```{r, echo=FALSE}
lynx.aar2b <- aar(lynx[1:80],2)
out.of.sample <- design.matrix.from.ts(lynx[-(1:78)],2)
lynx.preds <- predict(lynx.aar2b,newdata=out.of.sample)
plot(lynx, main="Fitting an AAR(2) to the First 80 Years")
lines(1823:1900,fitted(lynx.aar2b),lty="dashed")
lines(1901:1934,lynx.preds,col="grey")
legend("topright", legend=c("Data", "In-sample", "Out-of-sample"),
       lty=c("solid", "dashed", "solid"), col=c("black", "black", "grey"),
       cex=0.8)
```



## What About the Distribution?

- Suppose we want to find $p(y|z)$ from our $(z_i, y_i)$ data
- Use **kernel conditional density estimation**
\[
\hat{p}(y|z) = \frac{\hat{p}(z,y)}{\hat{p}(z)} = \frac{\frac{1}{nh_{Z} h_{Y}}\sum_{i=1}^{n}{K(|y-y_i|/h_Y) K(\|z-z_i\|/h_Z)}}{\frac{1}{nh_{Z}}\sum_{j=1}^{n}{K(\|z-z_j\|/h)}}
\]
- Pick both $h_Z$ and $h_Y$ by cross-validation [@Hall-Racine-Li-CV]
- Again, for us $Y$ is played by $X(t+1)$, and $Z$ by the history or time-delay vector $S(t)$
- We did this for the lynxes back in [Lecture 19](http://www.stat.cmu.edu/~cshalizi/dst/20/lectures/19/lecture-19.html#(24))


## Choosing the Memory Order $p$

- Ideally, we'd let $p \rightarrow\infty$ as $n\rightarrow\infty$
    + An exception: if $X(t)=f(U(t))$ and $U(t+1) = g(U(t))$ (we're observing a function of _deterministic_ dynamics), and $U(t)$ has dimension $m$, then we only need $p=2m+1$ [@Geometry-from-a-time-series;@Takens-embedding;@Kantz-Schreiber-2nd]
- In practice, picking a big $p$ leads to curse-of-dimensionality trouble
- Some information theory suggests using $p=O(\log{n})$, at least if we want to predict the conditional density [@Marton-Shields]
- Recommendation: treat $p$ as an additional control setting, like the $k$ for nearest neighbors or the $h$ for kernel smoothing




## When Will These Tricks Work?

- Nonparametric regression converges on $\Expect{Y|Z}$ when we get IID $(Z,Y)$ pairs as $n\rightarrow\infty$
- So nonparametric regression for time series works better the closer we get to independent and identically distributed data
- _Independent_: Replace with **weak dependence**, $X(t)$ and $X(s)$ should approach independence as $|t-s|\rightarrow\infty$
    + Often really want $X(t:t+k)$ and $X(s:s+k)$ to approach independence for abritrary $k$
- _Identically distributed_: Replace with **strict stationarity** or **strong stationarity**, $\Prob{X(t:t+k) = x(t:t+k)} = \Prob{X(s:s+k) = x(t:t+k)}$ for all $k$
- So these methods will work best for weakly-dependent, stationary processes [@Bosq-nonparametric]
    + They'll work _eventually_ for stationary, ergodic processes which aren't weakly-dependent, provided $n\rightarrow\infty$ [@Wiener-prediction]
- Non-stationary or strongly-dependent processes really need explicit models (next time)


## References

