---
title: Linear Prediction for Time Series
date: 20 September 2018
author: 36-467/36-667
output: slidy_presentation
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,
               tidy=TRUE,  tidy.opts=list(comment=FALSE),
               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}}
\]

## In our last episode

- General approach to optimal linear prediction
    + Predict $Y$ from $\vec{Z} = [Z_1, Z_2, \ldots Z_p]$
	+ Prediction is $\alpha + \vec{\beta} \cdot \vec{Z}$
	+ Best $\alpha = \Expect{Y} - \vec{\beta} \cdot \Expect{\vec{Z}}$
	+ Best $\beta = \Var{\vec{Z}}^{-1} \Cov{\vec{Z}, Y}$
- Today: application to time series
    + What can this do for us?
	+ How do we find the covariances?

## Optimal linear prediction for time series {.smaller}

- Given: $X(t_1), X(t_2), \ldots X(t_n)$
- Desired: prediction of $X(t_0)$

\begin{eqnarray}
\EstRegFunc(t_0) & = & \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

- 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/18/data/kyoto.csv")
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)
```

- A lot of this 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)
```

## When _did_ the cherries flower in 1015?

> - We need $\Cov{X(1015), X(t_i)}$ for every year $t_i$ where we _do_ 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 was $X(2016)$?
- What will $X(2019)$ 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

```{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)))
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

> - A time series is **weakly stationary** when $\Expect{X(t)} = \Expect{X(s)}$ (constant in time) and $\Cov{X(t), X(s)} = \gamma(|t-s|)$
>     + $\gamma(h)$ is the **autocovariance function** at **lag** $h$
>     + a.k.a. "second-order" or "wide-sense" stationarity
> - Assumes all the expectations are the same
> - Assumes $\Cov{X(t), X(t+h)} = \Cov{X(s), X(s+h)}$ so lots of covariances are the same
> - 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 & \gamma(h)
\end{eqnarray}

## The autocovariance function

- $\gamma(0) = \Var{X(t)}$ (constant in $t$)
- $\rho(h) \equiv \frac{\gamma(h)}{\gamma(0)} =$ **autocorrelation** function

## Assuming stationarity...

\begin{eqnarray}
\EstRegFunc(t_0) & = & \overline{x} + \left[\begin{array}{c} X(t_1) -\overline{x}\\
X(t_2) - \overline{x} \\ \vdots \\ X(t_n) -\overline{x}\end{array}\right]
{\left[\begin{array}{ccc} \gamma(0) & \gamma(|t_1 - t_2|) & \ldots & \gamma(|t_1 - t_n|)\\
\gamma(|t_2-t_1|) & \gamma(0) & \ldots & \gamma(|t_2-t_n|)\\
\vdots & \vdots & \ldots & \vdots\\
\gamma(|t_n-t_1|) & \gamma(|t_n-t_2|) & \ldots & \gamma(0)\end{array}\right]}^{-1}
\left[\begin{array}{c} \gamma(|t_0 - t_1|)\\ \gamma(|t_0-t_2|)\\ \vdots\\ \gamma(|t_0 -t_n|)\end{array}\right]\\
& = & \overline{x} + \left[\begin{array}{c} X(t_1) -\overline{x}\\
X(t_2) - \overline{x} \\ \vdots \\ X(t_n) -\overline{x}\end{array}\right] \mathbf{v}^{-1} \mathbf{c}\\
\Var{X(t_0) - \EstRegFunc(t_0)} & = & \gamma(0) - \mathbf{c}^T \mathbf{v}^{-1} \mathbf{c}
\end{eqnarray}

## In-class exercise

\begin{eqnarray}
\EstRegFunc(t_0) & = & \overline{x} + \left[\begin{array}{c} X(t_1) -\overline{x}\\
X(t_2) - \overline{x} \\ \vdots \\ X(t_n) -\overline{x}\end{array}\right]
{\left[\begin{array}{ccc} \gamma(0) & \gamma(|t_1 - t_2|) & \ldots & \gamma(|t_1 - t_n|)\\
\gamma(|t_2-t_1|) & \gamma(0) & \ldots & \gamma(|t_2-t_n|)\\
\vdots & \vdots & \ldots & \vdots\\
\gamma(|t_n-t_1|) & \gamma(|t_n-t_2|) & \ldots & \gamma(0)\end{array}\right]}^{-1}
\left[\begin{array}{c} \gamma(|t_0 - t_1|)\\ \gamma(|t_0-t_2|)\\ \vdots\\ \gamma(|t_0 -t_n|)\end{array}\right]\\
\gamma(h) & = & \gamma(0) \rho(h)
\end{eqnarray}

Re-write the predictor in terms of $\rho$ --- do we need $\gamma(0)$?

## In-class exercise --- solution

\begin{eqnarray}
\EstRegFunc(t_0) & = & \overline{x} + \left[\begin{array}{c} X(t_1) -\overline{x}\\
X(t_2) - \overline{x} \\ \vdots \\ X(t_n) -\overline{x}\end{array}\right]
{\left[\begin{array}{ccc} \gamma(0) & \gamma(|t_1 - t_2|) & \ldots & \gamma(|t_1 - t_n|)\\
\gamma(|t_2-t_1|) & \gamma(0) & \ldots & \gamma(|t_2-t_n|)\\
\vdots & \vdots & \ldots & \vdots\\
\gamma(|t_n-t_1|) & \gamma(|t_n-t_2|) & \ldots & \gamma(0)\end{array}\right]}^{-1}
\left[\begin{array}{c} \gamma(|t_0 - t_1|)\\ \gamma(|t_0-t_2|)\\ \vdots\\ \gamma(|t_0 -t_n|)\end{array}\right]\\
\gamma(h) & = & \gamma(0) \rho(h)\\
\EstRegFunc(t_0) & = & \overline{x} + \left[\begin{array}{c} X(t_1) -\overline{x}\\
X(t_2) - \overline{x} \\ \vdots \\ X(t_n) -\overline{x}\end{array}\right]
{\left[\begin{array}{ccc} 1 & \rho(|t_1 - t_2|) & \ldots & \rho(|t_1 - t_n|)\\
\rho(|t_2-t_1|) & 1 & \ldots & \rho(|t_2-t_n|)\\
\vdots & \vdots & \ldots & \vdots\\
\rho(|t_n-t_1|) & \rho(|t_n-t_2|) & \ldots & 1\end{array}\right]}^{-1} \frac{1}{\gamma(0)} \gamma(0)
\left[\begin{array}{c} \rho(|t_0 - t_1|)\\ \rho(|t_0-t_2|)\\ \vdots\\ \rho(|t_0 -t_n|)\end{array}\right]\\
\end{eqnarray}







## In R

```{r, eval=FALSE, echo=TRUE}
acf(x, lag.max, type, plot, na.action, ...)
```

> - `x` = time series _or_ data frame of multiple time series
> - `lag.max` = maximum value of lag $h$
> - `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}

