---
title: "STAT481/581: Introduction to Time Series Analysis"
author: "Ch3. The forecasters' toolbox"
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 = TRUE, cache=TRUE)

library(tidyverse)
library(fable)
library(tsibble)
library(feasts)
library(lubridate)
library(tsibbledata)
library(TSA)
source("nicefigs.R")
source("box_cox.R")
options(width=60)
```

# A tidy forecasting workflow

<!-- TODO: Restructure -->
<!-- 3a: toolbox -->
<!-- 1. Tidy forecasting workflow (with progress update on stage in workflow) -->
<!-- 2. Simple forecasting methods (using model() and forecast()) -->
<!-- 3. The model definition: Transformations, adjustments and model specials -->
<!-- 4. Distributional forecasts and the prediction intervals -->

## A tidy forecasting workflow

The process of producing forecasts can be split up into a few fundamental steps.

1. Preparing data
2. Data visualisation
3. Specifying a model
4. Model estimation
5. Accuracy \& performance evaluation
6. Producing forecasts

## A tidy forecasting workflow

```{r workflow, echo = FALSE}
line_curve <- function(x, y, xend, yend, ...){
  geom_curve(
    aes(x = x, y = y, xend = xend, yend = yend),
    arrow = arrow(type = "closed", length = unit(0.03, "npc")),
    ...
  )
}

ggplot() +
  geom_text(
    aes(x = x, y = y, label = label),
    data = tribble(
      ~ x, ~ y, ~ label,
      1, 0, "Tidy",
      7/3, 0, "Visualise",
      3, 0.5, "Specify",
      11/3, 0, "Estimate",
      3, -0.5, "Evaluate",
      5, 0, "Forecast"
    ),
    size = 5
  ) +
  geom_segment(
    aes(x = x, y = y, xend = xend, yend = yend),
    data = tribble(
      ~ x, ~ y, ~ xend, ~ yend,
      1.3, 0, 1.9, 0,
      4.1, 0, 4.6, 0
    ),
    arrow = arrow(type = "closed", length = unit(0.03, "npc"))
  ) +
  line_curve(7/3, 0.1, 8/3, 0.5, angle = 250, curvature = -0.3) +
  line_curve(10/3, 0.5, 11/3, 0.1, angle = 250, curvature = -0.3) +
  line_curve(8/3, -0.5, 7/3, -0.1, angle = 250, curvature = -0.3) +
  line_curve(11/3, -0.1, 10/3, -0.5, angle = 250, curvature = -0.3) +
  theme_void() +
  xlim(0.8, 5.2) +
  ylim(-0.6, 0.6) +
  coord_equal(ratio = 1)
```

# Some simple forecasting methods

## Some simple forecasting methods

```{r ausbeer, fig.height=4.6, echo=FALSE}
new_production <- aus_production %>%
  filter(year(Quarter) >= 1992)
new_production %>% autoplot(Beer) +
  xlab("Year") + ylab("megalitres") +
    ggtitle("Australian quarterly beer production")
```

\begin{textblock}{7}(0.2,8.7)
\begin{alertblock}{}
\small{How would you forecast these series?}
\end{alertblock}
\end{textblock}

## Some simple forecasting methods

```{r pigs, fig.height=4.6, echo=FALSE}
aus_livestock %>% filter(between(year(Month), 1990, 1995),
  Animal == "Pigs", State == "Victoria") %>%
  autoplot(Count) +
  xlab("Year") + ylab("thousands") +
  ggtitle("Number of pigs slaughtered in Victoria, 1990-1995")
```

\begin{textblock}{7}(0.2,8.7)
\begin{alertblock}{}
\small{How would you forecast these series?}
\end{alertblock}
\end{textblock}

## Some simple forecasting methods

```{r dj, fig.height=4.6, echo=FALSE}
gafa_stock %>%
  filter(Symbol == "FB", Date >= ymd("2018-01-01")) %>%
  autoplot(Close) +
  labs(title = "Facebook closing stock price in 2018",
       x = NULL, y = "Closing price ($USD)")
