---
title: "STAT481/581: Introduction to Time Series Analysis"
author: "Ch9. ARIMA models"
date: "OTexts.org/fpp3/"
toc: true
colortheme: monashwhite
output:
  binb::monash:
    fig_width: 7
    fig_height: 3.5
    includes:
      in_header: header.tex
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE, cache=TRUE, warning=FALSE, message=FALSE)
options(digits=4, width=60)
library(tidyverse)
library(fable)
library(tsibble)
library(tsibbledata)
library(feasts)
library(lubridate)

usmelec <- as_tsibble(fpp2::usmelec) %>%
  rename(Month = index, Generation = value)

us_change <- read_csv("data/uschange.csv") %>%
  mutate(Time = yearquarter(Time)) %>%
  as_tsibble(index = Time)

eu_retail <- as_tsibble(fpp2::euretail)

h02 <- tsibbledata::PBS %>%
  filter(ATC2 == "H02") %>%
  summarise(Cost = sum(Cost))
```

# Stationarity and differencing

## Stationarity

\begin{block}{Definition}
If $\{y_t\}$ is a strong stationary time series, then for all $s$, the distribution of $(y_t,\dots,y_{t+s})$ does not depend on $t$.
\end{block} 

A **strong stationary series** is:

* roughly horizontal
* constant variance
* no patterns predictable in the long-term

## Stationarity

A **stationary series or weak stationary series** is:

* roughly horizontal
* constant variance
* no patterns predictable in the long-term

## Stationary?

```{r}
gafa_stock %>%
  filter(Symbol == "GOOG", year(Date) == 2018) %>%
  autoplot(Close) + ylab("Google closing stock price") + xlab("Day")
```

## Stationary?

```{r}
gafa_stock %>%
  filter(Symbol == "GOOG", year(Date) == 2018) %>%
  autoplot(difference(Close)) + ylab("Google closing stock price") + xlab("Day")
```

## Stationary?

```{r}
as_tsibble(fma::strikes) %>% autoplot(value) +
  ylab("Number of strikes") + xlab("Year")
```

## Stationary?

```{r}
as_tsibble(fma::hsales) %>% autoplot(value) + xlab("Year") + ylab("Total sales") +
  ggtitle("Sales of new one-family houses, USA")
```

## Stationary?
```{r}
as_tsibble(fma::eggs) %>% autoplot(value) + xlab("Year") + ylab("$") +
  ggtitle("Price of a dozen eggs in 1993 dollars")
```

## Stationary?
```{r}
aus_livestock %>%
  filter(
    Animal == "Pigs",
    State == "Victoria",
    year(Month) >= 2010
  ) %>%
  autoplot(Count/1e3) + xlab("Year") + ylab("thousands") +
  ggtitle("Number of pigs slaughtered in Victoria")
```

## Stationary?

```{r}
pelt %>% autoplot(Lynx) + xlab("Year") + ylab("Number trapped") +
  ggtitle("Annual Canadian Lynx Trappings")
```

## Stationary?

```{r}
aus_production %>% filter_index("1992" ~ .) %>%
  autoplot(Beer) + xlab("Year") + ylab("megalitres") +
  ggtitle("Australian quarterly beer production")
```

## Stationarity

Transformations help to **stabilize the variance**.

For ARIMA modelling, we also need to **stabilize the mean**.

## Non-stationarity in the mean
\alert{Identifying non-stationary series}

* time plot.
* The ACF of stationary data drops to zero relatively quickly
* The ACF of non-stationary data decreases slowly.
* For non-stationary data, the value of $r_1$ is often large and positive.

## Example: Google stock price

```{r}
gafa_stock %>%
  filter(Symbol == "GOOG", year(Date) == 2018) %>%
  autoplot(Close) + ylab("Closing stock price ($USD)") + xlab("Day")
```

## Example: Google stock price

```{r}
google_2018 <- gafa_stock %>%
  filter(Symbol == "GOOG", year(Date) == 2018) %>%
  mutate(trading_day = row_number()) %>%
  update_tsibble(index = trading_day, regular = TRUE)
google_2018 %>%
  ACF(Close) %>% autoplot()
```

## Example: Google stock price

```{r}
gafa_stock %>%
  filter(Symbol == "GOOG", year(Date) == 2018) %>%
  autoplot(difference(Close)) +
  ylab("Change in Google closing stock price ($USD)") + xlab("Day")
```

## Example: Google stock price

```{r}
google_2018 %>% ACF(difference(Close)) %>% autoplot()
```

## Differencing

* Differencing helps to **stabilize the mean**.
* The differenced series is the *change* between each observation in the original series: $y'_t = y_t - y_{t-1}$.
* The differenced series will have only $T-1$ values since it is not possible to calculate a difference $y_1'$ for the first observation.

## Second-order differencing

Occasionally the differenced data will not appear stationary and it may be necessary to difference the data a second time: 
\begin{align*}
y''_{t} & = y'_{t} - y'_{t - 1} \\
        & = (y_t - y_{t-1}) - (y_{t-1}-y_{t-2})\\
        & = y_t - 2y_{t-1} +y_{t-2}.
\end{align*} 

* $y_t''$ will have $T-2$ values.
* In practice, it is almost never necessary to go beyond second-order differences.

## Seasonal differencing

A seasonal difference is the difference between an observation and the corresponding observation from the previous year. 
$$
 y'_t = y_t - y_{t-m}
$$
where $m=$ number of seasons. 

* For monthly data $m=12$.
* For quarterly data $m=4$.

## Electricity production
\fontsize{11}{13}\sf

```{r, echo=TRUE, fig.height=3.5}
usmelec %>% autoplot(
  Generation
)
```

## Electricity production
\fontsize{11}{13}\sf

```{r, echo=TRUE, fig.height=3.5}
usmelec %>% autoplot(
  log(Generation)
)
```

## Electricity production
\fontsize{11}{13}\sf

```{r, echo=TRUE, fig.height=3.5}
usmelec %>% autoplot(
  log(Generation) %>% difference(12)
)
```

## Electricity production
\fontsize{11}{13}\sf

```{r, echo=TRUE, fig.height=3.5}
usmelec %>% autoplot(
  log(Generation) %>% difference(12) %>% difference()
)
```

## Electricity production

* Seasonally differenced series is closer to being stationary.
* Remaining non-stationarity can be removed with further first difference.

If $y'_t = y_t - y_{t-12}$ denotes seasonally differenced series, then twice-differenced series is

\begin{block}{}
\begin{align*}
y^*_t &= y'_t - y'_{t-1} \\
      &= (y_t - y_{t-12}) - (y_{t-1} - y_{t-13}) \\
      &= y_t - y_{t-1} - y_{t-12} + y_{t-13}\: .
\end{align*}
\end{block}\vspace*{10cm}

## Seasonal differencing

When both seasonal and first differences are applied\dots 

* it makes no difference which is done first---the result will be the same.
* If seasonality is strong, we recommend that seasonal differencing be done first because sometimes the resulting series will be stationary and there will be no need for further first difference. 

It is important that if differencing is used, the differences are interpretable.

## Interpretation of differencing

* first differences are the change between **one observation and the next**;
* seasonal differences are the change between **one year to the next**.
 

But taking lag 3 differences for yearly data, for example, results in a model which cannot be sensibly interpreted.

## Unit root tests

\alert{Statistical tests to determine the required order of differencing.}

  1. Augmented Dickey Fuller test: null hypothesis is that the data are non-stationary and non-seasonal.
  2. Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test: null hypothesis is that the data are stationary and non-seasonal.
  3. Other tests available for seasonal data.

## KPSS test
\fontsize{9}{10}\sf

```{r, echo=TRUE}
google_2018 %>%
  features(Close, unitroot_kpss)
```

 

```{r, echo=TRUE}
google_2018 %>%
  features(Close, unitroot_ndiffs)
