---
title: "STAT481/581: Introduction to Time Series Analysis"
author: "Exponential smoothing"
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)

library(tidyverse)
library(fable)
library(tsibble)
library(feasts)
library(lubridate)
library(tsibbledata)

options(digits=4, width=55)
```

# Exponential smoothing

## Historical perspective

 * Developed in the 1950s and 1960s as methods (algorithms) to produce point forecasts.
 * Combine a "level", "trend" (slope) and "seasonal" component to describe a time series.
 * The rate of change of the components are controlled by "smoothing parameters":\newline $\alpha$, $\beta$ and $\gamma$ respectively.
  * Need to choose best values for the smoothing parameters (and initial states).
  * Equivalent ETS state space models developed in the 1990s and 2000s.

## Big idea: control the rate of change (smoothing)

$\alpha$ controls the flexibility of the **level**

* If $\alpha = 0$, the level never updates (mean)
* If $\alpha = 1$, the level updates completely (naive)

$\beta$ controls the flexibility of the **trend**

* If $\beta = 0$, the trend is linear (regression trend)
* If $\beta = 1$, the trend updates every observation

## Big idea: control the rate of change (smoothing)

$\gamma$ controls the flexibility of the **seasonality**

* If $\gamma = 0$, the seasonality is fixed (seasonal means)
* If $\gamma = 1$, the seasonality updates completely (seasonal naive)

## A model for levels, trends, and seasonalities

We want a model that captures the level ($\ell_t$), trend ($b_t$) and seasonality ($s_t$).

How do we combine these elements?

 

### Additively?

$y_t = \ell_{t-1} + b_{t-1} + s_{t-m} + \varepsilon_t$

 

### Multiplicatively?

$y_t = \ell_{t-1}b_{t-1}s_{t-m}(1 + \varepsilon_t)$

 

### Perhaps a mix of both?

$y_t = (\ell_{t-1} + b_{t-1}) s_{t-m} + \varepsilon_t$

## ETS models

\begin{block}{}
\hspace*{-0.25cm}\begin{tabular}{l@{}p{2.3cm}@{}c@{}l}
\structure{General n\rlap{otation}}
    &       & ~E T S~  & ~:\hspace*{0.3cm}\textbf{E}xponen\textbf{T}ial \textbf{S}moothing               \\ [-0.2cm]
    & \hfill{$\nearrow$\hspace*{-0.1cm}}        & {$\uparrow$} & {\hspace*{-0.2cm}$\nwarrow$} \\
    & \hfill{\textbf{E}rror\hspace*{0.2cm}} & {\textbf{T}rend}      & {\hspace*{0.2cm}\textbf{S}eason}
\end{tabular}
\end{block}

\alert{\textbf{E}rror:} Additive (`"A"`) or multiplicative (`"M"`)
 

\alert{\textbf{T}rend:} None (`"N"`), additive (`"A"`), multiplicative (`"M"`), or damped (`"Ad"` or `"Md"`).
 

\alert{\textbf{S}easonality:} None (`"N"`), additive (`"A"`) or multiplicative (`"M"`)

# Simple exponential smoothing

## Simple methods
\fontsize{11}{12}\sf

Time series $y_1,y_2,\dots,y_T$.

\begin{block}{Random walk forecasts}
  \centerline{$\pred{y}{T+h}{T} = y_T$}
\end{block} 

\begin{block}{Average forecasts}
  \centerline{$\displaystyle\pred{y}{T+h}{T} = \frac1T\sum_{t=1}^T y_t$}
\end{block} \vspace*{-0.2cm}

* Want something in between these methods.
* Most recent data should have more weight.

<!-- * Simple exponential smoothing uses a weighted moving average with weights that decrease exponentially. -->

## Simple Exponential Smoothing

\begin{block}{Forecast equation}
$\pred{y}{T+1}{T} = \alpha y_T + \alpha(1-\alpha) y_{T-1} + \alpha(1-\alpha)^2 y_{T-2}+ \cdots$
\end{block}
where $0 \le \alpha \le1$. \vspace*{0.2cm}

\small\begin{tabular}{lllll}
\toprule
& \multicolumn{4}{l}{Weights assigned to observations for:}\\
Observation  &   $\alpha = 0.2$   &   $\alpha = 0.4$  &   $\alpha = 0.6$  & $\alpha = 0.8$ \\
\midrule
$y_{T}$      & 0.2         & 0.4          & 0.6         & 0.8\\
$y_{T-1}$    & 0.16        & 0.24         & 0.24        & 0.16\\
$y_{T-2}$    & 0.128       & 0.144        & 0.096       & 0.032\\
$y_{T-3}$    & 0.1024      & 0.0864       & 0.0384      & 0.0064\\
$y_{T-4}$    & $(0.2)(0.8)^4$  & $(0.4)(0.6)^4$   & $(0.6)(0.4)^4$  & $(0.8)(0.2)^4$\\
$y_{T-5}$    & $(0.2)(0.8)^5$  & $(0.4)(0.6)^5$   & $(0.6)(0.4)^5$  & $(0.8)(0.2)^5$\\
\bottomrule
\end{tabular}

## Simple Exponential Smoothing

```{r alpha-static2, fig.height=5, fig.width=8}
algeria_economy <- tsibbledata::global_economy %>%
  filter(Country == "Algeria")
fit <- algeria_economy %>%
  model(
    ETS(Exports ~ error("A") + trend("N") + season("N"))
  )
alpha_static <- map_dfr(list(0, as.numeric(tidy(fit)[1,4]), 1), function(alpha){
  fit <- algeria_economy %>%
    model(ETS(Exports ~ error("A") + trend("N", alpha = alpha, alpha_range = c(-0.01,1),
          beta_range = c(-1,1)) + season("N", gamma_range = c(-1,1)), bounds = "admissible"))
  fit %>%
    augment() %>%
    mutate(alpha = tidy(fit)$estimate[tidy(fit)$term == "alpha"]) %>%
    as_tibble()
}) %>%
  mutate(alpha = factor(format(alpha)))