```

\begin{textblock}{7}(0.2,8.7)
\begin{alertblock}{}
\small{How would you forecast these series?}
\end{alertblock}
\end{textblock}

## Some simple forecasting methods
\fontsize{13}{14}\sf

### `MEAN(y)`: Average method

  * Forecast of all future values is equal to mean of historical data $\{y_1,\dots,y_T\}$.
  * Forecasts: $\hat{y}_{T+h|T} = \bar{y} = (y_1+\dots+y_T)/T$

```{r mean-method-explained, echo=FALSE, message=FALSE, warning=FALSE, fig.height = 3.3}
bricks <- aus_production %>%
  filter(!is.na(Bricks)) %>%
  mutate(average = mean(Bricks))

bricks %>%
  filter(!is.na(Bricks)) %>%
  model(MEAN(Bricks)) %>%
  forecast(h = "5 years") %>%
  autoplot(filter(bricks, year(Quarter) > 1990), level = NULL) +
  geom_point(data = slice(bricks, n()), colour = "blue") +
  ggtitle("Clay brick production in Australia")
```

## Some simple forecasting methods
\fontsize{13}{14}\sf

### `NAIVE(y)`: Naïve method

  * Forecasts equal to last observed value.
  * Forecasts: $\hat{y}_{T+h|T} =y_T$.

```{r naive-method-explained, echo = FALSE, warning = FALSE, fig.height = 3.1}
bricks %>%
  filter(!is.na(Bricks)) %>%
  model(NAIVE(Bricks)) %>%
  forecast(h = "5 years") %>%
  autoplot(filter(bricks, year(Quarter) > 1990), level = NULL) +
  geom_point(data = slice(bricks, n()), colour = "blue") +
  ggtitle("Clay brick production in Australia")
```

## Some simple forecasting methods
\fontsize{13}{14}\sf

### `SNAIVE(y ~ lag(m))`: Seasonal naïve method

  * Forecasts equal to last value from same season.
  * Forecasts: $\hat{y}_{T+h|T} =y_{T+h-m(k+1)}$, where $m=$ seasonal period and $k$ is the integer part of $(h-1)/m$.

```{r snaive-method-explained, echo = FALSE, warning = FALSE, fig.height = 2.8}
bricks %>%
  model(SNAIVE(Bricks ~ lag("year"))) %>%
  forecast(h = "5 years") %>%
  autoplot(filter(bricks, year(Quarter) > 1990), level = NULL) +
  geom_point(data = slice(bricks, (n()-3):n()), colour = "blue") +
  ggtitle("Clay brick production in Australia")
```

## Some simple forecasting methods
\fontsize{13}{14}\sf

### `RW(y ~ drift())`: Drift method

 * Forecasts equal to last value plus average change.
 * Forecasts:\vspace*{-.7cm}

 \begin{align*}
 \hat{y}_{T+h|T} & =  y_{T} + \frac{h}{T-1}\sum_{t=2}^T (y_t-y_{t-1})\\
                 & = y_T + \frac{h}{T-1}(y_T -y_1).
 \end{align*}\vspace*{-0.2cm}

   * Equivalent to extrapolating a line drawn between first and last observations.

## Some simple forecasting methods

### Drift method

```{r drift-method-explained, echo = FALSE, warning = FALSE}
aus_production %>%
  filter(!is.na(Bricks)) %>%
  model(RW(Bricks ~ drift())) %>%
  forecast(h = "5 years") %>%
  autoplot(aus_production, level = NULL) +
  geom_line(data = slice(aus_production, range(cumsum(!is.na(Bricks)))),
            linetype = "dashed", colour = "blue") +
  ggtitle("Clay brick production in Australia")
```

# The workflow in action

## Data preparation and visualisation
\fontsize{10}{13}\sf


```{r beer-plot, fig.height = 3.2}
# Set training data from 1992 to 2007
train <- aus_production %>%
  filter(between(year(Quarter), 1992, 2007))
