# TSPP: A Unified Benchmarking Tool for Time-series Forecasting

Jan Bączek
&Dmytro Zhylko
&Gilberto Titericz
&Sajad Darabi
&Jean-Francois Puget
&Izzy Putterman
&Dawid Majchrowski
&Anmol Gupta
&Kyle Kranen
&Pawel Morkisz

###### Abstract

While machine learning has witnessed significant advancements, the emphasis
has largely been on data acquisition and model creation. However, achieving a
comprehensive assessment of machine learning solutions in real-world settings
necessitates standardization throughout the entire pipeline. This need is
particularly acute in time series forecasting, where diverse settings impede
meaningful comparisons between various methods. To bridge this gap, we propose
a unified benchmarking framework that exposes the crucial modelling and
machine learning decisions involved in developing time series forecasting
models. This framework fosters seamless integration of models and datasets,
aiding both practitioners and researchers in their development efforts. We
benchmark recently proposed models within this framework, demonstrating that
carefully implemented deep learning models with minimal effort can rival
gradient-boosting decision trees requiring extensive feature engineering and
expert knowledge. The framework is accessible on here.

NVIDIA

##  1 Introduction

Time series forecasting plays a crucial role in diverse fields like finance,
weather prediction, and demand forecasting. Accurate forecasts empower
businesses to optimize decision-making, enhance operations, and improve
overall efficiency. However, the inherent complexities within time series
data, including trends, noise, missing values, and evolving relationships
between variables, present significant challenges to achieving accurate
forecasts.

To address the challenges of time series forecasting, a diverse landscape of
methods has emerged. Gradient boosting machines (GBM) [9] have become popular
in Kaggle competitions due to their effectiveness, but require substantial
feature engineering effort. While promising, deep learning models are less
frequently used independently, primarily due to data limitations. The NN3
competition aimed to showcase the viability of neural networks (NNs) for time
series forecasting, but initial results published in 2011 [7] were
disappointing, highlighting the dominance of statistical methods. This was
further confirmed by the 2017 Web Traffic Time Series Forecasting competition
[12], where machine learning methods performed better with abundant and
homogeneous data. Similarly, the renowned M4 competition [26] demonstrated the
superiority of statistical approaches on heterogeneous data. Despite this
trend, an ensemble combining an RNN and a statistical model won the
competition [32], although it was unclear whether this signaled a shift
towards deep learning. The successor M5 competition [27], utilizing real-world
sales data from Walmart, witnessed the dominance of gradient boosting
machines, with the second-place solution leveraging an ensemble of GBM and a
neural network [35, 28]. Suggesting a potential shift towards deep learning
architectures becoming an integral part of the time series forecasting
pipeline.

The rapid evolution of diverse time series forecasting approaches has made it
difficult to compare and establish consensus on what constitutes a high-
performing deep learning model. While efforts exist to standardize datasets
for benchmarking these models [11], comparisons often overlook crucial aspects
across the entire machine learning lifecycle, leading to ambiguous results
that offer limited guidance for model improvement. For instance, established
techniques for training neural networks, such as curriculum learning,
exponential moving average, and drift estimation, are readily available but
largely ignored in the current literature. These studies frequently adopt a
model-centric perspective, focusing solely on the forecasting algorithm
itself. Conversely, our work aims to address this gap by comprehensively
considering all moving parts within the machine learning lifecycle. To achieve
this, we propose a unified benchmarking framework that standardizes the
evaluation of forecasting models. This framework facilitates the
implementation and comparison of the most influential time series forecasting
models, enabling a more comprehensive and objective assessment of their
performance.

Our contribution are as follows: 1) We develop and open-source a standardized
framework for time series forecasting benchmarking. The framework has modular
components that enable fast and easy integration of datasets, models, and
training techniques. 2) Using the proposed framework, recent popular deep
learning models for time-series forecasting are re-implemented and improved
upon, assesing their preformance on commonly used datasets across these
methods. 3) Various observations are outlined from the extensive computational
resources allocated to tuning these models and comparing them. The complete
code and detailed reports for our experiments here.

##  2 Related work

The majority of research in time series forecasting has revolved around
competitions [18], where participants utilize diverse methodologies. These
methods fall into three main categories: 1) statistical methods, such as ARIMA
[3], 2) classical machine learning algorithms, including support vector
machines [17], gradient boosting machines [6], and 3) neural network-based
approaches, exemplified by TFT [24] (a transformer-based model), NHITS &
NBEATS [5, 28] (which utilize a hierarchical approach with multi-layer
perceptrons), and MTGNN [36] (another model leveraging the power of Graph
Neural Networks (GNNs) by incorporating graph structure). However, the use of
datasets with diverse settings in the methods published hinders accurate
performance evaluation and makes it difficult to identify a clear frontrunner.

Despite efforts to standardize benchmarking practices, the variability of
datasets across methods hinders meaningful comparisons. In [14], authors
highlight the limitations of dataset comparisons due to diverse features like
size, number of features, and noise level. Similar observations are echoed by
benchmarks like Monash Time Series Forecasting Benchmark[11], which curates
datasets specifically for model comparison. However, initiatives like [20]
focus on classifying model operating modes (univariate/multivariate,
local/global) as the basis for comparison. In contrast, [22] takes a model-
centric approach, exploring architectural changes in neural networks. Their
findings suggest the absence of a universal best architecture, highlighting
the data-dependence of model performance. While these studies acknowledge the
lack of standardization, their focus is either model-centric or data-centric
which hinders a comprehensive understanding of the interplay between the two.
Establishing a standardized framework that considers both aspects
simultaneously is crucial for fostering meaningful comparisons and
advancements in time series forecasting.

While deep learning models have gained popularity, recent works challenge
their effectiveness against simpler, more interpretable methods [8, 39, 25].
However, we argue that deep learning models haven’t received their due credit.
With careful consideration throughout the entire machine learning life-cycle
and appropriate optimization techniques, we believe these models can achieve
further improvements beyond their initially claimed capabilities.

##  3 Time-series Forecasting Framework

Figure 1: An overview of the pipeline standardization considered to guide
model development and comparison. The framework is comprised of three main
components, data, training, inference, which a tuner has access to for
adequate exploration to ensure reliable model training.

Time series encompass a sequence of data points indexed by time, and
forecasting is the process of predicting future values from past historical
data. Albeit the task can be defined simply in words, the machine learning
life cycle for such models are quite complex and an end-to-end machine
learning pipeline covering this from data curation to deployment is highly
needed. An overview of the framework proposed is depicted in Figure 1 which
comprises of the following components:

  * •

Data: Decisions made for improving the data and shaping it for model training.
This includes data selection, curation and cleaning.

  * •

Training: Decisions made with respect to the input data and model in order to
improve downstream performance. This includes model design tailored to input
data, optimization, and selection criterion.

  * •

Inference: Decisions made with respect to deployment and making predictions on
unseen data inputs.

  * •

Tuner: Once the above components have been defined, this component selects the
top configuration and uses it for post-deployment model monitoring,
retraining, and uncertainty quantification.

###  3.1 Task formulation

In time series forecasting, we are provided with a time series
$X=(x_{1},\dots,x_{T})$, $x_{i}\in\mathbb{R}$ and a pre-specified prediction
horizon $h\in\mathbb{N}$ as the objective. Then a quality metric
$\mathcal{L}:\mathbb{R}^{h}\times\mathbb{R}^{h}\to\mathbb{R}$ is defined to
evaluate the prediction quality of the forecasting model
$M:\mathbb{R}^{T}\to\mathbb{R}^{h}$. The objective is to minimize the
expectation of this metric $\mathbb{E}[\mathcal{L}(M(X),X^{+})]$, where
$X^{+}=(x_{T+1},\dots,x_{T+h})$ is a sequence representing the ground truth
future values.