algeria_economy %>%
  ggplot(aes(x = Year, y = Exports)) +
  geom_line() +
  geom_line(aes(y = .fitted, colour = alpha), data = alpha_static) +
  ylab("Exports (% of GDP)") +
  ggtitle("Algerian exports of goods and services: level")
```

## Simple Exponential Smoothing
### Methods

  * Algorithms that return point forecasts.

### Models

  * Generate same point forecasts but can also generate forecast distributions.
  * A stochastic (or random) data generating process that can generate an entire forecast distribution.
  * Allow for "proper" model selection.


## ETS(A,N,N)

State space form:
\begin{block}{}\vspace*{-0.4cm}
\begin{align*}
\text{Measurement equation}&& y_t &= \ell_{t-1} + \varepsilon_t\\
\text{State equation}&& \ell_t&=\ell_{t-1}+\alpha \varepsilon_t
\end{align*}
\end{block}
where $\varepsilon_t\sim\text{NID}(0,\sigma^2)$.

  * "innovations" or "single source of error" because equations have the same error process, $\varepsilon_t$.
  * Measurement equation: relationship between observations and states.
  * Transition/state equation(s): evolution of the state(s) through time.

## ETS(A,N,N)

## ETS(A,N,N)
\fontsize{11}{12}\sf

\begin{block}{Component form}\vspace*{-0.4cm}
\begin{align*}
\text{Forecast equation}&&\pred{y}{t+h}{t} &= \ell_{t}\\
\text{Smoothing equation}&&\ell_{t} &= \alpha y_{t} + (1 - \alpha)\ell_{t-1}
\end{align*}
\end{block}\vspace*{-0.2cm}

* $\ell_t$ is the level (or the smoothed value) of the series at time t.
* $\pred{y}{t+1}{t} = \alpha y_t + (1-\alpha) \pred{y}{t}{t-1}$\newline
  Iterate to get exponentially weighted moving average form.

\begin{block}{Weighted average form}
$\displaystyle\pred{y}{T+1}{T}=\sum_{j=0}^{T-1} \alpha(1-\alpha)^j y_{T-j}+(1-\alpha)^T \ell_{0}$
\end{block}

## Optimising smoothing parameters

  * Need to choose best values for $\alpha$ and $\ell_0$.
  * Similarly to regression, choose optimal parameters by minimising SSE:
$$
  \text{SSE}=\sum_{t=1}^T(y_t - \pred{y}{t}{t-1})^2.
$$
  * Unlike regression there is no closed form solution --- use numerical optimization.


## ETS(A,N,N): Specifying the model

\fontsize{13}{15}\sf

```{r ann-spec, echo = TRUE, results = "hide"}
ETS(y ~ error("A") + trend("N") + season("N"))
```

\fontsize{11}{12}\sf

By default, an optimal value for $\alpha$ and $\ell_0$ is used.

$\alpha$ can be chosen manually in `trend()`.

```{r alpha-spec, echo = TRUE, eval = FALSE}
trend("N", alpha = 0.5)
trend("N", alpha_range = c(0.2, 0.8))
```

## Algeria economy data 

\fontsize{11}{12}\sf

```{r, echo=TRUE}
algeria_economy <- tsibbledata::global_economy %>%
  filter(Country == "Algeria")
fit <- algeria_economy %>%
  model(ETS(Exports ~ error("A") + trend("N") 
              +season("N"),opt_crit = "mse"))
tidy(fit)$term
tidy(fit)$estimate
```

## Algeria economy data 

  * For Algerian Exports example:
    - $\hat\alpha = `r sprintf("%4.4f",tidy(fit)[1,4])`$
    - $\hat\ell_0 = `r sprintf("%4.2f",tidy(fit)[2,4])`$
  * Here $l_0$



## Example: Algerian Exports

\fontsize{9}{10}\sf

```{r ses-fit, echo=TRUE, cache=TRUE}
report(fit)
```

## Example: Algerian Exports

\fontsize{10}{11}\sf\vspace*{-0.2cm}

```{r ses-cmp, echo = TRUE}
components(fit) %>%
  left_join(fitted(fit), by = c("Country", ".model", "Year"))
```


## Simple Exponential Smoothing is the initial state. 
\fontsize{11}{12}\sf

```{r alpha-static,echo=TRUE, fig.height=3}
fc <- fit %>% forecast(h = 5)
fc %>%autoplot(algeria_economy) + geom_line(aes(y = .fitted, 
colour = "Fitted"), data = augment(fit)) + 
  ylab("Exports (% of GDP)") + xlab("Year")
```

## ETS(M,N,N)

SES with multiplicative errors.
\begin{block}{}\vspace*{-0.4cm}
\begin{align*}
\text{Measurement equation}&& y_t &= \ell_{t-1}(1 + \varepsilon_t)\\
\text{State equation}&& \ell_t&=\ell_{t-1}(1+\alpha \varepsilon_t)
\end{align*}
\end{block}


  * Models with additive and multiplicative errors with the same parameters generate the same point forecasts but different prediction intervals.

## ETS(M,N,N)


# Models with trend

## ETS(A,A,N)

Holt's linear method with additive errors.

  * $b_t$: slope
  * State space form:
  \begin{align*}
      y_t&=\ell_{t-1}+b_{t-1}+\varepsilon_t\\
      \ell_t&=\ell_{t-1}+b_{t-1}+\alpha \varepsilon_t\\
      b_t&=b_{t-1}+\alpha\beta^* \varepsilon_t
  \end{align*}
  * For simplicity, set $\beta=\alpha \beta^*$.


## Holt's linear trend

\begin{block}{Component form}\vspace*{-.4cm}
\begin{align*}
\text{Forecast }&& \pred{y}{t+h}{t} &= \ell_{t} + hb_{t} \\
\text{Level }&& \ell_{t} &= \alpha y_{t} + (1 - \alpha)(\ell_{t-1} + b_{t-1})\\
\text{Trend }&& b_{t} &= \beta^*(\ell_{t} - \ell_{t-1}) + (1 -\beta^*)b_{t-1},
\end{align*}
\end{block}
 \vspace*{-0.2cm}

  * Two smoothing parameters $\alpha$ and $\beta^*$ ($0\le\alpha,\beta^*\le1$).
  
## Holt's linear trend
   * $\ell_t$ level: weighted average between $y_t$ and one-step ahead forecast for time $t$, $(\ell_{t-1} + b_{t-1}=\pred{y}{t}{t-1})$.
  * $b_t$ slope: weighted average of $(\ell_{t} - \ell_{t-1})$ and $b_{t-1}$, current and previous estimate of slope.
  * Choose $\alpha, \beta^*, \ell_0, b_0$ to minimise SSE.

## Holt's linear trend


## Exponential smoothing: trend/slope

```{r beta-static, fig.height=5, fig.width=8}
apple_stock <- gafa_stock %>%
  filter(Symbol == "AAPL") %>%
  mutate(trading_day = row_number()) %>%
  update_tsibble(index = trading_day, regular = TRUE)

