Optimal Linear Prediction, Especially over Time
36-740
8 September 2026 (Lecture 5)
\[
\newcommand{\Expect}[1]{\mathbb{E}\left[ #1 \right]}
\newcommand{\Var}[1]{\mathrm{Var}\left[ #1 \right]}
\newcommand{\SampleVar}[1]{\widehat{\mathrm{Var}}\left[ #1 \right]}
\newcommand{\Cov}[1]{\mathrm{Cov}\left[ #1 \right]}
\newcommand{\TrueRegFunc}{\mu}
\newcommand{\EstRegFunc}{\widehat{\TrueRegFunc}}
\DeclareMathOperator{\tr}{tr}
\DeclareMathOperator*{\argmin}{argmin}
\DeclareMathOperator{\det}{det}
\newcommand{\TrueNoise}{\epsilon}
\newcommand{\EstNoise}{\widehat{\TrueNoise}}
\newcommand{\AutoCov}{\gamma}
\newcommand{\CovKernel}{\kappa}
\newcommand{\LinearExpect}[1]{\mathbb{L}\left[ #1 \right]}
\newcommand{\HilbertInner}[3]{\langle #1, #2\rangle_{#3}}
\newcommand{RKHSInner}[2]{\HilbertInner{#1}{#2}{\mathcal{K}}}
\]
In our previous episodes
- Linear smoothers
- Predictions are linear combinations of the data
- How to choose the weights?
- PCA
- Use covariance to break the data into additive components
- Fourier analysis
- Break the data into sinusoidal components
- Closely related to the covariance function
Today: use covariance to do prediction
Optimal prediction in general
What’s the best constant guess for a random variable \(Y\)?
\[\begin{eqnarray}
\TrueRegFunc & = & \argmin_{m \in \mathbb{R}}{\Expect{(Y-m)^2}}\\
& = & \argmin_{m \in \mathbb{R}}{\left(\Var{(Y-m)} + (\Expect{Y-m})^2\right)}\\
& = & \argmin_{m \in \mathbb{R}}{\left(\Var{Y} + (\Expect{Y} - m)^2\right)}\\
& = & \argmin_{m \in \mathbb{R}}{ (\Expect{Y} - m)^2}\\
& = & \Expect{Y}
\end{eqnarray}\]
Optimal prediction in general
- Now we get a covariate \(Z\) which takes values in some arbitrary space \(\mathcal{Z}\)
- What’s the best function of \(Z\) to guess for \(Y\)? \[\begin{eqnarray}
\TrueRegFunc & = & \argmin_{m: ~ \mathcal{Z} \mapsto \mathbb{R}}{\Expect{(Y-m(Z))^2}}\\
& = & \argmin_{m: ~ \mathcal{Z} \mapsto \mathbb{R}}{\Expect{\Expect{(Y-m(Z))^2|Z}}}\\
& = & \argmin_{m: ~ \mathcal{Z} \mapsto \mathbb{R}}{\left( \int_{\mathcal{Z}}{\Expect{(Y-m(z))^2|Z=z} p(z) dz}\right)}
\end{eqnarray}\]
For each \(z \in \mathcal{Z}\), best \(m(z)\) is \(\Expect{Y|Z=z}\) (by previous slide), so \[
\TrueRegFunc(z) = \Expect{Y|Z=z}
\]
Optimal prediction in general
- Learning arbitrary functions is hard!
- Who knows what the right function might be?
- What if we decide to make our predictions linear?
Optimal linear prediction with univariate predictor
- Our prediction will be of the form \[
m(z) = a + b z
\] and we want the best \(a, b\)
- When we want to emphasize that this is a linear approximation to the expectation, write \[
\LinearExpect{Y|Z}
\]
Optimal linear prediction, univariate case
\[
(\alpha, \beta) = \argmin_{a \in \mathbb{R}, b \in \mathbb{R}}{\Expect{(Y-(a+bZ))^2}}
\]
Remember: \(\Expect{W^2} = (\Expect{W})^2 + \Var{W}\)
\[\begin{eqnarray}
\Expect{(Y-(a+bZ))^2} & = & (\Expect{Y - (a+bZ)})^2 + \Var{Y - (a + bZ)}\\
& = & ((\Expect{Y} - b\Expect{Z}) - a)^2 + \Var{Y - bZ}\\
& = & ((\Expect{Y} - b\Expect{Z}) - a)^2 + \Var{Y} + b^2\Var{Z} - 2b\Cov{Y,Z}
\end{eqnarray}\]
The intercept
\[
\Expect{(Y-(a+bZ))^2} = ((\Expect{Y} - b\Expect{Z}) - a)^2 + \Var{Y} + b^2\Var{Z} - 2b\Cov{Y,Z}
\]
- The only place \(a\) shows up is inside the square
- We can always make the square zero with \(a = \Expect{Y} - b\Expect{Z}\)
- \(\therefore\) Once we find \(b\) we can always minimize the square
- \(\therefore\) \(\alpha = \Expect{Y} - \beta \Expect{Z}\)
- \(\therefore\) \(\therefore\) optimal linear predictor looks like \[
\Expect{Y} + \beta(Z-\Expect{Z})
\] \(\Rightarrow\) centering \(Z\) and/or \(Y\) won’t change the slope
The slope
\[\begin{eqnarray}
\left. \frac{\partial}{\partial b}\Expect{(Y-(a+bZ))^2}\right|_{a=\alpha, b=\beta} & = & 2(\Expect{Y} - (\alpha + \beta \Expect{Z})(-\Expect{Z})\\
& & + 2\beta\Var{Z} - 2\Cov{Y,Z}\\
0 & = & 2(\beta\Var{Z} - \Cov{Y,Z})\\
\beta \Var{Z} & = & \Cov{Y,Z}\\
\beta & = & \frac{\Cov{Y,Z}}{\Var{Z}}
\end{eqnarray}\]
The optimal linear predictor of \(Y\) from \(Z\)
The optimal linear predictor of \(Y\) from a single \(Z\) is always
\[
\LinearExpect{Y|Z} = \alpha + \beta Z = \Expect{Y} + \left(\frac{\Cov{Z,Y}}{\Var{Z}}\right) (Z - \Expect{Z})
\]
What did we not assume?
- That the true relationship between \(Y\) and \(Z\) is linear
- That anything is Gaussian
- That anything has constant variance
- That anything is independent or even uncorrelated
- NONE OF THAT MATTERS for the optimal linear predictor
A little worked example (I)
- We see \(X(r_1, t_1)\), for short \(X_1\)
- We want to guess \(X(r_0, t_0)\), for short \(X_0\)
- Assume: \(\Expect{X_0} = \Expect{X_1} = \mu\)
- Assume: \(\Var{X_0} = \Var{X_1} = \sigma^2\)
- Assume: \(\Cov{X_0, X_1} = \gamma\)
We know: \[
\LinearExpect{Y|Z} = \Expect{Y} + \frac{\Cov{Y, Z}}{\Var{Z}}(Z - \Expect{Z})
\]
Substituting in: \[
\LinearExpect{X_1} = \mu + \frac{\gamma}{\sigma^2}(X_1 - \mu)
\]
(Apparently first done by Karl Pearson’s former student F. E. Cave-Browne-Cave for weather forecasting (Cave-Browne-Cave 1905), see discussion in Klein (1997), pp. 262–264)
Some general properties of the optimal linear predictor
- The prediction errors average out to zero
- The prediction errors are uncorrelated with \(Z\)
- The variance of the prediction errors \(\leq\) the variance of \(Y\)
- \(\Cov{Y,Z} = \Cov{\LinearExpect{Y|Z}, Z}\)
The prediction errors average out to zero
\[\begin{eqnarray}
\Expect{Y-\LinearExpect{Y|Z}} & = & \Expect{Y - (\Expect{Y} + \beta(Z-\Expect{Z}))}\\
& = & \Expect{Y} - \Expect{Y} - \beta(\Expect{Z} - \Expect{Z}) = 0
\end{eqnarray}\]
- If they didn’t average to zero, we’d adjust the coefficients until they did
- Important: In general, \(\Expect{Y-\LinearExpect{Y|Z}|Z} \neq 0\)
How big are the prediction errors?
\[\begin{eqnarray}
\Var{Y-\LinearExpect{Y|Z}} & = & \Var{Y - \alpha - \beta Z}\\
& = & \Var{Y - \beta Z}\\
& = & \Var{Y} + \beta^2\Var{Z} - 2\beta\Cov{Y,Z}
\end{eqnarray}\]
but \(\beta = \Cov{Y,Z}/\Var{Z}\) so
\[\begin{eqnarray}
\Var{Y-\LinearExpect{Y|Z}} & = & \Var{Y} + \frac{(\Cov{Y,Z})^2}{\Var{Z}} - 2\frac{(\Cov{Y,Z})^2}{\Var{Z}}\\
& = & \Var{Y} - \frac{(\Cov{Y,Z})^2}{\Var{Z}}\\
& < & \Var{Y} ~ \text{unless}\ \Cov{Y,Z} = 0
\end{eqnarray}\]
\(\Rightarrow\) Optimal linear predictor is (almost) always better than nothing…
The predictions reproduce the covariances
\[\begin{eqnarray}
\Cov{Z, \LinearExpect{Y|Z}} & = & \Cov{Z, \alpha + \beta Z}\\
& = & \Cov{Z, \beta Z}\\
& = & \beta \Var{Z}\\
& = & \frac{\Cov{Y,Z}}{\Var{Z}} \Var{Z} = \Cov{Y,Z}
\end{eqnarray}\]
Multivariate case
We try to predict \(Y\) from a whole bunch of variables
Bundle those predictor variables into \(\vec{Z}\)
Solution:
\[
\LinearExpect{Y|Z} = \alpha+\vec{\beta}\cdot \vec{Z} = \Expect{Y} + (\Var{\vec{Z}})^{-1} \Cov{\vec{Z},Y} \cdot (\vec{Z} - \Expect{\vec{Z}})
\]
and
\[
\Expect{(Y - \LinearExpect{Y|Z})^2} = \Var{Y} - \Cov{Y,\vec{Z}}^T (\Var{\vec{Z}})^{-1} \Cov{Y,\vec{Z}}
\]
(Gory algebraic details in backup slides)
The PCA view
The factor of \((\Var{\vec{Z}})^{-1}\) rotates and scales \(\vec{Z}\) to uncorrelated, unit-variance variables
- Start with the eigendecomposition of \(\Var{\vec{Z}}\): \[\begin{eqnarray}
\Var{\vec{Z}} & = & \mathbf{w} \mathbf{\Lambda} \mathbf{w}^T\\
(\Var{\vec{Z}})^{-1} & = & \mathbf{w} \mathbf{\Lambda}^{-1} \mathbf{w}^T\\
(\Var{\vec{Z}})^{-1} & = & (\mathbf{w} \mathbf{\Lambda}^{-1/2}) (\mathbf{w} \mathbf{\Lambda}^{-1/2})^T\\
& = & (\Var{\vec{Z}})^{-1/2} \left((\Var{\vec{Z}})^{-1/2}\right)^T\\
\end{eqnarray}\]
- Use this to motivate defining new, uncorrelated, unit-variance variables: \[\begin{eqnarray}
\vec{U} & \equiv & \vec{Z} \Var{\vec{Z}}^{-1/2}\\
\Var{\vec{U}} & = & \mathbf{I}\\
\end{eqnarray}\]
Now replace \(\vec{Z}\) in the linear predictor with \(\vec{U}\) and see how simple it is: \[\begin{eqnarray}
\vec{Z}\cdot\vec{\beta} & = & \vec{Z} \cdot (\Var{\vec{Z}})^{-1} \Cov{\vec{Z}, Y}\\
& = & \vec{Z} \Var{\vec{Z}}^{-1/2} \cdot \left(\Var{\vec{Z}}^{-1/2}\right)^T \Cov{\vec{Z}, Y}\\
& = & \vec{U} \cdot \Cov{\vec{U}, Y}\\
\end{eqnarray}\]
What we don’t assume, again
- That anything is Gaussian
- Anything else about the distributions of \(Y\) or \(\vec{Z}\)
- That the linear predictor is correct
Some possible contexts
- Interpolating or extrapolating one variable over space and/or time
- Predicting one variable from another
- Predicting one variable from 2+ others
Prediction for \(X_0\) is a linear combination of \(X_1, \ldots X_n\) \[\begin{eqnarray}
\LinearExpect{X_0|X_1, X_2, \ldots X_n} & = & \alpha + \vec{\beta} \cdot \left[\begin{array}{c} X_1 \\ X_2 \\ \vdots \\X_n \end{array}\right]\\
\alpha & = & \Expect{X_0} - \vec{\beta} \cdot \left[\begin{array}{c} \Expect{X_1}\\ \Expect{X_2} \\ \vdots \\ \Expect{X_n}\end{array}\right] ~ \text{(goes away if everything's centered)}\\
\vec{\beta} & = & {\left[\begin{array}{cccc} \Var{X_1} & \Cov{X_1,X_2} & \ldots & \Cov{X_1,X_n}\\
\Cov{X_1,X_2} & \Var{X_2} & \ldots & \Cov{X_2,X_n}\\
\vdots & \vdots & \ldots & \vdots\\
\Cov{X_1,X_n} & \Cov{X_2,X_n} & \ldots & \Var{X_n}\end{array}\right]}^{-1} \left[\begin{array}{c} \Cov{X_0, X_1}\\
\Cov{X_0,X_2}\\ \vdots \\ \Cov{X_0,X_n}\end{array}\right]
\end{eqnarray}\]
- looks a lot like a linear smoother
- specifically \(\vec{\beta}\) looks like one row of \(\mathbf{w}\)
- best choice of weights \(\mathbf{w}\) comes from variances and covariances
A little worked example (II)
- We see \(X(r_1, t_1)\) and \(X(r_2, t_2)\), for short \(X_1\) and \(X_2\)
- We want to predict \(X(r_0, t_0)\), for short \(X_0\)
- Assume: \(\Expect{X(r, t)} = \mu\) for all \(r,t\)
- Assume: \(\Var{X(r,t)} = \sigma^2\) ditto
- Assume: \(\Cov{X_1, X_0} = \Cov{X_2, X_0} = \gamma\)
- Assume: \(\Cov{X_1, X_2} = \rho\)
Work out \(\vec{\beta}\) (off-line!) and get \[\begin{eqnarray}
\LinearExpect{X_0|X_1, X_@} = \mu + \frac{\gamma}{\sigma^2 + \rho}\left( (x_1 - \mu) + (x_2 - \mu)\right)
\end{eqnarray}\] vs. with one predictor \[
\mu + \frac{\gamma}{\sigma^2}(x_1 - \mu)
\]
Predicting one variable from another
- Given: values of variable \(U\) at many points, \(U(r_1, t_1), \ldots U(r_n, t_n)\)
- Desired: estimate of \(X\) at point \((r_0, t_0)\), \(X\neq U\)
\[\begin{eqnarray}
Y & = & X(r_0, t_0)\\
\vec{Z} & = & [U(r_1, t_1), U(r_2, t_2), \ldots U(r_n, t_n)]\\
\end{eqnarray}\]
- Need to find covariances of the \(U\)s with each other, and their covariances with \(X\)
Predicting one variable from 2+ others
- Given: values of two variables \(U\), \(V\) at many points
- Desired: estimate of \(X\) at one point
\[\begin{eqnarray}
Y & = & X(r_0, t_0)\\
\vec{Z} & = & [U(r_1, t_1), V(r_1, t_1), U(r_2, t_2), V(r_2, t_2), \ldots U(r_n, t_n), V(r_n, t_n)]
\end{eqnarray}\]
- Need to find covariances of \(U\)s and \(V\)s with each other, and with \(X\)
Covariance functions, kernels, etc.
- Karhunen-Loeve: expand the stochastic process in basis functions: \[
X(t) = \sum_{i=1}^{\infty}{S_i \phi_i(t)}
\] with \[
\int{\phi_i(t) \phi_j(t) dt} = \delta_{ij}
\] and \[
S_i = \int{X(t) \phi_i(t) dt}
\]
- Good (minimal MSE) basis functions solve an integral-equation eigenproblem: \[
\int{\Cov{X(t), X(s)} \phi_i(s) ds} = \lambda_i \phi_i(t)
\] since then \(\Cov{S_i, S_j} = \lambda_i \delta_{ij}\) and \[
\Cov{X(t), X(s)} = \sum_{i}{\lambda_i \phi_i(t) \phi_i(s)}
\]
- Abbreviate: \(\CovKernel(t, s) \equiv \Cov{X(t), X(s)}\)
- \(\CovKernel\) is the kernel of the integral operator
The function space implicit in the covariance
- Define \(\mathcal{K}\) as functions of the form \[
v(t) = \sum_{i=1}^{\infty}{v_i \phi_i(t)}
\]
- The realizations of our stochastic process are \(\in \mathcal{K}\)
- This is a Hilbert space, with inner product \[
\RKHSInner{u}{v} = \sum_{i=1}^{\infty}{\frac{1}{\lambda_i} u_i v_i}
\]
- Remember \(v_i = \int{v(t) \phi_i(t) dt}\)
- This makes the space into a reproducing-kernel Hilbert space (RKHS) \[
\RKHSInner{v}{\CovKernel(\cdot, t)} = \sum_{i=1}^{\infty}{\frac{1}{\lambda_i} v_i \lambda_i \phi_i(t)} = \sum_{i=1}^{\infty}{v_i \phi_i(t)} = v(t)
\]
- \(\CovKernel(\cdot, t)\) acts like the Dirac \(\delta(t)\), but is an actual function
After Parzen (1963); Parzen (1961)
\[
Z(t) = \sum_{i=1}^{\infty}{S_i \phi_i(t)}
\] with \(\Cov{S_i, S_j} = \lambda_i \delta_{ij}\)
\[\begin{eqnarray}
\Cov{\RKHSInner{h}{Z}, \RKHSInner{g}{Z}} & = & \Cov{\sum_{i=1}^{\infty}{\frac{1}{\lambda_i}{h_i S_i}}, \sum_{j=1}^{\infty}{\frac{1}{\lambda_j}{g_j S_j}}}\\
& = & \sum_{i, j}{\frac{1}{\lambda_i \lambda_j} h_i g_j \Cov{S_i, S_j}}\\
& = & \sum_{i, j}{\frac{1}{\lambda_i \lambda_j} h_i g_j \lambda_i \delta_{ij}}\\
& = & \sum_{i=1}^{\infty}{\frac{1}{\lambda_i} h_i g_i}\\
& = & \RKHSInner{h}{g}
\end{eqnarray}\]
Similarly \[\begin{eqnarray}
\Cov{Y, \RKHSInner{h}{Z}} & = & \sum_{i=1}^{\infty}{\frac{1}{\lambda_i} h_i \Cov{Y, S_i}}\\
& = & \RKHSInner{h}{\Cov{Y,Z}}
\end{eqnarray}\]
After Parzen (cont’d).
\[\begin{eqnarray}
\Expect{(Y - (a + \RKHSInner{b}{Z}))^2} & = & ((\Expect{Y} - \Expect{\RKHSInner{b}{Z}}) - a)^2 + \Var{Y} + \Var{\RKHSInner{b}{Z}} - 2\Cov{Y, \RKHSInner{b}{Z}}\\
& = & ((\Expect{Y} - \Expect{\RKHSInner{b}{Z}}) - a)^2 + \Var{Y} + \RKHSInner{b}{b} - 2 \RKHSInner{b}{\Cov{Y,Z}}\\
& = & ((\Expect{Y} - \Expect{\RKHSInner{b}{Z}}) - a)^2 + \Var{Y} + \RKHSInner{b- \Cov{Y, Z}}{b-\Cov{Y,Z}} - \RKHSInner{\Cov{Y,Z}}{\Cov{Y,Z}}
\end{eqnarray}\]
- \(a\) only appears in the square \(\Rightarrow\) send square to zero
- \(b\) only shows up in the norm (3rd term), minimized at \(b = \Cov{Y,Z}\)
The optimal linear predictor is therefore \[
\LinearExpect{Y|Z} = (\Expect{Y} - \RKHSInner{\Cov{Y,Z}}{\Expect{Z}}) + \RKHSInner{\Cov{Y,Z}}{Z}
\] and the minimized expected MSE is \[
\Var{Y} - \RKHSInner{\Cov{Y,Z}}{\Cov{Y,Z}}
\]
- Everything we did before is a special case
- Q: Wait, where did the inverse variance matrix go?
- A: It’s implicit in the inner product, because \(\RKHSInner{u}{v} = \sum_{i=1}^{n}{\lambda_i^{-1} u_i v_i}\)
Summing up on linear prediction in general
- We can always decide to use a linear predictor, \(\LinearExpect{Y|\vec{Z}} = \alpha + \vec{\beta} \cdot \vec{Z}\)
- The optimal linear predictor of \(Y\) from \(\vec{Z}\) always takes the same form: \[
\LinearExpect{Y|\vec{Z}} = \Expect{Y} + (\Var{\vec{Z}})^{-1} \Cov{Y,\vec{Z}} \cdot (\vec{Z} - \Expect{\vec{Z}})
\]
- Doing linear prediction requires finding the covariances
Optimal linear prediction for time series
- Given: \(X(t_1), X(t_2), \ldots X(t_n)\)
- Not necessarily equally spaced in time
- Desired: prediction of \(X(t_0)\)
\[\begin{eqnarray}
\LinearExpect{X(t_0)|X(t_1), \ldots X(t_n)} & = & \alpha + \vec{\beta} \cdot \left[\begin{array}{c} X(t_1) \\ X(t_2) \\ \vdots \\ X(t_n) \end{array}\right]\\
\alpha & = & \Expect{X(t_0)} - \vec{\beta} \cdot \left[\begin{array}{c} \Expect{X(t_1)}\\ \Expect{X(t_2)} \\ \vdots \\ \Expect{X(t_n)}\end{array}\right] ~ \text{(goes away if everything's centered)}\\
\vec{\beta} & = & {\left[\begin{array}{cccc} \Var{X(t_1)} & \Cov{X(t_1), X(t_2)} & \ldots & \Cov{X(t_1), X(t_n)}\\
\Cov{X(t_1), X(t_2)} & \Var{X(t_2)} & \ldots & \Cov{X(t_2), X(t_n)}\\
\vdots & \vdots & \ldots & \vdots\\
\Cov{X(t_1), X(t_n)} & \Cov{X(t_2), X(t_n)} & \ldots & \Var{X(t_n)}\end{array}\right]}^{-1} \left[\begin{array}{c} \Cov{X(t_0), X(t_1)}\\
\Cov{X(t_0), X(t_2)}\\ \vdots \\ \Cov{X(t_0), X(t_n)}\end{array}\right]
\end{eqnarray}\]
- What is this good for?
- Interpolation
- Extrapolation / prediction
Interpolation
- Time series often have gaps
- Instruments fail, people mess up, circumstances…
- What happened at the times in between the observations?
Back to Kyoto
- A lot of what we see in this plot is just made up
What we didn’t tell R to make up
When did the cherries flower in 1015?
- We need \(\Cov{X(1015), X(t_i)}\) for every year \(t_i\) where we have data
- We need \(\Expect{X(t_i)}\), \(\Var{X(t_i)}\) and \(\Cov{X(t_i), X(t_j)}\) ditto
- We need \(\Expect{X(1015)}\)
Getting the expectations and covariances
- We only see each \(X(t_i)\) once
- Maybe gives us an idea of \(\Expect{X(t_i)}\)
- but not \(\Var{X(t_i)}\) or \(\Cov{X(t_i), X(t_j)}\)
- let alone \(\Cov{X(t_0), X(t_i)}\)
- We could repeat the experiment many times
- We could make assumptions
Repeating the experiment
\[\begin{eqnarray}
\overline{x}(t) \equiv \frac{1}{n}\sum_{i=1}^{n}{x^{(i)}(t)} & \rightarrow & \Expect{X(t)}\\
\frac{1}{n}\sum_{i=1}^{n}{(x^{(i)}(t) - \overline{x}(t)) (x^{(i)}(s) - \overline{x}(s))} & \rightarrow & \Cov{X(t), X(s)}
\end{eqnarray}\]
Making assumptions
- Assume some covariances (and expectations) are equal
- Weaker: assume some of them are similar
Stationarity (assuming covariances are equal)
- A time series is weakly stationary when \(\Expect{X(t)} = \Expect{X(s)}\) (constant in time) and \(\Cov{X(t), X(s)} = \AutoCov(|t-s|)\)
- \(\AutoCov(h)\) is the autocovariance function at lag \(h\)
- this is also called “second-order” or “wide-sense” stationarity
- We’ll cover strong / full / strict stationarity later
- Weak stationarity implies all the expectations are the same
- Weak stationarity implies \(\Cov{X(t), X(t+h)} = \Cov{X(s), X(s+h)}\) so lots of covariances are the same
- Weak stationarity lets us pool information: \[\begin{eqnarray}
\overline{x} \equiv \frac{1}{n}\sum_{i=1}^{n}{X(t_i)} & \rightarrow & \Expect{X(0)}\\
\frac{1}{n}\sum_{i=1}^{n}{(x(t_i) - \overline{x})(x(t_i + h) - \overline{x})} & \rightarrow & \AutoCov(h)
\end{eqnarray}\]
The autocovariance function
- \(\AutoCov(0) = \Var{X(t)}\) (constant in \(t\))
- \(\rho(h) \equiv \frac{\AutoCov(h)}{\AutoCov(0)} =\) autocorrelation function
Assuming stationarity…
\[\begin{eqnarray}
\mathbf{v} & \equiv & \left[\begin{array}{ccc} \AutoCov(0) & \AutoCov(|t_1 - t_2|) & \ldots & \AutoCov(|t_1 - t_n|)\\
\AutoCov(|t_2-t_1|) & \AutoCov(0) & \ldots & \AutoCov(|t_2-t_n|)\\
\vdots & \vdots & \ldots & \vdots\\
\AutoCov(|t_n-t_1|) & \AutoCov(|t_n-t_2|) & \ldots & \AutoCov(0)\end{array}\right]\\
\mathbf{c} & = & \left[\begin{array}{c} \AutoCov(|t_0 - t_1|)\\ \AutoCov(|t_0-t_2|) \\ \vdots \\ \AutoCov(|t_0 -t_n|)\end{array} \right]\\
\vec{\beta} & = & \mathbf{v}^{-1}\mathbf{c}\\
\LinearExpect{X(t_0)|X(t_1), \ldots X(t_n)} & = & \overline{x} + \vec{\beta} \cdot \left[\begin{array}{c} X(t_1) -\overline{x}\\
X(t_2) - \overline{x} \\ \vdots \\ X(t_n) -\overline{x}\end{array}\right] = \overline{x} + \mathbf{c}^{T}\mathbf{v}^{-1} \left[\begin{array}{c} X(t_1) -\overline{x}\\
X(t_2) - \overline{x} \\ \vdots \\ X(t_n) -\overline{x}\end{array}\right]\\
\Expect{(X(t_0) - \LinearExpect{X(t_0)|X(t_1), \ldots X(t_n)})^2} & = & \AutoCov(0) - \mathbf{c}^T\mathbf{v}^{-1} \mathbf{c}
\end{eqnarray}\]
In R
acf(x, lag.max, type, plot, na.action, ...)
x = a time series (e.g., a vector) or a data frame of multiple time serieslag.max = maximum value of lag \(h\)
- Calculates \(\AutoCov(0), \AutoCov(1), \ldots \AutoCov(h)\) and stops there
type = correlation (default) or covariance?plot = make a plot? (default)na.action = how to handle NAs? default is give up, na.pass will use complete pairs
In R
kyoto.acf <- acf(kyoto$Flowering.DOY, lag.max=100, type="covariance", na.action=na.pass)
In R
kyoto.acf[0:5] # Why 0?
##
## Autocovariances of series 'kyoto$Flowering.DOY', by lag
##
## 0 1 2 3 4 5
## 42.93 10.43 10.02 10.31 8.17 8.35
kyoto.acf[0:5]$acf
## , , 1
##
## [,1]
## [1,] 42.932054
## [2,] 10.430414
## [3,] 10.023484
## [4,] 10.314843
## [5,] 8.172978
## [6,] 8.350650
When was the flowering in 1015?
# Find years within +- 49 of 1015
# Why +- 49 when we have covariances out to lag 100?
available.years <- with(na.omit(kyoto), Year.AD[Year.AD > 1015-49 & Year.AD < 1015+49])
historical.mean <- mean(kyoto$Flowering.DOY, na.rm=TRUE)
CovYZ <- matrix(kyoto.acf[abs(available.years-1015)]$acf, ncol=1)
year.lags <- outer(available.years, available.years, "-")
year.lags <- abs(year.lags)
VarZ <- kyoto.acf[year.lags]$acf
VarZ <- matrix(VarZ, ncol=length(available.years), byrow=FALSE)
Z <- with(kyoto, Flowering.DOY[Year.AD %in% available.years])
(fitted.value <- historical.mean + (Z-historical.mean) %*% solve(VarZ) %*% CovYZ)
## [,1]
## [1,] 106.3017
(fitted.value.se <- sqrt(as.matrix(kyoto.acf[0]$acf) - t(CovYZ) %*% solve(VarZ) %*% CovYZ))
## [,1]
## [1,] 5.862828
When was the flowering in 1015?
This is a lot of work…
…and we’d need to re-do most of it for every other year
… so we should write a function (comments in .Rmd)
kyoto.prediction <- function(times) {
historical.mean <- mean(kyoto$Flowering.DOY, na.rm = TRUE)
historical.variance <- as.matrix(kyoto.acf[0]$acf)
fits <- matrix(0, nrow = length(times), ncol = 3)
colnames(fits) <- c("time", "fit", "se")
fits[, "time"] <- times
rownames(fits) <- paste(times)
for (t in times) {
available.years <- with(na.omit(kyoto), Year.AD[Year.AD > t - 49 & Year.AD <
t + 49])
available.years <- setdiff(available.years, t)
CovYZ <- matrix(kyoto.acf[abs(available.years - t)]$acf, ncol = 1)
year.lags <- outer(available.years, available.years, "-")
year.lags <- abs(year.lags)
VarZ <- kyoto.acf[year.lags]$acf
VarZ <- matrix(VarZ, ncol = length(available.years), byrow = FALSE)
Z <- with(kyoto, Flowering.DOY[Year.AD %in% available.years])
fits[paste(t), "fit"] <- historical.mean + (Z - historical.mean) %*% solve(VarZ) %*%
CovYZ
fits[paste(t), "se"] <- sqrt(max(0, historical.variance - t(CovYZ) %*% solve(VarZ) %*%
CovYZ))
}
return(data.frame(fits))
}
Finally, some interpolation: close up
Finally, some interpolation: zoom out
Stationarity and trend-plus-fluctuation
- We can always say \(X(t) = \TrueRegFunc(t) + \TrueNoise(t)\), with \(\Expect{X(t)} =\TrueRegFunc(t)\) and \(\Expect{\TrueNoise(t)} = 0\)
- If \(X\) is (weakly) stationary, then
- \(\TrueRegFunc(t) = \mu(0)\), a constant
- \(\Cov{\TrueNoise(t), \TrueNoise(s)} = \Cov{X(t), X(s)} = \AutoCov(|t-s|)\)
- Assuming a flat trend is pretty strong…
Stationary fluctuations around a trend
- If \(\TrueNoise\) is weakly stationary, then \(X - \TrueRegFunc\) is weakly stationary
- so \(X - \EstRegFunc\) should be (approximately) weakly stationary
- Recipe:
- Use smoothing to estimate \(\Expect{X(t)}\) as \(\EstRegFunc(t)\)
- Find ACF from \(\EstNoise(t) = X(t) - \EstRegFunc(t)\)
- Warning: Yule-Slutsky says there’s bound to be some correlations in \(\EstNoise\)…
- Now get the coefficients and predict
Detrending the cherry blossoms with a spline
# Get the spline
kyoto.spline <- with(na.omit(kyoto), smooth.spline(x=Year.AD, y=Flowering.DOY))
# Calculate residuals but pad them out with NAs for easier plotting
# Start with a vector of all NAs
residuals.with.NAs <- rep(NA, times=nrow(kyoto))
# Replace the entries where we can calculate a sensible residual
residuals.with.NAs[!is.na(kyoto$Flowering.DOY)] <- na.omit(kyoto$Flowering.DOY) - kyoto.spline$y
# Add it to the data frame as a new column
kyoto$fluctuation <- residuals.with.NAs
Detrending the cherry blossoms with a spline
Detrending the cherry blossoms with a spline
Detrending the cherry blossoms with a spline
Covariance estimation
acf uses the sample covariance- Estimates covariance at each lag separately
- Sometimes inefficient
- Covariances might have known form, e.g., \(\AutoCov(h) = \AutoCov(0)e^{-h/\tau}\)
- We’ll look more at that next time, in the spatial context
- Even if we don’t believe something like that, we might believe that \(\AutoCov(h)\) should be close to \(\AutoCov(h-1)\) and \(\AutoCov(h+1)\) \(\Rightarrow\) use smoothing (again)
- Estimating the power spectrum, by smoothing the periodogram and then inverse-Fourier-transforming, imposes smoothness on the estimated ACF
Checking stationarity
- Divide the data into intervals
- Re-estimate ACF on each interval
- ACF shouldn’t change much
Checking stationarity: raw data
Checking stationarity: after spline detrending
Checking stationarity
- Divide the data into intervals
- Re-estimate ACF on each interval
- Shouldn’t change much
- How much difference is too much? (Will answer later with simulation)
- Some inherent pitfalls:
- Very slight non-stationarity looks pretty stationary
- Stationary but with \(\AutoCov(h) \rightarrow 0\) very slowly looks pretty non-stationary
- “Long-range correlations”, “slowly-decaying correlations”
- Trends show up as slowly-decaying correlations
Multiple time series
- \(X\) and \(U\) are jointly weakly stationary when
- \(X\) is weakly stationary, with ACF \(\AutoCov_X\)
- \(U\) is weakly stationary, with ACF \(\AutoCov_U\)
- \(\Cov{X(t), U(s)} = \AutoCov_{UX}(|t-s|)\)
ccf in R will calculate the cross-covariance (or cross-correlation) function
- So will
acf if it’s given a matrix or data frame
- Plug and chug: if predicting \(X\) from \(U\)
- \(\AutoCov_U\) goes into the variance matrix \(\mathbf{v}\)
- \(\AutoCov_{UX}\) goes into the covariance matrix \(\mathbf{c}\)
Where did all this come from?
- Norbert Wiener (1894–1964)
- One of the great mathematicians of the 20th century
- If you think you’ve got issues with parents pushing you to succeed in school, read Wiener (1953)
- Worked very closely with engineers
- His ideas are why we talk about “feedback”, “information”, “cyber-”, etc. (Wiener 1961, 1954)
- (Photo from Kline (2017))
Wiener’s research in 1942
- Anti-aircraft fire control
- Aim at where the target aircraft will be
- The target moves erratically
- \(\Rightarrow\) Need to extrapolate a random process
- Wiener’s solution: basically what we’ve just done
- Plus the continuous-time version
- Plus implementation using 1940s electronics
- Declassified after the war as Wiener (1949)
- Parallel work by Kolmogorov (1941)
- I’ve never read any explanation for what Kolmogorov was up to
- Wiener and Kolmogorov’s two big ideas:
- Optimal linear prediction didn’t need the usual regression assumptions (Gaussian noise, independent variables, etc.)
- For continuously-observed data, you can do everything in the Fourier domain
Summing up
- Optimal linear prediction needs to know the trend and the autocovariance function
- Estimating them from one time series needs assumptions
- Remove trend
- Find ACF from residuals / estimated fluctuations
- Once we have the ACF, finding coefficients is just linear algebra
Backup: Gory details for multivariate predictors
\[\begin{eqnarray}
m(\vec{Z}) & = & a + \vec{b} \cdot \vec{Z}\\
(\alpha, \vec{\beta}) & = & \argmin_{a \in \mathbb{R}, \vec{b} \in \mathbb{R}^n}{\Expect{(Y-(a + \vec{b} \cdot \vec{Z}))^2}}\\
\Expect{(Y-(a+\vec{b} \cdot \vec{Z}))^2} & = & \left(\Expect{Y} - a - \vec{b} \cdot \Expect{\vec{Z}}\right)^2 + \Var{Y - \vec{b} \cdot \vec{Z}}\\
& = & \left(\left(\Expect{Y} - \vec{b} \cdot\Expect{\vec{Z}}\right) - a\right)^2\\
& & + \Var{Y} + \vec{b} \cdot \Var{\vec{Z}} \vec{b} - 2\vec{b} \cdot \Cov{Y, \vec{Z}}
\end{eqnarray}\]
\(\Var{\vec{Z}}\) is a square matrix, \(\Cov{Y, \vec{Z}}\) is a vector
Backup: Gory details: the intercept
\(a\) only shows up in the square, which we can always zero out by setting \[
a = \Expect{Y} - \vec{b} \cdot \vec{Z}
\]
More formally, taake derivative w.r.t. \(a\), set to 0 at \(a=\alpha\), \(\vec{b}=\vec{\beta}\):
\[\begin{eqnarray}
0 & = & -2\Expect{Y} + 2\vec{\beta} \cdot \Expect{\vec{Z}} + 2\alpha \\
\alpha & = & \Expect{Y} - \vec{\beta} \cdot \Expect{\vec{Z}}\\
\end{eqnarray}\]
just like when \(Z\) was univariate
Backup: Gory details: the slopes
Take derivative and set to zero at \(\vec{b}=\vec{\beta}\):
\[\begin{eqnarray}
2 \Var{\vec{Z}} \vec{\beta} - 2 \Cov{Y, \vec{Z}} & = & 0\\
\Var{\vec{Z}}\vec{\beta} & = & \Cov{Y, \vec{Z}}\\
\vec{\beta} & = & \Var{\vec{Z}}^{-1} \Cov{Y, \vec{Z}}
\end{eqnarray}\]
Reduces to \(\Cov{Y,Z}/\Var{Z}\) when \(Z\) is univariate
Backup: Estimation I: “plug-in”
- We don’t see the true expectations, variances, covariances
- But we can have sample/empirical values
- One estimate of the optimal linear predictor: plug in the sample values
so for univariate \(Z\), \[
\hat{m}(z) = \overline{y} - \frac{\widehat{\Cov{Y,Z}}}{\widehat{\Var{Z}}}(z-\overline{z})
\]
Backup: Estimation II: ordinary least squares
- We don’t see the true expected squared error, but we do have the sample mean squared error
- Minimize that
- Leads to exactly the same results as plug-in approach!
Backup: Estimation: When does OLS/plug-in work?
- Jointly sufficient conditions:
- Sample means converge on expectation values
- Sample covariances converge on true covariance
- Sample variances converge on true, invertible variance
- Then by continuity OLS coefficients converge on true \(\beta\)
- This can all happen even when everything is dependent on everything else!
Backup: Square roots of a matrix
- A square matrix \(\mathbf{d}\) is a square root of \(\mathbf{c}\) when \(\mathbf{c} = \mathbf{d} \mathbf{d}^T\)
- If there are any square roots, there are many square roots
- Pick any orthogonal matrix \(\mathbf{o}^T = \mathbf{o}^{-1}\)
- \((\mathbf{d}\mathbf{o})(\mathbf{d}\mathbf{o})^T = \mathbf{d}\mathbf{d}^T\)
- Just like every real number has two square roots…
- If \(\mathbf{c}\) is diagonal, define \(\mathbf{c}^{1/2}\) as the diagonal matrix of square roots
- If \(\mathbf{c} = \mathbf{w}\mathbf{\Lambda}\mathbf{w}^T\), one square root is \(\mathbf{w}\mathbf{\Lambda}^{1/2}\)
Backup: Not inverting the variance matrix
- \(\vec{\beta} = \Var{\vec{Z}}^{-1} \Cov{\vec{Z}, Y}\)
- But inverting an \(n\times n\) matrix takes \(O(n^3)\) operations
- Then matrix multiplying \(\Var{\vec{Z}}^{-1}\) by \(\Cov{\vec{Z}, Y}\) is \(O(n^2)\), and multiplying \(\vec{\beta}\) by \(Z\) is \(O(n)\), and we need to do everything \(n\) times so over-all time complexity of getting \(n\) fits is \(n(O(n^3) + O(n^2) + O(n)) = O(n^4)\)
- Strictly speaking we don’t need to invert \(\Var{\vec{Z}}\); we just need to find the \(\vec{\beta}\) which solves the equation \(\Var{\vec{Z}} \vec{\beta} = \Cov{\vec{Z}, Y}\)
- Solving a linear system of equations is also \(O(n^3)\), but with better constants than doing the full inversion
- Details of how to solve without inverting are left to numerical linear algebra texts
- Programming exercise: re-do my code for the optimal linear predictor so it doesn’t invert the variance matrix but does
solve() for the optimal coefficients
Backup: How good is the optimal linear predictor?
all.fits$std.residuals <- (kyoto$Flowering.DOY - all.fits$fit)/all.fits$se
mean(all.fits$std.residuals, na.rm=TRUE)
## [1] -0.007427299
sd(all.fits$std.residuals, na.rm=TRUE)
## [1] 1.012488
Ideally: mean 0, standard deviation 1
Backup: How good is the optimal linear predictor?
This is a pretty good random scatter of points
Backup: fancier covariance estimation
- Two (equivalent) definitions of the covariance: \[
\Cov{X(t), X(s)} = \Expect{X(t) X(s)} - \Expect{X(t)}\Expect{X(s)} = \Expect{(X(t) - \Expect{X(t)}) (X(s) - \Expect{X(s)})}
\]
- Use the 2nd form: define \[
\Gamma(t,s) \equiv (X(t) - \Expect{X(t)})(X(s) - \Expect{X(s)})
\] so \(\Cov{X(t), X(s)} = \Expect{\Gamma(t,s)}\)
- Now assume constant mean so we can estimate \[
\widehat{\Gamma}(t,s) = (X(t) - \overline{x}) (X(s) - \overline{x})
\]
- Now assume stationarity so \(\Cov{X(t), X(s)} = \Expect{\Gamma(t,s)} = \AutoCov(|t-s|)\)
- Finally assume \(\AutoCov(h)\) changes slowly in \(h\)
- We can estimate \(\AutoCov(h)\) by finding pairs \(t, s\) with \(|t-s| \approx h\) and averaging \(\widehat{\Gamma}(t,s)\)
- Or kernel smoothing or spline smoothing or…
Backup: stationary linear predictor vs. spline
Backup: A little bit about how Wiener approached the problem
- \(X(t)\) observed continuously on some domain \(\mathbb{D}\), we desire a prediction for \(X(t_0)\), generally \(t_0 \not in \mathbb{D}\)
- The prediction should be a linear operation on \(X(t)\), say \(\beta\)
- This is another continuous function
- Minimizing expected squared error as usual leads to \[
\int_{\mathbb{D}}{\beta(s) \CovKernel(s, t) ds} = \Cov{X(t_0), X(t)}
\]
- A.k.a. a Wiener-Hopf integral equation
- Assuming stationarity, \[
\int_{\mathbb{D}}{\beta(s) \AutoCov{s-t} ds} = \AutoCov{t_0 - t}
\]
- If \(\mathbb{D}= (-\infty, \infty)\) then this is a convolution, so the solution would be easy in the frequency domain
- In practice \(\mathbb{D}\) is at best semi-infinite and finding explicit solutions by frequency-domain methods is very tricky
- Parzen drastically simplified the problem by the RKHS approach (Parzen 1963, 1961 )
- Pretty nearly the first use of RKHS in statistics
Backup: Hilbert space
A Hilbert space \(\mathcal{H}\) is a space which
- Is a vector space: \(u, v \in H\) and \(a, b \in \mathbb{R}\) \(\Rightarrow\) \(au+bv \in H\), \(\exists 0 \in H\) s.t. \(0+v = v\), etc.
- Or \(a, b \in \mathbb{C}\) or…
- Has an inner product \(\HilbertInner{\cdot}{\cdot}{\mathcal{H}}\), and a norm \(\|v\|_{\mathcal{H}} = \sqrt{\HilbertInner{v}{v}{\mathcal{H}}}\)
- Inner products are symmetric and bilinear, 0$, \(\HilbertInner{v}{v}{\mathcal{H}} = 0\) iff \(v=0\)
- Contains the limit of all Cauchy-convergent sequences.
- Is separable: you can name any \(\epsilon > 0\), and I can find a finite set \(\mathcal{S}_n\) of points in \(\mathcal{H}\) where radius-\(\epsilon\) balls around the \(\mathcal{S}_n\) points contain the whole of \(\mathcal{H}\)
Most (interesting) Hilbert spaces are spaces of functions on some domain \(\mathbb{D}\). A reproducing kernel for \(\mathcal{H}\) is a function \(K: \mathbb{D} \times \mathbb{D} \mapsto \mathbb{R}\) with the following properties:
- For each \(t \in \mathbb{D}\), \(K(\cdot, t) \in \mathcal{H}\)
- For each \(t\), \(\HilbertInner{v}{K(\cdot, t)}{\mathcal{H}} = v(t)\)
- Compare \(\int{v(s) \delta(s-t) ds} = v(t)\)
A reproducing kernel implies that \(K(s,t) = \sum_{i=1}^{\infty}{\lambda_i \phi_i(s) \phi_i(t)}\) for eigenfunctions \(\phi_i\) with eigenvalues \(\lambda_i\)
This means we can (usually) avoid having to explicitly work with the eigenfunctions, we can just evaluate the kernel
References
Cave-Browne-Cave, F. E. 1905. “On the Influence of the Time Factor on the Correlation Between the Barometric Heights at Stations More Than 1000 Miles Apart.” Proceedings of the Royal Society of London 74:403–13. https://www.jstor.org/stable/116708.
Klein, Judy L. 1997. Statistical Visions in Time: A History of Time Series Analysis, 1662–1938. Cambridge, England: Cambridge University Press.
Kline, Ronald R. 2017. The Cybernetics Moment: Or Why We Call Our Age the Information Age. Baltimore, Maryland: Johns Hopkins University Press.
———. 1963. “A New Approach to the Synthesis of Optimal Smoothing and Prediction Systems.” In Mathematical Optimization Techniques, edited by Richard Bellman, 75–108. Berkeley: University of California Press. https://doi.org/10.1525/9780520319875-007.
Wiener, Norbert. 1949. Extrapolation, Interpolation, and Smoothing of Stationary Time Series: With Engineering Applications. Cambridge, Massachusetts: The Technology Press of the Massachusetts Institute of Technology. https://doi.org/10.7551/mitpress/2946.001.0001.
———. 1953. Ex-Prodigy: My Childhood and Youth. New York: Simon; Schuster.
———. 1954. The Human Use of Human Beings: Cybernetics and Society. 2nd ed. Garden City, New York: Doubleday.