###  3.2 Data Aspects

Data is a critical component machine learning applications and shapes most
decisions made within the pipeline. Real-world time series datasets often
contain missing values, noise, and do not follow regularities seen in many
benchmark datasets. Most prior works consider the data to be static and
improvements are solely sought at the model/algorithmic level. We instead
consider data modifications as part of the overall modelling effort when
making comparisons.

Extending the time series definition, the dataset $\mathcal{D}=(X,C)$ consists
of a set of $n$ series $X=\\{X^{1},\dots,X^{n}\\}$, with
$X^{i}=(x_{1}^{i},x_{2}^{i},\dots,x_{T}^{i})$, where $T$ is the length of the
time series. Accompanying $X$ is a set of covariates
$C=\\{C^{1},\dots,C^{n}\\}$, $C^{i}=(c_{1}^{i},c_{2}^{i},\dots,c_{T}^{i})$.
Here, the time series and covariates can be vectors of arbitrary dimension,
i.e. $x_{t}^{i}\in\mathbb{R}^{d}$ and $c_{t}^{i}\in\mathbb{R}^{k}$, where $d$
and $k$ are the dimensions of the time series and covariates, respectively.
For the case when $d=1$ the time series is considered univariate and when
$d>1$ it is considered multivariate.

In this way of defining the dataset multiple assumptions can be violated 1)
Series in $X$ are independent of each other \- in this case a static data $S$
may be present within the dataset $\mathcal{D}=(X,C,S)$ which can contain
correlation information between the series (e.g. spatial location of the
sensors which were used for data collection) 2) Series in $X$ are of the same
length - in this case the length $T$ can vary between series 3) Prior
knowledge about the future is lacking where in fact it could indeed be
contained within the covariates or a variable. \- in this case covariates $C$
are split into two disjoint sets: $C_{o}$ which contains covariates that we
have no future knowledge on and $C_{k}$ where we do have such knowledge. For
example we may know an event occurring in the future in advance or we have
full control over a variables.

Using the dataset $\mathcal{D}$, the data-split module generates the splits
$(D_{\textrm{train}},D_{\textrm{val}},D_{\textrm{test}})$ corresponding to the
train, validation and test sets, respectively, and is fixed amongst methods.
Beyond this point, how the data is shaped can dramatically effect how it can
be consumed by downstream models. For example, an underlying method may be
sensitive to the dynamic range and may require normalization, or on other hand
cannot accommodate covariates as additional input, hence these must be
excluded. The pre-processing function $p$ implements a variety of
functionalities, such as filtering, outlier removal, imputation, label
correction, feature selection, and so on. It may also be the case that at each
pre-processing step multiple viable alternatives can be applied to the data to
improve its quality for a particular model and should be exposed to a tuner
when evaluating the overall performance. In general $p$ is set to satisfy
certain modeling conditions and downstream assumption.

###  3.3 Training Aspects

Once the dataset is prepared, we must determine how much exploration is
necessary for model training and how to ensure reliable training. In general,
a forecasting model $M$ is a function which takes a series of observed values
as well as the covariates as inputs, and predicts a series of forecast
outputs:

$\displaystyle M(S,X,C)\mapsto\hat{X}^{+}\in\mathbb{R}^{d\times h}$ (1)

Due to finite computational resources, it’s common to limit the data used for
training. This is achieved through a look-back window, which varies across
models. A standard assumption is that the most recent values of a series
contain the most relevant prior for future predictions. A fixed-length look-
back window $l$ is used, starting at the index where future predictions are
desired. More concretely, let $X^{\prime}$ be the set of series from $X$
trimmed to the last $l$ observations, i.e.
$X^{\prime}:=\\{(x^{i}_{T-l+1},\dots,x^{i}_{T})|(x^{i}_{1},\dots x^{i}_{T})\in
X\\}$ and analogously define $C_{k}^{\prime}$ and $C_{o}^{\prime}$. Following
this assumption, the history is approximated using the look back window:

$\displaystyle M(S,X,C_{k},C_{o})\approx
M(S,X^{\prime},C_{k}^{\prime},C_{o}^{\prime})$ (2)

Here, the model $M$ can either operate at a global level in which model
parameters are estimated using all available time series in the dataset or on
the other extreme at a local level where a single time series in the dataset
is used to estimate model parameters [11]. The choice of $M$ may also result
in different training parameters which can include optimization parameters,
initialization, objectives, whether to apply curriculum learning, exponential
moving averaging (EMA) and so on. These sets of parameters must also be
considered jointly with model exploration as will be discussed in Sec, 3.5.
Thus it is important to consider parameter ranges that control model
complexity and training complexity.

The quality of the model $M$ is measured by evaluating model’s predictions
given ground truth values using a loss function

$\displaystyle\mathcal{L}:\mathbb{R}^{d\times h}\times\mathbb{R}^{d\times
h}\to\mathbb{R}$ (3)

Here, often we assume there is only one metric that is useful to serve for
both as a training objective and for model selection, which may not be the
case. The stopping criterion used to terminate training early can also effect
final model selection. Having pre-specified dataset decisions, in the training
portion the objective is then to find the best data-informed model that has
the lowest loss on a held out dataset.

###  3.4 Inference Aspects

Successful models often rely on certain assumptions about the series they aim
to predict, such as stationarity or no gaps in the data. Generally the
forecast window predicted by the model is set by the downstream task,
resulting in predictions $X^{i}_{+}=(x^{i}_{T+1},\dots,x^{i}_{T+h})$. Similar
to how a preprocessing function $p$ was defined in Sec. 3.2 to transform the
raw series, a post-processing function $p^{\prime}$ operates on the forecast
window to reverse certain aspects of the preprocessing decisions steps (for
example, normalization) and incorporate domain expert knowledge about how the
predicted window should be interpreted (for example, product sales cannot be
negative and setting a lower bound).

Successful models often rely on certain assumptions about the data they aim to
predict, such as stationarity and no gaps in the series. These assumptions
simplify the modeling process and allow for more accurate predictions.
However, real-world data rarely adheres perfectly to such assumptions. To
compensate for this, similar to preprocessing techniques are often employed to
transform the raw data into a format that satisfies the model’s output
requirements.

The forecast window predicted by the model is typically determined by the
downstream task it serves. For example, a sales forecasting model might be
designed to predict product demand for the next month, resulting in a forecast
window $X^{i}_{+}=(x^{i}_{T+1},\dots,x^{i}_{T+h})$, where $T$ is the current
time and $h$ is the prediction horizon.

While preprocessing prepares the data for modeling, a post-processing step
ensures the predicted values are interpretable and aligned with domain
knowledge. This post-processing function, denoted as $p^{\prime}$, operates on
the forecast window and performs several key tasks. It reverses certain
aspects of the preprocessing steps, such as denormalization, to return the
predicted values to their original units. Additionally, it incorporates
domain-specific knowledge to ensure the predictions are realistic and
consistent with real-world constraints. For instance, a post-processing
function for product sales might enforce a non-negativity constraint,
preventing the model from predicting negative sales figures.

By combining preprocessing and post-processing steps, we can leverage the
strengths of machine learning models while ensuring the predictions are
interpretable and applicable to the specific domain.

###  3.5 Tuner Hyperparameter Selection

In the aforementioned sections we have considered a series of hyperparameter
configurations spaces which must be explored to ensure reliable model training
and comparison. There is a little incentive for authors to reproduce other
paper’s results in their training setup, and this in turn may put the
reference models at a disadvantage. To ensure a fair comparison models should
be implemented within the same framework with access to the same modifications
across the pipeline as well as hyperparameter optimization (HPO) algorithms.
Let $\Lambda=\Lambda_{1}\times\Lambda_{2}\times\cdots\Lambda_{n}$ define the
space of all hyperparameters, and a vector in this hyperparameter space be
denoted by $\lambda\in\Lambda$. Within this set of hyperparameters a series of
ranges for each are considered, and is denoted by $\mathcal{A}$. The
hyperparameters consists of both model specific $\lambda_{M}$ and dataset or
training specific parameters $\lambda_{P}$, that is
$\lambda=[\lambda_{M};\lambda_{P}]$. Then the objective of the tuner is to
find

