---
title: "Fitting ARIMA Models"
subtitle: "Lecture 12"
author: "Dr. Colin Rundel"
footer: "Sta 344 - Fall 2022"
format:
  revealjs:
    theme: slides.scss
    transition: fade
    slide-number: true
execute:
  echo: true
  warning: false
  collapse: true
---

```{r setup, include=FALSE}
library(tidyverse)
library(forecast)
library(brms)

knitr::opts_chunk$set(
  fig.align = "center"
)

ggplot2::theme_set(ggplot2::theme_bw())

source("util-crps.R")
source("fix_pred_draws.R")

set.seed(20221010)
```



# Model Fitting

## Fitting ARIMA

For an $ARIMA(p,d,q)$ model, 

* Assumes that the data is stationary after differencing

* Handling $d$ is straight forward, just difference the original data $d$ times (leaving $n-d$ observations)
$$ y'_t = \Delta^d \, y_t $$

* After differencing, fit an $ARMA(p,q)$ model to $y'_t$.

* To keep things simple we'll assume $w_t \overset{iid}{\sim} \mathcal{N}(0,\sigma^2_w)$


## MLE - Stationarity & iid normal errors

If both of these assumptions are met, then the time series $y_t$ will also be normal.

. . .

In general, the vector $\boldsymbol{y} = (y_1, y_2, \ldots, y_t)'$ will have a multivariate normal distribution with mean $\{\boldsymbol\mu\}_i = E(y_i) = E(y_t)$  and covariance $\boldsymbol\Sigma$ where $\{\boldsymbol{\Sigma}\}_{ij} = \gamma(i-j)$.

The joint density of $\boldsymbol y$ is given by