```

## Automatically selecting differences

STL decomposition: $y_t = T_t+S_t+R_t$

Seasonal strength $F_s = \max\big(0, 1-\frac{\text{Var}(R_t)}{\text{Var}(S_t+R_t)}\big)$

If $F_s > 0.64$, do one seasonal difference.

\fontsize{9}{10}\sf

```{r, echo=TRUE}
usmelec %>% mutate(log_gen = log(Generation)) %>%
  features(log_gen, list(unitroot_nsdiffs, feat_stl))
```

## Automatically selecting differences

\fontsize{9}{10}\sf

```{r, echo=TRUE}
usmelec %>% mutate(log_gen = log(Generation)) %>%
  features(log_gen, unitroot_nsdiffs)
usmelec %>% mutate(d_log_gen = difference(log(Generation), 12)) %>%
  features(d_log_gen, unitroot_ndiffs)
```

## Your turn

For the `tourism` dataset, compute the total number of trips and find an appropriate differencing (after transformation if necessary) to obtain stationary data.

## Backshift notation

A very useful notational device is the backward shift operator, $B$, which is used as follows:
$$
  B y_{t} = y_{t - 1}
$$ 
In other words, $B$, operating on $y_{t}$, has the effect of **shifting the data back one period**.  
Two applications of $B$ to $y_{t}$ **shifts the data back two periods**:
$$
  B(By_{t}) = B^{2}y_{t} = y_{t-2}
$$ 
For monthly data, if we wish to shift attention to "the same month last year", then $B^{12}$ is used, and the notation is $B^{12}y_{t} = y_{t-12}$.

## Backshift notation

The backward shift operator is convenient for describing the process of *differencing*.  
A first difference can be written as
$$
  y'_{t} = y_{t} - y_{t-1} = y_t - By_{t} = (1 - B)y_{t}
$$ 
Note that a first difference is represented by $(1 - B)$. 

Similarly, if second-order differences (i.e., first differences of first differences) have to be computed, then:
\[
  y''_{t} = y_{t} - 2y_{t - 1} + y_{t - 2} = (1 - B)^{2} y_{t}
\]

## Backshift notation

* Second-order difference is denoted $(1- B)^{2}$.
* *Second-order difference* is not the same as a *second difference*, which would be denoted $1- B^{2}$;
* In general, a $d$th-order difference can be written as
$$
  (1 - B)^{d} y_{t}
$$
* A seasonal difference followed by a first difference can be written as
$$
  (1-B)(1-B^m)y_t
$$

## Backshift notation

The "backshift" notation is convenient because the terms can be multiplied together to see the combined effect.
\begin{align*}
  (1-B)(1-B^m)y_t & = (1 - B - B^m + B^{m+1})y_t \\
                  & = y_t-y_{t-1}-y_{t-m}+y_{t-m-1}.
\end{align*} 
For monthly data, $m=12$ and we obtain the same result as earlier.

# Non-seasonal ARIMA models

## Autoregressive models

\begin{block}{Autoregressive (AR) models:}
$$
  y_{t} = c + \phi_{1}y_{t - 1} + \phi_{2}y_{t - 2} + \cdots + \phi_{p}y_{t - p} + \varepsilon_{t},
$$
where $\varepsilon_t$ is white noise. This is a multiple regression with \textbf{lagged values} of $y_t$ as predictors.
\end{block}

```{r arp, echo=FALSE, fig.height=2.5}
set.seed(1)
p1 <- tsibble(idx = seq_len(100), sim = 10 + arima.sim(list(ar = -0.8), n = 100), index = idx) %>%
  autoplot(sim) + ylab("") + ggtitle("AR(1)")
p2 <- tsibble(idx = seq_len(100), sim = 20 + arima.sim(list(ar = c(1.3, -0.7)), n = 100), index = idx) %>%
  autoplot(sim) + ylab("") + ggtitle("AR(2)")
gridExtra::grid.arrange(p1,p2,nrow=1)
```

## AR(1) model

\begin{block}{}
\centerline{$y_{t} = 2 -0.8 y_{t - 1} + \varepsilon_{t}$}
\end{block}
\rightline{$\varepsilon_t\sim N(0,1)$,\quad $T=100$.}

```{r, echo=FALSE}
p1
```

## AR(1) model

\begin{block}{}
  \centerline{$y_{t} = c + \phi_1 y_{t - 1} + \varepsilon_{t}$}
\end{block}

* When $\phi_1=0$, $y_t$ is **equivalent to WN**
* When $\phi_1=1$ and $c=0$, $y_t$ is **equivalent to a RW**
* When $\phi_1=1$ and $c\ne0$, $y_t$ is **equivalent to a RW with drift**
* When $\phi_1<0$, $y_t$ tends to **oscillate between positive and negative values**.

## AR(2) model

\begin{block}{}
  \centerline{$y_t = 8 + 1.3y_{t-1} - 0.7 y_{t-2} + \varepsilon_t$}
\end{block}
\rightline{$\varepsilon_t\sim N(0,1)$, \qquad $T=100$.}

```{r}
p2
```

## Stationarity conditions
\fontsize{14}{15}\sf

We normally restrict autoregressive models to stationary data, and then some constraints on the values of the parameters are required.

\begin{block}{General condition for stationarity}
  Complex roots of $1-\phi_1 z - \phi_2 z^2 - \dots - \phi_pz^p$ lie outside the unit circle on the complex plane.
\end{block} 

* For $p=1$: $-1<\phi_1<1$.
* For $p=2$:\newline $-1<\phi_2<1\qquad \phi_2+\phi_1 < 1 \qquad \phi_2 -\phi_1 < 1$.
* More complicated conditions hold for $p\ge3$.
* Estimation software takes care of this.

## Moving Average (MA) models

\begin{block}{Moving Average (MA) models:}
$$
  y_{t} = c + \varepsilon_t + \theta_{1}\varepsilon_{t - 1} + \theta_{2}\varepsilon_{t - 2} + \cdots + \theta_{q}\varepsilon_{t - q},
$$
where $\varepsilon_t$ is white noise.
This is a multiple regression with \textbf{past \emph{errors}} as predictors. \emph{Don't confuse this with moving average smoothing!}
\end{block}

## Moving Average (MA) models

```{r maq, fig.height=2.5}
set.seed(2)
p1 <- tsibble(idx = seq_len(100), sim = 20 + arima.sim(list(ma = 0.8), n = 100), index = idx) %>%
  autoplot(sim) + ylab("") + ggtitle("MA(1)")
p2 <- tsibble(idx = seq_len(100), sim = arima.sim(list(ma = c(-1, +0.8)), n = 100), index = idx) %>%
  autoplot(sim) + ylab("") + ggtitle("MA(2)")

