11-advanced.qmd

---
title: "ETC3550/ETC5550 Applied forecasting"
author: "Ch11. Advanced methods"
institute: "OTexts.org/fpp3/"
pdf-engine: pdflatex
fig-width: 7.5
fig-height: 3.5
format:
  beamer:
    theme: monash
    aspectratio: 169
    fontsize: 14pt
    section-titles: false
    knitr:
      opts_chunk:
        dev: "cairo_pdf"
include-in-header: header.tex
execute:
  echo: false
  message: false
  warning: false
---


```{r setup, include=FALSE}
source("setup.R")
library(readr)
```

# Complex seasonality

## Examples

```{r gasolinedata}
us_gasoline |> autoplot(Barrels) +
  labs(
    x = "Year", y = "Thousands of barrels per day",
    title = "Weekly US finished motor gasoline products"
  )
```

## Examples

```{r callsdata, message=FALSE, warning=FALSE}
calls <- read_tsv("http://robjhyndman.com/data/callcenter.txt") |>
  rename(time = `...1`) |>
  pivot_longer(-time, names_to = "date", values_to = "volume") |>
  mutate(
    date = as.Date(date, format = "%d/%m/%Y"),
    datetime = as_datetime(date) + time
  ) |>
  as_tsibble(index = datetime)
calls |>
  fill_gaps() |>
  autoplot(volume) +
  labs(
    x = "Weeks", y = "Call volume",
    title = "5 minute call volume at North American bank"
  )
```

## Examples

```{r, fig.height=5}
library(sugrrants)
calls |>
  filter(yearmonth(date) == yearmonth("2003 August")) |>
  ggplot(aes(x = time, y = volume)) +
  geom_line() +
  facet_calendar(date) +
  labs(
    x = "Weeks", y = "Call volume",
    title = "5 minute call volume at North American bank"
  )
```

## Examples

```{r turk}
turkey_elec <- read_csv("data/turkey_elec.csv", col_names = "Demand") |>
  mutate(Date = seq(ymd("2000-01-01"), ymd("2008-12-31"), by = "day")) |>
  as_tsibble(index = Date)
turkey_elec |> autoplot(Demand) +
  labs(
    title = "Turkish daily electricity demand",
    x = "Year", y = "Electricity Demand (GW)"
  )
```

## TBATS model

\begin{alertblock}{\Large TBATS}
\textbf{\Large T}rigonometric terms for seasonality\\
\textbf{\Large B}ox-Cox transformations for heterogeneity\\
\textbf{\Large A}RMA errors for short-term dynamics\\
\textbf{\Large T}rend (possibly damped)\\
\textbf{\Large S}easonal (including multiple and\\\hfill non-integer periods)
\end{alertblock}

## TBATS model
\vspace*{-0.8cm}\fontsize{13}{14}\sf

\begin{align*}
y_t&= \text{observation at time $t$}\\
y_t^{(\omega)} &= \begin{cases} (y_t^\omega-1)/\omega & \text{if $\omega\ne0$};\\
\log y_t & \text{if $\omega=0$}.\end{cases}\hspace*{10cm}\\
y_t^{(\omega)} &= \ell_{t-1} + \phi b_{t-1} + \sum_{i=1}^M s_{t-m_i}^{(i)} + d_t\\
\ell_t &= \ell_{t-1} + \phi b_{t-1} + \alpha d_t\\
b_t &= (1-\phi) b + \phi b_{t-1} + \beta d_{t}\\
d_t &= \sum_{i=1}^p \phi_i d_{t-i} + \sum_{j=1}^q \theta_j \varepsilon_{t-j} + \varepsilon_t\\
s_t^{(i)} &= \sum_{j=1}^{k_i} s_{j,t}^{(i)}
\end{align*}
\begin{textblock}{8}(4.3,7.5)
\begin{align*}
s_{j,t}^{(i)} &= \phantom{-}s_{j,t-1}^{(i)}\cos \lambda_j^{(i)} + s_{j,t-1}^{*(i)}\sin \lambda_j^{(i)} + \gamma_1^{(i)} d_t\\
s_{j,t}^{(i)} &= -s_{j,t-1}^{(i)}\sin \lambda_j^{(i)} + s_{j,t-1}^{*(i)}\cos \lambda_j^{(i)} + \gamma_2^{(i)} d_t
\end{align*}
\end{textblock}