$$ f_{\boldsymbol y}(\boldsymbol y) = \frac{1}{(2\pi)^{t/2}\,\det(\boldsymbol\Sigma)^{1/2}} \times \exp\left( -\frac{1}{2}(\boldsymbol y - \boldsymbol \mu)' \, \Sigma^{-1} \, (\boldsymbol y - \boldsymbol \mu) \right) $$

# AR

## Fitting $AR(1)$

$$ y_t = \delta + \phi \, y_{t-1} + w_t $$

We need to estimate three parameters: $\delta$, $\phi$, and $\sigma_w^2$, we know

$$ 
\begin{aligned}
E(y_t) = \frac{\delta}{1-\phi} \quad&\quad Var(y_t) = \frac{\sigma_w^2}{1-\phi^2} \\
\gamma(h) &= \frac{\sigma_w^2}{1-\phi^2} \phi^{|h|}
\end{aligned} 
$$

Using these properties it is possible to write the distribution of $\boldsymbol{y}$ as a MVN but that does not make it easy to write down a (simplified) closed form for the MLE estimate for $\delta$, $\phi$, and $\sigma_w^2$.


## Conditional Density

We can also rewrite the density as follows,

$$
\begin{aligned}
f(\boldsymbol y)
  &= f(y_t,\,y_{t-1},\,\ldots,\,y_2,\,y_1) \\
  &= f(y_t|\,y_{t-1},\,\ldots,\,y_2,\,y_1) f(y_{t-1}|y_{t-2},\,\ldots,\,y_2,\,y_1) \cdots f(y_2|y_1) f(y_1) \\
  &= f(y_t|\,y_{t-1}) f(y_{t-1}|y_{t-2}) \cdots f(y_2|y_1) f(y_1)
\end{aligned}
$$

where, 

$$
\begin{aligned}
y_1 &\sim \mathcal{N}\left(\delta, \, \frac{\sigma^2_w}{1-\phi^2} \right) \\
y_{t}|y_{t-1} &\sim \mathcal{N}\left(\delta+\phi\, y_{t-1}, \, \sigma^2_w \right) \\
f(y_{t}|y_{t-1}) &= \frac{1}{\sqrt{2\pi \, \sigma^2_w}} \exp \left( -\frac{1}{2}\frac{(y_t -\delta+\phi\, y_{t-1})^2 }{\sigma^2_w} \right)
\end{aligned}
$$

## Log likelihood of AR(1)

$$
\log f(y_{t} | y_{t-1}) = -\frac{1}{2}\left( \log 2\pi + \log \sigma^2_w + \frac{1}{\sigma_w^2} (y_t -\delta+\phi\, y_{t-1})^2 \right)
$$

. . .

::: {.medium}
$$
\begin{aligned}
\ell(\delta, \phi, \sigma^2_w) 
  &= \log f(\boldsymbol{y}) = \log f(y_1) + \sum_{i=2}^t \log f(y_{i}|y_{i-1}) \\
  &= - \frac{1}{2} \bigg(\log 2\pi + \log \sigma_w^2 - \log (1-\phi^2) + \frac{(1-\phi^2)}{\sigma_w^2 }(y_1-\delta)^2 \bigg) \\
  & ~~~~ - \frac{1}{2} \bigg( (n-1) \log 2\pi + (n-1) \log \sigma_w^2 + \frac{1}{\sigma_w^2} \sum_{i=2}^n (y_i -\delta+\phi\, y_{i-1})^2 \bigg) \\
  &= - \frac{1}{2} \bigg( n \log 2\pi + n \log \sigma_w^2 - \log (1-\phi^2) \\
  &~~~~~~~~~~~~~~~+ \frac{1}{\sigma_w^2} \bigg( (1-\phi^2)(y_1-\delta)^2 + \sum_{i=2}^n (y_i -\delta+\phi\, y_{i-1})^2 \bigg) \bigg)
\end{aligned}
$$
:::

## AR(1) Example

with $\phi = 0.75$, $\delta=0.5$, and $\sigma_w^2=1$,

```{r echo=FALSE}
ar1 = arima.sim(n=500, model = list(order=c(1,0,0), ar=0.75), mean=0.5)
forecast::ggtsdisplay(ar1)
```

## ARIMA

```{r}
( ar1_arima = forecast::Arima(ar1, order = c(1,0,0)) )
```

## mean vs $\delta$?

The reported mean value from the ARIMA model is $E(y_t)$ and not $\delta$ - for an ARIMA(1,0,0) 
$$
E(y_t) = \frac{\delta}{1-\phi} ~~\Rightarrow~~ \delta = E(y_t) * (1-\phi)
$$

. . .

:::: {.columns}
::: {.column width='50%'}
True $E(y_t)$:
```{r}
0.5 / (1-0.75)
```
:::

::: {.column width='50%'}
Sample $\delta$:
```{r}
ar1_arima$coef[2] * 
  (1 - ar1_arima$model$phi)
```
:::
::::




## lm

```{r}
d = tsibble::as_tsibble(ar1) %>%
  as_tibble() %>%
  rename(y = value)
summary({ ar1_lm = lm(y~lag(y), data=d) })
```

## Bayesian AR(1) Model

::: {.small}
```{r cache=TRUE, message=FALSE}
library(brms) # must be loaded for arma to work
( ar1_brms = brm(y ~ arma(p = 1, q = 0), data=d, refresh=0) )
```
:::


## Chains

```{r fig.align="center"}
plot(ar1_brms)
```

## PP Checks

```{r}
pp_check(ar1_brms, ndraws=100)
```


## Posteriors

```{r include=FALSE}
ar1_brms_post = tidybayes::spread_draws(
  ar1_brms, b_Intercept, `ar[1]`, sigma
) %>%
  rename(mean = b_Intercept, phi = `ar[1]`, sigma2_w = sigma) %>%
  mutate(
    delta = mean * (1 - phi)
  ) %>%
  pivot_longer(mean:delta, names_to = ".variable", values_to = ".value")

ar1_ests = bind_rows(
  tibble(
    model = "truth",
    .variable = c("delta", "phi", "sigma2_w"), 
    estimate = c(0.5, 0.75, 1)
  ),
  tibble(
    model = "lm",
    .variable = c("delta", "phi", "sigma2_w"), 
    estimate = c(coef(ar1_lm), var(ar1_lm$residuals))
  ),
  tibble(
    model = "ARIMA",
    .variable = c("delta", "phi", "sigma2_w"), 
    estimate = c((1-ar1_arima$model$phi)*ar1_arima$coef[2], ar1_arima$model$phi, ar1_arima$sigma2)
  )
)
```


```{r echo=FALSE, fig.height=4.5}  
ar1_brms_post %>%
  filter(.chain == 1, .variable != "mean") %>%
  ggplot(aes(x=.value)) +
    geom_density(fill="lightgrey") +
    geom_vline(data=ar1_ests, aes(xintercept = estimate, linetype=model, color=model), size=1.5, alpha=0.75) +
    facet_wrap(~.variable, ncol=3, scales = "free_x") +
    labs(x = "")
```

## Predictions

```{r echo=FALSE}
ar1_brms %>%
  predicted_draws_fix(newdata = d) %>%
  filter(.chain == 1) %>%
  ggplot(aes(y=y, x=index)) +
    geom_point() +
    tidybayes::stat_lineribbon(
      aes(y=.prediction), alpha=0.25
    )
```

## Forecasting

:::{.small}
```{r}
ar1_brms_fc = ar1_brms %>%
  predicted_draws_fix(
    newdata = tibble(index=501:550, y=NA)
  ) %>%
  filter(.chain == 1)
```

```{r echo=FALSE}
ar1_brms %>%
  predicted_draws_fix(newdata = d) %>%
  filter(.chain == 1) %>%
  ggplot(aes(y=y, x=index)) +
    geom_point() +
    geom_line() +
    tidybayes::stat_lineribbon(
      data = ar1_brms_fc,
      aes(y=.prediction), alpha=0.25
    )
```
:::
:::

# Fitting AR(p)



## Lagged Regression

As with the AR(1), we can rewrite the density using conditioning,
$$
\begin{aligned}
f(\boldsymbol y)
  &= f(y_t, \,y_{t-1}, \,\ldots, \,y_{2}, \,y_{1}) \\
  &= f(y_{n}|y_{n-1},\ldots,y_{n-p}) \cdots  f(y_{p+1}|y_p,\ldots,y_1)  f(y_p, \,\ldots, y_1)
\end{aligned}
$$

. . .

Regressing $y_t$ on $y_{t-1}, \ldots, y_{t-p}$ gets us an approximate solution, but it ignores the $f(y_1, \, y_2, \,\ldots, y_p)$ part of the likelihood. 

How much does this matter (vs. using the full likelihood)?

* If $p$ is near to $n$ then probably a lot

* If $p << n$ then probably not much


## Method of Moments

Recall for an AR(p) process,

$$
\begin{aligned}
\gamma(0) &= \sigma^2_w + \phi_1 \gamma(1) + \phi_2 \gamma(2) + \ldots + \phi_p \gamma(p) \\
\gamma(h) &= \phi_1 \gamma(h-1) + \phi_2 \gamma(h-2) + \ldots \phi_p \gamma(h-p)
\end{aligned}
$$
We can rewrite the first equation in terms of $\sigma^2_w$,
$$
\sigma^2_w =  \gamma(0) - \phi_1 \gamma(1) - \phi_2 \gamma(2) - \ldots - \phi_p \gamma(p)
$$
these are called the Yule-Walker equations.

## Yule-Walker

These equations can be rewritten into matrix notation as follows

$$
\underset{p \times p}{\boldsymbol\Gamma_p} \underset{p \times 1}{\boldsymbol\phi} = \underset{p \times 1}{\boldsymbol\gamma_p}
\qquad\qquad
\underset{1 \times 1}{\sigma^2_w} = \underset{1 \times 1}{\gamma(0)} - \underset{1 \times p}{\boldsymbol{\phi'}}\underset{p \times 1}{\boldsymbol{\gamma_p}}
$$
where

::: {.small}
$$ 
\begin{aligned}
\underset{p \times p}{\boldsymbol{\Gamma_p}} &= \{\gamma(j-k)\}_{j,k} \\
\underset{p \times 1}{\boldsymbol\phi} &= (\phi_1, \phi_2, \ldots, \phi_p)' \\
\underset{p \times 1}{\boldsymbol\gamma_p} &= (\gamma(1), \gamma(2), \ldots, \gamma(p))'
\end{aligned}
$$
:::

. . .

If we estimate the covariance structure from the data we obtain $\hat{\boldsymbol\gamma_p}$ and $\hat{\boldsymbol\Gamma_p}$ which we can plug in and solve for $\boldsymbol{\phi}$ and $\sigma^2_w$,
$$
\hat{\boldsymbol\phi} =\hat{\boldsymbol{\Gamma}_p}^{-1}\hat{\boldsymbol{\gamma}_p}
\qquad\qquad
\hat{\sigma}^2_w = \gamma(0) - \hat{\boldsymbol{\gamma}_p}' \hat{\boldsymbol{\Gamma}_p^{-1}} \hat{\boldsymbol{\gamma}_p}
$$

# ARMA

## Fitting $ARMA(2,2)$

$$ y_t = \delta + \phi_1 \, y_{t-1} + \phi_2 \, y_{t-2} + \theta_1 w_{t-1} + \theta_2 w_{t-2} + w_t $$

We now need to estimate six parameters: $\delta$, $\phi_1$, $\phi_2$, $\theta_1$, $\theta_2$ and $\sigma_w^2$.

. . .

$~$

We could figure out $E(y_t)$, $Var(y_t)$, and $Cov(y_t, y_{t+h})$, but the last two are going to be pretty nasty and the full MVN likehood is similarly going to be unpleasant to work with.

. . .

$~$

Like the AR(1) and AR(p) processes we want to use conditioning to simplify things.
\tinyoutput
$$
\begin{aligned}
y_t | \delta, &y_{t-1}, y_{t-2}, w_{t-1}, w_{t-2} \\ 
&\sim \mathcal{N}(\delta + \phi_1 \, y_{t-1} + \phi_2 \, y_{t-2} + \theta_1 w_{t-1} + \theta_2 w_{t-2},~\sigma_w^2) 
\end{aligned}
$$

## ARMA(2,2) Example

with $\phi = (0.75,-0.5)$, $\theta = (0.5,0.2)$, $\delta=0$, and $\sigma_w^2=1$ using the same models 

```{r echo=FALSE}
set.seed(202210122)
y = arima.sim(n=500, model=list(ar=c(0.75,-0.5), ma=c(0.5,0.2))) %>% unclass()
forecast::ggtsdisplay(y, points = FALSE)
```

## ARIMA

```{r}
forecast::Arima(y, order = c(2,0,2), include.mean = FALSE) %>% summary()
```

## AR only lm

::: {.small}
```{r}
lm(y ~ lag(y,1) + lag(y,2)) %>% summary()
```
:::


## Hannan-Rissanen Algorithm

1. Estimate a high order AR (remember AR $\Leftrightarrow$ MA when stationary + invertible)

2. Use AR to estimate values for unobserved $w_t$ via `lm` with `lag`s

3. Regress $y_t$ onto $y_{t-1}, \ldots, y_{t-p}, \hat{w}_{t-1}, \ldots \hat{w}_{t-q}$

4. Update $\hat{w}_{t-1}, \ldots \hat{w}_{t-q}$ based on current model, 

5. Goto 3, repeat until convergence


## Hannan-Rissanen - Step 1 & 2

```{r include=FALSE}
options(width=50)
```

:::: {.columns .small}
::: {.column width='50%'}
```{r}
(ar = ar.mle(y, order.max = 10))
```
:::

::: {.column width='50%'}
```{r}
(ar = forecast::Arima(y, order = c(10,0,0)))
```
:::
::::


```{r include=FALSE}
options(width=80)
```


## Residuals

```{r}
forecast::ggtsdisplay(ar$resid)
```

## Hannan-Rissanen - Step 3

::: {.medium}
```{r}
d = tibble(
  y = y %>% strip_attr(),
  index = seq_along(y),
  w_hat1 = ar$resid %>% strip_attr()
)

(lm1 = lm(y ~ lag(y,1) + lag(y,2) + lag(w_hat1,1) + lag(w_hat1,2), data=d)) %>%
  summary()
```
:::

## Hannan-Rissanen - Step 4

::: {.medium}
```{r}
d = modelr::add_residuals(d,lm1,"w_hat2")

(lm2 = lm(y ~ lag(y,1) + lag(y,2) + lag(w_hat2,1) + lag(w_hat2,2), data=d)) %>%
  summary()
```
:::

## Hannan-Rissanen - Step 3.2 + 4.2

::: {.medium}
```{r}
d = modelr::add_residuals(d,lm2,"w_hat3")

(lm3 = lm(y ~ lag(y,1) + lag(y,2) + lag(w_hat3,1) + lag(w_hat3,2), data=d)) %>%
  summary()
```
:::

## Hannan-Rissanen - Step 3.3 + 4.3


::: {.medium}
```{r}
d = modelr::add_residuals(d,lm3,"w_hat4")

(lm4 = lm(y ~ lag(y,1) + lag(y,2) + lag(w_hat4,1) + lag(w_hat4,2), data=d)) %>%
  summary()
```
:::

## Hannan-Rissanen - Step 3.4 + 4.4

::: {.medium}
```{r}
d = modelr::add_residuals(d,lm4,"w_hat5")

(lm5 = lm(y ~ lag(y,1) + lag(y,2) + lag(w_hat5,1) + lag(w_hat5,2), data=d)) %>%
  summary()
```
:::

## BRMS

::: {.medium}
```{r eval=FALSE}
( arma22_brms = brm(
    y~arma(p=2,q=2)-1, data=d, 
    chains=2, refresh=0, iter = 5000, cores = 4
) )
```
```{r eval=FALSE, include=FALSE}
saveRDS(arma22_brms,"arma22_brms.rds")
```

```{r echo=FALSE}
(arma22_brms = readRDS("arma22_brms.rds"))
```

:::

## Chains

::: {.small}
```{r}
plot(arma22_brms)
```
:::

## Comparison

```{r echo=FALSE}
arma22_ests = bind_rows(
  tibble(
    model = "Truth",
    .variable = c("ar[1]", "ar[2]", "ma[1]", "ma[2]"), 
    .value = c(0.75, -0.5, 0.5, 0.2)
  ),
  tibble(
    model = "Arima",
    .variable = c("ar[1]", "ar[2]", "ma[1]", "ma[2]"), 
    .value =  c(0.7290,  -0.4967,  0.4896,  0.2543)
  ),
  tibble(
    model = "HR",
    .variable = c("ar[1]", "ar[2]", "ma[1]", "ma[2]"), 
    .value =  c(0.75159, -0.50072, 0.46345, 0.22949)
  )
)

arma22_brms_post = tidybayes::gather_draws(
  arma22_brms, ar[i], ma[i]
) %>%
  ungroup() %>%
  mutate(.variable = paste0(.variable,"[",i,"]"))
```

```{r echo=FALSE}
arma22_brms_post %>%
  group_by(.variable) %>%
  filter(.chain == 1) %>%
  ggplot(aes(x=.value)) +
    geom_density(fill="lightgrey") +
    geom_vline(
      data=arma22_ests, 
      aes(xintercept = .value, linetype=model, color=model), 
      size=1.5, alpha=0.75
    ) +
    facet_wrap(~.variable, ncol=2, scales = "free_x") +
    labs(x = "")
```

## Predictions

```{r echo=FALSE}
arma22_brms %>%
  predicted_draws_fix(newdata = d) %>%
  filter(.chain == 1) %>%
  ggplot(aes(y=y, x=index)) +
    geom_point() +
    tidybayes::stat_lineribbon(
      aes(y=.prediction), alpha=0.25
    )
```

## Forecasting

:::{.small}
```{r}
arma22_brms_fc = arma22_brms %>%
  predicted_draws_fix(
    newdata = tibble(index=501:550, y=NA)
  ) %>%
  filter(.chain == 1)
```

```{r echo=FALSE}
arma22_brms %>%
  predicted_draws_fix(newdata = d) %>%
  filter(.chain == 1) %>%
  ggplot(aes(y=y, x=index)) +
    geom_point() +
    geom_line() +
    tidybayes::stat_lineribbon(
      data = arma22_brms_fc,
      aes(y=.prediction), alpha=0.25
    )
```
:::


## Stan Code

::: {.small}
```{r}
arma22_brms %>% stancode()
```
:::