∎

11institutetext: V. Cerqueira (✉) 22institutetext: Dalhousie University,
Halifax, Canada
22email<EMAIL_ADDRESS>33institutetext: Luís Torgo 44institutetext:
Dalhousie University, Halifax, Canada
44email<EMAIL_ADDRESS>
# Multi-output Ensembles for Multi-step Forecasting

Vitor Cerqueira  Luis Torgo

(Received: date / Accepted: date)

###### Abstract

This paper studies the application of ensembles composed of multi-output
models for multi-step ahead forecasting problems. Dynamic ensembles have been
commonly used for forecasting. However, these are typically designed for one-
step ahead tasks. On the other hand, the literature regarding the application
of dynamic ensembles for multi-step ahead forecasting is scarce. Moreover, it
is not clear how the combination rule is applied across the forecasting
horizon. We carried out extensive experiments to analyze the application of
dynamic ensembles for multi-step forecasting. We resorted to a case study with
3568 time series and an ensemble of 30 multi-output models. We discovered that
dynamic ensembles based on arbitrating and windowing present the best
performance according to average rank. Moreover, as the horizon increases,
most approaches struggle to outperform a static ensemble that assigns equal
weights to all models. The experiments are publicly available in a repository.

###### Keywords:

Ensemble methods Time series forecasting Multi-output models Time series

##  1 Introduction

Ensemble methods combine the output of several models to make aggregated
predictions. These approaches have been shown to improve predictive
performance over single models in many tasks brown2005managing , including
forecasting cerqueira2017dets . One of the reasons to use an ensemble is that
it reduces the chance of selecting the wrong model hibon2005combine .

Ensembles have attracted much attention in the literature. Bates and Granger
bates1969combination  pioneered the combination of different models for time
series forecasting. Since then, hundreds of works have been devoted to this
topic. The simplest approach to forecast combination is to weigh each model
equally by averaging the prediction with the arithmetic mean. Weighted
averages can also be used, where the weights reflect how well we expect each
model to perform.

Time series are often non-stationary. This means that the data distribution
changes over time. Therefore, the combination rule for forecasting problems is
typically dynamic. The weights change over time to adapt to changes in the
time series.

There are many different forecast combination strategies in the literature.
Most approaches are based on regret minimization Cesa-Bianchi:2006:PLG:1137817
, windowing Blast , or meta-learning cerqueira2019arbitrage . However, dynamic
approaches are typically designed for one-step ahead forecasting problems.
They assume that there is immediate feedback concerning the actual value for
computing the error and updating the combination rule.

Forecasting multiple observations in advance is important for many application
domains. It helps to reduce the long-term uncertainty of time series.
Consequently, this results in better planning within organizations. Approaches
for multi-step forecasting include the recursive, direct, or multi-output
methods taieb2012review . Multi-output methods have been shown to provide
better forecasting performance in multi-step ahead prediction settings
taieb2012review . In effect, our work focuses on this type of approach.

However, since dynamic combination methods are usually designed for one-step
ahead settings, it is not clear how they should be applied for multi-step
forecasting. For example, one could compute different weights for each
horizon. Or estimate the weights for all forecasting horizons jointly.

In this work, we investigate the performance of ensembles composed of multi-
output models for multi-step forecasting problems. Our goal is two-fold.
First, we aim at studying which dynamic combination rule is the best for
multi-step ahead problems. Previous works have carried out such analysis for
one-step ahead tasks, e.g. cerqueira2019arbitrage . Conversely, the literature
on dynamic ensembles for multi-step ahead forecasting is scarce. The second
objective is to analyze which is the best approach for computing the ensemble
weights along the forecasting horizon. To our knowledge, no prior work has
addressed this topic.

The analysis presented in this work is based on an empirical study. The
experiments included a set of 3568 univariate time series from multiple
domains. We test different combination rules using an ensemble composed of 30
individual multi-output models. These models are standard supervised learning
regression algorithms.

The results suggest that the dynamic combination rules based on arbitrating
cerqueira2019arbitrage  and windowing newbold1974experience  present the best
average rank. Yet, most approaches struggle to outperform the combination rule
which assigns equal weights to all models. We discovered that the forecasting
horizon has a significant impact on these results. Dynamic combination rules
tend to perform better for short-term forecasting. On the other hand, these
approaches decrease their performance as the forecasting horizon increases. In
terms of computing the ensemble weights along the horizon, no considerable
differences were found between the tested approaches. The data sets and
experiments are publicly available in an online
repository111https://github.com/vcerqueira/mo_ensembles.

The rest of the paper is organized as follows. In the next section (Section
2), we provide a background to our work. We start by defining the time series
predictive task. We also overview the literature related to our work. We focus
on the topics of dynamic ensembles and multi-step forecasting. In Section 3,
we present the materials and methods used in this paper. We describe the case
study which contains time series from several application domains. We also
describe the methods used in detail. Moreover, we explain the experimental
design used in the experiments. Then, we present the experiments in Section 4.
In that section, we outline the research questions and provide an empirical
answer to them. We discuss the results obtained in Section 5, pointing out
some future directions. Finally, the paper is concluded in Section 6.

##  2 Background

We split the literature related to our work into five parts: time series
forecasting (Section 2.1), ensemble methods (Section 2.2), dynamic ensembles
for forecasting (Section 2.3), ensemble pruning (Section 2.4), and multi-step
ahead forecasting (Section 2.5).

