\[ \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)