gridExtra::grid.arrange(p1,p2,nrow=1)
```

## MA(1) model

\begin{block}{}
  \centerline{$y_t = 20 + \varepsilon_t + 0.8 \varepsilon_{t-1}$}
\end{block}
\rightline{$\varepsilon_t\sim N(0,1)$,\quad $T=100$.}

```{r}
p1
```

## MA(2) model

\begin{block}{}
  \centerline{$y_t = \varepsilon_t -\varepsilon_{t-1} + 0.8 \varepsilon_{t-2}$}
\end{block}
\rightline{$\varepsilon_t\sim N(0,1)$,\quad $T=100$.}

```{r}
p2
```

## MA($\infty$) models

It is possible to write any stationary AR($p$) process as an MA($\infty$) process.

**Example: AR(1)**
\begin{align*}
y_t &= \phi_1y_{t-1} + \varepsilon_t\\
&= \phi_1(\phi_1y_{t-2} + \varepsilon_{t-1}) + \varepsilon_t\\
&= \phi_1^2y_{t-2} + \phi_1 \varepsilon_{t-1} + \varepsilon_t\\
&= \phi_1^3y_{t-3} + \phi_1^2\varepsilon_{t-2} + \phi_1 \varepsilon_{t-1} + \varepsilon_t\\
&\dots
\end{align*} 
Provided $-1 < \phi_1 < 1$:
\[
  y_t = \varepsilon_t + \phi_1 \varepsilon_{t-1} + \phi_1^2 \varepsilon_{t-2} + \phi_1^3 \varepsilon_{t-3} + \cdots
\]

## Invertibility

* Any MA($q$) process can be written as an AR($\infty$) process if we impose some constraints on the MA parameters.
* Then the MA model is called "invertible".
* Invertible models have some mathematical properties that make them easier to use in practice.
* Invertibility of an ARIMA model is equivalent to forecastability of an ETS model.

## Invertibility

\begin{block}{General condition for invertibility}
  Complex roots of $1+\theta_1 z + \theta_2 z^2 + \dots + \theta_qz^q$ lie outside the unit circle on the complex plane.
\end{block} 

* For $q=1$: $-1<\theta_1<1$.
* For $q=2$:\newline $-1<\theta_2<1\qquad \theta_2+\theta_1 >-1 \qquad \theta_1 -\theta_2 < 1$.
* More complicated conditions hold for $q\ge3$.
* Estimation software takes care of this.

## Inverting MA($1$)

## ARIMA models

\begin{block}{Autoregressive Moving Average models:}
\begin{align*}
  y_{t} &= c + \phi_{1}y_{t - 1} + \cdots + \phi_{p}y_{t - p} \\
        & \hspace*{2.4cm}\text{} + \theta_{1}\varepsilon_{t - 1} + \cdots + \theta_{q}\varepsilon_{t - q} + \varepsilon_{t}.
\end{align*}
\end{block} 

* Predictors include both **lagged values of $y_t$ and lagged errors.**
* Conditions on coefficients ensure stationarity.
* Conditions on coefficients ensure invertibility.
 

## Autoregressive Integrated Moving Average models
* Combine ARMA model with **differencing**.
* $(1-B)^d y_t$ follows an ARMA model.

## ARIMA models

\alert{Autoregressive Integrated Moving Average models}
\begin{block}{ARIMA($p, d, q$) model}
\begin{tabular}{rl}
AR:& $p =$ order of the autoregressive part\\
I: & $d =$ degree of first differencing involved\\
MA:& $q =$ order of the moving average part.
\end{tabular}
\end{block}

* White noise model: ARIMA(0,0,0)
* Random walk: ARIMA(0,1,0) with no constant
* Random walk with drift: ARIMA(0,1,0) with \rlap{const.}
* AR($p$): ARIMA($p$,0,0)
* MA($q$): ARIMA(0,0,$q$)

## Backshift notation for ARIMA

* ARMA model:\vspace*{-1cm}\newline
\parbox{12cm}{\small\begin{align*}
\hspace*{-1cm}
y_{t} &= c + \phi_{1}By_{t} + \cdots + \phi_pB^py_{t}
           + \varepsilon_{t} + \theta_{1}B\varepsilon_{t} + \cdots + \theta_qB^q\varepsilon_{t} \\
\hspace*{-1cm} \text{or}\quad
      & (1-\phi_1B - \cdots - \phi_p B^p) y_t = c + (1 + \theta_1 B + \cdots + \theta_q B^q)\varepsilon_t
\end{align*}}

* ARIMA(1,1,1) model:\vspace*{-0.5cm}

\[
\begin{array}{c c c c}
(1 - \phi_{1} B) & (1 - B) y_{t} &= &c + (1 + \theta_{1} B) \varepsilon_{t}\\
{\uparrow} & {\uparrow} & &{\uparrow}\\
{\text{AR(1)}} & {\text{First}} &  &{\text{MA(1)}}\\
& {\hbox to 0cm{\hss\text{difference}\hss}}\\
\end{array}
\] 
Written out:
$$
  y_t = c + y_{t-1} + \phi_1 y_{t-1}- \phi_1 y_{t-2} + \theta_1\varepsilon_{t-1} + \varepsilon_t
$$

## Two forms

\fontsize{13}{16}\sf

\begin{block}{Intercept form}
  \centerline{$(1-\phi_1B - \cdots - \phi_p B^p) y_t' = c + (1 + \theta_1 B + \cdots + \theta_q B^q)\varepsilon_t$}
\end{block}

\begin{block}{Mean form}
  \centerline{$(1-\phi_1B - \cdots - \phi_p B^p)(y_t' - \mu) = (1 + \theta_1 B + \cdots + \theta_q B^q)\varepsilon_t$}
\end{block}

 * $y_t' = (1-B)^d y_t$
 * $\mu$ is the mean of $y_t'$.
 * $c = \mu(1-\phi_1 - \cdots - \phi_p )$.
 
 
## Two forms
 * R uses mean form
 * fable uses intercept form

## Parameter redundancy in ARIMA
* Consider the following lines of thought:
\begin{eqnarray*}
x_t & = &  \varepsilon_t\\
x_t-.5x_{t-1}&  = &  \varepsilon_t-.5x_{t-1}\\
x_t -.5x_{t-1}& = &  \varepsilon_t-.5 \varepsilon_{t-1}\\
x_t & = & .5x_{t-1} +  \varepsilon_t-.5 \varepsilon_{t-1}
\end{eqnarray*}
This looks like an ARMA$(1,1)$, but we know that $x_t$ is just white noise. 

* Solution: remove the common factor of $\phi(z)$ and $\theta(z)$.

## Example on redundancy, stationarity, and invertibility.

Consider the process $x_t=0.4x_{t-1}+0.45x_{t-2}+ \varepsilon_t+ \varepsilon_{t-1}+0.25 \varepsilon_{t-2}$ for parameter redundancy, stationarity, and invertibility. In some cases, use 'polyroot' function in R to calculate roots. If complex roots are obtained, use 'Mod' function in R to find the module of complex numbers.

## Australian household expenditure

\fontsize{10}{10}\sf

```{r, eval = FALSE, echo = TRUE}
us_change <- read_csv(
  "https://otexts.com/fpp3/extrafiles/us_change.csv") %>%
  mutate(Time = yearquarter(Time)) %>%
  as_tsibble(index = Time)
```

```{r}
us_change %>% autoplot(Consumption) +
  xlab("Year") +
  ylab("Quarterly percentage change") +
  ggtitle("US consumption")
```

## US personal consumption
\fontsize{9}{10}\sf

```{r, echo=TRUE}
fit <- us_change %>% model(arima = ARIMA(Consumption ~ PDQ(0,0,0)))
report(fit)
```

 \vfill

```{r}
if(!grepl("ARIMA\\(1,0,3\\)", format(fit$arima)))
  warning("Needs fixing")
```

### ARIMA(1,0,3) model:
\centerline{$
  y_t = `r format(tidy(fit)$estimate[5], nsmall=3, digits=3)`
      + `r format(tidy(fit)$estimate[1], nsmall=3, digits=3)`y_{t-1}
      + `r format(tidy(fit)$estimate[2], nsmall=3, digits=3)` \varepsilon_{t-1}
      + `r format(tidy(fit)$estimate[3], nsmall=3, digits=3)` \varepsilon_{t-2}
      + `r format(tidy(fit)$estimate[4], nsmall=3, digits=3)` \varepsilon_{t-2}
      + \varepsilon_{t},
$}
where $\varepsilon_t$ is white noise with a standard deviation of $`r format(sqrt(glance(fit)$sigma2), nsmall=3, digits=3)` = \sqrt{`r format(glance(fit)$sigma2, nsmall=3, digits=3)`}$.

## US personal consumption
\fontsize{12}{15}\sf

```{r, echo=TRUE, fig.height=3.8}
fit %>% forecast(h=10) %>%
  autoplot(slice(us_change, (n()-80):n()))