We start by formalizing the problem of time series forecasting in Section 2.1.
Section 2.2 outlines the motivation for using ensemble methods. There, we
refer to seminal works on this topic. Section 2.3 overviews the application of
ensemble methods to forecasting problems. We focus on a particular type of
ensemble, namely dynamic and heterogeneous ensembles. These are commonly used
to tackle forecasting tasks. Within Section 2.3, we also list several methods
that are split into three categories: windowing (Section 2.3.1), regret
minimization (Section 2.3.2), and meta-learning (Section 2.3.3). Then, we
describe ensemble pruning in Section 2.4, which is an important step for
improving forecasting performance. Finally, we list state-of-the-art
approaches to multi-step ahead forecasting in Section 2.5. These are split
into two types: single-output methods and multi-output methods. We explain the
difference between the two and our choice for the second type.

###  2.1 Time Series Forecasting

We define a univariate time series as a sequence of values
$Y=\\{y_{1},y_{2},\dots,$ $y_{n}\\}$, where
$y_{i}\in\mathcal{Y}\subset\mathbb{R}$ is the numeric value of the time series
at time $i$ and $n$ is the length of $Y$. We assume that the observations of
the series are captured at regular periods (e.g. every hour). This work
addresses multi-step ahead forecasting (also referred to as long-term
forecasting taieb2012review ) problems in univariate time series. In practice,
we aim at predicting the value of the $H$ upcoming observations of the time
series, $y_{n+1},\ldots,y_{n+H}$, where $H$ denotes the forecasting horizon.

The predictive task is formalized as an auto-regressive problem. Thus, each
observation is represented based on the past and recent values before it. This
is accomplished by reconstructing the time series using time delay embedding
Takens1981 . This method is used to transform the time series from a sequence
to a tabular format.

The transformation based on time delay embedding leads to a data set
$\mathcal{D}=\\{$($X$, $y$)$\\}$. As described before, each observation
$y_{i}$ is modeled according to past $q$ values before it:
$X_{i}=\\{y_{i-1},y_{i-2},\dots,y_{i-q}\\}$, where
$y_{i}\in\mathcal{Y}\subset\mathbb{R}$. In effect, $y$ represents the target
variable which represents the observation we want to predict, and
$X_{i}\in\mathcal{X}\subset\mathbb{R}^{q}$ represents the $i$-th embedding
vector. Then, we train a multiple regression model $f$ that can be written as
$y_{i}=f(X_{i})$. This definition is designed for one-step ahead forecasting.
In Section 2.5, we describe how this formalization is extended to the multi-
step ahead case.

###  2.2 Ensemble Methods

The No Free Lunch theorem for supervised learning wolpert1996lack  postulates
that no learning algorithm is the most appropriate for all problems. All
methods have strengths and limitations. This is the key motivation for
ensemble methods brown2005managing .

Ensemble methods aim at combining the output of several different predictive
models. These methods have been shown to perform better than individual models
in many tasks and domains of applications.

A key aspect of developing accurate ensembles is the diversity among
individual models. The models composing the ensemble should be accurate but
different from each other. Brown et al. brown2005diversity  survey several
methods to achieve this. In this paper, we focus on heterogeneous ensembles
kuncheva2004classifier . These ensembles are composed of individual models
which are trained with distinct learning algorithms, which is a typical way of
encouraging diversity.

###  2.3 Dynamic Ensembles for Forecasting

Time series forecasting is one of the tasks in which ensembles have been shown
to provide state-of-the-art performance bates1969combination ;
cerqueira2019arbitrage ; oreshkin2019n . Hibon and Evgeniou hibon2005combine
argue that combining multiple different forecasting models reduces the chance
of picking the wrong individual model.

Time series data is amenable to change due to sources of non-stationarity.
Consequently, different forecasting models usually show varying relative
performance aiolfi2006persistence . Aiolfi and Timmermann study this
phenomenon. They discovered that some models perform better than others in
some periods over a time series. This is the main motivation for using dynamic
ensembles for forecasting.

A dynamic ensemble combines the predictions of individual models using
weighted averages, where the weights change over time. The weights change to
adapt to the current process generating the time series. The main problem when
using dynamic ensembles is defining how to compute these weights at each time-
step. We overview several approaches to accomplish this. These are split into
three types: windowing approaches (Section 2.3.1), regret minimization
approaches (Section 2.3.2), and meta-learning approaches (Section 2.3.3).

####  2.3.1 Windowing Approaches

The static forecast combination approach which assigns equal weights to all
available models is a robust combination method clemen1986combining  (Simple).
Another static combination approach is a weighted average, in which the
weights are set according to the performance of models in the training data
(LossTrain). A variant of this approach is to select the model with the best
performance on the training data (Best). Notwithstanding, dynamic combination
rules are also typically used in forecasting problems.

Using past recent performance is a typical way of dynamically combining
ensembles. The idea is to give a higher weight to models that performed better
in recent observations (Window). This approach has shown promising results in
different works newbold1974experience ; cerqueira2017dynamic . A variant of
this is the method called Blast, which selects the best model in recent
observations for forecasting Blast .

An alternative to computing performance in a window of recent data is to use
forgetting factors. These take into account all past observations but give
more relevance to more recent ones. An example method that follows this
approach is AEC sanchez2008adaptive .