\only<2->{\begin{textblock}{5}(7.3,1.8)\begin{block}{}Box-Cox transformation\end{block}\end{textblock}}
\only<3->{\begin{textblock}{5}(7.3,3.4)\begin{block}{}$M$ seasonal periods\end{block}\end{textblock}}
\only<4->{\begin{textblock}{5}(7.3,4.7)\begin{block}{}global and local trend\end{block}\end{textblock}}
\only<5->{\begin{textblock}{5}(7.3,6.0)\begin{block}{}ARMA error\end{block}\end{textblock}}
\only<6->{\begin{textblock}{5}(7.3,7.3)\begin{block}{}Fourier-like seasonal terms\end{block}\end{textblock}}

\only<7>{\begin{textblock}{4}(2.5,2.5)\large\begin{alertblock}{TBATS}
\textbf{T}rigonometric\\
\textbf{B}ox-Cox\\
\textbf{A}RMA\\
\textbf{T}rend\\
\textbf{S}easonal
\end{alertblock}\end{textblock}}

## Complex seasonality
\fontsize{12}{14}\sf

```{r gasoline, echo=TRUE, fig.height=4, eval = FALSE}
gasoline |>
  tbats() |>
  forecast() |>
  autoplot()
```

## Complex seasonality
\fontsize{12}{14}\sf

```{r callcentref, echo=TRUE, fig.height=4, eval = FALSE}
calls |>
  tbats() |>
  forecast() |>
  autoplot()
```

## Complex seasonality
\fontsize{12}{14}\sf

```{r telecf, echo=TRUE, fig.height=4, eval = FALSE}
telec |>
  tbats() |>
  forecast() |>
  autoplot()
```

## TBATS model

\begin{alertblock}{\Large TBATS}\fontsize{13.5}{19}\sf
\textbf{\Large T}rigonometric terms for seasonality\\
\textbf{\Large B}ox-Cox transformations for heterogeneity\\
\textbf{\Large A}RMA errors for short-term dynamics\\
\textbf{\Large T}rend (possibly damped)\\
\textbf{\Large S}easonal (including multiple and non-integer periods)
\end{alertblock}\vspace*{0.1cm}\fontsize{13}{15}\sf

* Handles non-integer seasonality, multiple seasonal periods.
* Entirely automated
* Prediction intervals often too wide
* Very slow on long series

\vspace*{10cm}

# Vector autoregression

# Neural network models

## Neural network models

\alert{Simplest version: linear regression}