```

## Understanding ARIMA models
\fontsize{14}{16}\sf

* If $c=0$ and $d=0$, the long-term forecasts will go to zero.
* If $c=0$ and $d=1$, the long-term forecasts will go to a non-zero constant.
* If $c=0$ and $d=2$, the long-term forecasts will follow a straight line.

* If $c\ne0$ and $d=0$, the long-term forecasts will go to the mean of the data.
* If $c\ne0$ and $d=1$, the long-term forecasts will follow a straight line.
* If $c\ne0$ and $d=2$, the long-term forecasts will follow a quadratic trend.

## Understanding ARIMA models
\fontsize{14}{15.5}\sf

### Forecast variance and $d$
  * The higher the value of $d$, the more rapidly the prediction intervals increase in size.
  * For $d=0$, the long-term forecast standard deviation will go to the standard deviation of the historical data.

### Cyclic behaviour
  * For cyclic forecasts, $p\ge2$ and some restrictions on coefficients are required.
  * If $p=2$, we need $\phi_1^2+4\phi_2<0$. Then average cycle of length
\[
  (2\pi)/\left[\text{arc cos}(-\phi_1(1-\phi_2)/(4\phi_2))\right].
\]

# Mean, Variance, ACF, PACF

## Mean of the time series

Mean of the time series, denoted by $\mu_t$ is defined as:

\[\mu_t = E(y_t)\]

* For stationary time series, mean is a constant.

## Mean of ARMA $(p, q)$

* For AR$(p)$, \[\mu_t = \frac{c}{1-\phi_1-\phi_2-\dots-\phi_p}\]
* For MA $(q)$, \[\mu_t = c\]
* For ARMA $(p,q)$, \[\mu_t = \frac{c}{1-\phi_1-\phi_2-\dots-\phi_p}\]

## Variance of the time series

Variance of the time series, denoted by $\sigma^2_t$ is defined as:

\[\sigma^2_t = Var(y_t)\]

* For stationary time series, variance is a constant.

## Variance of ARMA $(p, q)$

* For AR(1), \[\sigma^2_t  = \frac{\sigma_{ \varepsilon}^2}{1-\phi_1^2}\]
where $\sigma_{ \varepsilon}^2$ is the variance of the error term.
* For AR(2), \[\sigma^2_t = \frac{1-\phi_2}{1+\phi_2}*\frac{\sigma_{ \varepsilon}^2}{(1-\phi_2)^2-\phi_1^2}\] 
* For MA $(q)$, \[\sigma^2_t = \sigma_{\varepsilon}^2(1+\theta_1^2+\theta_2^2+\dots+\theta_q^2)\]

## Procedure of finding variance for ARMA $(p, q)$

* Assume that $\phi(B)y_t = c+\theta(B)\varepsilon_t$ is stationary where the roots of $\phi(z)$ are outside the unit cycle. 
* Write \[y_t = \frac{c}{\phi(B)}+\frac{\theta(B)}{\phi(B)} \varepsilon_t = \sum_{j=0}^\infty \psi_j \varepsilon_{t-j}+c^*\]

* Variance of the time series is 
\[\sigma^2_t = \sigma_{\varepsilon}^2 \sum_{j=0}^\infty \psi_j^2\]

## ACF

Autocorrelation of stationary time series, denoted by $\rho_h$ is defined as:

\[\rho_h = cor(y_t,y_{t-h})\]

* For nonstationary time series, ACF cannot be defined.

## ACF 

* For AR $(1)$: \[\rho_h = \phi_1^h.\]
* For MA $(1)$: \[\rho_1 = \theta_1/(1+\theta_1^2)\] and $\rho_h = 0$ for $h \geq 2$.
* For ARMA $(p, q)$: \[\rho_h = \frac{ \sum_{j=0}^\infty \psi_j \psi_{j+h}}{ \sum_{j=0}^\infty \psi_j^2}.\]

## Sample ACF $r_h$

* White noise: For large n,
\[r_h \sim AN(0, \frac{1}{n})\]
for $h=1,2,\dots$. This explains why $\pm 2/\sqrt{n}$ serve as approximate margin of error bounds for $r_h$. Values of $r_h$ outside these bounds would be ``unusual'' under the white noise model assumption.

## Sample ACF

* AR $(1)$: For large $n$, \[r_h \sim AN(\rho_h, \sigma^2_{r_h}) \] where
 \[\sigma^2_{r_h}=\frac{1}{n} \left[ \frac{(1+\phi_1^2)(1-\phi_1^{2h})}{1-\phi_1^2}-2h\phi_1^{2h}\right].\]
* MA $(q)$:  For large $n$, \[r_{q+k} \sim AN \left [0, \frac{1}{n} \left ( 1+2 \sum_{j=1}^q \rho_j^2\right )\right] \] for $k=1,2,\dots$. 

## Sample ACF

* For other cases, results are usually much more complicated. General results can be found in Shumway and Stoffer's book `Time Series Analysis and its application`.

\end{itemize} 

## ARIMA models

* Mean is not a constant.
* Variance is not a constant and may approach infinity as $t\rightarrow \infty$.
* ACF does not exist as ACF is defined only for stationary time series.

## PACF

* \alert{Partial autocorrelations} measure relationship\newline between $y_{t}$ and $y_{t - h}$, when the effects of other time lags --- $1, 2, 3, \dots, h - 1$ --- are removed. 

* For a time series with Normality assumption of the error process, the partial autocorrelation between $y_t$ and $y_{t-h}$ is defined as the conditional correlation between $y_t$ and $y_{t-h}$, conditional on $y_{t-h+1}, \dots , y_{t-1}$, the set of observations that come between the time points $t$ and $t-h$.



## PACF

* The $1^{st}$ order partial autocorrelation is defined to equal the $1^{st}$ order autocorrelation.

* The $2^{nd}$ order (lag) partial autocorrelation is 
\[\frac{\mbox{Covariance}(y_t,y_{t-2}|y_{t-1})}{\mbox{std.Deviation}(y_t|y_{t-1}) \mbox{std.Deviation}(y_{t-2}|y_{t-1})}\]

* The two variances in the denominator will equal each other in a stationary series.

## PACF
\fontsize{11}{13}\sf

* The $3^{rd}$ order (lag) partial autocorrelation is
\[\frac{\mbox{Covariance}(y_t,y_{t-3}|y_{t-1},y_{t-2})}{\mbox{std.Deviation}(y_t|y_{t-1},y_{t-2}) \mbox{std.Deviation}(y_{t-3}|y_{t-1},y_{t-2})}\]

## PACF of AR $(1)$ model

## PACF of MA $(1)$ model

## PACF

* Under the hypothesis that an AR($p$) is correct, the sample partial autocorrelations at lags greater than $p$ are approximately normally distributed with zero means and variances $1/n$ ($n$ is sample size). 
* For $h>p$, $\pm 2/\sqrt{n}$ can be used as critical limits on $\hat{\phi}_{hh}$ to test the null hypothesis that $\phi_{hh}=0$. 

# Estimation and order selection

## Maximum likelihood estimation

Having identified the model order, we need to estimate the parameters $c$, $\phi_1,\dots,\phi_p$, $\theta_1,\dots,\theta_q$. 

* MLE is very similar to least squares estimation obtained by minimizing
$$
  \sum_{t-1}^T e_t^2
$$
* The `ARIMA()` model allows CLS or MLE estimation.
* Non-linear optimization must be used in either case.
* Different software will give different estimates.

## Partial autocorrelations

\fontsize{13}{14}\sf

\alert{Partial autocorrelations} measure relationship\newline between $y_{t}$ and $y_{t - k}$, when the effects of other time lags --- $1, 2, 3, \dots, k - 1$ --- are removed. 
\begin{block}{}
\begin{align*}
\alpha_k &= \text{$k$th partial autocorrelation coefficient}\\
         &= \text{equal to the estimate of $\phi_k$ in regression:}\\
         & \hspace*{0.8cm} y_t = c + \phi_1 y_{t-1} + \phi_2 y_{t-2} + \dots + \phi_k y_{t-k}.
\end{align*}
\end{block} 

* Varying number of terms on RHS gives $\alpha_k$ for different values of $k$.
* There are more efficient ways of calculating $\alpha_k$.

## Example: Mink trapping

```{r}
mink <- as_tsibble(fma::mink)
mink %>% autoplot(value) +
  xlab("Year") +
  ylab("Minks trapped (thousands)") +
  ggtitle("Annual number of minks trapped")