beta_static <- map_dfr(list(0, 0.001),
                     function(beta){
    fit <- apple_stock %>%
      model(ETS(Close ~ error("A") + trend("A", alpha = 0.00001, alpha_range = c(-.1,1),
            beta = beta, beta_range = c(-0.1,0.005)) + season("N", gamma_range = c(-1,1)), bounds = "admissible"))
    fit %>%
      augment() %>%
      mutate(beta = tidy(fit)$estimate[tidy(fit)$term == "beta"]) %>%
      as_tibble()
  }) %>%
  mutate(beta = factor(format(beta))) %>%
  left_join(select(apple_stock, trading_day, Date), by = "trading_day")
apple_stock %>%
  ggplot(aes(x = Date, y = Close)) +
  geom_line() +
  geom_line(aes(y = .fitted, colour = beta), data = beta_static) +
  ylab("Closing price ($USD)") +
  ggtitle("Apple closing stock price: trend")
```

## ETS(A,A,N): Specifying the model

\fontsize{13}{15}\sf

```{r aan-spec, echo = TRUE, results = "hide"}
ETS(y ~ error("A") + trend("A") + season("N"))
```

\fontsize{11}{12}\sf

By default, optimal values for $\beta$ and $b_0$ are used.

$\beta$ can be chosen manually in `trend()`.

```{r beta-spec, echo = TRUE, eval = FALSE}
trend("A", beta = 0.004)
trend("A", beta_range = c(0, 0.1))
```


## AUS pupulation data
\fontsize{11}{12}\sf
```{r auseco,echo=TRUE, fig.height=5}
aus_economy <- global_economy %>%filter(Code == "AUS") %>%
  mutate(Pop = Population / 1e6)
fit <- aus_economy %>%
  model(AAN = ETS(Pop ~ error("A") + trend("A") + 
                    season("N")))
fc <- fit %>% forecast(h = 10)
```

## AUS pupulation data
\fontsize{11}{12}\sf
```{r auseco2,echo=TRUE}
report(fit)
```

## AUS pupulation data

\fontsize{10}{11}\sf\vspace*{-0.2cm}

```{r auseco3, echo = TRUE}
components(fit) %>%
  left_join(fitted(fit), by = c("Country", ".model", "Year"))
```


## AUS pupulation data
\fontsize{11}{12}\sf
```{r auseco4,echo=TRUE}
fc
```


## AUS pupulation data

\fontsize{12}{12}\sf

```{r holt-fc, echo=TRUE, cache=TRUE, dependson='holt-fit'}
fit %>%forecast(h = 10) %>%autoplot(aus_economy) +
  geom_line(aes(y=.fitted,colour="Fitted"),data = 
augment(fit))+ylab("Population") + xlab("Year")
```


## ETS(M,A,N)

Holt's linear method with multiplicative errors.

  * Following a similar approach as above, the innovations state space model underlying Holt's linear method with multiplicative errors is specified as\vspace*{-0.4cm}
  \begin{align*}
      y_t&=(\ell_{t-1}+b_{t-1})(1+\varepsilon_t)\\
      \ell_t&=(\ell_{t-1}+b_{t-1})(1+\alpha \varepsilon_t)\\
      b_t&=b_{t-1}+\beta(\ell_{t-1}+b_{t-1}) \varepsilon_t
  \end{align*}
  where again  $\beta=\alpha \beta^*$ and $\varepsilon_t \sim \text{NID}(0,\sigma^2)$.

## ETS(M,A,N) component form

## Damped trend method
\begin{block}{Component form}\vspace*{-0.4cm}
\begin{align*}
\pred{y}{t+h}{t} &= \ell_{t} + (\phi+\phi^2 + \dots + \phi^{h})b_{t} \\
\ell_{t} &= \alpha y_{t} + (1 - \alpha)(\ell_{t-1} + \phi b_{t-1})\\
b_{t} &= \beta^*(\ell_{t} - \ell_{t-1}) + (1 -\beta^*)\phi b_{t-1}.
\end{align*}
\end{block}
 

  * Damping parameter $0<\phi<1$.
  * If $\phi=1$, identical to Holt's linear trend.
  * As $h\rightarrow\infty$, $\pred{y}{T+h}{T}\rightarrow \ell_T+\phi b_T/(1-\phi)$.
  * Short-run forecasts trended, long-run forecasts constant.


## Your turn
\large

 * Write down the state space model for ETS(A,Ad,N)
 * Derive the forecast equation and the smoothing equations.
 
## Your turn

## Example: Australian population
\fontsize{10}{11}\sf

```{r, echo=TRUE, fig.height=3.6}
aus_economy %>%
  model(holt = ETS(Pop ~ error("A") + trend("Ad") + season("N"))) %>%
  forecast(h = 10) %>%
  autoplot(aus_economy)
```

## Example: Australian population
\fontsize{9}{12}\sf

```{r, echo=TRUE}
fit <- aus_economy %>%
  filter(Year <= 2010) %>%
  model(
    ses = ETS(Pop ~ error("A") + trend("N") + season("N")),
    holt = ETS(Pop ~ error("A") + trend("A") + season("N")),
    damped = ETS(Pop ~ error("A") + trend("Ad") + season("N"))
  )
