---
title: "STAT481/581: Introduction to Time Series Analysis"
author: "Advanced forecasting methods"
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(tsibbledata)
library(feasts)
library(lubridate)
library(forecast)
```

# Complex seasonality

## Examples

```{r gasolinedata}
gasoline <- as_tsibble(fpp2::gasoline)
gasoline %>% autoplot(value) +
  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") %>%
  gather("date", "volume", -X1) %>% transmute(
    time = X1, date = as.Date(date, format = "%d/%m/%Y"),
    datetime = as_datetime(date) + time, volume) %>%
  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 = TRUE}
fit <- tbats(USAccDeaths)
plot(forecast(fit))
```



## 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}


# 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.

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

```{r sunspot-nnetar, echo=TRUE, fig.height=4}
sunspots <- as_tsibble(fpp2::sunspotarea)
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, fig.height=4}
fit %>% forecast(h=20) %>% 
  autoplot(sunspots)
```