```

## Example: Mink trapping

```{r, echo = TRUE,  fig.height=3}
p1 <- mink %>% ACF(value) %>% autoplot()
p2 <- mink %>% PACF(value) %>% autoplot()
gridExtra::grid.arrange(p1,p2,nrow=1)
```


## Example: Mink trapping

```{r, echo=TRUE}
mink %>% gg_tsdisplay(value)
```


## ACF and PACF interpretation

**AR(1)**
\begin{align*}
\hspace*{1cm}
  \rho_k &= \phi_1^k\qquad\text{for $k=1,2,\dots$};\\
\alpha_1 &= \phi_1 \qquad\alpha_k = 0\qquad\text{for $k=2,3,\dots$}.
\end{align*}

So we have an AR(1) model when

  * autocorrelations exponentially decay
  * there is a single significant partial autocorrelation.

## ACF and PACF interpretation

**AR($p$)**

  * ACF dies out in an exponential or damped sine-wave manner
  * PACF has all zero spikes beyond the $p$th spike

So we have an AR($p$) model when

  * the ACF is exponentially decaying or sinusoidal
  * there is a significant spike at lag $p$ in PACF, but none beyond $p$

## ACF and PACF interpretation

**MA(1)**
\begin{align*}
\hspace*{1cm}\rho_1 &= \theta_1\qquad \rho_k = 0\qquad\text{for $k=2,3,\dots$};\\
\alpha_k &= -(-\theta_1)^k
\end{align*}

So we have an MA(1) model when

 * the PACF is exponentially decaying and
 * there is a single significant spike in ACF

## ACF and PACF interpretation

**MA($q$)**

 * PACF dies out in an exponential or damped sine-wave manner
 * ACF has all zero spikes beyond the $q$th spike

So we have an MA($q$) model when

  * the PACF is exponentially decaying or sinusoidal
  * there is a significant spike at lag $q$ in ACF, but none beyond $q$


## Information criteria

\alert{Akaike's Information Criterion (AIC):}
\centerline{$\text{AIC} = -2 \log(L) + 2(p+q+k+1),$}
where $L$ is the likelihood of the data,\newline
$k=1$ if $c\ne0$ and $k=0$ if $c=0$. \vspace*{0.2cm}

\alert{Corrected AIC:}
\centerline{$\text{AICc} = \text{AIC} + \frac{2(p+q+k+1)(p+q+k+2)}{T-p-q-k-2}.$} \vspace*{0.2cm}

## Information criteria

\alert{Bayesian Information Criterion:}
\centerline{$\text{BIC} = \text{AIC} + [\log(T)-2](p+q+k-1).$}
 \vspace*{-0.2cm}
\begin{block}{}Good models are obtained by minimizing either the AIC, \text{AICc}\ or BIC\@. Our preference is to use the \text{AICc}.\end{block}

# ARIMA modelling in R



## Modelling procedure with `ARIMA`-choose your own model 
\fontsize{12}{13}\sf

1. Plot the data. Identify any unusual observations.
2. If necessary, transform the data (using a Box-Cox transformation) to stabilize the variance.
3. If the data are non-stationary: take first differences of the data until the data are stationary.
4. Examine the ACF/PACF: Is an AR($p$) or MA($q$) model appropriate?
5. Try your chosen model(s), and use the \text{AICc} to search for a better model.

## Modelling procedure with `ARIMA`-choose your own model 
\fontsize{12}{13}\sf

6. Check the residuals from your chosen model by plotting the ACF of the residuals, and doing a portmanteau test of the residuals. If they do not look like white noise, try a modified model.
7. Once the residuals look like white noise, calculate forecasts.

## Automated ARIMA() 

\begin{block}{A non-seasonal ARIMA process}
\[
\phi(B)(1-B)^dy_{t} = c + \theta(B)\varepsilon_t
\]
Need to select appropriate orders: \alert{$p,q, d$}
\end{block}

\alert{Hyndman and Khandakar (JSS, 2008) algorithm:}

  * Select no.\ differences \alert{$d$} and \alert{$D$} via KPSS test and seasonal strength measure.
  * Select \alert{$p,q$} by minimising AICc.
  * Use stepwise search to traverse model space.

## How does automated ARIMA() work?
\fontsize{12}{12}\sf

\begin{block}{}
\centerline{$\text{AICc} = -2 \log(L) + 2(p+q+k+1)\left[1 +
\frac{(p+q+k+2)}{T-p-q-k-2}\right].$}
where $L$ is the maximised likelihood fitted to the \textit{differenced} data,
$k=1$ if $c\neq 0$ and $k=0$ otherwise.
\end{block} 

Step1:
: Select current model (with smallest AICc) from:\newline
ARIMA$(2,d,2)$\newline
ARIMA$(0,d,0)$\newline
ARIMA$(1,d,0)$\newline
ARIMA$(0,d,1)$
 \vspace*{-0.1cm}

Step 2:
: Consider variations of current model:

    * vary one of $p,q,$ from current model by $\pm1$;
    * $p,q$ both vary from current model by $\pm1$;
    * Include/exclude $c$ from current model.

  Model with lowest AICc becomes current model.

\alert{Repeat Step 2 until no lower AICc can be found.}




## Automatic modelling procedure with `ARIMA`
\fontsize{12}{13}\sf

1. Plot the data. Identify any unusual observations.
2. If necessary, transform the data (using a Box-Cox transformation) to stabilize the variance.

\vspace*{1.15cm}

3. Use `ARIMA` to automatically select a model.

\vspace*{1.15cm}

6. Check the residuals from your chosen model by plotting the ACF of the residuals, and doing a portmanteau test of the residuals. If they do not look like white noise, try a modified model.
7. Once the residuals look like white noise, calculate forecasts.

## Modelling procedure

\centerline{\includegraphics[height=8.cm]{Figure-8-10}}

## \large Seasonally adjusted electrical equipment
\fontsize{11}{12}\sf

```{r ee1, fig.height=3.3, echo=TRUE}
elec_equip <- as_tsibble(fpp2::elecequip)
elec_dcmp <- elec_equip %>%
  model(STL(value ~ season(window="periodic"))) %>%
  components() %>% select(-.model) %>% as_tsibble()
elec_dcmp %>%autoplot(season_adjust)
```

## \large Seasonally adjusted electrical equipment

1. Time plot shows sudden changes, particularly big drop in 2008/2009 due to global economic environment. Otherwise nothing unusual and no need for data adjustments.
2. No evidence of changing variance, so no Box-Cox transformation.
3. Data are clearly non-stationary, so we take first differences.

## \large Seasonally adjusted electrical equipment
\fontsize{11}{12}\sf

```{r ee2, echo=TRUE, fig.height=4}
elec_dcmp %>%
  gg_tsdisplay(difference(season_adjust), plot_type='partial')
```


## \large Seasonally adjusted electrical equipment

4. PACF is suggestive of AR(3). So initial candidate model is ARIMA(3,1,0). No other obvious candidates.
5. Fit ARIMA(3,1,0) model along with variations: ARIMA(4,1,0), ARIMA(2,1,0), ARIMA(3,1,1), etc. ARIMA(3,1,1) has smallest \text{AICc} value.

## \large Seasonally adjusted electrical equipment
\fontsize{11}{12}\sf

```{r, echo=TRUE}
fit <- elec_dcmp %>%
  model(
    arima = ARIMA(season_adjust ~ pdq(3,1,1) + PDQ(0,0,0))
  )
report(fit)
```

## \large Seasonally adjusted electrical equipment

6. ACF plot of residuals from ARIMA(3,1,1) model look like white noise.

\fontsize{11}{11}\sf

```{r, echo=FALSE, fig.height=3}
fit %>% gg_tsresiduals()
```

## \large Seasonally adjusted electrical equipment
\fontsize{11}{12}\sf

```{r}
augment(fit) %>%
  features(.resid, ljung_box, lag = 24, dof = 4)