```

```{r, echo = TRUE, results = 'hide'}
tidy(fit)
accuracy(fit)
```

## Example: Australian population
\fontsize{13}{15}\sf

```{r echo=FALSE}
fit_terms <- tidy(fit) %>%
  spread(.model, estimate) %>%
  mutate(term = factor(term, levels = c("alpha", "beta", "phi", "l", "b"), labels = c("$\\alpha$", "$\\beta^*$", "$\\phi$", "$\\ell_0$", "$b_0$"))) %>%
  arrange(term)

fit_accuracy <- accuracy(fit) %>%
  bind_rows(
    forecast(fit, h = 9) %>%
      accuracy(aus_economy)
  ) %>%
  gather(term, estimate, -Country, -.model, -.type) %>%
  spread(.model, estimate) %>%
  filter(term == "RMSE" | .type == "Test" & term %in% c("RMSE", "MAE", "MAPE", "MASE")) %>%
  arrange(desc(.type), desc(term)) %>%
  unite("term", .type, term, sep = " ")

bind_rows(fit_terms, fit_accuracy) %>%
  select(term, ses, holt, damped) %>%
  rename(SES = ses, `Linear trend` = holt, `Damped trend` = damped) %>%
  mutate_if(is.numeric, ~ ifelse(is.na(.), "", formatC(., format = "f", 2))) %>%
  knitr::kable(booktabs = TRUE, align='r')
```

## Your turn

`fma::eggs` contains the price of a dozen eggs in the United States from 1900–1993

 1. Use simple exponential smoothing (SES) and Holt’s method (with and without damping) to forecast “future” data.

     [Hint: use h=100 so you can clearly see the differences between the options when plotting the forecasts.]
 2. Which method gives the best training RMSE?
 3. Are these RMSE values comparable?
 4. Do the residuals from the best fitting method look like white noise?

# Models with seasonality

## ETS(A,A,A)

* Holt and Winters extended Holt's method to capture seasonality.

* Holt-Winters additive method in state space form:

\begin{block}{}\vspace*{-0.4cm}
\begin{align*}
\text{Observation equation}&& y_t&=\ell_{t-1}+b_{t-1}+s_{t-m} + \varepsilon_t\\
\text{State equations}&& \ell_t&=\ell_{t-1}+b_{t-1}+\alpha \varepsilon_t\\
&&        b_t&=b_{t-1}+\beta \varepsilon_t \\
&& s_t &= s_{t-m} + \gamma\varepsilon_t
\end{align*}
\end{block}

* $k$ is integer part of $(h-1)/m$.
* $s_t$: seasonal state.

## Holt-Winters additive method component form
\fontsize{13}{15}\sf

\begin{block}{Component form}\vspace*{-0.4cm}
\begin{align*}
\pred{y}{t+h}{t} &= \ell_{t} + hb _{t} + s_{t+h-m(k+1)} \\
\ell_{t} &= \alpha(y_{t} - s_{t-m}) + (1 - \alpha)(\ell_{t-1} + b_{t-1})\\
b_{t} &= \beta^*(\ell_{t} - \ell_{t-1}) + (1 - \beta^*)b_{t-1}\\
s_{t} &= \gamma (y_{t}-\ell_{t-1}-b_{t-1}) + (1-\gamma)s_{t-m}
\end{align*}
\end{block}\fontsize{12}{14}\sf

  * $k=$ integer part of $(h-1)/m$. Ensures estimates from the final year are used for forecasting.
  * Parameters:&nbsp; $0\le \alpha\le 1$,&nbsp; $0\le \beta^*\le 1$,&nbsp; $0\le \gamma\le 1-\alpha$&nbsp;  and $m=$  period of seasonality (e.g. $m=4$ for quarterly data).

## Holt-Winters additive method component form



## Holt-Winters additive method

  * Seasonal component is usually expressed as
        $s_{t} = \gamma^* (y_{t}-\ell_{t})+ (1-\gamma^*)s_{t-m}.$
  * Substitute in for $\ell_t$:
        $s_{t} = \gamma^*(1-\alpha) (y_{t}-\ell_{t-1}-b_{t-1})+ [1-\gamma^*(1-\alpha)]s_{t-m}$
  * We set $\gamma=\gamma^*(1-\alpha)$.
  * The usual parameter restriction is $0\le\gamma^*\le1$, which translates to $0\le\gamma\le(1-\alpha)$.

## Exponential smoothing: seasonality

```{r gamma-static, fig.height=5, fig.width=8}
j07 <- PBS %>%
  filter(ATC2 == "J07") %>%
  summarise(Cost = sum(Cost))
gamma_static <- map_dfr(list(0, NULL, 1), function(gamma){
  fit <- j07 %>%
    model(ETS(Cost ~ error("A") + trend("N", alpha = 0.001, alpha_range = c(-1,1),
              beta_range = c(-1,1)) + season("A", gamma = gamma, gamma_range = c(-1,1)), bounds = "admissible"))
  fit %>%
    augment() %>%
    mutate(gamma = tidy(fit)$estimate[tidy(fit)$term == "gamma"]) %>%
    as_tibble()
}) %>%
  mutate(gamma = factor(format(gamma)))
j07 %>%
  ggplot(aes(x = Month, y = Cost)) +
  geom_line() +
  geom_line(aes(y = .fitted, colour = gamma), data = gamma_static) +
  ylab("Cost of scripts ($AUD)") +
  ggtitle("Medicare Australia cost of vaccine scripts: seasonality")
