fforma.Rmd

---
title: "Feature-based forecasting algorithms for large collections of time series"
author: "Rob J Hyndman"
date: "25 January 2019"
abstract: "I will discuss two algorithms used in forecasting large collections of diverse time series. Each of these algorithms uses a meta-learning approach with vectors of features computed from  time series to guide the way the forecasts are computed. In FFORMS (Feature-based FORecast Model Selection), we use a random forest classifier to identify the best forecasting method using only time series features. A key advantage of our proposed framework is that the time-consuming process of building a classifier is handled in advance of the forecasting task at hand, and only the selected forecasting model needs to be computed in real time. In FFORMA (Feature-based FORecast Model Averaging), we use gradient boosting to obtain the weights for forecast combinations using as inputs only a vector of time series features. This is slower than FFORMS (because forecasts from all candidate models must be computed), but it provides substantially more accurate forecasts. Both approaches perform very well compared to competitive methods in large forecasting competitions, with FFORMA achieving 2nd place in the recent M4 forecasting competition."
fontsize: 14pt
titlefontsize: 22pt
output:
  binb::monash:
    fig_height: 5
    fig_width: 8
    highlight: tango
    incremental: no
    keep_tex: no
    includes:
      in_header: preamble.tex
colortheme: monashblue
toc: yes
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  echo = FALSE, message = FALSE, warning = FALSE, cache = TRUE,
  dev.args = list(bg = grey(0.9), pointsize = 11)
)
library(tidyverse)
library(Mcomp)
library(fpp2)
library(GGally)
library(tsfeatures)
source("functions.R")
set.seed(20180605)
options(digits = 3, width = 63)
```


# Makridakis forecasting competitions

## M competition: 1982

\placefig{0.1}{1.4}{height=8.2cm,width=10cm}{M1}

\only<2->{\placefig{1}{4.3}{height=4cm,width=10cm,keepaspectratio}{SMakridakis}}

\only<3>{\begin{textblock}{5.5}(6.5,2)
  \begin{block}{M-competition}
  \begin{itemize}
  \item 1001 series from demography, industry, economics.
  \item Annual, quarterly, monthly data.
  \item Anyone could submit forecasts.
  \item Multiple forecast measures used.
  \end{itemize}
  \end{block}\end{textblock}}
}

## M3 competition: 2000

\full{M3paper}

## M3 competition: 2000
\fontsize{13}{14}\sf
\begin{block}{}
``The M3-Competition is a final attempt by the authors to settle the accuracy issue of
various time series methods\dots\ The extension involves the inclusion of more methods/ researchers (in particular in the areas of neural networks and expert systems) and more series.''
\end{block}

  * 3003 series
  * All data from business, demography, finance and economics.
  * Series length between 14 and 126.
  * Either non-seasonal, monthly or quarterly.
  * All time series positive.

## M4 competition: 2018

\full{m4}

## M4 competition: 2018

 * January -- May 2018
 * 100,000 time series: yearly, quarterly, monthly, weekly, daily, hourly.
 * Point forecast and prediction intervals assessed.
 * Code must be public
 * 248 registrations, 50 submissions.

\pause

### Winning methods
 1. Hybrid of Recurrent Neural Network and Exponential Smoothing models
 2. FFORMA: Feature-based forecast combinations using xgboost to find weights


# Time series features

```{r scalem3}
scalem3 <- list()
for (i in 1:3003)
{
  scalem3[[i]] <- M3[[i]]$x - min(M3[[i]]$x)
  scalem3[[i]] <- as.numeric(scalem3[[i]] / max(scalem3[[i]]))
}
```

## Key idea
\placefig{9.1}{.5}{width=3.6cm}{tukey}
\begin{textblock}{3}(9.7,5.4)\small\textit{John W Tukey}\end{textblock}
\begin{textblock}{8}(0.7,1.2)
\begin{alertblock}{Cognostics}
Computer-produced diagnostics\\ (Tukey and Tukey, 1985).
\end{alertblock}
\end{textblock}\pause
\vspace*{2.5cm}

\alert{Examples for time series}

  * lag correlation
  * size and direction of trend
  * strength of seasonality
  * timing of peak seasonality
  * spectral entropy

\vspace*{0.3cm}
\begin{block}{}
Called ``features'' in the machine learning literature.
\end{block}



```{r M3data, include=FALSE}
M3data <- purrr::map(
  Mcomp::M3,
  function(x) {
    tspx <- tsp(x$x)
    ts(c(x$x, x$xx), start = tspx[1], frequency = tspx[3])
  }
)
```

\fontsize{9}{10}\sf

```{r M3Features, include=FALSE, dependson="M3data"}
lambda_stl <- function(x, ...) {
  lambda <- forecast::BoxCox.lambda(x, lower = 0, upper = 1, method = "loglik")
  y <- forecast::BoxCox(x, lambda)
  c(stl_features(y, s.window = "periodic", robust = TRUE, ...),
    lambda = lambda
  )
}
M3Features <- bind_cols(
  tsfeatures(M3data, c("frequency", "entropy")),
  tsfeatures(M3data, "lambda_stl", scale = FALSE)
) %>%
  select(frequency, entropy, trend, seasonal_strength, e_acf1, lambda) %>%
  replace_na(list(seasonal_strength = 0)) %>%
  dplyr::rename(
    Frequency = frequency,
    Entropy = entropy,
    Trend = trend,
    Season = seasonal_strength,
    ACF1 = e_acf1,
    Lambda = lambda
  ) %>%
  mutate(Period = as.factor(Frequency))