```

## \large Seasonally adjusted electrical equipment
\fontsize{11}{12}\sf

```{r, echo=TRUE}
fit %>% forecast() %>% autoplot(elec_dcmp)
```


## Roots checking
\fontsize{11}{12}\sf
* The three red dots correspond to the roots of the polynomials $\phi(B)$ (left) and $\theta(B)$ (right)of the ARIMA(3,1,1) model.

```{r, echo=TRUE,fig.height=3}
gg_arma(fit)
```

## Roots checking

* They are all inside the unit circle, as we would expect because R ensures the fitted model is both stationary and invertible.
* Any roots close to the unit circle may be numerically unstable, and the corresponding model will not be good for forecasting.
* The `ARIMA()` function will never return a model with inverse roots outside the unit circle.


# Forecasting

## Point forecasts

1. Rearrange ARIMA equation so $y_t$ is on LHS.
2. Rewrite equation by replacing $t$ by $T+h$.
3. On RHS, replace future observations by their forecasts, future errors by zero, and past errors by corresponding residuals.

Start with $h=1$. Repeat for $h=2,3,\dots$.

## Point forecasts
\fontsize{14}{14}\sf

\alert{ARIMA(3,1,1) forecasts: Step 1}
\begin{block}{}
\centerline{$(1-\phi_1B -\phi_2B^2-\phi_3B^3)(1-B) y_t = (1+\theta_1B)\varepsilon_{t},$}
\end{block}
 \vspace*{-0.4cm}
\begin{align*}
\left[1-(1+\phi_1)B +(\phi_1-\phi_2)B^2 + (\phi_2-\phi_3)B^3 +\phi_3B^4\right] y_t\\ = (1+\theta_1B)\varepsilon_{t},
\end{align*} \vspace*{-0.4cm}
\begin{align*}
y_t - (1+\phi_1)y_{t-1} +(\phi_1-\phi_2)y_{t-2} + (\phi_2-\phi_3)y_{t-3}\\ \mbox{}+\phi_3y_{t-4} = \varepsilon_t+\theta_1\varepsilon_{t-1}.
\end{align*} \vspace*{-0.4cm}
\begin{align*}
y_t = (1+\phi_1)y_{t-1} -(\phi_1-\phi_2)y_{t-2} - (\phi_2-\phi_3)y_{t-3}\\\mbox{} -\phi_3y_{t-4} + \varepsilon_t+\theta_1\varepsilon_{t-1}.
\end{align*}

## Point forecasts (h=1)
\fontsize{14}{14}\sf

\begin{block}{}
\begin{align*}
y_t = (1+\phi_1)y_{t-1} -(\phi_1-\phi_2)y_{t-2} - (\phi_2-\phi_3)y_{t-3}\\\mbox{} -\phi_3y_{t-4} + \varepsilon_t+\theta_1\varepsilon_{t-1}.
\end{align*}
\end{block} 
\alert{ARIMA(3,1,1) forecasts: Step 2}
\begin{align*}
y_{T+1} = (1+\phi_1)y_{T} -(\phi_1-\phi_2)y_{T-1} - (\phi_2-\phi_3)y_{T-2}\\\mbox{} -\phi_3y_{T-3} + \varepsilon_{T+1}+\theta_1\varepsilon_{T}.
\end{align*} 
\alert{ARIMA(3,1,1) forecasts: Step 3}
\begin{align*}
\hat{y}_{T+1|T} = (1+\phi_1)y_{T} -(\phi_1-\phi_2)y_{T-1} - (\phi_2-\phi_3)y_{T-2}\\\mbox{} -\phi_3y_{T-3} + \theta_1 e_{T}.
\end{align*}

## Point forecasts (h=2)
\fontsize{14}{14}\sf

\begin{block}{}
\begin{align*}
y_t = (1+\phi_1)y_{t-1} -(\phi_1-\phi_2)y_{t-2} - (\phi_2-\phi_3)y_{t-3}\\\mbox{} -\phi_3y_{t-4} + \varepsilon_t+\theta_1\varepsilon_{t-1}.
\end{align*}
\end{block} 
\alert{ARIMA(3,1,1) forecasts: Step 2}
\begin{align*}
y_{T+2} = (1+\phi_1)y_{T+1} -(\phi_1-\phi_2)y_{T} - (\phi_2-\phi_3)y_{T-1}\\\mbox{} -\phi_3y_{T-2} + \varepsilon_{T+2}+\theta_1\varepsilon_{T+1}.
\end{align*} 
\alert{ARIMA(3,1,1) forecasts: Step 3}
\begin{align*}
\hat{y}_{T+2|T} = (1+\phi_1)\hat{y}_{T+1|T} -(\phi_1-\phi_2)y_{T} - (\phi_2-\phi_3)y_{T-1}\\\mbox{} -\phi_3y_{T-2}.
\end{align*}


## Prediction intervals

\begin{block}{95\% prediction interval}
$$\hat{y}_{T+h|T} \pm 1.96\sqrt{v_{T+h|T}}$$
where $v_{T+h|T}$ is estimated forecast variance.
\end{block}

* Multi-step prediction intervals for ARIMA(0,0,$q$):
\centerline{$\displaystyle y_t = \varepsilon_t + \sum_{i=1}^q \theta_i \varepsilon_{t-i}.$}
\centerline{$\displaystyle
v_{T|T+h} = \hat{\sigma}^2 \left[ 1 + \sum_{i=1}^{h-1} \theta_i^2\right], \qquad\text{for~} h=2,3,\dots.$}

 

* AR(1): Rewrite as MA($\infty$) and use above result.
* Other models beyond scope of this subject.

## Prediction intervals

* Prediction intervals **increase in size with forecast horizon**.
* Prediction intervals can be difficult to calculate by hand
* Calculations assume residuals are **uncorrelated** and **normally distributed**.
* Prediction intervals tend to be too narrow.
    * the uncertainty in the parameter estimates has not been accounted for.
    * the ARIMA model assumes historical patterns will not change during the forecast period.
    * the ARIMA model assumes uncorrelated future \rlap{errors}

## Your turn

For the United States GDP data (from `global_economy`):

 * if necessary, find a suitable Box-Cox transformation for the data;
 * fit a suitable ARIMA model to the transformed data;
 * check the residual diagnostics;
 * produce forecasts of your fitted model. Do the forecasts look reasonable?

# Seasonal ARIMA models

## Seasonal ARIMA models

| ARIMA | $~\underbrace{(p, d, q)}$ | $\underbrace{(P, D, Q)_{m}}$ |
| ----: | :-----------------------: | :--------------------------: |
|       | ${\uparrow}$              | ${\uparrow}$                 |
|       | Non-seasonal part         | Seasonal part of             |
|       | of the model              | of the model                 |

where $m =$ number of observations per year.

## Seasonal ARIMA models
\fontsize{11}{12}\sf

E.g., ARIMA$(1, 1, 1)(1, 1, 1)_{4}$ model (without constant) 
$$(1 - \phi_{1}B)(1 - \Phi_{1}B^{4}) (1 - B) (1 - B^{4})y_{t} ~= ~
(1 + \theta_{1}B) (1 + \Theta_{1}B^{4})\varepsilon_{t}.
$$ 

\setlength{\unitlength}{1mm}
\begin{footnotesize}
\begin{picture}(100,25)(-5,0)
\thinlines
{\put(5,22){\vector(0,1){6}}}
{\put(22,10){\vector(0,1){18}}}
{\put(38,22){\vector(0,1){6}}}
{\put(52,10){\vector(0,1){18}}}
{\put(77,22){\vector(0,1){6}}}
{\put(95,10){\vector(0,1){18}}}
{\put(-10,17){$\left(\begin{array}{@{}c@{}} \text{Non-seasonal} \\ \text{AR(1)}
                    \end{array}\right)$}}
{\put(12,5){$\left(\begin{array}{@{}c@{}} \text{Seasonal} \\ \text{AR(1)}
                    \end{array}\right)$}}
{\put(25,17){$\left(\begin{array}{@{}c@{}} \text{Non-seasonal} \\ \text{difference}
                    \end{array}\right)$}}
{\put(40,5){$\left(\begin{array}{@{}c@{}} \text{Seasonal} \\ \text{difference}
                    \end{array}\right)$}}
{\put(65,17){$\left(\begin{array}{@{}c@{}} \text{Non-seasonal} \\ \text{MA(1)}
                    \end{array}\right)$}}
{\put(85,5){$\left(\begin{array}{@{}c@{}} \text{Seasonal} \\ \text{MA(1)}
                    \end{array}\right)$}}
\end{picture}
\end{footnotesize}

\vspace*{10cm}

## Seasonal ARIMA models
\fontsize{11}{12}\sf
E.g., ARIMA$(1, 1, 1)(1, 1, 1)_{4}$ model (without constant)
$$(1 - \phi_{1}B)(1 - \Phi_{1}B^{4}) (1 - B) (1 - B^{4})y_{t} ~= ~
(1 + \theta_{1}B) (1 + \Theta_{1}B^{4})\varepsilon_{t}.
$$

All the factors can be multiplied out and the general model
written as follows:
\begin{align*}
y_{t} &= (1 + \phi_{1})y_{t - 1} - \phi_1y_{t-2} + (1 + \Phi_{1})y_{t - 4}\\
&\text{}
 - (1 + \phi_{1} + \Phi_{1} + \phi_{1}\Phi_{1})y_{t - 5}
 + (\phi_{1} + \phi_{1} \Phi_{1}) y_{t - 6} \\
& \text{} - \Phi_{1} y_{t - 8} + (\Phi_{1} + \phi_{1} \Phi_{1}) y_{t - 9}
  - \phi_{1} \Phi_{1} y_{t - 10}\\
  &\text{}
  + \varepsilon_{t} + \theta_{1}\varepsilon_{t - 1} + \Theta_{1}\varepsilon_{t - 4} + \theta_{1}\Theta_{1}\varepsilon_{t - 5}.
\end{align*}
\vspace*{10cm}

## Common ARIMA models

The US Census Bureau uses the following models most often:\vspace*{0.5cm}

\begin{tabular}{|ll|}
\hline
ARIMA(0,1,1)(0,1,1)$_m$& with log transformation\\
ARIMA(0,1,2)(0,1,1)$_m$& with log transformation\\
ARIMA(2,1,0)(0,1,1)$_m$& with log transformation\\
ARIMA(0,2,2)(0,1,1)$_m$& with log transformation\\
ARIMA(2,1,2)(0,1,1)$_m$& with no transformation\\
\hline
\end{tabular}

## Seasonal ARIMA models
\fontsize{11}{12}\sf

The seasonal part of an AR or MA model will be seen in the seasonal lags of
the PACF and ACF.

\alert{ARIMA(0,0,0)(0,0,1)$_{12}$ will show:}

  * a spike at lag 12 in the ACF but no other significant spikes.
  * The PACF will show exponential decay in the seasonal lags; that is, at lags 12, 24, 36, \dots.

\alert{ARIMA(0,0,0)(1,0,0)$_{12}$ will show:}

  * exponential decay in the seasonal lags of the ACF
  * a single significant spike at lag 12 in the PACF.

## European quarterly retail trade

```{r, echo=TRUE, fig.height=3.6}
eu_retail %>% autoplot(value) +
  xlab("Year") + ylab("Retail index")
