---
title: "Seasonal Arima"
subtitle: "Lecture 11"
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)

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

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


## ARIMA - General Guidance

::: {.medium}
1. Positive autocorrelations out to a large number of lags usually indicates a need for differencing

2. Slightly too much or slightly too little differencing can be corrected by adding AR or MA terms respectively.

3. A model with no differencing usually includes a constant term, a model with two or more orders (rare) differencing usually does not include a constant term. 

4. After differencing, if the PACF has a sharp cutoff then consider adding AR terms to the model.

5. After differencing, if the ACF has a sharp cutoff then consider adding an MA term to the model.

6. It is possible for an AR term and an MA term to cancel each other's effects, so try models with fewer AR terms and fewer MA terms.
:::

::: {.aside}
Based on guidance in Robert Nau's notes on [regression and time series analysis](https://people.duke.edu/~rnau/411home.htm) - Chapter 5
:::

# Electrical Equipment Sales 

## Data

```{r echo=FALSE}
data(elecequip, package = "fpp")
elec_sales = forecast::seasadj(stl(elecequip, s.window="periodic"))

par(mar=c(2,2,2,2))
forecast::tsdisplay(elec_sales, points=FALSE)
```

## 1st order differencing

```{r echo=FALSE}
forecast::ggtsdisplay(diff(elec_sales,1), points=FALSE)
```

## 2nd order differencing

```{r echo=FALSE}
forecast::ggtsdisplay(diff(elec_sales,2), points=FALSE)
```

## Model

```{r}
forecast::Arima(elec_sales, order = c(3,1,0))
```

## Residuals

::: {.small}
```{r}
forecast::Arima(elec_sales, order = c(3,1,0)) %>% residuals() %>% 
  forecast::ggtsdisplay(points=FALSE)
```
:::

## Model Comparison

Model choices:
```{r}
forecast::Arima(elec_sales, order = c(3,1,0))$aicc

forecast::Arima(elec_sales, order = c(3,1,1))$aicc

forecast::Arima(elec_sales, order = c(4,1,0))$aicc

forecast::Arima(elec_sales, order = c(2,1,0))$aicc
```

## Automatic selection (AICc)

```{r}
forecast::auto.arima(elec_sales)
```

## Automatic selection (BIC)

```{r}
forecast::auto.arima(elec_sales, ic = "bic")
```

## Model fit

```{r}
plot(elec_sales, lwd=2, col=adjustcolor("black", alpha.f=0.75))
forecast::Arima(elec_sales, order = c(3,1,0)) %>% fitted() %>% 
  lines(col=adjustcolor('red',alpha.f=0.75),lwd=2)
```

## Model forecast

```{r}
forecast::Arima(elec_sales, order = c(3,1,0)) %>% 
  forecast::forecast() %>% autoplot()
```

## Model forecast - Zoom

```{r message=FALSE}
forecast::Arima(elec_sales, order = c(3,1,0)) %>% 
  forecast::forecast() %>% autoplot() + xlim(2009,2014)
```


# Seasonal Models


## Australian Wine Sales Example

Australian total wine sales by wine makers in bottles <= 1 litre. Jan 1980 – Aug 1994.

```{r echo=FALSE}
data(wineind, package = "forecast")
forecast::ggtsdisplay(wineind, points = FALSE)
```


## Differencing

```{r echo=FALSE}
forecast::ggtsdisplay(diff(wineind), points = FALSE)
```


## Seasonal Arima


We can extend the existing ARIMA model to handle these higher order lags (without having to include all of the intervening lags).


Seasonal $\text{ARIMA}\,(p,d,q) \times (P,D,Q)_s$:
$$ \Phi_P(L^s) \, \phi_p(L) \, \Delta_s^D \, \Delta^d \, y_t = \delta + \Theta_Q(L^s) \, \theta_q(L) \, w_t$$
. . .

where

::: {.small}
$$
\begin{aligned}
\phi_p(L) &= 1-\phi_1 L - \phi_2 L^2 - \ldots - \phi_p L^p\\
\theta_q(L) &= 1+\theta_1 L + \theta_2 L^2 + \ldots + \theta_p L^q \\
\Delta^d &= (1-L)^d\\
\\
\Phi_P(L^s) &= 1-\Phi_1 L^s - \Phi_2 L^{2s} - \ldots - \Phi_P L^{Ps} \\
\Theta_Q(L^s) &= 1+\Theta_1 L + \Theta_2 L^{2s} + \ldots + \theta_p L^{Qs} \\
\Delta_s^D &= (1-L^s)^D\\
\end{aligned}
$$
:::