$\displaystyle\lambda_{top}=\operatorname*{arg\,min}_{\lambda\in\mathcal{A}}\mathbb{E}_{(D_{train}\cup
D_{val},D_{test})\sim\mathcal{D}}V(\mathcal{L},\lambda,D_{train}\cup
D_{val},D_{test}),$ (4)

where $V$ measures the loss of the framework set with hyperparameters
$\mathcal{A}_{\lambda}$. When running this tuner for each model, a fixed
budget $B$ for finding the best configuration defined in terms of compute
hours is used when running the pipeline on the same system and environment. In
order to obtain the final model $M$ a two-fold optimization process is used:
1) a hyperparameter $\lambda$ is sampled from $\mathcal{A}$ 2) an optimization
procedure is applied on model weights using $D_{train}\cup D_{val}$. The model
optimization procedure terminates when it meets a stopping condition evaluated
on a validation set. 2) step 1 is repeated until the HPO budget $B$ is
exhausted. In the next section we explain how models protability and
generalizability is assessed. The process is depicted with more detail in the
flow chart depicted in Figure 2 and Section B of the Appendix.

###  3.6 Model Selection & Evaluation

In order to compare models, simply evaluating a single trained model is not
adequate due to the inherent variance in the performance for such tasks. For
example, simply changing random seeds will result in an accuracy spread which
may be deemed statistically significant as done in [30]. Instead, the top
configuration found by the tuner is used to initialize the family of
parametric models $q(\theta|\mathcal{D},\lambda_{top})$, from which model are
sampled and trained. The trained model is denoted as
$p(y|\mathcal{D},\theta))$. We sample $K$ such models $\theta_{m}\sim
q(\theta|\mathcal{D},\lambda_{top})$, and evaluate using the loss:

$\displaystyle\mathbb{E}[\mathcal{L}(p(y|\mathcal{D};\theta),y)]\approx\frac{1}{K}\sum_{m=1}^{K}\mathcal{L}(p(y|\mathcal{D},\theta_{m}),y)$
(5)

Similarly, the variance of the losses is computed and represented as
$\sigma^{2}=\text{Var}(\mathcal{L}(p(y|\mathcal{D};\theta),y))$. To compare
two models a statistical test is formulated to determine if model $M_{1}$ has
a lower expected loss than the other model $M_{2}$. Let $q_{1}$ be the
parametric family of models for $M_{1}$, along with its expected loss $E_{1}$,
and variance $\sigma_{1}^{2}$. Similarily, $q_{2}$ be the parametric model for
$M_{2}$, $E_{2}$ the expected loss and $\sigma_{2}^{2}$ its variance. The null
hypothesis posits $h_{0}$: there is no significant difference between the
expected losses of the models. Then the $t$ value is given by

$\displaystyle
t=\frac{E_{1}-E_{2}}{\frac{1}{K}\sqrt{\sigma_{1}^{2}+\sigma_{2}^{2}}}$ (6)

By specifying a desired significance level $\alpha$, the null hypothesis is
rejected if $t>t_{crit}$ where $t_{crit}$ is taken from the $t$-distribution
table. In the case where there is a significance difference it can easily be
determined which model is better by picking the one with lower expected loss.
This method of comparing is robust to random initialization and training
settings, which also better reflects the portability of the model across
datasets.

##  4 Experiments

###  4.1 Dataset Details

Table 1: Summary of datasets. | Electricity | PEMS-BAY | M5 | Wikitraffic
---|---|---|---|---
Num. Series ($n$) | 370 | 325 | 30490 | 145063
Sample Rate | 1h | 5 min | 1 day | 1 day
Look-back Window ($l$) | 168 | 12 | 28 | 64
Horizon Length ($h$) | 24 | 12 | 28 | 64
Covariates ($C$) | ✗ | ✓ | ✗ | ✗
Static ($S$) | ✗ | ✓ | ✓ | ✓

Time series model benchmarking can be challenging due to the diversity of data
characteristics. Unlike in natural language processing (NLP) or computer
vision (CV), where data tends to have more consistent characteristics, time
series data can vary significantly in terms of its structure and
complexity[15]. This makes it difficult to accurately compare the performance
of different models and methodologies, as there may not be a one-size-fits-all
approach that works well across all datasets. The data-centric approach taken
considers the variations in how the data is transformed as an additional
tunable parameter.

In our experiments, four widely used datasets in time series forecasting
benchmarking are considered and summarized in Table. 1: the Electricity, PEMS-
BAY, Wiki-Traffic Prediction, and the M5 dataset. For more details on each
dataset refer to Appendix C

For these datasets we consider various normalization preprocessing steps $p$:
$\textrm{Z}(x):=\frac{x-\mu}{\sigma}$ where $\mu$ is the mean and $\sigma$ the
standard deviation, log transformation $\textrm{log}(x):=\text{log}(x+1)$, and
the composite of these
$\textrm{Z}(\textrm{log}(x)):=\textrm{Z}(\textrm{log}(x+1))$.

Feature-engineering is also considered as a preprocessing step on the input
which is detailed in Appendix C and depends on the dataset properties. For
preprocessing functions that require fitting on an underlying data, they are
done so on the $\mathcal{D}_{\text{train}}$. The validation
$\mathcal{D}_{\text{val}}$ and test $\mathcal{D}_{\text{test}}$ data splits
are standardized based on the training split of each dataset.

###  4.2 Models Details

For our experiments, the considered models contain well-established reference
implementations and demonstrate state-of-the-art performance on the datasets
listed in the previous section. These models are re-implemented and optimized
within the proposed framework, achieving performance that either matches or
surpasses the reference implementations. The summary of these models is
presented in Table 2.

Table 2: Table summarizing models used in experiments, and their operating mode. In the table, RNN $:=$ Recurrent Neural Networks, CNN $:=$ Convolutional Neural Networks, GNN $:=$ Graph Neural Networks. Model | Summary | Univariate | Multivariate | Inputs
---|---|---|---|---
|  |  |  | covariates | static
ARIMA [3] | Statistical Model | ✓ | ✗ | Yes | No
DeepAR [31] | RNN | ✓ | ✗ | Yes | Yes
TFT [24] | Transformer+RNN | ✓ | ✓ | Yes | Yes
N-BEATS [28] | MLP | ✓ | ✗ | No | No
N-HITS [5] | MLP | ✓ | ✗ | No | No
MTGNN [36] | CNN + GNN | ✗ | ✓ | Yes | Yes
DCRNN [23] | GNN + RNN | ✗ | ✓ | Yes | Yes
XGBoost [6] | Boosting Decision Trees | ✓ | ✓ | Yes | Yes

Each model has a corresponding set of hyperparameters $\lambda_{m}$ which
pertain to the hidden dimension sizes, number of layers, etc, the complete
list and the ranges considered for each dataset and model is specified in
Appendix D.1. Note that amongst these models is also a strong baseline
XGBoost, often used to compete in competitions with additional feature
engineering pre-processing steps.

###  4.3 Evaluation Metrics

The following metrics are considered when evaluating each model:

  * •

Mean Absolute Error (MAE): $\frac{1}{n}\sum^{n}|X^{+}-\hat{X}^{+}|$

  * •

Root Mean Squared Error (RMSE):
$\sqrt{\frac{1}{n}\sum^{n}(X^{+}-\hat{X}^{+})^{2}}$

  * •