\begin{tikzpicture}[shorten >=1pt,->,draw=black!50, node distance=2.5cm]
    \tikzstyle{every pin edge}=[<-,shorten <=1pt]
    \tikzstyle{neuron}=[circle,fill=black!25,minimum size=17pt,inner sep=0pt]
    \tikzstyle{input neuron}=[neuron, fill=green!50];
    \tikzstyle{output neuron}=[neuron, fill=red!50];
    \tikzstyle{hidden neuron}=[neuron, fill=blue!50];
    \tikzstyle{annot} = [text width=4em, text centered]
    % Draw the input layer nodes
    \foreach \name / \y in {1,...,4}
    % This is the same as writing \foreach \name / \y in {1/1,2/2,3/3,4/4}
        \node[input neuron, pin=left:Input \#\y] (I-\name) at (0,-\y) {};
    % Draw the output layer node
    \node[output neuron,pin={[pin edge={->}]right:Output}, right of=I-2, yshift=-0.5cm] (O) {};
    % Connect every node in the  layer with the output layer
    \foreach \source in {1,...,4}
        \path (I-\source) edge (O);
    % Annotate the layers
    \node[annot] (input) {Input layer};
    \node[annot,right of=input] {Output layer};
\end{tikzpicture}\pause\fontsize{12}{14}\sf\vspace*{-0.2cm}

* Coefficients attached to predictors are called "weights".
* Forecasts are obtained by a linear combination of inputs.
* Weights selected using a "learning algorithm" that minimises a "cost function".

## Neural network models

\alert{Nonlinear model with one hidden layer}

\begin{tikzpicture}[shorten >=1pt,->,draw=black!50, node distance=2.5cm]
    \tikzstyle{every pin edge}=[<-,shorten <=1pt]
    \tikzstyle{neuron}=[circle,fill=black!25,minimum size=17pt,inner sep=0pt]
    \tikzstyle{input neuron}=[neuron, fill=green!50];
    \tikzstyle{output neuron}=[neuron, fill=red!50];
    \tikzstyle{hidden neuron}=[neuron, fill=blue!50];
    \tikzstyle{annot} = [text width=4em, text centered]
    % Draw the input layer nodes
    \foreach \name / \y in {1,...,4}
    % This is the same as writing \foreach \name / \y in {1/1,2/2,3/3,4/4}
        \node[input neuron, pin=left:Input \#\y] (I-\name) at (0,-\y) {};
    % Draw the hidden layer nodes
    \foreach \name / \y in {1,...,3}
        \path[yshift=-.5cm]
            node[hidden neuron] (H-\name) at (2.5cm,-\y cm) {};
    % Draw the output layer node
    \node[output neuron,pin={[pin edge={->}]right:Output}, right of=H-2] (O) {};
    % Connect every node in the input layer with every node in the
    % hidden layer.
    \foreach \source in {1,...,4}
        \foreach \dest in {1,...,3}
            \path (I-\source) edge (H-\dest);
    % Connect every node in the hidden layer with the output layer
    \foreach \source in {1,...,3}
        \path (H-\source) edge (O);
    % Annotate the layers
    \node[annot,above of=H-2, node distance=2.5cm] (hl) {Hidden layer};
    \node[annot,left of=hl] {Input layer};
    \node[annot,right of=hl] {Output layer};
\end{tikzpicture}
\pause\vspace*{-0.1cm}\fontsize{13}{13}\sf
* A **multilayer feed-forward network** where each layer of nodes receives inputs from the previous layers.
* Inputs to each node combined using linear combination.
* Result modified by nonlinear function before being output.

## Neural network models

Inputs to hidden neuron $j$ linearly combined:
$$
z_j = b_j + \sum_{i=1}^4 w_{i,j} x_i.
$$
Modified using nonlinear function such as a sigmoid:
$$
s(z) = \frac{1}{1+e^{-z}},
$$
This tends to reduce the effect of extreme input values, thus making the network somewhat robust to outliers.

## Neural network models

* Weights take random  values to begin with, which are then updated using the observed data.
* There is an element of randomness in the predictions. So the network is usually trained several times using different random starting points, and the results are averaged.
* Number of hidden layers, and the number of nodes in each hidden layer, must be specified in advance.

## NNAR models
\fontsize{14}{15}\sf

* Lagged values of the time series can be used as inputs to a neural network.
* NNAR($p,k$): $p$ lagged inputs and $k$ nodes in the single hidden layer.
* NNAR($p,0$) model is equivalent to an ARIMA($p,0,0$) model but without stationarity restrictions.
* Seasonal NNAR($p,P,k$): inputs $(y_{t-1},y_{t-2},\dots,y_{t-p},y_{t-m},y_{t-2m},y_{t-Pm})$ and $k$ neurons in the hidden layer.
* NNAR($p,P,0$)$_m$ model is equivalent to an ARIMA($p,0,0$)($P$,0,0)$_m$ model but without stationarity restrictions.

## NNAR models in R
\fontsize{14}{16}\sf

* The `nnetar()` function fits an NNAR($p,P,k$)$_m$ model.
* If $p$ and $P$ are not specified, they are automatically selected.
* For non-seasonal time series, default $p=$ optimal number of lags (according to the AIC) for a linear AR($p$) model.
* For seasonal time series, defaults are $P=1$ and $p$ is chosen from the optimal linear model fitted to the seasonally adjusted data.
* Default $k=(p+P+1)/2$ (rounded to the nearest integer).

## Sunspots

* Surface of the sun contains magnetic regions that appear as dark spots.
* These affect the propagation of radio waves and so telecommunication companies like to predict sunspot activity in order to plan for any future difficulties.
* Sunspots follow a cycle of length between 9 and 14 years.

## Sunspots
\fontsize{11}{11}\sf

```{r sunspots, echo=TRUE}
sunspots <- sunspot.year |> as_tsibble()
sunspots |> autoplot(value)
```

## NNAR(9,5) model for sunspots
\fontsize{11}{11}\sf

```{r sunspot-nnetar, echo=TRUE}
sunspots <- sunspot.year |> as_tsibble()
fit <- sunspots |> model(NNETAR(value))
fit |>
  forecast(h = 20, times = 1) |>
  autoplot(sunspots, level = NULL)
```

## Prediction intervals by simulation
\fontsize{11}{11}\sf

```{r sunspot-nnetar-pi, echo=TRUE}
fit |>
  forecast(h = 20) |>
  autoplot(sunspots)
```

# Bootstrapping and bagging