####  2.3.2 Regret Minimization

Combining the output of multiple models is a well-studied topic in the online
learning literature Cesa-Bianchi:2006:PLG:1137817 . Example methods include
the exponentially weighted average (EWA), the polynomially weighted average
(MLpol), and the fixed share aggregation (FS).

These approaches are designed to minimize regret. Regret is the average error
suffered relative to the best error we could have obtained. We refer the
reader to the second chapter of the seminal work by Cesa-Bianchi and Lugosi
Cesa-Bianchi:2006:PLG:1137817 .

####  2.3.3 Meta-learning

Meta-learning is another type of strategy that can be used to combine the
output of multiple models. Arguably, the most popular meta-learning approach
is Stacking, which was proposed by Wolpert wolpert1992stacked .

Cerqueira et al. cerqueira2019arbitrage  proposed a meta-learning approach
called Arbitrated Dynamic Ensemble (ADE) for dynamic forecast combinations.
ADE works by building a meta-model for each model in the ensemble. Each meta-
model is designed to model and forecast the error of the corresponding (base)
model. Then, the models in the ensemble are weighted according to the error
forecasts provided by meta-models.

###  2.4 Ensemble Pruning

Pruning is an important stage in the creation of ensembles mendes2012ensemble
. The goal of this process is to remove poor-performing models from the
ensemble.

Pruning has been shown to improve the performance of forecasting ensembles
jose2008simple ; cerqueira2019arbitrage . Jose and Winkler jose2008simple
trim a percentage of the worst models based on past performance. We also take
this approach in this work.

###  2.5 Multi-step Forecasting

Multi-step ahead forecasting denotes the process of predicting multiple
instances of a time series. This task reduces the long-term uncertainty of
time series, which is desirable across many application domains. Several works
have been devoted to multi-step ahead forecasting problems taieb2012review ;
venkatraman2015improving ; de2019batch .

In this section, we review several state-of-the-art approaches for multi-step
ahead forecasting. We split these into two types: single-output methods
(Section 2.5.1), and multi-output methods (Section 2.5.2). Taieb et al.
taieb2012review  are followed closely to describe some of these methods.

####  2.5.1 Single Output Methods

One of the most popular approaches to multi-step ahead forecasting is the
Recursive method (also known as the Iterative). Recursive works by fitting a
model $f$ for one-step ahead forecasting (c.f. Section 2.1):

$y_{i+1}=f(\\{y_{i},y_{i-1},\dots,y_{i-q+1}\\})$

To get predictions for the next $H$ observations, the model $f$ is iterated
$H$. To be more precise, in the $i$-th time-step,
$\\{y_{i},y_{i-1},\dots,y_{i-q+1}\\}$ is the input used to forecast the value
of $y_{i+1}$. Let $\hat{y}_{i+1}$ denote that prediction. Then,
$\\{\hat{y}_{i+1},y_{i},\dots,y_{i-q}\\}$ is the input vector for forecasting
$\hat{y}_{i+2}$ in the same time-step.

The Recursive strategy is known to propagate errors along the forecasting
horizon. Thus, this approach requires an accurate model specification to work
well.

Some methods have been proposed to mitigate the error propagation problem of
Recursive. The method RECNOISY, proposed by Taieb and Bontempi
taieb2011recursive , adds noise to the input data at each forecasting horizon.
The idea behind this approach is to make the forecasting model robust to error
propagation.

The method Data as Demonstrator (DaD) by Venkatraman et al.
venkatraman2015improving ; venkatraman2016improved  is also specifically
designed to deal with the error propagation problem. The idea is to use the
training set to correct errors that occur during multi-step prediction. This
is accomplished by iteratively augmenting the training set with these
corrections. After the augmentation process terminates, the Recursive approach
is applied. The augmented training data is used to train the forecasting
model.

Another popular multi-step ahead forecasting method is the Direct approach.
This strategy works by building a forecasting model for each horizon:

$y_{i+h}=f_{h}(\\{y_{i},y_{i-1},\dots,y_{i-q+1}\\}),$

where $h\in\\{1,\dots,H\\}$. Since no model is iterated along the horizon, the
Direct approach does not suffer from error propagation. Notwithstanding, this
method leads to greater computational costs because it trains a model for each
horizon. Besides, it assumes that each horizon is independent, which is not
generally true.

The method DirRec attempts to bridge the best aspects of Recursive and Direct.
Similarly to Direct, this approach builds one model for each forecasting
horizon. Moreover, for each successive horizon, the set of inputs is augmented
with the predictions from previous steps (following Recursive). This approach
is also referred to as classifier (or regression) chains in the machine
learning literature read2021classifier .

####  2.5.2 Multiple Output Methods

As the name implied, single output approaches model one horizon at a time.
Thus, they do not account for the stochastic dependency along the forecasting
horizon. This aspect is taken into account by multi-input multi-output (MIMO)
models taieb2012review . In this work, we also refer to these approaches by
multi-output models.

A multi-output model is trained on the complete forecasting horizon jointly:

$[y_{i+1},y_{i+2},\dots,y_{i+H}]=F(\\{y_{i},y_{i-1},\dots,y_{i-q+1}\\}),$

where $F$ denotes the multi-output model. Essentially, there are $H$ target
variables instead of one, and these are modeled jointly by a single model.