train %>% autoplot(Beer)
```

## Model estimation

The `model()` function trains models to data.

\fontsize{10}{13}\sf


```{r beer-model}
# Fit the models
beer_fit <- train %>%
  model(
    Mean = MEAN(Beer),
    `Naïve` = NAIVE(Beer),
    `Seasonal naïve` = SNAIVE(Beer), #or SNAIVE(Beer~lag(4))
    # or SNAIVE(Beer~lag("year"))
    Drift = RW(Beer ~ drift())
  )
```

## Model estimation

\fontsize{10}{13}\sf

```{r beer-mable, echo = TRUE, dependson='beer-model'}
beer_fit
```

A `mable` is a model table, each cell corresponds to a fitted model.

## Producing forecasts

\fontsize{10}{13}\sf

```{r beer-fc, echo = TRUE, dependson='beer-model'}
beer_fc <- beer_fit %>%
  forecast(h = 11)
```

```{r beer-fbl, echo = FALSE, dependson='beer-fc'}
print(beer_fc, n = 4)
```

A `fable` is a forecast table with point forecasts and distributions.

## Visualising forecasts

\footnotesize

```{r beer-fc-plot, warning=FALSE, message=FALSE, fig.height=3}
beer_fc %>%
  autoplot(train, level = NULL) + 
  #level=NULL means no prediction interval
  ggtitle("Forecasts for quarterly beer production") +
  xlab("Year") + ylab("Megalitres") +
  guides(colour=guide_legend(title="Forecast"))
#guide_legend: let the legend value determine the colour
```

## Facebook closing stock price

\fontsize{9}{10}\sf

```{r fbf, eval=FALSE}
# Extract training data
fb_stock <- gafa_stock %>%
  group_by(Symbol) %>%
  mutate(trading_day = row_number()) %>%
  update_tsibble(index=trading_day, regular=TRUE) %>%
  filter(Symbol == "FB",
         between(Date, ymd("2018-01-01"), ymd("2018-09-01")))
# Specify, estimate and forecast
fb_stock %>%
  model(
    Mean = MEAN(Close),
    `Naïve` = NAIVE(Close),
    Drift = RW(Close ~ drift())
  ) %>%
  forecast(h=42) %>%
  autoplot(fb_stock, level = NULL) +
  ggtitle("Facebook closing stock price (daily ending Sep 2018)") +
  xlab("Day") + ylab("") +
  guides(colour=guide_legend(title="Forecast"))
```

## Facebook closing stock price

```{r fb-fc-plot, echo=FALSE, fig.height=4.6}
# Extract training data
fb_stock <- gafa_stock %>%
  mutate(trading_day = row_number()) %>%
  update_tsibble(index=trading_day, regular=TRUE) %>%
  filter(Symbol == "FB",
         between(Date, ymd("2018-01-01"), ymd("2018-09-01")))

# Specify, estimate and forecast
fb_stock %>%
  model(
    Mean = MEAN(Close),
    `Naïve` = NAIVE(Close),
    Drift = RW(Close ~ drift())
  ) %>%
  forecast(h=42) %>%
  autoplot(fb_stock, level = NULL) +
  ggtitle("Facebook closing stock price (daily ending Sep 2018)") +
  xlab("Day") + ylab("") +
  guides(colour=guide_legend(title="Forecast"))