```


## Your turn
\large

 * Write down the state space model and component form for ETS(A,N,A).
 
## Your turn


## ETS(M,A,M)

* Holt-Winters multiplicative method with multiplicative errors for when seasonal variations are changing proportional to the level of the series.

\begin{block}{}\vspace*{-0.4cm}
\begin{align*}
\text{Observation equation}&& y_t&= (\ell_{t-1}+b_{t-1})s_{t-m}(1 + \varepsilon_t)\\
\text{State equations}&& \ell_t&=(\ell_{t-1}+b_{t-1})(1+\alpha \varepsilon_t)\\
&&        b_t&=b_{t-1}(1+\beta \varepsilon_t) \\
&&s_t &= s_{t-m}(1 + \gamma\varepsilon_t)
\end{align*}
\end{block}

* $k$ is integer part of $(h-1)/m$.

## Holt-Winters multiplicative method
\fontsize{13}{14}\sf\vspace*{-0.1cm}


\begin{block}{Component form}\vspace*{-0.3cm}
    \begin{align*}
        \pred{y}{t+h}{t} &= (\ell_{t} + hb_{t})s_{t+h-m(k+1)} \\
        \ell_{t} &= \alpha \frac{y_{t}}{s_{t-m}} + (1 - \alpha)(\ell_{t-1} + b_{t-1})\\
        b_{t} &= \beta^*(\ell_{t}-\ell_{t-1}) + (1 - \beta^*)b_{t-1}        \\
        s_{t} &= \gamma \frac{y_{t}}{(\ell_{t-1} + b_{t-1})} + (1 - \gamma)s_{t-m}
    \end{align*}
\end{block}\vspace*{-0.2cm}\fontsize{11}{12}\sf

  * $k$ is integer part of $(h-1)/m$.
  * With additive method $s_t$ is in absolute terms:\newline within each year $\sum_i s_i \approx 0$.
  * With multiplicative method $s_t$ is in relative terms:\newline within each year $\sum_i s_i \approx m$.

## Holt-Winters multiplicative method component form

## Example: Australian holiday tourism

\fontsize{9}{10}\sf

```{r 7-HW, echo=TRUE}
aus_holidays <- tourism %>% filter(Purpose == "Holiday") %>%
  summarise(Trips = sum(Trips))
fit <- aus_holidays %>%
  model(
    additive = ETS(Trips ~ error("A") + trend("A") + season("A")),
    multiplicative = ETS(Trips ~ error("M") + trend("A") + season("M"))
  )
fc <- fit %>% forecast()
```

```{r, fig.height=2.6}
fc %>%
  autoplot(aus_holidays, level = NULL) + xlab("Year") +
  ylab("Overnight trips (thousands)") +
  scale_color_brewer(type = "qual", palette = "Dark2")
```

## Estimated components

```{r, echo = TRUE, results = 'hide'}
components(fit)
```

```{r fig-7-LevelTrendSeas}
components(fit) %>%
  gather("state", "value", -.model, -Quarter, factor_key = TRUE) %>%
  group_by(.model) %>%
  group_split() %>%
  map(
    ~ ggplot(., aes(x = Quarter, y = value)) +
      geom_line() +
      facet_grid(state ~ ., scales = "free") +
      xlab("Year") + ylab(NULL) +
      ggtitle(str_to_title(unique(.$.model)) %>% paste("states"))
  ) %>%
  invoke(gridExtra::grid.arrange, ., ncol = 2)