```

## European quarterly retail trade
\fontsize{11}{12}\sf

```{r, echo=TRUE, fig.height=4}
eu_retail %>% gg_tsdisplay(
  value %>% difference(4), plot_type='partial')
```

## European quarterly retail trade
\fontsize{11}{12}\sf

```{r, echo=TRUE, fig.height=3.8}
eu_retail %>% gg_tsdisplay(
  value %>% difference(4) %>% difference(1),plot_type='partial')
```

## European quarterly retail trade

  * $d=1$ and $D=1$ seems necessary.
  * Significant spike at lag 1 in ACF suggests non-seasonal MA(1) component.
  * Significant spike at lag 4 in ACF suggests seasonal MA(1) component.
  * Initial candidate model: ARIMA(0,1,1)(0,1,1)$_4$.
  * We could also have started with ARIMA(1,1,0)(1,1,0)$_4$.

## European quarterly retail trade
\fontsize{11}{12}

```{r, echo=TRUE, fig.height=3.5}
fit <- eu_retail %>%
  model(arima = ARIMA(value ~ pdq(0,1,1) + PDQ(0,1,1)))
augment(fit) %>% gg_tsdisplay(.resid, plot_type = "hist")
```

## European quarterly retail trade
\fontsize{11}{12}\sf

```{r, echo = TRUE}
augment(fit) %>%
  features(.resid, ljung_box, lag = 8, dof = 2)
```


## European quarterly retail trade

```{r, echo=FALSE}
aicc <- eu_retail %>%
  model(
    mdl_1 = ARIMA(value ~ pdq(0,1,1) + PDQ(0,1,1)),
    mdl_2 = ARIMA(value ~ pdq(0,1,2) + PDQ(0,1,1)),
    mdl_3 = ARIMA(value ~ pdq(0,1,3) + PDQ(0,1,1)),
    mdl_4 = ARIMA(value ~ pdq(0,1,4) + PDQ(0,1,1))
  ) %>%
  glance %>%
  pull(AICc)
```

  * ACF and PACF of residuals show significant spikes at lag 2, and maybe lag 3.
  * AICc of ARIMA(0,1,1)(0,1,1)$_4$ model is `r round(aicc[1],2)`
  * AICc of ARIMA(0,1,2)(0,1,1)$_4$ model is `r round(aicc[2],2)`.
  * AICc of ARIMA(0,1,3)(0,1,1)$_4$ model is `r round(aicc[3],2)`.
  * AICc of ARIMA(0,1,4)(0,1,1)$_4$ model is `r round(aicc[4],2)`.

## European quarterly retail trade

\fontsize{11}{12}\sf

```{r, echo=TRUE}
fit <- eu_retail %>%
  model(
    arima013011 = ARIMA(value ~ pdq(0,1,3) + PDQ(0,1,1))
  )
report(fit)
```

## European quarterly retail trade
\fontsize{13}{15}\sf

```{r, echo=TRUE, fig.height=4}
augment(fit) %>%
  gg_tsdisplay(.resid, plot_type = "hist")
```

## European quarterly retail trade
\fontsize{13}{15}\sf

```{r, echo=TRUE}
augment(fit) %>%
  features(.resid, ljung_box, lag = 8, dof = 4)
```

## European quarterly retail trade

```{r, echo=TRUE, fig.height=3.8}
fit %>% forecast(h = "3 years") %>%
  autoplot(eu_retail)
```

## European quarterly retail trade
\fontsize{11}{14}\sf

```{r, echo=TRUE}
eu_retail %>% model(ARIMA(value)) %>% report()
```

## European quarterly retail trade
\fontsize{11}{14}\sf

```{r eu_retailtryharder, echo=TRUE}
eu_retail %>% model(ARIMA(value, stepwise = FALSE,
  approximation = FALSE)) %>% report()
```


## Cortecosteroid drug sales

```{r h02}
h02 %>%
  mutate(log(Cost)) %>%
  gather() %>%
  ggplot(aes(x = Month, y = value)) +
  geom_line() +
  facet_grid(key ~ ., scales = "free_y") +
  xlab("Year") + ylab("") +
  ggtitle("Cortecosteroid drug scripts (H02)")
```

## Cortecosteroid drug sales
```{r h02b}
h02 %>% gg_tsdisplay(difference(log(Cost),12), lag_max = 36,plot_type='partial')
```

## Cortecosteroid drug sales

  * Choose $D=1$ and $d=0$.
  * Spikes in PACF at lags 12 and 24 suggest seasonal AR(2) term.
  * Spikes in PACF suggests possible non-seasonal AR(3) term.
  * Initial candidate model: ARIMA(3,0,0)(2,1,0)$_{12}$.

## Cortecosteroid drug sales

```{r h02aicc, echo=FALSE}
models <- list(
  c(3,0,0,2,1,0),
  c(3,0,1,2,1,0),
  c(3,0,2,2,1,0),
  c(3,0,1,1,1,0),
  c(3,0,1,0,1,1),
  c(3,0,1,0,1,2),
  c(3,0,1,1,1,1)
)
model_defs <- map(models, ~ ARIMA(log(Cost) ~ 0 + pdq(!!.[1], !!.[2], !!.[3]) + PDQ(!!.[4], !!.[5], !!.[6])))
model_defs <- set_names(model_defs, map_chr(models,
  ~ sprintf("ARIMA(%i,%i,%i)(%i,%i,%i)[12]", .[1], .[2], .[3], .[4], .[5], .[6])))