Taieb et al. taieb2012review  compared several approaches for multi-step ahead
forecasting, including the ones described above. They used the 111 time series
from the NN5 forecasting competition. The main conclusion is that multi-output
strategies perform better than single-output ones. For this reason, we focus
on multi-output methods in this work. Notwithstanding, the dynamic ensembles
used in this paper can also be applied with other multi-step ahead forecasting
strategies.

##  3 Materials and Methods

This section details the materials and methods used in this work. We describe
the data sets used in the experiments in Section 3.1. Then, we detail the
methodology carried out to evaluate each approach (Section 3.2). We list all
approaches in Section 3.3, including learning algorithms, dynamic combination
rules, and different strategies for weighting models along the horizon.
Finally, we detail the experimental design in Section 3.4, including the
cross-validation procedure and evaluation metric.

###  3.1 Data

The experiments encompassed 3568 time series. These were collected from the
following 7 popular databases:

  * •

electricity_nips: 370 hourly time series concerning the electricity
consumption of 370 customers yu2015high ;

  * •

nn5_daily_without_missing: 111 daily time series about cash withdrawal amounts
taieb2012review ;

  * •

solar-energy 137 hourly time series representing power produced from solar
energy lai2018modeling :

  * •

traffic_nips: 963 hourly time series concerning the occupancy rates (between
zero and one) of car lanes in San Francisco bay area freeways yu2015high ;

  * •

taxi_30min: 1214 time series related to the number of taxi rides in different
locations around New York in 30-minute windows borchert2022multi ;

  * •

m4_hourly: 414 hourly time series from the M4 time series database
makridakis2020m4  which cover different application domains;

  * •

m4_weekly: 414 weekly time series from the M4 time series database
makridakis2020m4  which cover different application domains.

Table 1 presents a summary of the data sets, including the number of time
series and their average length. The name of the data sets refers to their
identifier in the Python library gluonts.

Table 1: Summary of the data sets Name | Frequency | # Time Series | Avg. Length
---|---|---|---
electricity_nips | Hourly | 370 | 5833
nn5_daily_without_missing | Daily | 111 | 735
solar-energy | Hourly | 137 | 7009
traffic_nips | Hourly | 963 | 4001
taxi_30min | Half-hourly | 1214 | 1488
m4_hourly | Hourly | 414 | 960
m4_weekly | Weekly | 359 | 934

###  3.2 Methodology

The goal of this paper is to analyze how several dynamic ensembles perform for
multi-step ahead forecasting problems. This is accomplished by carrying out a
set of experiments.

For each time series in the case study (c.f. Section 3.1), we apply the
following procedure. The available time series is split into training and test
sets. As we describe below in Section 3.4, this process is repeated using a
cross-validation procedure.

The models in the ensemble are fit using the training data. As mentioned in
Section 2.3.2, pruning the ensemble usually leads to better forecasting
performance. In effect, we prune the ensemble and keep only the best 75% of
models in the available pool. The pruning process is carried out using a
nested validation procedure. The training set is further split into two parts:
an inner training set and a validation set. The inner training set contains
70% of the observations of the complete (outer) training set. The validation
set contains the final 30% of the complete training set. All models are fit
using the inner training set and evaluated in the validation set. The 25% of
models with the worst forecasting performance are discarded. The remaining
ones are re-fit using the complete training set.

After the pruning and re-fitting process, the ensemble is applied to the test
set. A combination rule is applied to aggregate the predictions of all
available models. The list of all combination rules applied in the experiments
is described in Section 3.3.2. Note that we test these combination rules with
different approaches concerning the estimation of weights across the horizon.
We test four strategies that are listed in Section 3.3.3. Finally, the
performance of each method is evaluated according to its performance in the
test set. The evaluation metric is described in Section 3.4.

###  3.3 Methods

This section details the methods used in the experiments. First, we describe
the regression learning algorithms used to create multi-output forecasting
models (Section 3.3.1). Then, we describe the 9 methods used to combine the
predictions of those models (Section 3.3.2). Finally, we detail four different
approaches to compute the ensemble weights across the forecasting horizon
(Section 3.3.3).

####  3.3.1 Learning Algorithms

The ensembles used in the experiments are composed of 40 individual multi-
output models. These were created using the following learning algorithms:
random forest regression, extra trees regression, bagging of decision trees,
projection pursuit regression, LASSO regression, ridge regression, elastic-net
regression, k-nearest neighbors regression, principal components regression,
and partial least squares regression. We used the implementation in the
scikit-learn pedregosa2011scikit  Python library to use these methods. Table 2
describes the different parameters used for each learning algorithm. In total,
there are 40 different learning approaches.

k-nearest neighbors regression, principal components regression, and partial
least squares regression.

Table 2: Summary of the parameters of the learning algorithms