```

## Your turn

 * Produce forecasts from the appropriate method for Amazon closing price (`gafa_stock`) and Australian takeaway food turnover (`aus_retail`).
 * Plot the results using `autoplot()`.

# Transformations

## Variance stabilization

\fontsize{13}{15}\sf

If the data show different variation at different levels of the series, then a transformation can be useful.
\pause

Denote original observations as $y_1,\dots,y_n$ and transformed
observations as $w_1, \dots, w_n$.
\pause

\begin{block}{\footnotesize Mathematical transformations for stabilizing
variation}
\begin{tabular}{llc}
Square root & $w_t = \sqrt{y_t}$ & $\downarrow$ \\[0.2cm]
Cube root & $w_t = \sqrt[3]{y_t}$ & Increasing \\[0.2cm]
Logarithm & $w_t = \log(y_t)$  & strength
\end{tabular}
\end{block}
\pause

Logarithms, in particular, are useful because they are more interpretable: changes in a log value are **relative (percent) changes on the original scale**.

## Variance stabilization

```{r food, echo=TRUE}
food <- aus_retail %>%
  filter(Industry == "Food retailing") %>%
  summarise(Turnover = sum(Turnover))
```

```{r food-plot, echo = FALSE, fig.height=3.4}
food %>% autoplot(Turnover) +
  labs(y = "Turnover ($AUD)")
```

## Variance stabilization

```{r food-sqrt1, echo=TRUE, fig.height=4}
food %>% autoplot(sqrt(Turnover)) +
  labs(y = "Square root turnover")
```

## Variance stabilization

```{r food-cbrt, echo=TRUE, fig.height=4}
food %>% autoplot(Turnover^(1/3)) +
  labs(y = "Cube root turnover")
```

## Variance stabilization

```{r food-log, echo=TRUE, fig.height=4}
food %>% autoplot(log(Turnover)) +
  labs(y = "Log turnover")
```

## Variance stabilization

```{r food-inverse, echo=TRUE, fig.height=4}
food %>% autoplot(-1/Turnover) +
  labs(y = "Inverse turnover")
```

## Box-Cox transformations

Each of these transformations is close to a member of the
family of \textbf{Box-Cox transformations}:
$$w_t = \left\{\begin{array}{ll}
        \log(y_t),      & \quad \lambda = 0; \\
        (y_t^\lambda-1)/\lambda ,         & \quad \lambda \ne 0.
\end{array}\right.
$$\pause

* $\lambda=1$: (No substantive transformation)
* $\lambda=\frac12$: (Square root plus linear transformation)
* $\lambda=0$: (Natural logarithm)
* $\lambda=-1$: (Inverse plus 1)


## Box-Cox transformations
\fontsize{13}{14}\sf

```{r food-bc, echo=TRUE,fig.height=4}
food %>% autoplot(box_cox(Turnover, 1/3)) +
  labs(y = "Box-Cox transformed turnover")
```

## Box-Cox transformations

* $y_t^\lambda$ for $\lambda$ close to zero behaves like logs.
* If some $y_t=0$, then must have $\lambda>0$
* if some $y_t<0$, no power transformation is possible unless all $y_t$ adjusted by \textbf{adding a constant to all values}.
* Simple values of $\lambda$ are easier to explain.
* Results are  relatively insensitive to  $\lambda$.
* Often no transformation ($\lambda=1$) needed.
* Transformation can have very large effect on PI.
* Choosing $\lambda=0$ is a simple way to force forecasts to be positive

## Box-Cox transformations

\fontsize{10}{13}\sf


```{r food-lambda, echo=TRUE}
food %>%
  features(Turnover, features = guerrero)
