--- title: Optimal Linear Prediction, Especially over Time author: 36-740 date: 8 September 2026 (Lecture 5) output: slidy_presentation: math_method: engine: mathjax url: https://bactra.org/mathjax/tex-svg.js bibliography: locusts.bib --- ```{r, include=FALSE} # General set-up options library(knitr) opts_chunk$set(size="small",background="white", highlight=FALSE, cache=TRUE, autodep=TRUE, warning=FALSE, message=FALSE, echo=FALSE) ``` \[ \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](https://mathshistory.st-andrews.ac.uk/Biographies/Cave-Browne-Cave_Evelyn/) for weather forecasting [@Cave-Browne-Cave-invents-spatio-temporal-prediction], see discussion in @Klein-statistical-visions, pp. 262--264) ## Some general properties of the optimal linear predictor 1. The prediction errors average out to zero 2. The prediction errors are uncorrelated with $Z$ 3. The variance of the prediction errors $\leq$ the variance of $Y$ 4. $\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$ ## The prediction errors are uncorrelated with $Z$ \begin{eqnarray} \Cov{Z, Y- \LinearExpect{Y|Z}} & = & \Expect{Z(Y-\LinearExpect{Y|Z})} ~\text{(by previous slide)}\\ & = & \Expect{Z(Y - \Expect{Y} - \frac{\Cov{Y,Z}}{\Var{Z}}(Z-\Expect{Z}))}\\ & = & \Expect{ZY - Z\Expect{Y} - \frac{\Cov{Y,Z}}{\Var{Z}}(Z^2) + \frac{\Cov{Y,Z}}{\Var{Z}} (Z \Expect{Z})}\\ & = & \Expect{ZY} - \Expect{Z}\Expect{Y} - \frac{\Cov{Y,Z}}{\Var{Z}}\Expect{Z^2} + \frac{\Cov{Y,Z}}{\Var{Z}} (\Expect{Z})^2\\ & = & \Cov{Z,Y} - \frac{\Cov{Y,Z}}{\Var{Z}}(\Var{Z})\\ & = & 0 \end{eqnarray} > - If they weren't uncorrelated, we'd adjust the coefficients until they were ## The prediction errors are uncorrelated with $Z$ Alternate take: \begin{eqnarray} \Cov{Z, Y-\LinearExpect{Y|Z}} & = & \Cov{Z, Y} - \Cov{Z, \alpha + \beta Z}\\ & = & \Cov{Y,Z} - \Cov{Z, \beta Z}\\ & = & \Cov{Y,Z} - \beta\Cov{Z,Z}\\ & = & \Cov{Y,Z} - \beta\Var{Z}\\ & = & \Cov{Y,Z} - \Cov{Y,Z} = 0 \end{eqnarray} ## 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 {.smaller} - 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 ## Interpolating or extrapolating a single variable - Given: $X(r_1, t_1), X(r_2, t_2), \ldots X(r_n, t_n)$, for short $X_1, \ldots X_n$ - Desired: estimate/guess at $X(r_0, t_0)$, for short $X_0$ \begin{eqnarray} Y & = & X_0\\ \vec{Z} & = & [X_1, X_2, \ldots X_n] \end{eqnarray} ## {.smaller} - 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) {.smaller} - 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. {.smaller} - 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 {.smaller} - _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-synthesis-of-linear-smoothers; @Parzen-time-series-analysis {.smaller} \[ 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). {.smaller} \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 {.smaller} - 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 ```{r} kyoto <- read.csv("http://www.stat.cmu.edu/~cshalizi/dst/26/data/kyoto-2026.csv") plot(Flowering.DOY ~ Year.AD, data=kyoto, type="o", pch=16, cex=0.5, ylab="Day in year of full flowering", xlab="Year (AD)", main="Cherry blossoms in Kyoto") rug(side=1, x=na.omit(kyoto)$Year.AD, ticksize=-0.02) ``` - A lot of what we see in this plot is just made up ## What we didn't tell R to make up ```{r} plot(Flowering.DOY ~ Year.AD, data=kyoto, type="p", pch=16, cex=0.3, ylab="Day in year of full flowering", xlab="Year (AD)", main="Cherry blossoms in Kyoto") rug(side=1, x=na.omit(kyoto)$Year.AD, ticksize=-0.02) ``` ## 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)}$ ## Similarly for extrapolation - What was $X(800)$? ("retrodiction") - What will $X(2027)$ be? ("prediction" or "forecast" in the strictest sense) ## 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 {.smaller} ```{r} # A simulation example wiener.bridge <- function(n, sigma2=1) { # Divide the time interval [0,1] into n equal chunks; on each chunk, # take an independent Gaussian step of variance $\sigma^2/n$; # then subtract the straight-line path from the origin to the end. W <- cumsum(rnorm(n, mean=0, sd=sqrt(sigma2/n))) B <- W - W[n]*(1:n)/n B <- c(0, B) # The process always starts at 0 return(B) } # Make an empty plot of the right size plot(0, type="n", xlim=c(0,1), ylim=c(-1,1), xlab="t", ylab=expression(X(t))) # Now fill it in with multiple runs of the process invisible(replicate(20, lines((0:100)/100, wiener.bridge(100), col="grey"))) ``` \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 ```{r, eval=FALSE, echo=TRUE} acf(x, lag.max, type, plot, na.action, ...) ``` > - `x` = a time series (e.g., a vector) _or_ a data frame of multiple time series > - `lag.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 ```{r, echo=TRUE} kyoto.acf <- acf(kyoto$Flowering.DOY, lag.max=100, type="covariance", na.action=na.pass) ``` ## In R ```{r, echo=TRUE} kyoto.acf[0:5] # Why 0? kyoto.acf[0:5]$acf ``` ## When was the flowering in 1015? {.smaller} ```{r, echo=TRUE} # 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) (fitted.value.se <- sqrt(as.matrix(kyoto.acf[0]$acf) - t(CovYZ) %*% solve(VarZ) %*% CovYZ)) ``` ## When was the flowering in 1015? {.smaller} ```{r} plot(available.years, solve(VarZ) %*% CovYZ, xlab="Year", ylab=expression(hat(beta)), main="Coefficient for estimating X(1015)") abline(h=0) abline(v=1015,lty="dashed") ``` ## This is a lot of work... {.smaller} ...and we'd need to re-do most of it for every other year ... so we should write a function (comments in .Rmd) ```{r, echo=TRUE, tidy=TRUE, tidy.opts=list(comment=FALSE)} # Assuming stationarity, calculate predictions from the Kyoto flowering # data # Uses observations up to 49 years on either side of the target time # Inputs: vector of years at which predictions are desired (times) # Outputs: 3-column matrix giving times, fitted values, and standard errors # Presumes: # kyoto is in memory (and is a data frame with the right column names) # kyoto.acf is in memory (and has covariances out to lag 100 years) # Every year in times has at least one data point within 49 years kyoto.prediction <- function(times) { # NOTE: Some bad programming practices here: presumes the existence of kyoto, # kyoto.acf, etc. the magic number 49 is written in to the code, etc. # Can you make those arguments to the function, so this # could be used on other, similar time series? # Find the global mean historical.mean <- mean(kyoto$Flowering.DOY, na.rm=TRUE) # Find the global variance historical.variance <- as.matrix(kyoto.acf[0]$acf) # Prepare a 3-column matrix for storing times, fits, and standard errors fits <- matrix(0, nrow=length(times), ncol=3) colnames(fits) <- c("time", "fit", "se") # Filling the first column is easy... fits[,"time"] <- times # ... but it'll also help to access rows by name rownames(fits) <- paste(times) # Do the calculation for each time for (t in times) { # What years have data to serve as predictors? available.years <- with(na.omit(kyoto), Year.AD[Year.AD > t-49 & Year.AD < t +49]) # Don't include the year itself # Though you could! available.years <- setdiff(available.years, t) # Covariances between the predictors and the predictand CovYZ <- matrix(kyoto.acf[abs(available.years-t)]$acf, ncol=1) # What are all the different lags between the predictors? year.lags <- outer(available.years, available.years, "-") year.lags <- abs(year.lags) # What are the implied covariances between the predictors? VarZ <- kyoto.acf[year.lags]$acf VarZ <- matrix(VarZ, ncol=length(available.years), byrow=FALSE) # What are the actual values of the predictors? Z <- with(kyoto, Flowering.DOY[Year.AD %in% available.years]) # What's the predicted value? fits[paste(t), "fit"] <- historical.mean + (Z-historical.mean) %*% solve(VarZ) %*% CovYZ # What's the standard error = sqrt of variance of error? # This should be 0 at every point which was in the data used for # estimation, but round-off error sometimes gives small negative # values inside the square root --- rectify those to 0 to avoid # annoying warning messages fits[paste(t), "se"] <- sqrt(max(0,historical.variance - t(CovYZ) %*% solve(VarZ) %*% CovYZ)) } return(data.frame(fits)) } ``` ```{r} # Calculate fits for all years # Avoid some magic numbers... first.year <- min(kyoto$Year.AD) last.year <- max(kyoto$Year.AD) all.fits <- kyoto.prediction(first.year:last.year) ``` ## Finally, some interpolation: close up ```{r} plot(Flowering.DOY ~ Year.AD, data=kyoto, type="p", pch=16, cex=0.5, ylab="Day in year of full flowering", xlab="Year (AD)", xlim=c(1000,1050)) lines(fit ~ time, data=all.fits, col="blue") lines(fit-se ~ time, data=all.fits, col="blue", lty="dashed") lines(fit+se ~ time, data=all.fits, col="blue", lty="dashed") legend("bottomleft", legend=c("Observation X(t)", expression(paste("Interpolation ", hat(m)(t))), expression(hat(m)(t) %+-% sqrt(Var(X(t) - hat(m(t)))))), lty=c(NA, "solid", "dashed"), pch=c(16, NA, NA), col=c("black", "blue", "blue"), cex=0.75) # lines(all.fits$time, all.fits$fit+all.fits$se, col="blue", lty="dashed") ``` ## Finally, some interpolation: zoom out ```{r} plot(Flowering.DOY ~ Year.AD, data=kyoto, type="p", ylab="Day in year of full flowering", xlab="Year (AD)", pch=16, cex=0.5) lines(first.year:last.year, all.fits$fit, col="blue") lines(first.year:last.year, all.fits$fit-all.fits$se, col="blue", lty="dashed") lines(first.year:last.year, all.fits$fit+all.fits$se, col="blue", lty="dashed") ``` ## Some extrapolation ```{r, echo=TRUE} plot(Flowering.DOY ~ Year.AD, data=kyoto, type="p", pch=16, cex=0.5, xlim=c(1970, 2030)) new.fits <- kyoto.prediction(1970:2030) lines(fit ~ time, data=new.fits, col="blue") lines(fit-se ~ time, data=new.fits, col="blue", lty="dashed") lines(fit+se ~ time, data=new.fits, col="blue", lty="dashed") ``` ## Some further extrapolation ```{r} plot(Flowering.DOY ~ Year.AD, data=kyoto, type="p", pch=16, cex=0.5, xlim=c(1970, 2060)) new.fits <- kyoto.prediction(1970:2060) lines(fit ~ time, data=new.fits, col="blue") lines(fit-se ~ time, data=new.fits, col="blue", lty="dashed") lines(fit+se ~ time, data=new.fits, col="blue", lty="dashed") abline(h=mean(kyoto$Flowering.DOY, na.rm=TRUE), col="grey") ``` What's going on here? ## 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 ```{r, echo=TRUE} # 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 ```{r} plot(Flowering.DOY ~ Year.AD, data=kyoto, type="l", ylab="Day in year of full flowering", xlab="Year (AD)", main="Cherry blossoms in Kyoto") rug(side=1, x=na.omit(kyoto)$Year.AD) lines(kyoto.spline, col="red", lwd=2) ``` ## Detrending the cherry blossoms with a spline ```{r} plot(fluctuation ~ Year.AD, data=kyoto, type="l", ylab="Deviation from spline fit", xlab="Year (AD)", main="Detrended cherry blossoms in Kyoto") rug(side=1, x=na.omit(kyoto)$Year.AD) ``` ## Detrending the cherry blossoms with a spline ```{r} acf(kyoto$fluctuation, type="covariance", na.action=na.pass, lag.max=100) ``` ## 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 ```{r} par(mfrow=c(1,2)) with(kyoto, acf(Flowering.DOY[Year.AD < 1410], lag.max=100, type="covariance", na.action=na.pass, main="ACF from years < 1410")) with(kyoto, acf(Flowering.DOY[Year.AD >= 141], lag.max=100, type="covariance", na.action=na.pass, main="ACF from years >= 1410")) par(mfrow=c(1,1)) ``` ## Checking stationarity: after spline detrending ```{r} par(mfrow=c(1,2)) with(kyoto, acf(fluctuation[Year.AD < 1410], lag.max=100, type="covariance", na.action=na.pass, main="ACF from years < 1410")) with(kyoto, acf(fluctuation[Year.AD >= 141], lag.max=100, type="covariance", na.action=na.pass, main="ACF from years >= 1410")) par(mfrow=c(1,1)) ``` ## 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? ```{r, fig.retina=NULL, out.width=400, echo=FALSE} knitr::include_graphics("wiener.png") ``` - 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-ex-prodigy + Worked _very_ closely with engineers + His ideas are why we talk about "feedback", "information", "cyber-", etc. [@Wiener-cybernetics;@Wiener-human] + (Photo from @Kline-cybernetics-moment) ## 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-time-series + Parallel work by @Kolmogorov-interpolation-extrapolation * I've never read any explanation for what Kolmogorov was up to - Wiener and Kolmogorov's two big ideas: 1. Optimal linear prediction didn't need the usual regression assumptions (Gaussian noise, independent variables, etc.) 2. 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? ```{r, echo=TRUE} all.fits$std.residuals <- (kyoto$Flowering.DOY - all.fits$fit)/all.fits$se mean(all.fits$std.residuals, na.rm=TRUE) sd(all.fits$std.residuals, na.rm=TRUE) ``` Ideally: mean 0, standard deviation 1 ## Backup: How good is the optimal linear predictor? ```{r} par(mfrow=c(1,2)) plot(kyoto$Year.AD, all.fits$std.residuals, type="p", pch=16, cex=0.5, xlab="Year (AD)", ylab="Standardized residuals from Wiener interpolation", main="Actual residuals") plot(kyoto$Year.AD, sample(all.fits$std.residuals), type="p", pch=16, cex=0.5, col="red", ylab="", main="Randomly shuffled residuals", xlab="") par(mfrow=c(1,1)) ``` 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 ```{r} par(pty="s") plot(all.fits$fit, predict(kyoto.spline, x=first.year:last.year)$y, xlab="Fitted value from Wiener interpolation", ylab="Fitted value from spline smoothing") abline(a=0,b=1, col="grey") par(pty="m") ``` ## 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-synthesis-of-linear-smoothers; @Parzen-time-series-analysis ] * Pretty nearly the first use of RKHS in statistics ## Backup: Hilbert space A **Hilbert space** $\mathcal{H}$ is a space which 1. 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... 2. 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, \HilbertInner{v}{v}{\mathcal{H}} \geq 0$, $\HilbertInner{v}{v}{\mathcal{H}} = 0$ iff $v=0$ 3. Contains the limit of all Cauchy-convergent sequences[^Cauchy]. 4. 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: 1. For each $t \in \mathbb{D}$, $K(\cdot, t) \in \mathcal{H}$ 2. 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 [^Cauchy]: A sequence $v_n$ is **Cauchy** if, no matter how small an $\epsilon > 0$ you name, I can always find an $n(\epsilon)$ such that all points in the sequence after $n(\epsilon)$ are within $\epsilon$ of each other, i.e., $n, m \geq n(\epsilon)$ $\Rightarrow$ $\|v_n - v_m\| \leq \epsilon$. ## References