ID | Algorithm | Parameter(s) | Value(s)
---|---|---|---
BAGGING_1 | Bagging of decision trees | No. trees | 50
BAGGING_2 | 100
RF_1 | Random Forest | {No. trees, max depth} | {50, default}
RF_2 | {50, 3}
RF_3 | {50, 5}
RF_4 | {100, default}
RF_5 | {100, 3}
RF_6 | {100, 5}
ET_1 | Extra trees regression | {No. trees, max depth} | {50, default}
ET_2 | {50, 3}
ET_3 | {50, 5}
ET_4 | {100, default}
ET_5 | {100, 3}
ET_6 | {100, 5}
KNN_1 | K-nearest neighbors | {K, weight} | {1, uniform}
KNN_2 | {5, uniform}
KNN_3 | {10, uniform}
KNN_4 | {20, uniform}
KNN_5 | {50, uniform}
KNN_6 | {1, distance}
KNN_7 | {5, distance}
KNN_8 | {10, distance}
KNN_9 | {20, distance}
KNN_10 | {50, distance}
PPR | Projection pursuit regression |  | default
LASSO_1 | LASSO regression | {regularization} | {1}
LASSO_2 | {0.75}
LASSO_3 | {0.5}
LASSO_4 | {0.25}
RIDGE_1 | Ridge regression | {regularization} | {1}
RIDGE_2 | {0.75}
RIDGE_3 | {0.5}
RIDGE_4 | {0.25}
EN | Elastic-net regression |  | default
PLS_1 | Partial least squares regression | No. components | 2
PLS_2 | 3
PLS_3 | 5
PCR_1 | Principal components regression | No. components | 2
PCR_2 | 3
PCR_3 | 5

####  3.3.2 Forecast Combination Methods

In terms of forecast combination methods, we focus on the following
approaches.

  * •

Simple: Combination rule which assigns equal weights to all models. In
practice, the predictions of the available models are combined using the
arithmetic mean;

  * •

Window: Dynamic weighted average of the predictions of the available models
newbold1974experience . The weights are computed according to the forecasting
performance in the last $\lambda$ observations;

  * •

Blast: A variant of the Window approach Blast . Instead of using past recent
performance to weigh the available models, the idea is to select the model
with the best performance in the last $\lambda$ observations;

  * •

ADE: A dynamic combination approach based on a meta-learning strategy called
arbitrating cerqueira2019arbitrage . The idea is to build a meta model (a
Random Forest) for each (base) model in the ensemble. Each meta model is
designed to predict the error of the corresponding base model. Then, the
models in the ensemble are weighted according to the error forecasts. We refer
to the work by Cerqueira et al. cerqueira2019arbitrage  for a complete read on
this method;

  * •

EWA: A dynamic combination rule based on an exponentially weighted average.
This method follows the popular weighted majority algorithm Cesa-
Bianchi:2006:PLG:1137817 ;

  * •

FS: The fixed share dynamic combination approach. This method is designed to
handle non-stationary time series;

  * •

MLpol: A dynamic combination method based on a polynomially weighted average;

  * •

Best: A baseline which selects the individual model in the ensemble with the
best performance in the training data to predict all the test instances;

  * •

LossTrain: Another baseline which weights the available models based on the
error on the training set. The weights are static and fixed for all testing
observations;

Most of these combination approaches are dynamic to cope with the non-
stationarities present in the time series. The exceptions are LossTrain and
Simple. We followed the study by Cerqueira et al. cerqueira2019arbitrage  to
set the value of the $\lambda$ parameter to 50 observations.

####  3.3.3 Weighting Approaches Over the Horizon

Dynamic ensemble methods typically assume immediate feedback. They are
designed only for one-step ahead forecasting. Thus, it is not clear how the
ensemble weights should be computed along the forecasting horizon. We study
the following approaches to estimate the weights at each time-step:

  * •

Complete Horizon (CH): The weights of individual models are estimated using
their average performance over the complete forecasting horizon;

  * •

Individual Horizon (IH): The ensemble estimates different weights for each
horizon;

  * •

First Horizon Forward (FHF): The weights computed for the first horizon are
propagated over the rest of the horizon;

  * •

Last Horizon Backward (LHB): For completeness, we include the inverse approach
to FHF. According to LHB, the weights computed for the last horizon are
propagated backward to all horizons before it.

We test all these variants with the combination rules presented above. The
exception is Simple, whose weights are static and not dependent on forecasting
performance. This leads to a total of 33 variants for analysis.

###  3.4 Experimental Design

We estimate the forecasting performance of models using a Monte Carlo cross-
validation procedure picard1984cross , which is also referred to as repeated
holdout cerqueira2020evaluating . This estimation method is applied with 10
folds. The training and test sizes of each fold are set to 60% and 10% of the
size of the input time series, respectively. Monte Carlo cross-validation
provides competitive performance estimates relative to other approaches
cerqueira2020evaluating .

We preprocess each time series as follows. We take the first differences to
remove the trend. Then, we apply time delay embedding to transform the series
for auto-regression. We set the parameter $q$ (the number of lags) to 5. This
means that the future values of a time series are modeled based on the
previous 5 observations. Finally, we set the maximum forecasting horizon to 18
observations ($H=18$).

The mean absolute error (MAE) is used as the evaluation metric. This metric
has the limitation of being scale-dependent. However, we will focus on ranks
and percentage differences to compare results across multiple time series.

##  4 Experiments

This section presents the experiments carried out to analyze the performance
of different dynamic ensembles for multi-step ahead forecasting. First, we
outline the research questions (Section 4.1). Then, we answer these questions
in turn (Section 4.2).

###  4.1 Research Questions

The experiments are designed to address the following research questions:

  * •

RQ1: How do the dynamic ensemble methods compare with each other in terms of
their rank across all data sets?

  * •

RQ2: What is the best approach for computing the weights along the forecasting
horizon?

  * •