Symmetric mean absolute percentage error (SMAPE):
$\frac{1}{n}\sum^{n}\frac{|X^{+}-\hat{X}^{+}|}{(|X^{+}|+|\hat{X}^{+}|)/2}$ *
100

Each metric is computed pointwise regardless of the forecast window $h$. Many
publications use MAPE[40, 33, 36, 37] or q-risk[24, 31] metrics, and since
MAPE is unbounded (if the ground truth approaches 0) we chose to use SMAPE
instead.

###  4.4 Training Loop & Optimizers

Irrespective of models and datasets, there are techniques that can be
leveraged when training time series forecasting models to improve performance.
In particular the choice of optimizer, for example Adam [21] and its
hyperparameters such as initial learning rate, momentum, etc. Techniques that
directly modify the training loop like exponential move averaging (EMA), where
an offline model keeps track of the model-weights using a exponential moving
average similar to [13, 19] or curriculum learning (CL) [2] are also
considered to be tunable.

###  4.5 Tuner & Compute Details

The hyperparameters considered for the tuner to optimize are summarized in
Table 3 and their corresponding ranges setting the subspace of valid
hyperparameters $\mathcal{A}$.

All our experiments are conducted on DGX1-32G system with 8xV100 GPUs. The
computational budget $B$ is set to one DGX week i.e. 1344 GPU hours for each
model-dataset combination. For evaluating models, the top configuration is
trained $K=256$ times for Electricity and PEMS and $K=64$ for M5 and
Wikitraffic datasets, and then evaluated on the test set using the metrics
described in the section 4.3.

Table 3: Summary of hyperparameters considered specific hyper-parameters for various models and a more extensive list is provided in the Appendix D.1.  Component | Parameters | Values
---|---|---
|  | Electricity | PEMS-BAY | WikiTraffic | M5
Dataset | Norm | Z | Z | [Z, log1p, Z(log1p)] | [log1p, Z(log1p)]
Covariates | F | T/F | F | F
Static | F | T/F | T/F | T/F
Training | CL | T/F | T/F | T/F | T/F
EMA | T/F | T/F | T/F | T/F
EMA decay | (0.9, 0.9999) | (0.9, 0.9999) | (0.9, 0.9999) | (0.9, 0.9999)
Optimizer | [Adam] | [Adam] | [Adam] | [Adam]
| Criterion | [L1, GLL, MSE] | [L1, GLL, MSE] | [L1, GLL, MSE] | [L1, GLL, MSE, Tweedie]
| Learning Rate | (1e-5, 1e-2) | (1e-5, 1e-2) | (1e-5, 1e-2) | (1e-5, 1e-2)
| Dropout | (0, 0.5) | (0, 0.5) | (0, 0.99) | (0, 0.9)
Inference | post-processing | - | - | clip to zero and round | clip to zero and round

##  5 Results

Table 4 summarizes the empirical performance of the different models on the
target datasets within the proposed framework and experimental setup. Notably,
deep learning models exhibit strong competitive performance against XGBoost
models, which leverage selective feature engineering and expert knowledge, on
several individual datasets. This finding challenges the assumption that
traditional, feature-engineered models consistently outperform deep learning
approaches in time series forecasting tasks. Moreover, the results showcase
the absence of a single model that reigns supreme across all datasets,
highlighting the importance of considering data-specific characteristics when
selecting a forecasting model.

Building upon the individual model performance, we explored ensemble learning,
a widely recognized strategy to boost forecasting accuracy in time series
tasks [10]. To this end, we trained up to 256 instances of each model with
optimized hyper-parameters derived from different random seeds, resulting in a
diverse ensemble of predictors. The final prediction was obtained through
decision fusion via prediction averaging, leveraging the collective wisdom of
the ensemble; Table 5 presents the summary of results for the ensemble models.
Through the experiments conducted on these models a series of observations can
be made:

Metrics: Different models can be superior depending on the chosen metric,
emphasizing the importance of considering a collective of metrics for model
comparison. We excluded the TDI metric as the tuner was unable to effectively
utilize it for model discrimination. This highlights the need of including
multiple evaluation metrics that are both informative and readily usable by
automated optimization tools.

Effect of Batch Size: Contrary to domain-specific recommendations favoring
smaller or larger batch sizes, our findings suggest a tighter coupling between
batch size and the specific dataset and model. For example, DeepAR prefers
smaller batch sizes, while N-BEATS exhibits insensitivity to this
hyperparameter, and TFT requires a specific range.

CL & EMA: In contrast to previous claims, the investigated training techniques
did not consistently improve performance across datasets. Notably, curriculum
learning (CL) provided negligible improvement despite efforts to encourage the
tuner to enable it. Exponential moving averaging (EMA) produced mixed results,
potentially due to the noise inherent in the forecasting task and the
detrimental effects of EMA smoothing.

Context Length (l): The optimal context length exhibits a strong dependency on
the specific model and dataset combination. This implies that a single,
universal context length recommendation is unlikely to be effective across all
models and datasets. N-BEATS and NHITS, for example, demonstrate poor
performance when provided with limited context (less than 256 data points).
However, their performance significantly improves with access to a richer
context (greater than 256 data points).

Weak Models as Strong Baselines in Ensembles: DeepAR, despite its subpar
performance when trained individually (e.g., on the WikiTraffic dataset),
exhibits a remarkable transformation when incorporated into ensembles. As
shown in Table 5, ensembles utilizing DeepAR achieve significantly superior
performance compared to other models. This intriguing finding underscores the
potential of seemingly weak models to act as powerful baselines within
ensemble configurations. By leveraging the diverse strengths of individual
models, ensembles can collectively achieve superior performance, even when
some components under-perform in isolation.

Table 4: This table presents the mean error performance of various forecasting
models, along with their corresponding standard deviations. Lower values
indicate better performance.

