$$\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)$, and the expected squared error approaches $2\Var{X(t+1)|S(t)}$
  • 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 and Shalizi 2017)

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^{-0.038})$ 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

Fitting an AAR(2) to the Lynx Data

• Note: a high count two years ago predicts a low count this year $\Rightarrow$ cycling

Fitting an AAR(2) to the Lynx Data

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, and Li 2004)
• 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

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$ (Packard et al. 1980; Takens 1981; Kantz and Schreiber 2004)
• 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 and Shields 1994)
• 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 1998)
  • They’ll work eventually for stationary, ergodic processes which aren’t weakly-dependent, provided $n\rightarrow\infty$ (Wiener 1956)
• Non-stationary or strongly-dependent processes really need explicit models (next time)