RQ3: How does each dynamic ensemble method compare with a static ensemble that
assigns equal weights to all models?

  * •

RQ4: Does the forecasting horizon affect the results obtained?

###  4.2 Results

This section answers each research question in turn.

####  4.2.1 RQ1

Figure 1: Rank distribution of each method across all time series. Methods are
ordered by median rank.

As an exploratory analysis of the results, we start by studying the rank
distribution of each approach across the 3568 time series. An approach gets a
rank of 1 in a given data set if it shows the lowest error. The average rank
quantifies the average position of a method relative to the remaining ones.
This analysis is presented in Figure 1. Window_FHF, ADE_FHF, and Simple show
the top three best median rank. The distributions suggest that all methods
perform best in some data sets. There is not a single approach that is the
best one for all problems. This result aligns with the No Free Lunch theorem
for supervised learning put forth by Wolpert wolpert2002supervised .

Figure 2: Average rank, and respective standard deviation, of each method
across all time series.

Figure 2 shows the average (mean) rank, and respective standard deviation, of
each method across the 3568 time series. The best approaches are variants of
ADE, Window, and MLpol. On the other hand variants of Blast, Best, FS, and EWA
occupy the bottom positions. All methods show a considerable standard
deviation of rank. This corroborates the idea that no method dominates over
the rest.

####  4.2.2 RQ2

Figure 3: Average rank of each method, but aggregated by weighting strategy.

Figure 3 illustrates the average rank of all methods but aggregated by
weighting strategy. This analysis uncovers interesting outcomes. Varying the
weighting strategy (FHF, LHB, IH, CH) does not have a significant impact in
average rank. On the other hand, different combination approaches show
significantly different scores. For example, ADE shows better performance
relative to Best irrespective of the weighting strategy. The two best forecast
combination methods, ADE and Window, show similar behavior in terms of
relative scores for the different weighting strategies. Their best score is
achieved with FHF, followed by IH. Conversely, the worst combination
approaches, namely Blast, Best, EWA, and FS, have their best score when
coupled with CH.

In Figure 4, we repeated the average rank analysis but focused on the best
version of each combination method. The best weighting scheme for each
forecast combination method is detailed in Table 3.

Table 3: Best weighting scheme for each combination rule Combination Rule | Best Approach
---|---
ADE | FHF
Best | CH
Blast | CH
EWA | CH
FS | CH
LossTrain | FHF
MLpol | IH
Window | FHF
Figure 4: Average rank, and respective standard deviation, of each combination
rule applied with the respective best approach for computing the weights
across the forecasting horizon.

The best approach is Window_FHF followed by ADE_FHF. The results showed so far
suggest that most methods struggle to outperform the baseline Simple: the
approach which combines the individual models with equal weights in all
forecasting horizons. Indeed, in Figure 4 Simple shows the third-best average
rank.

####  4.2.3 RQ3

So far, the results were analyzed according to the average rank. However,
average rank ignores the magnitude of differences in predictive performance
benavoli2016should . To overcome this limitation, we also study the percentage
difference in performance between each method and a reference method. We set
Simple as the reference method, which assigns equal weights to the models in
the ensemble.

For each method m, the percentage difference is computed as follows.

$100\times\frac{\text{MAE}_{\texttt{m}}-\text{MAE}_{\texttt{Simple}}}{\text{MAE}_{\texttt{Simple}}}$

where $\text{MAE}_{\texttt{m}}$ and $\text{MAE}_{\texttt{Simple}}$ represent
the MAE of method m and Simple, respectively. Negative values denote better
performance by method m.

Figure 5: Distribution of percentage difference in MAE between each method and
Simple across all time series. Negative values denote better performance of
the respective method.

Figure 5 depicts the distribution of the percentage difference in MAE between
each method and Simple. The methods are ordered by decreasing median
percentage difference in MAE.

The order of the methods is similar to that obtained in the previous analyzes.
Only ADE_FHF and Window_FHF show a median performance difference below zero.
This indicates that Simple outperforms the other dynamic combination methods
more times than not.

Figure 5 indicates that the percentage difference is close to zero in many
problems. We analyzed the percentage difference in performance but considered
small differences to be negligible. More specifically, we consider percentage
differences below 1% to be a draw between the pair of methods under analysis.

Figure 6: Paired comparisons between each method and Simple. Each barplot
shows the probability of each method winning in blue (result below -1%),
drawing in light grey (result within [-1%, 1%]), or losing in red (result
above 1%) against Simple.

The outcome of this analysis is illustrated in Figure 6. This figure shows the
probability of each method winning in blue (percentage difference below -1%),
drawing in light grey (percentage difference below 1%), or losing in red
(percentage difference above 1%) against Simple. As an example, ADE_FHF wins
in around 40% of the datasets. A draw occurs in about 35% of the time series,
while Simple wins in around 25% of the cases. The scores obtained indicate
that only ADE_FHF, Window_FHF, MLpol_FHF, and LossTrain_FHF are competitive
with Simple. When compared with the remaining methods, the probability of
Simple winning is considerably higher than the opposite outcome.

####  4.2.4 RQ4

The analysis presented so far quantifies the performance of each method for
long-term forecasting. Specifically, predicting 18 step in advance. However,
long-term forecasts are typically less accurate than short-term ones. We
analyzed the impact of the forecasting horizon in the results obtained.