```

## Holt-Winters damped method
Often the single most accurate forecasting method for seasonal data:
\begin{block}{}\vspace*{-0.4cm}
\begin{align*}
\pred{y}{t+h}{t} &= [\ell_{t} + (\phi+\phi^2 + \dots + \phi^{h})b_{t}]s_{t+h-m(k+1)} \\
\ell_{t} &= \alpha(y_{t} / s_{t-m}) + (1 - \alpha)(\ell_{t-1} + \phi b_{t-1})\\
b_{t} &= \beta^*(\ell_{t} - \ell_{t-1}) + (1 - \beta^*)\phi b_{t-1}       \\
s_{t} &= \gamma \frac{y_{t}}{(\ell_{t-1} + \phi b_{t-1})} + (1 - \gamma)s_{t-m}
\end{align*}
\end{block}

## Your turn

Apply Holt-Winters’ multiplicative method to the Gas data from `aus_production`.

 1. Why is multiplicative seasonality necessary here?
 2. Experiment with making the trend damped.
 3. Check that the residuals from the best method look like white noise.

# Innovations state space models

## Exponential smoothing methods
\fontsize{12}{14}\sf

\begin{block}{}
\begin{tabular}{ll|ccc}
& &\multicolumn{3}{c}{\bf Seasonal Component} \\
\multicolumn{2}{c|}{\bf Trend}& N & A & M\\
\multicolumn{2}{c|}{\bf Component}  & (None)    & (Additive)  & (Multiplicative)\\
\cline{3-5} &&&&\\[-0.4cm]
N & (None) & (N,N) & (N,A) & (N,M)\\
&&&&\\[-0.4cm]
A & (Additive) & (A,N) & (A,A) & (A,M)\\
&&&&\\[-0.4cm]
A\damped & (Additive damped) & (A\damped,N) & (A\damped,A) & (A\damped,M)
\end{tabular}
\end{block}\fontsize{12}{14}\sf

\begin{tabular}{lp{9.7cm}}
\textcolor[rgb]{0.90,0.,0.00}{(N,N)}:        &Simple exponential smoothing\\
\textcolor[rgb]{0.90,0.,0.00}{(A,N)}:        &Holt's linear method\\
\textcolor[rgb]{0.90,0.,0.00}{(A\damped,N)}: &Additive damped trend method\\
\textcolor[rgb]{0.90,0.,0.00}{(A,A)}:~~ &Additive Holt-Winters' method\\
\textcolor[rgb]{0.90,0.,0.00}{(A,M)}: &Multiplicative Holt-Winters' method\\
\textcolor[rgb]{0.90,0.,0.00}{(A\damped,M)}: &Damped multiplicative Holt-Winters' method
\end{tabular}

\begin{block}{}\fontsize{12}{14}\sf
There are also multiplicative trend methods (not recommended).
\end{block}

## ETS models
\fontsize{11}{12}\sf

\begin{block}{}
\begin{tabular}{ll|ccc}
  \multicolumn{2}{l}{\alert{\bf Additive Error}} &        \multicolumn{3}{c}{\bf Seasonal Component}         \\
          \multicolumn{2}{c|}{\bf Trend}         &         N         &         A         &         M         \\
        \multicolumn{2}{c|}{\bf Component}       &     ~(None)~      &    (Additive)     & (Multiplicative)  \\ \cline{3-5}
           &                                     &                   &                   &  \\[-0.3cm]
  N        & (None)                              &       A,N,N       &       A,N,A       &    A,N,M     \\
           &                                     &                   &                   &  \\[-0.3cm]
  A        & (Additive)                          &       A,A,N       &       A,A,A       &    A,A,M     \\
           &                                     &                   &                   &  \\[-0.3cm]
  A\damped & (Additive damped)                   &   A,A\damped,N    &   A,A\damped,A    & A,A\damped,M
\end{tabular}
\end{block}

\begin{block}{}
\begin{tabular}{ll|ccc}
  \multicolumn{2}{l}{\alert{\bf Multiplicative Error}} &     \multicolumn{3}{c}{\bf Seasonal Component}      \\
             \multicolumn{2}{c|}{\bf Trend}            &      N       &         A         &        M         \\
           \multicolumn{2}{c|}{\bf Component}          &   ~(None)~   &    (Additive)     & (Multiplicative) \\ \cline{3-5}
           &                                           &              &                   &  \\[-0.3cm]
  N        & (None)                                    &    M,N,N     &       M,N,A       &      M,N,M       \\
           &                                           &              &                   &  \\[-0.3cm]
  A        & (Additive)                                &    M,A,N     &       M,A,A       &      M,A,M       \\
           &                                           &              &                   &  \\[-0.3cm]
  A\damped & (Additive damped)                         & M,A\damped,N &   M,A\damped,A    &   M,A\damped,M
\end{tabular}
\end{block}

## Additive error models

\placefig{0}{1.5}{width=12.8cm,trim=0 120 0 0,clip=true}{fig_7_ets_add.pdf}

## Multiplicative error models

\placefig{0}{1.5}{width=12.8cm,trim=0 120 0 0,clip=true}{fig_7_ets_multi.pdf}

## Estimating ETS models

  * Smoothing parameters $\alpha$, $\beta$, $\gamma$ and $\phi$, and the initial states $\ell_0$, $b_0$, $s_0,s_{-1},\dots,s_{-m+1}$ are estimated by maximising the "likelihood" = the probability of the data arising from the specified model.
  * For models with additive errors equivalent to minimising SSE.
  * For models with multiplicative errors, \textbf{not} equivalent to minimising SSE.

## Innovations state space models
\fontsize{12}{14}\sf

Let $\bm{x}_t = (\ell_t, b_t, s_t, s_{t-1}, \dots, s_{t-m+1})$ and
$\varepsilon_t\stackrel{\mbox{\scriptsize iid}}{\sim}
\mbox{N}(0,\sigma^2)$.
\begin{block}{}
\begin{tabular}{lcl}
$y_t$ &=& $\underbrace{h(\bm{x}_{t-1})} +
\underbrace{k(\bm{x}_{t-1})\varepsilon_t}$\\
&& \hspace*{0.5cm}$\mu_t$ \hspace*{1.45cm} $e_t$ \\[0.2cm]
$\bm{x}_t$ &=& $f(\bm{x}_{t-1}) +
g(\bm{x}_{t-1})\varepsilon_t$\\
\end{tabular}
\end{block}

Additive errors
: \mbox{}\vspace*{-0.5cm}\newline
  $k(x)=1$.\qquad $y_t = \mu_{t} + \varepsilon_t$.

Multiplicative errors
: \mbox{}\vspace*{-0.5cm}\newline
  $k(\bm{x}_{t-1}) = \mu_{t}$.\qquad $y_t = \mu_{t}(1 + \varepsilon_t)$.\newline
  $\varepsilon_t = (y_t - \mu_t)/\mu_t$ is relative error.

## Innovations state space models

\structure{Estimation}\vspace*{0.5cm}

\begin{block}{}
\begin{align*}
L^*(\bm\theta,\bm{x}_0) &= n\log\!\bigg(\sum_{t=1}^n \varepsilon^2_t/k^2(\bm{x}_{t-1})\!\bigg) + 2\sum_{t=1}^n \log|k(\bm{x}_{t-1})|\\
&= -2\log(\text{Likelihood}) + \mbox{constant}
\end{align*}
\end{block}

* Estimate parameters $\bm\theta = (\alpha,\beta,\gamma,\phi)$ and
initial states $\bm{x}_0 = (\ell_0,b_0,s_0,s_{-1},\dots,s_{-m+1})$ by
minimizing $L^*$.

## Parameter restrictions
\fontsize{12}{14}\sf

### *Usual* region

  * Traditional restrictions in the methods $0< \alpha,\beta^*,\gamma^*,\phi<1$\newline (equations interpreted as weighted averages).
  * In models we set $\beta=\alpha\beta^*$ and $\gamma=(1-\alpha)\gamma^*$.
  * Therefore $0< \alpha <1$, &nbsp;&nbsp; $0 < \beta < \alpha$ &nbsp;&nbsp; and $0< \gamma < 1-\alpha$.
  * $0.8<\phi<0.98$ --- to prevent numerical difficulties.
  
## Parameter restrictions
\fontsize{12}{14}\sf

### *Admissible* region

  * To prevent observations in the distant past having a continuing effect on current forecasts.
  * Usually (but not always) less restrictive than the \textit{traditional} region.
  * For example for ETS(A,N,N): \newline \textit{traditional} $0< \alpha <1$ --- \textit{admissible} is $0< \alpha <2$.

## Model selection
\fontsize{13}{15}\sf

\begin{block}{Akaike's Information Criterion}
\[
\text{AIC} = -2\log(\text{L}) + 2k
\]
\end{block}\vspace*{-0.2cm}
where $L$ is the likelihood and $k$ is the number of parameters initial states estimated in the model. 

\begin{block}{Corrected AIC}
\[
\text{AIC}_{\text{c}} = \text{AIC} + \frac{2(k+1)(k+2)}{T-k}
\]
\end{block}
which is the AIC corrected (for small sample bias).
 
\begin{block}{Bayesian Information Criterion}
\[
\text{BIC} = \text{AIC} + k(\log(T)-2).
\]
\end{block}

## AIC and cross-validation

\Large

\begin{alertblock}{}
Minimizing the AIC assuming Gaussian residuals is asymptotically equivalent to minimizing one-step time series cross validation MSE.
\end{alertblock}

## Automatic forecasting

**From Hyndman et al.\ (IJF, 2002):**

* Apply each model that is appropriate to the data.
Optimize parameters and initial values using MLE (or some other
criterion).
* Select best method using AICc:
* Produce forecasts using best method.
* Obtain forecast intervals using underlying state space model.

Method performed very well in M3 competition.

## Some unstable models

* Some of the combinations of (Error, Trend, Seasonal) can lead to numerical difficulties; see equations with division by a state.
* These are: ETS(A,N,M), ETS(A,A,M), ETS(A,A\damped,M).
* Models with multiplicative errors are useful for strictly positive data, but are not numerically stable with data containing zeros or negative values. In that case only the six fully additive models will be applied.

## Exponential smoothing models
\fontsize{11}{12}\sf

\begin{block}{}
\begin{tabular}{ll|ccc}
  \multicolumn{2}{l}{\alert{\bf Additive Error}} &        \multicolumn{3}{c}{\bf Seasonal Component}         \\
          \multicolumn{2}{c|}{\bf Trend}         &         N         &         A         &         M         \\
        \multicolumn{2}{c|}{\bf Component}       &     ~(None)~      &    (Additive)     & (Multiplicative)  \\ \cline{3-5}
           &                                     &                   &                   &  \\[-0.3cm]
  N        & (None)                              &       A,N,N       &       A,N,A       &    \st{A,N,M}     \\
           &                                     &                   &                   &  \\[-0.3cm]
  A        & (Additive)                          &       A,A,N       &       A,A,A       &    \st{A,A,M}     \\
           &                                     &                   &                   &  \\[-0.3cm]
  A\damped & (Additive damped)                   &   A,A\damped,N    &   A,A\damped,A    & \st{A,A\damped,M}
\end{tabular}
\end{block}

\begin{block}{}
\begin{tabular}{ll|ccc}
  \multicolumn{2}{l}{\alert{\bf Multiplicative Error}} &     \multicolumn{3}{c}{\bf Seasonal Component}      \\
             \multicolumn{2}{c|}{\bf Trend}            &      N       &         A         &        M         \\
           \multicolumn{2}{c|}{\bf Component}          &   ~(None)~   &    (Additive)     & (Multiplicative) \\ \cline{3-5}
           &                                           &              &                   &  \\[-0.3cm]
  N        & (None)                                    &    M,N,N     &       M,N,A       &      M,N,M       \\
           &                                           &              &                   &  \\[-0.3cm]
  A        & (Additive)                                &    M,A,N     &       M,A,A       &      M,A,M       \\
           &                                           &              &                   &  \\[-0.3cm]
  A\damped & (Additive damped)                         & M,A\damped,N &   M,A\damped,A    &   M,A\damped,M
\end{tabular}
\end{block}

## Example: Australian holiday tourism
\fontsize{10}{10}\sf

```{r, echo=TRUE}
fit <- aus_holidays %>% model(ETS(Trips))
report(fit)
```

## Example: Australian holiday tourism

Model selected: ETS(M,N,M)
\begin{align*}
y_{t} &= \ell_{t-1}s_{t-m}(1 + \varepsilon_t)\\
\ell_t &= \ell_{t-1}(1 + \alpha \varepsilon_t)\\
s_t &=  s_{t-m}(1+ \gamma \varepsilon_t).
\end{align*}

$\hat\alpha=`r format(tidy(fit)$estimate[1],nsmall=4,digits=4)`$,
and $\hat\gamma=`r format(tidy(fit)$estimate[2],nsmall=3,digits=3, scientific=FALSE)`$.

## Example: Australian holiday tourism

```{r MAMstates, fig.height=3.5,fig.width=6, echo=TRUE, results = "hide"}
components(fit)
```

```{r r MAMstates, fig.height=3.5,fig.width=6}
components(fit) %>%
  autoplot() +
  ggtitle("ETS(M,N,M) components")