# it uses the `BoxCoxLambda` function in the `forecast package`
# use the guerrero function for an automated approach
# Guerrero, V.M. (1993) Time-series analysis supported by 
# power transformations. 
# Journal of Forecasting, 12, 37--48.
```

## Box-Cox transformations

\fontsize{10}{13}\sf

* The `guerrero` approach attempts to balance the seasonal fluctuations and random variation across the series.
* Always check the results.
* A low value of $\lambda$ can give extremely large prediction intervals.

## Box-Cox transformations

\fontsize{10}{13}\sf


```{r food-lambda2, warning=FALSE, fig.height=4}
#use library(TSA)
lambda.fit <- BoxCox.ar(food$Turnover)
```

## Box-Cox transformations

\fontsize{10}{13}\sf

`BoxCox.ar` is to estimate the power transformation so that the transformed time series is approximately a Gaussian AR process (no seasonal/cyclic pattern or trend is allowed). Hence the food turner over data is not appropriate here.

```{r food-lambda3, warning=FALSE}
lambda.fit$mle
lambda.fit$ci
```


## Back-transformation

We must reverse the transformation (or \textit{back-transform}) to obtain
forecasts on the original scale. The reverse Box-Cox transformations are given
by
$$ y_t = \left\{\begin{array}{ll}
        \exp(w_t),      & \quad \lambda = 0; \\
        (\lambda W_t+1)^{1/\lambda} ,   & \quad \lambda \ne 0.
\end{array}\right.$$

## Modelling with transformations

Transformations used in the left of the formula will be automatically back-transformed. To model log-transformed food retailing turnover, you could use:

\fontsize{13}{15}\sf

```{r food-bt-fit}
fit <- food %>%
  model(SNAIVE(log(Turnover) ~ lag("year")))
```
```{r food-bt-mbl, echo = FALSE}
fit
```

## Forecasting with transformations

```{r food-bt-fc}
fc <- fit %>%
  forecast(h = "3 years")
```

\fontsize{10}{13}\sf


```{r food-bt-fbl, echo = FALSE}
print(fc, n = 6)
```

## Forecasting with transformations
\fontsize{12}{13}\sf

```{r elec9,echo=TRUE,fig.height=4}
fc %>% autoplot(filter(food,year(Month)>2010))
```

## Your turn

Find a transformation that works for the Australian gas production (`aus_production`).

## Bias adjustment

  * Back-transformed median forecasts for $w_{T+h}$ are median forecasts for $y_{T+h}$.
  * Back-transformed PI have the correct coverage.
  * A forecast $\hat{y}_{T+h|T}$ is (usually) the mean of the conditional distribution $y_{T+h} \mid y_1, \dots, y_{T}$.
 
\pause

**Back-transformed means**

Let $X$ be have mean $\mu$ and variance $\sigma^2$.

Let $f(x)$ be back-transformation function, and $Y=f(X)$.



## Bias adjustment

 
Taylor series expansion about $\mu$:
$$f(X) = f(\mu) + (X-\mu)f'(\mu) + \frac{1}{2}(X-\mu)^2f''(\mu).$$\pause

\begin{alertblock}{}
\centerline{$\E[Y] = \E[f(X)] = f(\mu) + \frac12 \sigma^2 f''(\mu)$}
\end{alertblock}


## Bias adjustment

\fontsize{13}{15}\sf

**Box-Cox back-transformation:**
\begin{align*}
y_t &= \left\{\begin{array}{ll}
        \exp(w_t)      & \quad \lambda = 0; \\
        (\lambda W_t+1)^{1/\lambda}  & \quad \lambda \ne 0.
\end{array}\right. \\
f(x) &= \begin{cases}
                        e^x & \quad\lambda=0;\\
 (\lambda x + 1)^{1/\lambda} & \quad\lambda\ne0.
 \end{cases}\\
f''(x) &= \begin{cases}
                        e^x & \quad\lambda=0;\\
 (1-\lambda)(\lambda x + 1)^{1/\lambda-2} & \quad\lambda\ne0.
 \end{cases}
\end{align*}

## Bias adjustment

\fontsize{13}{15}\sf

\begin{alertblock}{}
\centerline{$\E[Y] = \begin{cases}
                        e^\mu\left[1+\frac{\sigma^2}{2}\right] & \quad\lambda=0;\\
 (\lambda \mu + 1)^{1/\lambda}\left[1+\frac{\sigma^2(1-\lambda)}{2(\lambda\mu+1)^2}\right] & \quad\lambda\ne0.
 \end{cases}$}
\end{alertblock}

## Bias adjustment

\fontsize{9}{9}\sf

```{r biasadj, warning=FALSE}
eggs <- as_tsibble(fma::eggs)
fit <- eggs %>% model(RW(log(value) ~ drift()))
fc <- fit %>% forecast(h=50, bias_adjust = TRUE)
fc_biased <- fit %>% forecast(h=50, bias_adjust = FALSE)