Figure 7 shows the average (median) percentage difference of each method
relative to Simple over the forecasting horizon. As before, the average is
computed across the 3568 time series.

Figure 7: Median percentage difference of each method relative to Simple in
each forecasting horizon.

The figure shows a clear trend which indicates that the methods decrease their
performance relative to Simple as the forecasting horizon increases. For t+1
(one-step ahead forecasting), 5 out of 8 combination methods outperform
Simple. But, for long-term forecasting (t+18), only two methods achieve better
performance. It is also interesting to note that there are slight changes in
the relative performance of methods. For example, for t+1 Window_FHF is only
the fourth best approach (ADE_FHF shows the best average rank). However, by
t+18, Window_FHF shows the best score.

##  5 Discussion

This paper investigates the performance of several dynamic ensembles for
multi-step ahead forecasting problems. We focus on ensembles composed of
multi-output models. Notwithstanding, the combination rules tested in this
work can also be applicable to models following a different strategy regarding
multi-step ahead prediction (c.f. Section 2.3.3).

The main motivation for this work is that the literature concerning dynamic
ensembles for forecasting is focused on one-step ahead predictions. Most
approaches assume immediate feedback from the environment to compute the error
of each method. Previous works concerning dynamic ensembles for multi-step
ahead forecasting are scarce. Thus, it is not clear how these methods should
weigh the individual models across the forecasting horizon.

We carried out a set of experiments that suggest the following conclusions.
All methods show high variability in their rank over the time series used in
the experiments. Notwithstanding, variants of ADE, and Window show the best
average rank (RQ1).

Four different approaches were tested to compute the ensemble weights across
the forecasting horizon. We did not find a significant impact on relative
performance between the different approaches. Yet, the methods with the best
overall performance maximize their scores when applied with the FHF strategy
(RQ2).

The dynamic ensembles were compared with a static ensemble which assigns equal
weights to all models. We found that only two methods systematically
outperform this approach, namely ADE_FHF, and Window_FHF (RQ3).

We also discovered that the forecasting horizon significantly affects the
relative performance of dynamic ensembles. All methods decrease their
performance relative to Simple as the forecasting horizon increases (RQ4).

Future research opportunities include analyzing how dynamic ensembles can cope
with an increasing forecasting horizon, and remain competitive with static
combination rules.

##  6 Conclusions

This work addresses multi-step ahead forecasting problems using ensembles
composed of multi-output models. A substantial amount of work have studied how
dynamic ensembles perform in forecasting tasks. Yet, the literature is focused
on one-step ahead predictions, ignoring multi-step ahead problems.

The goal of this work was to bridge this gap in the literature. The main
contribution is an extensive empirical study that analyzes how different
dynamic ensembles perform for multi-step ahead forecasting. We also study what
is the best approach for weighting the models composing the ensemble along the
forecasting horizon.

The dynamic methods ADE and Window present the best average rank. On the other
hand, other approaches struggle to outperform a static ensemble which assigns
equal weights to all models (Simple). We found that different approaches for
computing the ensemble weights along the horizon did not affect performance
significantly. Yet, propagating the weights computed for one-step ahead
predictions along the horizon maximized the performance of the best dynamic
ensembles. Another important conclusion from this work is that the forecasting
horizon has a significant impact on the results obtained. Most dynamic
ensembles struggle to outperform Simple as the forecasting horizon increases.