| Electricity | PEMS-BAY | WikiTraffic | M5
---|---|---|---|---
Model | MAE$\downarrow$ | RMSE$\downarrow$ | SMAPE$\downarrow$ | MAE$\downarrow$ | RMSE$\downarrow$ | SMAPE $\downarrow$ | MAE$\downarrow$ | RMSE$\downarrow$ | SMAPE$\downarrow$ | MAE $\downarrow$ | RMSE$\downarrow$ | SMAPE$\downarrow$
TFT | 41.72 $\tiny\pm$ 0.714 | 352.85 $\tiny\pm$ 12.74 |  8.65 $\tiny\pm$ 0.0572 | 1.64 $\tiny\pm$ 0.024 | 3.83 $\tiny\pm$ 0.072 | 3.34 $\tiny\pm$ 0.0573 |  308.4 $\tiny\pm$ 6.683 | 12343.06 $\tiny\pm$ 817.28 |  38.83 $\tiny\pm$ 0.48 | 1.07 $\tiny\pm$ 0.01 | 2.60 $\tiny\pm$ 0.09 | 146.88 $\tiny\pm$ 0.55
N-BEATS | 43.15 $\tiny\pm$ 0.227 | 341.65 $\tiny\pm$ 3.29 | 9.45 $\tiny\pm$ 0.0527 | 1.97 $\tiny\pm$ 0.001 | 4.67 $\tiny\pm$ 0.007 | 3.99 $\tiny\pm$ 0.0027 | 322.66 $\tiny\pm$ 14.38 |  12191.58 $\tiny\pm$ 1504.04 | 40.38 $\tiny\pm$ 0.41 | 1.02 $\tiny\pm$ 6e-4 | 2.46 $\tiny\pm$ 4e-3 | 159.90 $\tiny\pm$ 0.25
N-HITS | 42.30 $\tiny\pm$ 0.093 | 340.56 $\tiny\pm$ 1.51 | 9.25 $\tiny\pm$ 0.0146 | 1.95 $\tiny\pm$ 0.001 | 4.62 $\tiny\pm$ 0.005 | 3.92 $\tiny\pm$ 0.0023 | 332.30 $\tiny\pm$ 8.44 | 13554.27 $\tiny\pm$ 855.14 | 39.86 $\tiny\pm$ 0.18 | 1.01 $\tiny\pm$ 5e-4 | 2.46 $\tiny\pm$ 3e-3 | 159.91 $\tiny\pm$ 0.10
MTGNN | 43.00 $\tiny\pm$ 0.537 | 342.64 $\tiny\pm$ 6.56 | 9.44 $\tiny\pm$ 0.1032 | 1.55 $\tiny\pm$ 0.005 | 3.61 $\tiny\pm$ 0.014 | 3.11 $\tiny\pm$ 0.0107 | $-^{*}$ | $-^{*}$ | $-^{*}$ | $-^{*}$ | $-^{*}$ | $-^{*}$
DCRNN | $-^{**}$ | $-^{**}$ | $-^{**}$ |  1.54 $\tiny\pm$ 0.007 |  3.57 $\tiny\pm$ 0.0187 |  3.08 $\tiny\pm$ 0.0138 | $-^{**}$ | $-^{**}$ | $-^{**}$ | $-^{**}$ | $-^{**}$ | $-^{**}$
DeepAR | 46.98 $\tiny\pm$ 1.53 | 367.54 $\tiny\pm$ 28.86 | 10.27 $\tiny\pm$ 0.1546 | 1.73 $\tiny\pm$ 0.011 | 3.88 $\tiny\pm$ 0.0279 | 3.51 $\tiny\pm$ 0.0224 | 9108.69 $\tiny\pm$ 47456.31 | 10633457.81 $\tiny\pm$ 55977165.11 | 39.47 $\tiny\pm$ 0.38 | 1.01 $\tiny\pm$ 3e-3 | 2.47 $\tiny\pm$ 0.02 | 158.88 $\tiny\pm$ 0.65
XGBoost |  40.67 $\tiny\pm$ 0.03 |  327.68 $\tiny\pm$ 0.59 | 8.78 $\tiny\pm$ 2e-3 | 1.66 $\tiny\pm$ 2e-4 | 3.87 $\tiny\pm$ 5e-4 | 3.35 $\tiny\pm$ 4e-4 | 317.58 $\tiny\pm$ 0.46 | 12801.7 $\tiny\pm$ 85.2 | 38.93 $\tiny\pm$ 0.01 |  1.002 $\tiny\pm$ 4e-4 |  2.45 $\tiny\pm$ 2e-3 |  127.91 $\tiny\pm$ 0.086
Notes: Experiments not performed due to memory requirements.
DCRNN requires static graph which is available only for PEMS-BAY |

Table 5: Error of different forecasting models using an ensemble of models the
mean performance is reported for each. Ensemble of 256 model was used for
Electricity and PEMS-BAY and ensemble of 64 models for Wikitraffic and M5.

| Electricity | PEMS-BAY | WikiTraffic | M5
---|---|---|---|---
Model | MAE$\downarrow$ | RMSE$\downarrow$ | SMAPE$\downarrow$ | MAE$\downarrow$ | RMSE$\downarrow$ | SMAPE$\downarrow$ | MAE$\downarrow$ | RMSE$\downarrow$ | SMAPE$\downarrow$ | MAE$\downarrow$ | RMSE$\downarrow$ | SMAPE$\downarrow$
TFT | 40.40 | 340.04 | 8.43 | 1.55 | 3.57 | 3.13 | 297.95 | 11896.13 | 37.86 | 1.07 | 2.58 | 146.83
N-BEATS | 42.57 | 338.70 | 9.27 | 2.03 | 4.46 | 4.11 | 303.87 | 10883.19 | 39.55 | 1.01 | 2.45 | 159.61
N-HITS | 41.71 | 338.17 | 9.06 | 1.94 | 4.60 | 3.91 | 320.81 | 12908.48 | 39.45 | 1.01 | 2.46 | 159.75
MTGNN | 41.71 | 331.11 | 9.20 | 1.52 | 3.54 | 3.05 | $-^{*}$ | $-^{*}$ | $-^{*}$ | $-^{*}$ | $-^{*}$ | $-^{*}$
DCRNN | $-^{**}$ | $-^{**}$ | $-^{**}$ | 1.47 | 3.40 | 2.95 | $-^{**}$ | $-^{**}$ | $-^{**}$ | $-^{**}$ | $-^{**}$ | $-^{**}$
DeepAR | 42.34 | 329.03 | 9.58 | 1.67 | 3.68 | 3.36 | 292.07 | 11412.32 | 37.94 | 1.00 | 2.43 | 158.42
XGBoost | 40.44 | 325.58 | 8.73 | 1.65 | 3.85 | 3.33 | 317.31 | 12793.09 | 38.91 | 1.002 | 2.44 | 127.85
Notes: *Experiments not performed due to memory requirements.
**DCRNN requires static graph which is available only for PEMS-BAY |

##  6 Conclusion

The lack of standardization within the time series forecasting domain presents
a significant impediment to progress. This issue manifests in inconsistent
reporting practices, even when using the same dataset, as evidenced by
discrepancies in metrics reported across various studies [24, 36, 38].
Additionally, the frequent reporting of results obtained with suboptimal
hyper-parameters further obscures the true performance potential of different
models, making it difficult to discerningly judge their relative strengths and
weaknesses [24]. This work addresses these challenges by proposing a
standardized workflow and framework that facilitates seamless integration of
diverse datasets and training/inference logic. Furthermore, this framework
transparently exposes decision-making processes throughout the workflow,
enabling optimization tools to guide researchers and developers in their
exploration. Our empirical evaluations reveal that deep learning methods
demonstrate competitive performance against XGBoost models, challenging the
traditional perception that feature-engineered approaches consistently
dominate this domain. By promoting standardization and transparency, this work
paves the way for more robust and reproducible research in the field of time
series forecasting.