```

```{r M3examples, include=FALSE, dependson="M3Features"}
# Consider only long series
n <- unlist(lapply(M3, function(x) {
  x$n
}))
M3Featureslong <- M3Features[n > 50, ]
M3long <- M3[names(M3)[n > 50]]
fnames <- c("M3Freq", "M3spec", "M3trend", "M3season", "M3acf", "M3lambda")
k <- NROW(M3Featureslong)
for (i in 1:6)
{
  j <- order(M3Featureslong[[i]])
  savepdf(paste(fnames[i], "Lo", sep = ""), width = 20, height = 7)
  print(autoplot(M3long[[j[1]]]$x) +
    ylab(M3long[[j[1]]]$sn) + xlab(""))
  endpdf()
  savepdf(paste(fnames[i], "Hi", sep = ""), width = 20, height = 7)
  print(autoplot(M3long[[j[k]]]$x) +
    ylab(M3long[[j[k]]]$sn) + xlab(""))
  endpdf()
}
```

## Distribution of Period for M3

```{r M3period, dependson="M3Features"}
ggally_barDiag(M3Features,
  mapping = aes(Period), width = 0.2,
  colour = "#cc5900", fill = "#cc5900"
)
```

## Distribution of Seasonality for M3

```{r M3season, dependson="M3Features"}
gghist(M3Features, aes_string("Season"))
```

\only<2->{
\begin{textblock}{6}(0.2,3)
  \begin{alertblock}{Low Seasonality}
    \includegraphics[width=6cm]{M3seasonLo.pdf}
  \end{alertblock}
\end{textblock}
}
\only<3>{
\begin{textblock}{6}(6.6,3)
  \begin{alertblock}{High Seasonality}
    \includegraphics[width=6cm]{M3seasonHi.pdf}
  \end{alertblock}
\end{textblock}
}

## Distribution of Trend for M3

```{r M3trend, dependson="M3Features"}
gghist(M3Features, aes_string("Trend"))
```

\only<2->{
\begin{textblock}{6}(0.2,3)
  \begin{alertblock}{Low Trend}
    \includegraphics[width=6cm]{M3trendLo.pdf}
  \end{alertblock}
\end{textblock}
}
\only<3>{
\begin{textblock}{6}(6.6,3)
  \begin{alertblock}{High Trend}
    \includegraphics[width=6cm]{M3trendHi.pdf}
  \end{alertblock}
\end{textblock}
}

## Distribution of Residual ACF1 for M3

```{r M3ACF1, dependson="M3Features"}
gghist(M3Features, aes_string("ACF1"))
```

\only<2->{
\begin{textblock}{6}(0.2,3)
  \begin{alertblock}{Low ACF1}
    \includegraphics[width=6cm]{M3acfLo.pdf}
  \end{alertblock}
\end{textblock}
}
\only<3>{
\begin{textblock}{6}(6.6,3)
  \begin{alertblock}{High ACF1}
    \includegraphics[width=6cm]{M3acfHi.pdf}
  \end{alertblock}
\end{textblock}
}

## Distribution of Spectral Entropy for M3

```{r M3entropy, dependson="M3Features"}
gghist(M3Features, aes_string("Entropy"))
```

\only<2->{
\begin{textblock}{6}(0.2,3)
  \begin{alertblock}{Low Entropy}
    \includegraphics[width=6cm]{M3specLo.pdf}
  \end{alertblock}
\end{textblock}
}
\only<3>{
\begin{textblock}{6}(6.6,3)
  \begin{alertblock}{High Entropy}
    \includegraphics[width=6cm]{M3specHi.pdf}
  \end{alertblock}
\end{textblock}
}

## Feature distributions

```{r ACF1SE, dependson="M3Features"}
ggplot(M3Features, aes(x = Entropy, y = ACF1)) + geom_point()
```

## Feature distributions

```{r TrendSE, dependson="M3Features"}
ggplot(M3Features, aes(x = Entropy, y = Trend)) + geom_point()
```

## Feature distributions

```{r M3pairs, dependson="M3Features"}
# Fig 1 of paper
yk_ggally_densityDiag <- wrap(gghist, adjust = 0.5)
yk_ggally_barDiag <- wrap(ggally_barDiag,
  colour = "#cc5900",
  fill = "#cc5900", width = 0.2
)
M3Features %>%
  select(Period, Entropy, Trend, Season, ACF1, Lambda) %>%
  ggpairs(
    diag = list(
      continuous = yk_ggally_densityDiag,
      discrete = yk_ggally_barDiag
    ),
    axisLabels = "none",
    lower = list(continuous = wrap("points", alpha = 0.5, size = 0.2))
  ) -> p