## References

  * (1) Aiolfi, M., Timmermann, A.: Persistence in forecasting performance and conditional combination strategies.  Journal of Econometrics 135(1), 31–53 (2006)
  * (2) Bates, J.M., Granger, C.W.: The combination of forecasts.  Journal of the Operational Research Society 20(4), 451–468 (1969)
  * (3) Benavoli, A., Corani, G., Mangili, F.: Should we really use post-hoc tests based on mean-ranks?  The Journal of Machine Learning Research 17(1), 152–161 (2016)
  * (4) Borchert, O., Salinas, D., Flunkert, V., Januschowski, T., Günnemann, S.: Multi-objective model selection for time series forecasting.  arXiv preprint arXiv:2202.08485 (2022)
  * (5) Brown, G., Wyatt, J., Harris, R., Yao, X.: Diversity creation methods: a survey and categorisation.  Information Fusion 6(1), 5–20 (2005)
  * (6) Brown, G., Wyatt, J.L., Tiňo, P.: Managing diversity in regression ensembles.  Journal of machine learning research 6(Sep), 1621–1650 (2005)
  * (7) Cerqueira, V., Torgo, L., Mozetič, I.: Evaluating time series forecasting models: An empirical study on performance estimation methods.  Machine Learning 109(11), 1997–2028 (2020)
  * (8) Cerqueira, V., Torgo, L., Oliveira, M., Pfahringer, B.: Dynamic and heterogeneous ensembles for time series forecasting.  In: 2017 IEEE International Conference on Data Science and Advanced Analytics (DSAA), pp. 242–251 (2017).  DOI 10.1109/DSAA.2017.26
  * (9) Cerqueira, V., Torgo, L., Oliveira, M., Pfahringer, B.: Dynamic and heterogeneous ensembles for time series forecasting.  In: 2017 IEEE International Conference on Data Science and Advanced Analytics (DSAA), pp. 242–251. IEEE (2017)
  * (10) Cerqueira, V., Torgo, L., Pinto, F., Soares, C.: Arbitrage of forecasting experts.  Machine Learning 108(6), 913–944 (2019)
  * (11) Cesa-Bianchi, N., Lugosi, G.: Prediction, Learning, and Games.  Cambridge University Press, New York, NY, USA (2006)
  * (12) Clemen, R.T., Winkler, R.L.: Combining economic forecasts.  Journal of Business & Economic Statistics 4(1), 39–46 (1986)
  * (13) De Stefani, J., Le Borgne, Y.A., Caelen, O., Hattab, D., Bontempi, G.: Batch and incremental dynamic factor machine learning for multivariate and multi-step-ahead forecasting.  International Journal of Data Science and Analytics 7(4), 311–329 (2019)
  * (14) Hibon, M., Evgeniou, T.: To combine or not to combine: selecting among forecasts and their combinations.  International journal of forecasting 21(1), 15–24 (2005)
  * (15) Jose, V.R.R., Winkler, R.L.: Simple robust averages of forecasts: Some empirical results.  International Journal of Forecasting 24(1), 163–169 (2008)
  * (16) Kuncheva, L.I.: Classifier ensembles for changing environments.  In: International Workshop on Multiple Classifier Systems, pp. 1–15. Springer (2004)
  * (17) Lai, G., Chang, W.C., Yang, Y., Liu, H.: Modeling long-and short-term temporal patterns with deep neural networks.  In: The 41st international ACM SIGIR conference on research & development in information retrieval, pp. 95–104 (2018)
  * (18) Makridakis, S., Spiliotis, E., Assimakopoulos, V.: The m4 competition: 100,000 time series and 61 forecasting methods.  International Journal of Forecasting 36(1), 54–74 (2020)
  * (19) Mendes-Moreira, J., Soares, C., Jorge, A.M., Sousa, J.F.D.: Ensemble approaches for regression: A survey.  Acm computing surveys (csur) 45(1), 1–40 (2012)
  * (20) Newbold, P., Granger, C.W.: Experience with forecasting univariate time series and the combination of forecasts.  Journal of the Royal Statistical Society. Series A (General) pp. 131–165 (1974)
  * (21) Oreshkin, B.N., Carpov, D., Chapados, N., Bengio, Y.: N-beats: Neural basis expansion analysis for interpretable time series forecasting.  arXiv preprint arXiv:1905.10437 (2019)
  * (22) Pedregosa, F., Varoquaux, G., Gramfort, A., Michel, V., Thirion, B., Grisel, O., Blondel, M., Prettenhofer, P., Weiss, R., Dubourg, V., et al.: Scikit-learn: Machine learning in python.  the Journal of machine Learning research 12, 2825–2830 (2011)
  * (23) Picard, R.R., Cook, R.D.: Cross-validation of regression models.  Journal of the American Statistical Association 79(387), 575–583 (1984)
  * (24) Read, J., Pfahringer, B., Holmes, G., Frank, E.: Classifier chains: a review and perspectives.  Journal of Artificial Intelligence Research 70, 683–718 (2021)
  * (25) van Rijn, J.N., Holmes, G., Pfahringer, B., Vanschoren, J.: Having a blast: Meta-learning and heterogeneous ensembles for data streams.  In: Data Mining (ICDM), 2015 IEEE International Conference on, pp. 1003–1008 (2015).  DOI 10.1109/ICDM.2015.55
  * (26) Sánchez, I.: Adaptive combination of forecasts with application to wind energy.  International Journal of Forecasting 24(4), 679–693 (2008)
  * (27) Taieb, S.B., Bontempi, G.: Recursive multi-step time series forecasting by perturbing data.  In: 2011 IEEE 11th International Conference on Data Mining, pp. 695–704. IEEE (2011)
  * (28) Taieb, S.B., Bontempi, G., Atiya, A.F., Sorjamaa, A.: A review and comparison of strategies for multi-step ahead time series forecasting based on the nn5 forecasting competition.  Expert systems with applications 39(8), 7067–7083 (2012)
  * (29) Takens, F.: Dynamical Systems and Turbulence, Warwick 1980: Proceedings of a Symposium Held at the University of Warwick 1979/80, chap. Detecting strange attractors in turbulence, pp. 366–381.  Springer Berlin Heidelberg, Berlin, Heidelberg (1981)
  * (30) Venkatraman, A., Capobianco, R., Pinto, L., Hebert, M., Nardi, D., Bagnell, J.A.: Improved learning of dynamics models for control.  In: International Symposium on Experimental Robotics, pp. 703–713. Springer (2016)
  * (31) Venkatraman, A., Hebert, M., Bagnell, J.A.: Improving multi-step prediction of learned time series models.  In: Twenty-Ninth AAAI Conference on Artificial Intelligence (2015)
  * (32) Wolpert, D.H.: Stacked generalization.  Neural networks 5(2), 241–259 (1992)
  * (33) Wolpert, D.H.: The lack of a priori distinctions between learning algorithms.  Neural computation 8(7), 1341–1390 (1996)
  * (34) Wolpert, D.H.: The supervised learning no-free-lunch theorems.  Soft computing and industry pp. 25–42 (2002)
  * (35) Yu, H.F., Rao, N., Dhillon, I.S.: High-dimensional time series prediction with missing values.  arXiv preprint arXiv:1509.08333 (2015)