## References

  * [1]
  * Bengio et al. [2009] Yoshua Bengio, Jérôme Louradour, Ronan Collobert, and Jason Weston. 2009\.  Curriculum learning. In _Proceedings of the 26th annual international conference on machine learning_. 41–48.
  * Box [1970] George EP Box. 1970\.  GM Jenkins Time Series Analysis: Forecasting and Control.  _San Francisco, Holdan-Day_ (1970).
  * Cao et al. [2020] Defu Cao, Yujing Wang, Juanyong Duan, Ce Zhang, Xia Zhu, Congrui Huang, Yunhai Tong, Bixiong Xu, Jing Bai, Jie Tong, et al. 2020\.  Spectral temporal graph neural network for multivariate time-series forecasting.  _Advances in neural information processing systems_ 33 (2020), 17766–17778.
  * Challu et al. [2022] Cristian Challu, Kin G Olivares, Boris N Oreshkin, Federico Garza, Max Mergenthaler-Canseco, and Artur Dubrawski. 2022.  N-hits: Neural hierarchical interpolation for time series forecasting.  _arXiv preprint arXiv:2201.12886_ (2022).
  * Chen et al. [2015] Tianqi Chen, Tong He, Michael Benesty, Vadim Khotilovich, Yuan Tang, Hyunsu Cho, Kailong Chen, Rory Mitchell, Ignacio Cano, Tianyi Zhou, et al. 2015\.  Xgboost: extreme gradient boosting.  _R package version 0.4-2_ 1, 4 (2015), 1–4.
  * Crone et al. [2011] Sven F Crone, Michele Hibon, and Konstantinos Nikolopoulos. 2011\.  Advances in forecasting with neural networks? Empirical evidence from the NN3 competition on time series prediction.  _International Journal of forecasting_ 27, 3 (2011), 635–660.
  * Elsayed et al. [2021] Shereen Elsayed, Daniela Thyssens, Ahmed Rashed, Hadi Samer Jomaa, and Lars Schmidt-Thieme. 2021\.  Do we really need deep learning models for time series forecasting?  _arXiv preprint arXiv:2101.02118_ (2021).
  * Friedman [2001] Jerome H Friedman. 2001\.  Greedy function approximation: a gradient boosting machine.  _Annals of statistics_ (2001), 1189–1232.
  * Ganaie et al. [2022] Mudasir A Ganaie, Minghui Hu, AK Malik, M Tanveer, and PN Suganthan. 2022.  Ensemble deep learning: A review.  _Engineering Applications of Artificial Intelligence_ 115 (2022), 105151\.
  * Godahewa et al. [2021] Rakshitha Godahewa, Christoph Bergmeir, Geoffrey I Webb, Rob J Hyndman, and Pablo Montero-Manso. 2021.  Monash time series forecasting archive.  _arXiv preprint arXiv:2105.06643_ (2021).
  * Google [2017] Google. 2017.  Web Traffic Time Series Forecasting | Kaggle.  https://www.kaggle.com/competitions/web-traffic-time-series-forecasting.
  * Grill et al. [2020] Jean-Bastien Grill, Florian Strub, Florent Altché, Corentin Tallec, Pierre Richemond, Elena Buchatskaya, Carl Doersch, Bernardo Avila Pires, Zhaohan Guo, Mohammad Gheshlaghi Azar, et al. 2020\.  Bootstrap your own latent-a new approach to self-supervised learning.  _Advances in neural information processing systems_ 33 (2020), 21271–21284.
  * Hahn et al. [2023a] Yannik Hahn, Tristan Langer, Richard Meyes, and Tobias Meisen. 2023a.  Time Series Dataset Survey for Forecasting with Deep Learning.  _Forecasting_ 5, 1 (2023), 315–335.  https://doi.org/10.3390/forecast5010017
  * Hahn et al. [2023b] Yannik Hahn, Tristan Langer, Richard Meyes, and Tobias Meisen. 2023b.  Time Series Dataset Survey for Forecasting with Deep Learning.  _Forecasting_ 5, 1 (2023), 315–335.
  * Harries et al. [1999] Michael Harries, New South Wales, et al. 1999\.  Splice-2 comparative evaluation: Electricity pricing.  _UNSW_ (1999).
  * Hearst et al. [1998] Marti A. Hearst, Susan T Dumais, Edgar Osuna, John Platt, and Bernhard Scholkopf. 1998.  Support vector machines.  _IEEE Intelligent Systems and their applications_ 13, 4 (1998), 18–28.
  * Hyndman [2020] Rob J Hyndman. 2020\.  A brief history of forecasting competitions.  _International Journal of Forecasting_ 36, 1 (2020), 7–14.
  * Izmailov et al. [2018] Pavel Izmailov, Dmitrii Podoprikhin, Timur Garipov, Dmitry Vetrov, and Andrew Gordon Wilson. 2018\.  Averaging weights leads to wider optima and better generalization.  _arXiv preprint arXiv:1803.05407_ (2018).
  * Januschowski et al. [2020] Tim Januschowski, Jan Gasthaus, Yuyang Wang, David Salinas, Valentin Flunkert, Michael Bohlke-Schneider, and Laurent Callot. 2020\.  Criteria for classifying forecasting methods.  _International Journal of Forecasting_ 36, 1 (2020), 167–177.
  * Kingma and Ba [2014] Diederik P Kingma and Jimmy Ba. 2014.  Adam: A method for stochastic optimization.  _arXiv preprint arXiv:1412.6980_ (2014).
  * Lara-Benítez et al. [2021] Pedro Lara-Benítez, Manuel Carranza-García, and José C Riquelme. 2021.  An experimental review on deep learning architectures for time series forecasting.  _International Journal of Neural Systems_ 31, 03 (2021), 2130001\.
  * Li et al. [2017] Yaguang Li, Rose Yu, Cyrus Shahabi, and Yan Liu. 2017\.  Diffusion convolutional recurrent neural network: Data-driven traffic forecasting.  _arXiv preprint arXiv:1707.01926_ (2017).
  * Lim et al. [2021] Bryan Lim, Sercan Ö Arık, Nicolas Loeff, and Tomas Pfister. 2021.  Temporal fusion transformers for interpretable multi-horizon time series forecasting.  _International Journal of Forecasting_ 37, 4 (2021), 1748–1764.
  * Liu et al. [2023] Xu Liu, Yuxuan Liang, Chao Huang, Hengchang Hu, Yushi Cao, Bryan Hooi, and Roger Zimmermann. 2023.  Do We Really Need Graph Neural Networks for Traffic Forecasting?  _arXiv preprint arXiv:2301.12603_ (2023).
  * Makridakis et al. [2020] Spyros Makridakis, Evangelos Spiliotis, and Vassilios Assimakopoulos. 2020.  The M4 Competition: 100,000 time series and 61 forecasting methods.  _International Journal of Forecasting_ 36, 1 (2020), 54–74.
  * Makridakis et al. [2022] Spyros Makridakis, Evangelos Spiliotis, and Vassilios Assimakopoulos. 2022.  M5 accuracy competition: Results, findings, and conclusions.  _International Journal of Forecasting_ 38, 4 (2022), 1346–1364.
  * Oreshkin et al. [2019] Boris N Oreshkin, Dmitri Carpov, Nicolas Chapados, and Yoshua Bengio. 2019.  N-BEATS: Neural basis expansion analysis for interpretable time series forecasting.  _arXiv preprint arXiv:1905.10437_ (2019).
  * PEMS [[n. d.]] PEMS. [n. d.].  Caltrans PeMS.  https://pems.dot.ca.gov/.
  * Picard [2021] David Picard. 2021\.  Torch. manual_seed (3407) is all you need: On the influence of random seeds in deep learning architectures for computer vision.  _arXiv preprint arXiv:2109.08203_ (2021).
  * Salinas et al. [2020] David Salinas, Valentin Flunkert, Jan Gasthaus, and Tim Januschowski. 2020.  DeepAR: Probabilistic forecasting with autoregressive recurrent networks.  _International Journal of Forecasting_ 36, 3 (2020), 1181–1191.
  * Smyl [2020] Slawek Smyl. 2020\.  A hybrid method of exponential smoothing and recurrent neural networks for time series forecasting.  _International Journal of Forecasting_ 36, 1 (2020), 75–85.
  * Tian and Chan [2021] Chenyu Tian and Wai Kin Chan. 2021.  Spatial-temporal attention wavenet: A deep learning framework for traffic prediction considering spatial-temporal dependencies.  _IET Intelligent Transport Systems_ 15, 4 (2021), 549–561.
  * UCI [[n. d.]] UCI. [n. d.].  UCI Machine Learning Repository: ElectricityLoadDiagrams20112014 Data Set.  https://archive.ics.uci.edu/ml/datasets/ElectricityLoadDiagrams20112014.
  * UNC [2020] UNC. 2020.  M5 Forecasting - Accuracy | Kaggle.  https://www.kaggle.com/competitions/m5-forecasting-accuracy/discussion/164599.
  * Wu et al. [2020] Zonghan Wu, Shirui Pan, Guodong Long, Jing Jiang, Xiaojun Chang, and Chengqi Zhang. 2020\.  Connecting the dots: Multivariate time series forecasting with graph neural networks. In _Proceedings of the 26th ACM SIGKDD international conference on knowledge discovery & data mining_. 753–763.
  * Wu et al. [2019] Zonghan Wu, Shirui Pan, Guodong Long, Jing Jiang, and Chengqi Zhang. 2019.  Graph wavenet for deep spatial-temporal graph modeling.  _arXiv preprint arXiv:1906.00121_ (2019).
  * Yu et al. [2016] Hsiang-Fu Yu, Nikhil Rao, and Inderjit S Dhillon. 2016.  Temporal regularized matrix factorization for high-dimensional time series prediction.  _Advances in neural information processing systems_ 29 (2016).
  * Zeng et al. [2022] Ailing Zeng, Muxi Chen, Lei Zhang, and Qiang Xu. 2022\.  Are transformers effective for time series forecasting?  _arXiv preprint arXiv:2205.13504_ (2022).
  * Zheng et al. [2020] Chuanpan Zheng, Xiaoliang Fan, Cheng Wang, and Jianzhong Qi. 2020\.  Gman: A graph multi-attention network for traffic prediction. In _Proceedings of the AAAI conference on artificial intelligence_ , Vol. 34. 1234–1241.