print(p)
savepdf("PairwisePlot")
print(p)
endpdf()
```



# Feature-based forecasting algorithms

## Features used to select a forecasting model

\begin{textblock}{12}(0.1,1.3)\small
\begin{multicols}{2}
  \begin{itemize}\tightlist
    \item length
    \item strength of seasonality
    \item strength of trend
    \item linearity
    \item curvature
    \item spikiness
    \item stability
    \item lumpiness
    \item parameter estimates of Holt's linear trend method
    \item spectral entropy
    \item Hurst exponent
    \item nonlinearity
    \item parameter estimates of Holt-Winters' additive method
    \item unit root test statistics
    \item crossing points, flat spots
    \item peaks, troughs
    \item ACF and PACF based features - calculated on raw, differenced, and remainder series.
    \item ARCH/GARCH statistics and ACF of squared series and residuals.
    \end{itemize}
\end{multicols}
\end{textblock}

## Features used to select a forecasting model

\alert{Why these features?}

 * Hyndman, Wang and Laptev. “Large scale unusual time series detection” (ICDM 2015).
 * Kang, Hyndman & Smith-Miles. “Visualising forecasting algorithm performance using time series instance spaces” (IJF 2017).
 * Talagala, Hyndman and Athanasopoulos. “Meta-learning how to forecast time series” (2018).
 * Implemented in the tsfeatures R package

## \fontsize{15}{15}\bf\sffamily FFORMS: Feature-based FORecast Model Selection

  * Using large collection of time series, each series split into training and test sets.
  * Features computed on training data.
  * All forecasting methods fitted to training data, and forecasts obtained for test data period.
  * Forecast accuracy for each method/series computed from test data.
  * Train a random forest to identify the most accurate forecasting method for a given time series using only a vector of features.

## \fontsize{15}{15}\bf\sffamily FFORMA: Feature-based FORecast Model Averaging

 * Like FFORMS but using gradient boosted trees (xgboost) rather than random forest.
 * Trained on temporal holdout version of M4 dataset, where size of test sets equal to required forecast horizons
 * Optimization criterion: forecast accuracy not classification accuracy.
 * Probability of each model being best is used to construct model weights for combination forecast.
 * 5 days computing time.


## \fontsize{15}{15}\bf\sffamily FFORMA: Feature-based FORecast Model Averaging

\placefig{0.2}{1.4}{width=12.4cm, trim=0 0 0 70, clip=true}{fforma_graphic}


\vspace*{5.2cm}

### M4 competition results (based on average OWA)

```{r, results='asis'}
tribble(
    ~Place, ~OWA, ~Method,
    "1st", 0.821, NA,
    "2nd", 0.838, "(FFORMA)",
    "3rd", 0.841, NA
  ) %>%
  baretable(digits=3)