```{r, echo=TRUE}
# Some bad programming practices here: presumes the existence of kyoto,
# kyoto.acf, etc. the magic number 49, etc.
# Can you make those arguments to the function, so this
# could be used on other, similar time series?
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=2)
    colnames(fits) <- c("fit", "se")
    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])
        # 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))
}
```

## Finally, some interpolation

```{r, echo=TRUE}
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))
all.fits <- kyoto.prediction(801:2015)
lines(801:2015, all.fits$fit, col="blue")
lines(801:2015, all.fits$fit-all.fits$se, col="blue", lty="dashed")
lines(801:2015, all.fits$fit+all.fits$se, col="blue", lty="dashed")
```

## Finally, some interpolation

```{r}
plot(Flowering.DOY ~ Year.AD, data=kyoto, type="l",
     ylab="Day in year of full flowering", xlab="Year (AD)")
all.fits <- kyoto.prediction(801:2015)
lines(801:2015, all.fits$fit, col="blue")
lines(801:2015, all.fits$fit-all.fits$se, col="blue", lty="dashed")
lines(801:2015, 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(1940,2020))
new.fits <- kyoto.prediction(1940:2020)
lines(1940:2020, new.fits$fit, col="blue")
lines(1940:2020, new.fits$fit-new.fits$se, col="blue", lty="dashed")
lines(1940:2020, new.fits$fit+new.fits$se, col="blue", lty="dashed")
```

## 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)} = \gamma(|t-s|)$
- Assuming a flat trend is pretty strong...

## Stationary fluctuations around a non-stationary 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)$
	+ Now get the coefficients and predict

## Detrending the cherry blossoms with a spline

```{r, echo=TRUE}
kyoto.spline <- with(na.omit(kyoto), smooth.spline(x=Year.AD, y=Flowering.DOY))
residuals.with.NAs <- rep(NA, times=nrow(kyoto))
residuals.with.NAs[!is.na(kyoto$Flowering.DOY)] <- na.omit(kyoto$Flowering.DOY) - kyoto.spline$y
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., $\gamma(h) = \gamma(0)e^{-h/\tau}$
	+ Fancier tricks (when we do Fourier methods)

## Checking stationarity

- Divide the data into intervals
- Re-estimate ACF on each interval
- Shouldn't change much
    + How much is too much? (Will answer later with simulation)
- Some inherent pitfalls:
    + Very slight non-stationarity looks pretty stationary
	+ Very slowly decaying correlations look pretty non-stationary

## Multiple time series

- $X$ and $U$ are jointly weakly stationary when
    + $X$ is weakly stationary, ACF $\gamma_X$
	+ $U$ is weakly stationary, ACF $\gamma_U$
	+ $\Cov{X(t), U(s)} = \gamma_{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$
    + $\gamma_U$ goes into the variance matrix $\mathbf{v}$
	+ $\gamma_{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("cover.png")
```

- Norbert Wiener (1894--1964)
    + One of the great mathematicians of the 20th century
	+ Worked _very_ closely with engineers
	+ Why we talk about "feedback", "information", "cyber-", etc.
	+ (Photo from @Kline-cybernetics-moment)


## Wiener's research in 1942

- Anti-aircraft fire control
    + Aim at where the target **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 in 1940s electronics
	+ Declassified after the war as @Wiener-time-series
    + Related work by @Kolmogorov-interpolation-extrapolation
- Wiener and Kolmogorov were the first to realize optimal linear prediction didn't need the usual regression assumptions (Gaussian noise, independent variables, etc.)

## 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: 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
- Solving the equation $\Var{\vec{Z}} \vec{\beta} = \Cov{\vec{Z}, Y}$ for _fixed_ right-hand side can be done in $O(n^2)$ operations
- Details left to numerical linear algebra texts

## Backup: How good is the prediction standard error?

```{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 prediction standard error?

```{r}
plot(kyoto$Year.AD, all.fits$std.residuals, type="l",
     xlab="Year (AD)", ylab="Standardized residuals from Wiener interpolation")
```



## References