##  Appendix A Additional Experiments

###  A.1 Impact of Lookback Window Length

In this section an additional study is conducted on the effect of the lookback
window length $l$ on the final model performance. In particular, MTGNN and
DCRNN are designed to work with large chunks of data comprising of many
series, i.e. they operate in a local multivariate setting ($n=1,d=325$), which
imposes a limit on either the batch size or lookback window. For the sake of
comparison, we also train TFT, DeepAR, N-BEATS and N-HITS on PEMS-BAY with the
same lookback window as MTGNN and DCRNN as in the results presented in Table
4. The ultimate goal of time series forecasting task should be predicting the
future values based on all the available history, regardless of the context
length. To this end, in addition to the lookback window length of 12 we sweep
it across 24, 48 and 288, then apply the same procedure we defined in section
3.5 for each model. The results for the best configuration is provided in
Table. 6.

Table 6: Summary of result for varying lookback window lengths. The top
configuration is selected based on the mean average error (MAE) as the
criterion.

| Context 24 | Context 48 | Context 288
---|---|---|---
Model | MAE$\downarrow$ | RMSE$\downarrow$ | SMAPE$\downarrow$ | MAE$\downarrow$ | RMSE$\downarrow$ | SMAPE $\downarrow$ | MAE$\downarrow$ | RMSE$\downarrow$ | SMAPE$\downarrow$
TFT | 1.566 | 3.68 | 3.165 | 1.583 | 3.693 | 3.202 | 1.64 | 3.795 | 3.355
N-BEATS | 1.937 | 4.632 | 3.912 | 1.909 | 4.571 | 3.848 | 1.74 | 4.145 | 3.545
N-HITS* | - | - | - | - | - | - | 1.824 | 3.994 | 3.721
MTGNN** | 1.534 | 3.566 | 3.074 | 1.549 | 3.617 | 3.113 | - | - | -
DCRNN** | 1.528 | 3.564 | 3.063 | 1.512 | 3.564 | 3.028 | - | - | -
DeepAR | 1.681 | 3.789 | 3.391 | 1.666 | 3.773 | 3.363 | 1.895 | 4.186 | 3.869
Notes: *Due similarity of N-BEATS and N-HITS architectures and their
performances on context length of 12 we skipped contexts 24 and 48.
**For DCRNN and MTGNN experiments for context length of 288 were not performed due the memory requirements to run these |

It is worth mentioning that the longer the lookback window generally results
in longer model training. This means that with the same computational budget
$B$, the HPO algorithm might not find an optimal configuration for longer
lookback windows hence there is an inherent trade-off. The results presented
in 6 are the results for the best configuration found by HPO. These are biased
towards lower errors due to the stochastic nature in training these models. A
stability test on the best configurations is conducted, summarized in Table. 7
using a sample $>$ 100 models. Additionally, the best configuration for the
XGBoost model which is manually crafted by a domain expert Kaggle Grandmaster,
which uses a much larger lookback window comparatively to other models.

Table 7: Error mean and variance for the best configurations found by HPO using a sample $>$ 100 models. Model | Lookback window length | MAE $\downarrow$ | RMSE $\downarrow$ | SMAPE $\downarrow$
---|---|---|---|---
TFT | 48 | 1.627 $\pm$ 0.0453 | 3.838 $\pm$ 0.1272 | 3.303 $\pm$ 0.1044
N-BEATS | 288 | 1.739 $\pm$ 0.0018 | 4.139 $\pm$ 0.0097 | 3.543 $\pm$ 0.0037
N-HITS | 288 | 1.723 $\pm$ 0.0016 | 4.109 $\pm$ 0.0079 | 3.511 $\pm$ 0.0036
MTGNN | 24 | 1.544 $\pm$ 0.0058 | 3.584 $\pm$ 0.0147 | 3.097 $\pm$ 0.0118
DCRNN | 48 |  1.53 $\pm$ 0.0045 |  3.562 $\pm$ 0.015 |  3.063 $\pm$ 0.0093
DeepAR | 48 | 1.685 $\pm$ 0.0186 | 3.825 $\pm$ 0.0419 | 3.398 $\pm$ 0.0369
XGBoost | 4043 | 1.577 $\pm$ 0.0012 | 3.685 $\pm$ 0.0043 | 3.216 $\pm$ 0.0026

###  A.2 Effect of Exponential Moving Average and Curriculum Learning

In our sweeps we consider EMA and CL as training techniques that can be turned
on or off when searching for the best model configuration. Even though TPE
algorithm holds an internal state containing information about how beneficial
particular hyperparameter are, it is hard to interpret this state in a
meaningful way. Instead we summarize which of the top configurations selected
this technique for the final model in Table. 8. In general EMA was found to be
quite useful for improving downstream performance.

Table 8: EMA usage amongst top configurations.

Model | Electricity | PEMS-BAY | PEMS-BAY 24 | PEMS-BAY 48 | PEMS-BAY 288 | M5 | Wiki Traffic
---|---|---|---|---|---|---|---
TFT | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓
N-BEATS | ✓ | ✓ | ✓ | ✗ | ✓ | ✗ | ✗
N-HITS | ✓ | ✓ | - | - | ✓ | ✓ | ✗
DeepAR | ✓ | ✓ | ✓ | ✓ | ✓ | ✓ | ✓
DCRNN | ✓ | ✓ | ✓ | ✓ | - | - | -
MTGNN | ✓ | ✗ | ✗ | ✗ | - | - | -
Notes: *Due similarity of N-BEATS and N-HITS architectures and their
performances on context length of 12 we skipped contexts 24 and 48.
**For DCRNN and MTGNN experiments for context length of 288 were not performed
due the memory requirements to run these

For the case of curriculum learning we considered this method for MTGNN which
was used explicitly by the authors in their proposed approach. As a result,
sweeps were conducted on this hyperparameter for the MTGNN model. From the
results summarized in 9, we found that on both electricity and PEMS-BAY HPO
algorithm tends to choose configurations with lower curriculum learning update
steps, which reflects that there is a tendency to for the HPO procedure to
turn off this feature.

Table 9: Best configurations found by HPO algorithm on MTGNN. Baseline is the mean error value for the best configuration without CL, reported in the main section of the paper.  | Electricity | PEMS-BAY
---|---|---
Configuration | MAE$\downarrow$ | RMSE$\downarrow$ | SMAPE$\downarrow$ | MAE$\downarrow$ | RMSE$\downarrow$ | SMAPE$\downarrow$
MTGNN w/o CL | 41.57 | 324.43 | 9.33 | 1.55 | 3.59 | 3.10
MTGNN + CL | 41.85 | 329.34 | 9.38 | 1.55 | 3.57 | 3.10

###  A.3 Performance analysis

In Table 10 the average time required to train each of the top performing
models is provided. The time is calculated with respect to the average number
of epochs required to converge to the best checkpoint. We observed that TFT
consistently requires less time to train than XGBoost, while achieving
comparable errors.