## Seasonal ARIMA - AR

Lets consider an $\text{ARIMA}(0,0,0) \times (1,0,0)_{12}$:
$$
\begin{aligned}
(1-\Phi_1 L^{12}) \, y_t = \delta + w_t \\
y_t = \Phi_1 y_{t-12} + \delta + w_t
\end{aligned}
$$

```{r}
(m1.1 = forecast::Arima(wineind, seasonal=list(order=c(1,0,0), period=12)))
```


## Fitted - Model 1.1

```{r echo=FALSE, fig.height=4}
cbPalette = c("#000000", "#ff3333", "#92c5de")

model = "Model 1.1 - Arima (0,0,0) x (1,0,0)[12]"
rmse = (wineind-m1.1$fitted)^2 %>% mean() %>% sqrt() %>% round(2) %>% paste0("      [RMSE: ", . ,"]")

m1.1 %>%
  {tibble(
    wineind = wineind %>% unclass(),
    model   = .$fitted %>% unclass(),
    time    = time(.$fitted) %>% unclass()
  )} %>% 
  tidyr::gather(var, sales, -time) %>%
  mutate(var = forcats::as_factor(var)) %>%
  ggplot(aes(x=time, y=sales, color=var)) + 
    geom_line(alpha=0.75, size=0.8) +
    scale_colour_manual(values=cbPalette, name="") +
    labs(title=paste(model, rmse))
```


## Seasonal Arima - Diff


Lets consider an $\text{ARIMA}(0,0,0) \times (0,1,0)_{12}$:
$$
\begin{aligned}
(1 - L^{12}) \, y_t = \delta + w_t \\
y_t = y_{t-12} + \delta + w_t
\end{aligned}
$$

```{r}
(m1.2 = forecast::Arima(
  wineind, seasonal=list(order=c(0,1,0), period=12)
))
```



## Fitted - Model 1.2

```{r echo=FALSE, fig.height=4}
cbPalette = c("#000000", "#ff3333", "#92c5de")

model = "Model 1.2 - Arima (0,0,0) x (0,1,0)[12]"
rmse = (wineind-m1.2$fitted)^2 %>% mean() %>% sqrt() %>% round(2) %>% paste0("      [RMSE: ", . ,"]")

m1.2 %>%
  {tibble(
    wineind = wineind %>% unclass(),
    model   = .$fitted %>% unclass(),
    time    = time(.$fitted) %>% unclass()
  )} %>% 
  tidyr::gather(var, sales, -time) %>%
  mutate(var = forcats::as_factor(var)) %>%
  ggplot(aes(x=time, y=sales, color=var)) + 
    geom_line(alpha=0.75, size=0.8) +
    scale_colour_manual(values=cbPalette, name="") +
    labs(title=paste(model, rmse))
```


## Residuals - Model 1.1 (SAR)

```{r echo=FALSE}
forecast::ggtsdisplay(m1.1$residuals)
```



## Residuals - Model 1.2 (SDiff)

```{r echo=FALSE}
forecast::ggtsdisplay(m1.2$residuals)
```



## Model 2


$\text{ARIMA}(0,0,0) \times (0,1,1)_{12}$:

::: {.small}
$$
\begin{aligned}
(1-L^{12}) y_t = \delta + (1+\Theta_1 L^{12})  w_t \\
y_t = \delta + y_{t-12} + w_t + \Theta_1 w_{t-12}
\end{aligned}
$$
:::


```{r}
(m2 = forecast::Arima(wineind, order=c(0,0,0), 
                      seasonal=list(order=c(0,1,1), period=12)))
```


## Fitted - Model 2

```{r echo=FALSE, fig.height=4}
model = "Model 2 - forecast::Arima (0,0,0) x (0,1,1)[12]"
rmse = (wineind-m2$fitted)^2 %>% mean() %>% sqrt() %>% round(2) %>% paste0("      [RMSE: ", . ,"]")

m2 %>%
  {tibble(
    wineind = wineind %>% unclass(),
    model   = .$fitted %>% unclass(),
    time    = time(.$fitted) %>% unclass()
  )} %>% 
  tidyr::gather(var, sales, -time) %>%
  mutate(var = forcats::as_factor(var)) %>%
  ggplot(aes(x=time, y=sales, color=var)) + 
    geom_line(alpha=0.75, size=0.8) +
    scale_colour_manual(values=cbPalette, name="") +
    labs(title=paste(model, rmse))
```