```

## \fontsize{15}{15}\bf\sffamily FFORMA: Feature-based FORecast Model Averaging

### Models included

1. Naive
1. Seasonal naive
1. Random walk with drift
1. Theta method
1. ARIMA
1. ETS
1. TBATS
1. STL decomposition with AR for seasonally adjusted series
1. Neural network autoregression


## \fontsize{15}{15}\bf\sffamily FFORMA: Feature-based FORecast Model Averaging

\only<1>{\placefig{0.2}{1.4}{width=12.4cm, trim=0 0 0 0, clip=true}{prototypes}}
\only<2>{\placefig{0.2}{1.4}{width=12.4cm, trim=0 0 0 0, clip=true}{prototype1}}
\only<3>{\placefig{0.2}{1.4}{width=12.4cm, trim=0 0 0 0, clip=true}{prototype2}}
\only<4>{\placefig{0.2}{1.4}{width=12.4cm, trim=0 0 0 0, clip=true}{prototype3}}
\only<5>{\placefig{0.2}{1.4}{width=12.4cm, trim=0 0 0 0, clip=true}{prototype4}}
\only<6>{\placefig{0.2}{1.4}{width=12.4cm, trim=0 0 0 0, clip=true}{prototype5}}



## Papers and packages
\fontsize{14}{17}\sf

\begin{block}{R packages}
\begin{itemize}\tightlist
 \item \alert{tsfeatures}: Calculating time series features. \newline\url{github.com/robjhyndman/tsfeatures}
 \item \alert{seer}: FFORMS --- selecting forecasting model using features. \newline\url{github.com/thiyangt/seer}
 \item \alert{M4metalearning}: FFORMA -- forecast combinations using features to choose weights. \newline\url{github.com/robjhyndman/M4metalearning}
\end{itemize}
\end{block}

\begin{alertblock}{Papers}
Available from \url{robjhyndman.com}
\end{alertblock}


## Acknowledgements

\begin{textblock}{12.5}(0.2,1.2)
\begin{block}{}\fontsize{9}{10}\sf
\centering\begin{tabular}{p{3.4cm}p{3.4cm}p{3.4cm}}
\includegraphics[height=3.4cm, width=10cm, keepaspectratio]{kate} &
\includegraphics[height=3.4cm, width=10cm, keepaspectratio]{yanfei}&
\includegraphics[height=3.4cm, width=10cm, keepaspectratio]{earowang} \\
Kate Smith-Miles & Yanfei Kang & Earo Wang \\
\includegraphics[height=3.4cm, width=10cm, keepaspectratio]{thiyanga} &
\includegraphics[height=3.4cm, width=10cm, keepaspectratio]{george} &
\includegraphics[height=3.4cm, width=10cm, keepaspectratio]{pablo} \\
Thiyanga Talagala & George Athanasopoulos & Pablo \rlap{Montero-Manso}
\end{tabular}
\end{block}
\end{textblock}