Table 10: Average time required to train top performing configurations. | Electricity | PEMS-BAY
---|---|---
Model | convergence epoch | training time[s] | convergence epoch | training time[s]
TFT | 2.98 | 700.37 | 9.43 | 1912.27
NBEATS | 15.78 | 570.73 | 19.27 | 1673.75
NHITS | 22.85 | 824.58 | 24.96 | 2227.26
MTGNN | 18.15 | 2386.161276 | 66.77 | 9841.90
DCRNN | -∗ | -∗ | 33.71 | 11078.38
XGBoost | - | 1989.19 | - | 2188.38
DeepAR | 1.07 | 1448.98 | 8.84 | 1303.71
Notes: *DCRNN requires static graph which is available only for PEMS-BAY. |

##  Appendix B Framework Flowchart

The high sensitivity of models to hyperparameter changes and the high
variability of accuracy require us to be rigorous in our experiments. To
ensure a fair comparison, certain parts of the experimental flow must be fixed
(train/val/test splits, metrics, hyperparameter optimization algorithm).
Others, however, must be sampled by the hyperparameter optimization algorithm
using the hyperparameter space $\mathbb{S}$: model architecture parameters,
training loop, pre- and post-processing algorithms.

Figure 2: A flowchart of the framework, where orange shows components that can
be changed by a hyperparameter optimizatoin algorithm.

##  Appendix C Datasets

Below we provide a more elaborate description for each dataset and how they
are processed/standardized for the experiments.

  1. 1.

Electricity

Electricity load diagrams[34] (not to be confused with Electricity[16]) is a
commonly [24, 31, 38] used dataset. It consists of 369 distinct univariate
time series without additional covariates. For our experiments the data is
processed as in [24] and [31], taking 500k randomly selected samples between
2014-01-01 to 2014-09-01 as a training set, from which the last week are used
as validation split. The test set spans one week immediately after validation
set, from 2014-09-01 00:00:00 to 2014-09-07 23:00:00. We compute a z-score on
the training set (excluding validation subsplit) on each series separately and
apply to each series separately.

  2. 2.

PEMS-BAY

The Caltrans Performance Measurement System (PeMS)[29] is a traffic
measurement system spanning freeway system across all major metropolitan areas
of the State of California. The system consists of more than 15000 sensors. In
the literature one can find models evaluated on different subsets of this
data. The most commonly used is PEMS-BAY as described [23].

For our experiments which require graph definition, we use the one precomputed
by DCRNN authors, to closely match other papers [40, 33, 4, 36]. Each data
point is a 5 minute aggregate of statistics of traffic collected by a single
sensor. The objective is to predict the next hour (12 time steps) of traffic
using previous hour as a prior. The dataset consists of the same 325 selected
time series, where the training split ranges from 2017-01-01 to 2017-05-08,
the validation split from 2017-05-08 to 2017-05-25 17:50:00, and the test
split from 2017-05-25 17:50:00 to 2017-07-01. The first 12 steps of the test
split are used as an encoder input to a model, thus labels range from
2017-05-25 18:50:00 to 2017-07-01.

  3. 3.

M5

Kaggle M5[27] competition’s objective was to accurately predict 30.5k values
of unit sales across different Walmart shops for the next 28 days. The dataset
used for the competition has a hierarchical structure, where each level
consists of entries aggregated by specific criterion, for example using same
location, product or product category, etc. In total there were 42,840 series
across 12 levels of aggregation. Evaluation was performed on all series,
including aggregates, with the RMSE metric. In our experiments exactly the
same dataset as in the M5 competition. We follow the standard train/valid/test
splits: data from 2011-01-29 to 2016-04-24 for training set, from 2016-04-25
to 2016-05-22 for validation set and from 2016-05-23 to 2016-06-19 for test
set. During our experiments an encoder length of 28 was used resulting in in
example total length of 56. Due to the nature of our problem formulation,
where we want to use one model to predict all the time series, we use bottom-
up approach, trying to predict only series at level 12(lowest level of
aggregation) and then aggregate predictions to produce higher levels without
any correction.

  4. 4.

Wiki traffic

Kaggle Web Traffic Time Series Forecasting [12] competition’s objective was to
predict next two months of daily ’hits’ (visits) for 145k pages collected
across 10 Wikimedia sites for different access and agent types. Authors choose
to use SMAPE as a target metric to evaluate models’ performance. The dataset
consists of time series with different levels of popularity, which introduces
an added complexity as models have to be robust to scale change, but even more
significant is the fact that data for this challenge was collected as it came
in and not all the pages from the original training set existed by the end of
the competition, those pages were excluded from test evaluation. Additionally,
original data contained 2 days gap between the validation and test set, due to
the final train and validation sets being published 2 days before test data
started being collected. That’s why we used modified train/validation/test
ranges, to account for the gap: for train set we used data points dated
2015-07-01 to 2017-07-08, for validation - 2017-05-06 to 2017-09-10 and for
test - 2017-07-09 to 2017-11-13. Which, after the exclusion of 64 encoder
inputs produces test set from 2017-09-11 to 2017-11-13 instead of from
2017-09-13 to 2017-11-13, as in the original competition data. We use longer
prediction horizon of 64 instead of 62 to make dataset continuous, but we
assign zero weight to the first 2 values in the test set to produce metrics
comparable to those used during the competition. We remove zero values from
the beginning of each series. Series that contain less than 64 observations in
the training set after removing the zero prefix, are excluded from the dataset
entirely.

##  Appendix D Tuner

In this section we provide more details on the specific ranges considered for
each hyperparamter (i.e. specifiying $\mathcal{A}$ and the tuner setup when
conducting the experiments & evaluation.

###  D.1 Model Hyper-parameters

A summary of the hyperparameters used for various models is provded in Table
11, where paranthesis "()" denote ranges, comma seperated values are explicit
options, and brackets enclose model specific configurations.

Table 11: Model specific hyper parameters and their ranges considered in this
work.

Model | Parameter name | Parameter range
---|---|---
TFT | #head | 1, 2, 4
| hidden size | 96, 128, 192, 256
| dropout | (0, 0.5)
N-BEATS | stack 0 type | trend, generic
| stack 1 type | seasonality, generic
| stack 0 # blocks | 2, 4, 8
| stack 1 # blocks | 2, 4, 8
| stack 0 theta dim | 2, 4, 8, 16
| stack 1 theta dim | 0, 2, 4, 8, 16
| stack 0 share weights | True, False
| stack 1 share weights | True, False
| stack 0 hidden size | 256, 512, 1024, 2048, 4096*
| stack 1 hidden size | 256, 512, 1024, 2048, 4096*
N-HITS | mlp layers | 2,3,4
| hidden size | 256,512,1024,2048
| activation | ReLU, Softplus, Tanh, SELU, LeakyReLU, PReLU, Sigmoid
| pooling mode | MaxPool1d, AvgPool1d
| #blocks | [1,1,1],[1,1,2],[1,2,1],[1,2,2],[2,1,1],[2,1,2],[2,2,1],[2,2,2]
| pool kernel size | [6,3,1],[6,2,1],[4,2,1],[3,3,1],[2,2,1]
| freq downsample | [6,3,1],[6,2,1],[4,2,1],[3,3,1],[2,2,1]
DeepAR | # layers | 2,3,4,5
| hidden size | 64,128, 256, 512, 1024, 2048**
| dropout | (0, 0.6)
MTGNN | gcn depth | 2,3,4
| dropout | (0,0.5)
| subgraph size | 10,15,20
| node dim | 32,40,64,128
| conv channels | 16,32,64
| residual channels | 16,32,64
| skip channels | 32,64
| end channels | 32,64
| #layers | 2,3,4
| propalpha | (0.01, 0.2)
| tanalpha | (2.0, 4.0)
| in dim | 16, 24, 32, 64
| embedding | True, False
DCRNN | max diffusion step | 1, 2
| #RNN layers | 2, 3
| RNN units | 32, 64, 128
| activation | tanh, ReLU