## Residuals

```{r echo=FALSE}
forecast::ggtsdisplay(m2$residuals)
```



## Model 3

$\text{ARIMA}(3,0,0) \times (0,1,1)_{12}$

::: {.small}
$$
\begin{aligned}
(1-\phi_1 L - \phi_2 L^2 - \phi_3 L^3) \, (1-L^{12}) y_t = \delta + (1 + \Theta_1 L)w_t \\
y_t = \delta + \sum_{i=1}^3 \phi_i y_{t-i}  + y_{t-12}  - \sum_{i=1}^3 \phi_i y_{t-12-i} + w_t + w_{t-12}
\end{aligned}
$$
:::

::: {.medium}
```{r}
(m3 = forecast::Arima(wineind, order=c(3,0,0), 
                      seasonal=list(order=c(0,1,1), period=12)))
```
:::

## Fitted model

```{r echo=FALSE, fig.height=4}
model = "Model 3 - forecast::Arima (3,0,0) x (0,1,1)[12]"
rmse = (wineind-m3$fitted)^2 %>% mean() %>% sqrt() %>% round(2) %>% paste0("      [RMSE: ", . ,"]")

m3 %>%
  {tibble(
    wineind = wineind %>% unclass(),
    model   = .$fitted %>% unclass(),
    time    = time(.$fitted) %>% unclass()
  )} %>% 
  tidyr::gather(var, sales, -time) %>%
  mutate(var = forcats::as_factor(var)) %>%
  ggplot(aes(x=time, y=sales, color=var)) + 
    geom_line(alpha=0.75, size=0.8) +
    scale_colour_manual(values=cbPalette, name="") +
    labs(title=paste(model, rmse))
```


## Model - Residuals

```{r echo=FALSE}
forecast::ggtsdisplay(m3$residuals)
```


# Federal Reserve Board Production Index


## `prodn` from the `astsa` package

Monthly Federal Reserve Board Production Index (1948-1978)

::: {.small}
```{r}
data(prodn, package="astsa")
forecast::ggtsdisplay(prodn, points = FALSE)
```
:::


## Differencing

Based on the ACF it seems like standard differencing may be required

. . .

```{r}
forecast::ggtsdisplay(diff(prodn))
```


## Differencing + Seasonal Differencing

Additional seasonal differencing also seems warranted

:::: {.columns .small}
::: {.column width='50%'}
```{r}
(fr_m1 = forecast::Arima(
  prodn, order = c(0,1,0), 
  seasonal = list(order=c(0,0,0), period=12)
))
```
:::

::: {.column width='50%'}
```{r}
(fr_m2 = forecast::Arima(
  prodn, order = c(0,1,0), 
  seasonal = list(order=c(0,1,0), period=12)
))
```
:::
::::

. . .

:::: {.columns .small}
::: {.column width='50%'}
```{r}
yardstick::rmse_vec(
  prodn %>% unclass(), 
  fr_m1$fitted %>% unclass()
)
```
:::

::: {.column width='50%'}
```{r}
yardstick::rmse_vec(
  prodn %>% unclass(), 
  fr_m2$fitted %>% unclass()
)
```
:::
::::



## Residuals - Model 2

```{r}
forecast::ggtsdisplay(fr_m2$residuals, points=FALSE, lag.max=36)
```


## Adding Seasonal MA

:::: {.columns .small}
::: {.column width='50%'}
```{r}
(fr_m3.1 = forecast::Arima(
  prodn, order = c(0,1,0), 
  seasonal = list(order=c(0,1,1), period=12)
))
```
:::

::: {.column width='50%'}
```{r}
(fr_m3.2 = forecast::Arima(
  prodn, order = c(0,1,0), 
  seasonal = list(order=c(0,1,2), period=12)
))
```
:::
::::

. . .

:::: {.columns .small}
::: {.column width='50%'}
```{r}
yardstick::rmse_vec(
  prodn %>% unclass(), 
  fr_m3.1$fitted %>% unclass()
)
```
:::

::: {.column width='50%'}
```{r}
yardstick::rmse_vec(
  prodn %>% unclass(), 
  fr_m3.2$fitted %>% unclass()
)
```
:::
::::


## Adding Seasonal MA (cont.)