```

## Bias adjustment

\fontsize{9}{9}\sf

```{r biasadj2, fig.height=3, warning=FALSE}
# the point forecasts are medians 
eggs %>% autoplot(value) +
  autolayer(fc_biased, series="Simple back transformation", level=80)
```

## Bias adjustment

\fontsize{9}{9}\sf

```{r biasadj3, fig.height=3, warning=FALSE}
# the point forecasts are bias adjusted means
eggs %>% autoplot(value) +
  autolayer(fc, series="Simple back transformation", level=80) 
```

# Distributional forecasts

## Forecast distributions

 * Most time series models produce normally distributed forecasts.
 * The forecast distribution describes the probability of observing any future value.

## Forecast distributions

\fontsize{14}{18}\sf

Assuming residuals are normal, uncorrelated, sd = $\hat\sigma$:

\begin{block}{}
\begin{tabular}{ll}
\bf Mean: & $\hat{y}_{T+h|T} \sim N(\bar{y}, (1 + 1/T)\hat{\sigma}^2)$\\[0.2cm]
\bf Naïve: & $\hat{y}_{T+h|T} \sim N(y_T, h\hat{\sigma}^2)$\\[0.2cm]
\bf Seasonal naïve: & $\hat{y}_{T+h|T} \sim N(y_{T+h-m(k+1)}, (k+1)\hat{\sigma}^2)$\\[0.2cm]
\bf Drift: & $\hat{y}_{T+h|T} \sim N(y_T + \frac{h}{T-1}(y_T - y_1),h\frac{T+h}{T}\hat{\sigma}^2)$
\end{tabular}
\end{block}

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

Note that when $h=1$ and $T$ is large, these all give the same approximate forecast variance: $\hat{\sigma}^2$.

## Prediction intervals

 * A prediction interval gives a region within which we expect $y_{T+h}$ to lie with a specified probability.
 * Assuming forecast errors are normally distributed, then a 95% PI is
 \begin{alertblock}{}
\centerline{$
  \hat{y}_{T+h|T} \pm 1.96 \hat\sigma_h
$}
\end{alertblock}
where $\hat\sigma_h$ is the st dev of the $h$-step distribution.

 * When $h=1$, $\hat\sigma_h$ can be estimated from the residuals.

## Prediction intervals

\fontsize{10}{13}\sf


```{r, echo = FALSE, cache = FALSE}
options(width = 50)
```

```{r fb-fc, echo=TRUE}
fit <- fb_stock %>% model(NAIVE(Close))
forecast(fit)
```

## Prediction intervals

\footnotesize

```{r, echo = FALSE, cache = FALSE}
options(width = 80)
```

```{r fb-pi-manual}
res_sd <- sqrt(mean(augment(fit)$.resid^2, na.rm = TRUE))
last(fb_stock$Close) + 1.96 * res_sd * c(-1,1)
```


## Prediction intervals

\footnotesize

```{r, echo = FALSE, cache = FALSE}
options(width = 80)
```

```{r fb-pi-hilo, echo=TRUE}
forecast(fit, h = 1) %>%
  transmute(interval = hilo(.distribution, level = 95))
# transmute: add a new variable to the tsibble object by forecast
# hilo: use the forecast object and compute its 95% PI
# .distribution is a column of the forecast object
```

## Prediction intervals

 * Point forecasts are often useless without a measure of uncertainty (such as prediction intervals).
 * Prediction intervals require a stochastic model (with random errors, etc).
 * Multi-step forecasts for time series require a more sophisticated approach (with PI getting wider as the forecast horizon increases).

## Prediction intervals

  * Computed automatically from the forecast distribution.
  * Use `level` argument to control coverage.
  * Check residual assumptions before believing them (we will see this next class).
  * Usually too narrow due to unaccounted uncertainty.