```

## Residuals
\fontsize{16}{18}\sf

### Response residuals
$$\hat{e}_t = y_t - \hat{y}_{t|t-1}$$

### Innovation residuals
Additive error model:
$$\hat\varepsilon_t = y_t - \hat{y}_{t|t-1}$$

Multiplicative error model:
$$\hat\varepsilon_t = \frac{y_t - \hat{y}_{t|t-1}}{\hat{y}_{t|t-1}}$$

## Example: Australian holiday tourism
\fontsize{9.5}{12}\sf

```{r, echo = TRUE, results = "hide"}
residuals(fit)
residuals(fit, type = "response")
```

```{r, echo=FALSE}
residuals(fit) %>%
  mutate(Type = "Innovation residuals") %>%
  bind_rows(residuals(fit, type = "response") %>% mutate(Type = "Response residuals")) %>%
  ggplot(aes(x = Quarter, y = .resid)) +
  geom_line() +
  facet_grid(Type ~ ., scales = "free_y") +
  ylab(NULL)
```

# Forecasting with exponential smoothing

## Forecasting with ETS models

\structure{Point forecasts:} iterate the equations for $t=T+1,T+2,\dots,T+h$ and set all $\varepsilon_t=0$ for $t>T$. 

* Not the same as $\text{E}(y_{t+h} | \bm{x}_t)$ unless trend and seasonality are both additive.
* Point forecasts for ETS(A,\*,\*) are identical to ETS(M,\*,\*) if the parameters are the same.

## Example: ETS(A,A,N)

\vspace*{-1.3cm}

\begin{align*}
y_{T+1} &= \ell_T + b_T  + \varepsilon_{T+1}\\
\hat{y}_{T+1|T} & = \ell_{T}+b_{T}\\
y_{T+2}         & = \ell_{T+1} + b_{T+1} + \varepsilon_{T+2}\\
                & =
                      (\ell_T + b_T + \alpha\varepsilon_{T+1}) +
                      (b_T + \beta \varepsilon_{T+1}) +
                      \varepsilon_{T+2} \\
\hat{y}_{T+2|T} &= \ell_{T}+2b_{T}
\end{align*}
etc.

## Example: ETS(M,A,N)
\fontsize{13}{16}\sf

\vspace*{-1.3cm}

\begin{align*}
y_{T+1} &= (\ell_T + b_T )(1+ \varepsilon_{T+1})\\
\hat{y}_{T+1|T} & = \ell_{T}+b_{T}.\\
y_{T+2}         & = (\ell_{T+1} + b_{T+1})(1 + \varepsilon_{T+2})\\
                & = \left\{
                    (\ell_T + b_T) (1+ \alpha\varepsilon_{T+1}) +
                    \left[b_T + \beta (\ell_T + b_T)\varepsilon_{T+1}\right]
                    \right\}
                   (1 + \varepsilon_{T+2}) \\
\hat{y}_{T+2|T} &= \ell_{T}+2b_{T}
\end{align*}
etc.

## Forecasting with ETS models

\structure{Prediction intervals:} can only generated using the models.

  * The prediction intervals will differ between models with additive and multiplicative errors.
  * Exact formulae for some models.
  * More general to simulate future sample paths, conditional on the last estimate of the states, and to obtain prediction intervals from the percentiles of these simulated future paths.

## Prediction intervals
\fontsize{12}{13}\sf\vspace*{-0.2cm}

PI for most ETS models: $\hat{y}_{T+h|T} \pm c \sigma_h$, where $c$ depends on coverage probability and $\sigma_h$ is forecast standard deviation.

\fontsize{10}{12}\sf\vspace*{0.2cm}

\hspace*{-0.8cm}\begin{tabular}{ll}
\hline
(A,N,N) & $\sigma_h = \sigma^2\big[1 + \alpha^2(h-1)\big]$\\
(A,A,N) & $\sigma_h = \sigma^2\Big[1 + (h-1)\big\{\alpha^2 + \alpha\beta h + \frac16\beta^2h(2h-1)\big\}\Big]$\\
(A,A$_d$,N) & $\sigma_h = \sigma^2\biggl[1 + \alpha^2(h-1) + \frac{\beta\phi h}{(1-\phi)^2} \left\{2\alpha(1-\phi) +\beta\phi\right\}$\\
      & \hspace*{1.5cm}$\mbox{} - \frac{\beta\phi(1-\phi^h)}{(1-\phi)^2(1-\phi^2)} \left\{ 2\alpha(1-\phi^2)+ \beta\phi(1+2\phi-\phi^h)\right\}\biggr]$\\
(A,N,A) &              $\sigma_h = \sigma^2\Big[1 + \alpha^2(h-1) + \gamma k(2\alpha+\gamma)\Big]$\\
(A,A,A) &              $\sigma_h = \sigma^2\Big[1 + (h-1)\big\{\alpha^2 + \alpha\beta h + \frac16\beta^2h(2h-1)\big\} + \gamma k \big\{2\alpha+ \gamma + \beta m (k+1)\big\} \Big]$\\
(A,A$_d$,A) &  $\sigma_h = \sigma^2\biggl[1 + \alpha^2(h-1) +\frac{\beta\phi h}{(1-\phi)^2} \left\{2\alpha(1-\phi)  + \beta\phi \right\}$\\
  & \hspace*{1.5cm}$\mbox{} - \frac{\beta\phi(1-\phi^h)}{(1-\phi)^2(1-\phi^2)} \left\{ 2\alpha(1-\phi^2)+ \beta\phi(1+2\phi-\phi^h)\right\}$ \\
  & \hspace*{1.5cm}$\mbox{} + \gamma k(2\alpha+\gamma)  + \frac{2\beta\gamma\phi}{(1-\phi)(1-\phi^m)}\left\{k(1-\phi^m) - \phi^m(1-\phi^{mk})\right\}\biggr]$
\end{tabular}

## Example: Corticosteroid drug sales

\fontsize{10}{10}\sf

```{r h02-plot, echo = TRUE}
h02 <- tsibbledata::PBS %>%
  filter(ATC2 == "H02") %>%
  summarise(Cost = sum(Cost))