:::: {.columns .small}
::: {.column width='50%'}
```{r}
(fr_m3.3 = forecast::Arima(
  prodn, order = c(0,1,0), 
  seasonal = list(order=c(0,1,3), period=12)
))
```
:::
::::

. . .

:::: {.columns .small}
::: {.column width='50%'}
```{r}
yardstick::rmse_vec(
  prodn %>% unclass(), 
  fr_m3.3$fitted %>% unclass()
)
```
:::
::::


## Residuals - Model 3.3

```{r echo=FALSE}
forecast::ggtsdisplay(fr_m3.3$residuals, points=FALSE, lag.max = 36)
```



## Adding AR

:::: {.columns .small}
::: {.column width='50%'}
```{r}
(fr_m4.1 = forecast::Arima(
  prodn, order = c(1,1,0), 
  seasonal = list(order=c(0,1,3), period=12)
))
```
:::

::: {.column width='50%'}
```{r}
(fr_m4.2 = forecast::Arima(
  prodn, order = c(2,1,0), 
  seasonal = list(order=c(0,1,3), period=12)
))
```
:::
::::

. . .

:::: {.columns .small}
::: {.column width='50%'}
```{r}
yardstick::rmse_vec(
  prodn %>% unclass(), 
  fr_m4.1$fitted %>% unclass()
)
```
:::

::: {.column width='50%'}
```{r}
yardstick::rmse_vec(
  prodn %>% unclass(), 
  fr_m4.2$fitted %>% unclass()
)
```
:::
::::





## Residuals - Model 4.1

```{r echo=FALSE}
forecast::ggtsdisplay(fr_m4.1$residuals, points=FALSE, lag.max = 36)
```


## Residuals - Model 4.2

```{r echo=FALSE}
forecast::ggtsdisplay(fr_m4.2$residuals, points=FALSE, lag.max = 36)
```


## Model Fit

```{r echo=FALSE, fig.height=4}
model = "Model 4.1 - forecast::Arima (1,1,0) x (0,1,3)[12]"
rmse = (prodn-fr_m4.1$fitted)^2 %>% mean() %>% sqrt() %>% round(3) %>% paste0("      [RMSE: ", . ,"]")

fr_m4.1 %>%
  {tibble(
    prodn = prodn %>% unclass(),
    model = .$fitted %>% unclass(),
    time  = time(.$fitted) %>% unclass()
  )} %>% 
  tidyr::gather(var, sales, -time) %>%
  mutate(var = forcats::as_factor(var)) %>%
  ggplot(aes(x=time, y=sales, color=var)) + 
    geom_line(alpha=0.75, size=0.8) +
    scale_colour_manual(values=cbPalette, name="") +
    labs(title=paste(model, rmse))
```


## Model Forecast

```{r}
forecast::forecast(fr_m4.1) %>% plot()
```


## Model Forecast (cont.)

```{r}
forecast::forecast(fr_m4.1) %>% plot(xlim=c(1975,1982))
```


## Model Forecast (cont.)

```{r}
forecast::forecast(fr_m4.1, 120) %>% plot()
```


## Model Forecast (cont.)

```{r}
forecast::forecast(fr_m4.1, 120) %>% plot(xlim=c(1975,1990))
```



## Auto ARIMA - Model Fit

```{r echo=FALSE, fig.height=4}
fr_aa = forecast::auto.arima(prodn)
model = "Model Auto ARIMA - forecast::auto.arima (2,0,1) x (0,1,1)[12]"
rmse = (prodn-fr_aa$fitted)^2 %>% mean() %>% sqrt() %>% round(3) %>% paste0("      [RMSE: ", . ,"]")

fr_aa %>%
  {tibble(
    prodn = prodn %>% unclass(),
    model = .$fitted %>% unclass(),
    time  = time(.$fitted) %>% unclass()
  )} %>% 
  tidyr::gather(var, sales, -time) %>%
  mutate(var = forcats::as_factor(var)) %>%
  ggplot(aes(x=time, y=sales, color=var)) + 
    geom_line(alpha=0.75, size=0.8) +
    scale_colour_manual(values=cbPalette, name="") +
    labs(title=paste(model, rmse))
```


## Exercise - Cortecosteroid Drug Sales

Monthly cortecosteroid drug sales in Australia from 1992 to 2008.

```{r fig.height=4}
data(h02, package="fpp")
forecast::ggtsdisplay(h02, points=FALSE)
```