fit <- h02 %>%
  model(!!!model_defs)

fit %>%
  glance %>%
  arrange(AICc) %>%
  select(.model, AICc) %>%
  knitr::kable(digits=2, row.names=FALSE, align='cc', booktabs=TRUE)
```


## Cortecosteroid drug sales
\fontsize{10}{12}\sf

```{r arimah02, echo=TRUE}
fit <- h02 %>%
  model(best = ARIMA(log(Cost) ~ 0 + pdq(3,0,1) + PDQ(0,1,2)))
report(fit)
```

## Cortecosteroid drug sales
\fontsize{10}{12}\sf

```{r h02res, echo=TRUE, fig.height=4}
augment(fit) %>%
  gg_tsdisplay(.resid, lag_max=36, plot_type = "hist")
```

## Cortecosteroid drug sales
\fontsize{11}{15}\sf

```{r h02resb, echo = TRUE, fig.height=4}
augment(fit) %>%
  features(.resid, ljung_box, lag = 36, dof = 6)
```

## Cortecosteroid drug sales
\fontsize{11}{13}\sf

```{r h02auto, echo=TRUE, fig.height=3.6}
fit <- h02 %>% model(auto = ARIMA(log(Cost)))
report(fit)
```

## Cortecosteroid drug sales
\fontsize{10}{12}\sf

```{r, echo=TRUE, fig.height=4}
augment(fit) %>%
  gg_tsdisplay(.resid, lag_max = 36, plot_type = "hist")
```

## Cortecosteroid drug sales
\fontsize{11}{15}\sf

```{r, echo = TRUE}
augment(fit) %>%
  features(.resid, ljung_box, lag = 36, dof = 5)
```


## Cortecosteroid drug sales
\fontsize{9}{9}\sf

```{r echo=TRUE, fig.height=3.6}
fit <- h02 %>%
  model(best = ARIMA(log(Cost), stepwise = FALSE,
                 approximation = FALSE,
                 order_constraint = p + q + P + Q <= 9))
report(fit)
```

## Cortecosteroid drug sales
\fontsize{10}{12}\sf

```{r, echo=TRUE, fig.height=4}
augment(fit) %>%
  gg_tsdisplay(.resid, lag_max = 36, plot_type = "hist")
```

## Cortecosteroid drug sales
\fontsize{11}{15}\sf

```{r, echo = TRUE}
augment(fit) %>%
  features(.resid, ljung_box, lag = 36, dof = 9)
```


## Cortecosteroid drug sales
\fontsize{10}{12}\sf

Training data: July 1991 to June 2006

Test data: July 2006--June 2008

```{r, echo = TRUE}
fit <- h02 %>%
  filter_index(~ "2006 Jun") %>%
  model(
    ARIMA(log(Cost) ~ pdq(3, 0, 0) + PDQ(2, 1, 0)),
    ARIMA(log(Cost) ~ pdq(3, 0, 1) + PDQ(2, 1, 0)),
    ARIMA(log(Cost) ~ pdq(3, 0, 2) + PDQ(2, 1, 0)),
    ARIMA(log(Cost) ~ pdq(3, 0, 1) + PDQ(1, 1, 0))
    # ... #
  )

fit %>%
  forecast(h = "2 years") %>%
  accuracy(h02 %>% filter_index("2006 Jul" ~ .))
```

## Cortecosteroid drug sales
\fontsize{10}{12}\sf

```{r, cache=TRUE}
models <- list(
  c(3,0,0,2,1,0),
  c(3,0,1,2,1,0),
  c(3,0,2,2,1,0),
  c(3,0,1,1,1,0),
  c(3,0,1,0,1,1),
  c(3,0,1,0,1,2),
  c(3,0,1,1,1,1),
  c(3,0,3,0,1,1),
  c(3,0,2,0,1,1),
  c(2,1,3,0,1,1),
  c(2,1,4,0,1,1),
  c(2,1,5,0,1,1),
  c(4,1,1,2,1,2))

model_defs <- map(models, ~ ARIMA(log(Cost) ~ 0 + pdq(!!.[1], !!.[2], !!.[3]) + PDQ(!!.[4], !!.[5], !!.[6])))
model_defs <- set_names(model_defs, map_chr(models,
  ~ sprintf("ARIMA(%i,%i,%i)(%i,%i,%i)[12]", .[1], .[2], .[3], .[4], .[5], .[6])))

fit <- h02 %>%
  filter_index(~ "2006 Jun") %>%
  model(!!!model_defs)
```

## Cortecosteroid drug sales
\fontsize{10}{12}\sf

```{r h02-rmse2, cache=TRUE}
fit %>%
  forecast(h = "2 years") %>%
  accuracy(h02 %>% filter_index("2006 Jul" ~ .)) %>%
  arrange(RMSE) %>%
  select(.model, RMSE) %>%
  knitr::kable(digits=4)
```

## Cortecosteroid drug sales

  * Models with lowest AICc values tend to give slightly better results than the other models.
  * AICc comparisons must have the same orders of differencing. But RMSE test set comparisons can involve any models.
  * Use the best model available, even if it does not pass all tests.

## Cortecosteroid drug sales
\fontsize{11}{14}\sf

```{r h02f, echo=TRUE, fig.height=3}
fit <- h02 %>%
  model(ARIMA(Cost ~ 0 + pdq(3,0,1) + PDQ(1,1,1)))
fit %>% forecast %>% autoplot(h02) +
  ylab("H02 Expenditure ($AUD)") + xlab("Year")
```

# ARIMA vs ETS

## ARIMA vs ETS

  * Myth that ARIMA models are more general than exponential smoothing.

  * Linear exponential smoothing models all special cases of ARIMA models.

  * Non-linear exponential smoothing models have no equivalent ARIMA counterparts.

  * Many ARIMA models have no exponential smoothing counterparts.

  * ETS models all non-stationary. Models with seasonality or non-damped trend (or both) have two unit roots; all other models have one unit \rlap{root.}

\vspace*{10cm}

## Equivalences

\fontsize{12}{14}\sf

|**ETS model**  | **ARIMA model**             | **Parameters**                       |
| :------------ | :-------------------------- | :----------------------------------- |
| ETS(A,N,N)    | ARIMA(0,1,1)                | $\theta_1 = \alpha-1$                |
| ETS(A,A,N)    | ARIMA(0,2,2)                | $\theta_1 = \alpha+\beta-2$          |
|               |                             | $\theta_2 = 1-\alpha$                |
| ETS(A,$A_d$,N)    | ARIMA(1,1,2)                | $\phi_1=\phi$                        |
|               |                             | $\theta_1 = \alpha+\phi\beta-1-\phi$ |
|               |                             | $\theta_2 = (1-\alpha)\phi$          |
| ETS(A,N,A)    | ARIMA(0,0,$m$)(0,1,0)$_m$   |                                      |
| ETS(A,A,A)    | ARIMA(0,1,$m+1$)(0,1,0)$_m$ |                                      |
| ETS(A,$A_d$,A)    | ARIMA(1,0,$m+1$)(0,1,0)$_m$ |                                      |


## Your turn
\fontsize{12}{14}\sf

For the `fma::condmilk` series:

  * Do the data need transforming? If so, find a suitable transformation.
  * Are the data stationary? If not, find an appropriate differencing which yields stationary data.
  * Identify a couple of ARIMA models that might be useful in describing the time series.
  * Which of your models is the best according to their AIC values?
  
## Your turn
\fontsize{12}{14}\sf

  * Estimate the parameters of your best model and do diagnostic testing on the residuals. Do the residuals resemble white noise? If not, try to find another ARIMA model which fits better.
  * Forecast the next 24 months of data using your preferred model.
  * Compare the forecasts obtained using `ets()`.