h02 %>%
  autoplot(Cost)
```

## Example: Corticosteroid drug sales
\fontsize{8}{8}\sf

```{r, echo=TRUE}
h02 %>% model(ETS(Cost)) %>% report
```

## Example: Corticosteroid drug sales
\fontsize{8}{8}\sf

```{r, echo=TRUE}
h02 %>% model(ETS(Cost ~ error("A") + trend("A") + season("A"))) %>% report
```

## Example: Corticosteroid drug sales

\fontsize{10}{10}\sf
```{r, echo=TRUE, fig.height=4}
h02 %>% model(ETS(Cost)) %>% forecast() %>% autoplot(h02)
```

## Example: Corticosteroid drug sales
\fontsize{11}{13}\sf

```{r, echo=TRUE, results = "hide"}
h02 %>%
  model(
    auto = ETS(Cost),
    AAA = ETS(Cost ~ error("A") + trend("A") + season("A"))
  ) %>%
  accuracy()
```

```{r}
h02 %>%
  model(
    auto = ETS(Cost),
    AAA = ETS(Cost ~ error("A") + trend("A") + season("A"))
  ) %>%
  accuracy() %>%
  transmute(Model = .model, ME, MAE, RMSE, MAPE, MASE) %>%
  knitr::kable(booktabs = TRUE)
```

## Your turn

* Use `ETS()` on some of these series:\vspace*{0.2cm}

> `tourism`, `gafa_stock`, `pelt`

* Does it always give good forecasts?

* Find an example where it does not work well. Can you figure out why?

