2022

[1]Rakshitha Godahewa

1]Department of Data Science and AI, Monash University, Melbourne, Victoria,
Australia

# SETAR-Tree: A Novel and Accurate Tree Algorithm for Global Time Series
Forecasting

<EMAIL_ADDRESS>Geoffrey I. Webb<EMAIL_ADDRESS>Daniel
Schmidt<EMAIL_ADDRESS>Christoph Bergmeir
<EMAIL_ADDRESS>[

###### Abstract

Threshold Autoregressive (TAR) models have been widely used by statisticians
for non-linear time series forecasting during the past few decades, due to
their simplicity and mathematical properties. On the other hand, in the
forecasting community, general-purpose tree-based regression algorithms
(forests, gradient-boosting) have become popular recently due to their ease of
use and accuracy. In this paper, we explore the close connections between TAR
models and regression trees. These enable us to use the rich methodology from
the literature on TAR models to define a hierarchical TAR model as a
regression tree that trains globally across series, which we call SETAR-Tree.
In contrast to the general-purpose tree-based models that do not primarily
focus on forecasting, and calculate averages at the leaf nodes, we introduce a
new forecasting-specific tree algorithm that trains global Pooled Regression
(PR) models in the leaves allowing the models to learn cross-series
information and also uses some time-series-specific splitting and stopping
procedures. The depth of the tree is controlled by conducting a statistical
linearity test commonly employed in TAR models, as well as measuring the error
reduction percentage at each node split. Thus, the proposed tree model
requires minimal external hyperparameter tuning and provides competitive
results under its default configuration. We also use this tree algorithm to
develop a forest where the forecasts provided by a collection of diverse
SETAR-Trees are combined during the forecasting process. In our evaluation on
eight publicly available datasets, the proposed tree and forest models are
able to achieve significantly higher accuracy than a set of state-of-the-art
tree-based algorithms and forecasting benchmarks across four evaluation
metrics.

###### keywords:

Threshold Autoregressive Models, Global Time Series Forecasting, Tree Models,
Forest Models

##  1 Introduction

Global time series forecasting models (GFM, Januschowski et al., 2020) have
shown their capability of providing accurate forecasts compared to traditional
univariate forecasting models such as Exponential Smoothing (ETS, Hyndman et
al., 2008) and Auto-Regressive Integrated Moving Average models (ARIMA, Box et
al., 2015), in particular with obtaining superior results in the recently held
M4 (Makridakis et al., 2018) and M5 (Makridakis et al., 2022) forecasting
competitions. In contrast to traditional univariate forecasting models that
build isolated models to forecast each series, GFMs are trained across
collections of time series allowing the models to learn cross-series
information, which also allows them to generalise better and control model
complexity on a global level (Montero-Manso and Hyndman, 2021).

While a lot of research in this space focuses on neural networks and deep-
learning models (e.g., Salinas et al., 2020; Hewamalage et al., 2021b;
Oreshkin et al., 2020), tree-based GFM algorithms have risen greatly in
popularity after dominating the winning methods of the M5 forecasting
competition (Makridakis et al., 2022). In particular, these approaches use
LightGBM (Ke et al., 2017) which is a highly efficient gradient boosted tree-
based algorithm, to train global models across time series. Also other tree-
based algorithms that can be trained as GFMs such as eXtreme Gradient Boosting
(XGBoost, Chen and Guestrin, 2016) and CatBoost (Prokhorenkova et al., 2018)
have been used in the forecasting space. The advantages of these methods
(especially compared with deep-learning methods) are their speed of training,
their ease to include external regressors (numerical and categorical), and
their relative ease and robustness of hyperparameter tuning. However, these
algorithms are general-purpose algorithms, and none of them is specifically
designed for forecasting. When performing regression, they use the average of
the training outputs at a leaf node as the prediction for an instance that
reaches that node during forecasting. Furthermore, while potentially easier to
handle than in other algorithms, their accuracy is still subject to tuning of
a number of hyperparameters such as the maximum tree depth, maximum number of
leaf nodes, minimum number of instances in a leaf node and learning rate.
These findings motivate us to develop a global tree-based algorithm
specifically for forecasting, achieving higher forecasting accuracy with
minimal hyperparameter tuning.

Piecewise linear models such as Threshold Autoregressive (TAR) models (Tong,
1978, 1993; Terasvirta, 1994) have greatly contributed to the econometric and
forecasting fields during the past few decades. These models partition the
instance space of a particular single time series to be forecast into several
subspaces where each subspace is modeled by a separate Autoregressive (AR)
model. There are many variations of TAR models such as Self Exciting Threshold
Autoregressive (SETAR, Tong, 1993), Smooth Transition Autoregressive (STAR,
Terasvirta, 1994), Exponential STAR (ESTAR), Logistic STAR (LSTAR) and Neuro-
Coefficient STAR (Medeiros and Veiga, 2000) models.

Similar to the works of Aznarte and Benítez (2010), Aznarte et al. (2010) and
Aznarte et al. (2007) that establish analogies between certain types of TAR
models and Fuzzy Rule Based Systems, we exploit the analogy of TAR models and
decision trees that fit a linear model per leaf, for global time series
forecasting, to introduce a new tree-based algorithm that uses the underlying
concept of SETAR models in defining the splits. Thus, we name this model as
SETAR-Tree. In particular, our method internally finds the optimal past lag
and the threshold that should be used to split the instances of a particular
tree node into child nodes. The instances whose corresponding values of the
optimal lag are greater than or equal to the selected optimal threshold and
less than the selected optimal threshold are grouped separately whereas each
instance group is considered as a child node. In contrast to the traditional
general purpose tree-based regression algorithms that compute the average of
the training outputs at a leaf node as the final prediction, our proposed
method is a forecasting-specific tree which builds a Pooled Regression model
(PR, Gelman and Hill, 2006; Montero-Manso and Hyndman, 2021) which is a global
linear AR model, at each leaf node allowing the model to learn cross-series
information. As our model is a forecasting-specific tree with global
regression models in the leaf nodes, it can be essentially identified as a
global (cross-series) hierarchical SETAR model. Furthermore, the SETAR-Tree
uses some time-series-specific stopping procedures. Exploiting the analogy of
TAR models and regression trees, our proposed SETAR-Tree internally tunes the
maximum tree depth using a statistical linearity test well established in the
literature for TAR models (Terasvirta, 1994), as well as measuring the error
reduction percentage at each node split. Thus, our proposed model requires
minimal external hyperparameter tuning whereas it provides competitive results
even with its default hyperparameters. The proposed method is also applicable
to time series forecasting problems that use external covariates during
modelling.

Proposing a novel tree algorithm for forecasting naturally lends itself to
extend the procedure to a forest, which is another main contribution of this
study. In line with the traditional forest-based algorithms such as Random
Forest (RF, Breiman, 2001), our forest model uses a collection of diverse
SETAR-Trees that are trained using a set of randomly selected time series. The
forecasts provided by all trees are averaged to obtain the final forecasts.
The trees are made diverse in terms of the significance of the statistical
linearity test and error reduction percentage that are used to split each
node. Thus, the forecasting accuracy of the forest model is even less
sensitive to the hyperparameters than the tree.

Our proposed SETAR-Tree and SETAR-Forest algorithms in particular also
outperform a set of state-of-the-art tree-based algorithms and forecasting
benchmarks with statistical significance across eight experimental datasets on
four error metrics. All implementations of this study are publicly available
at: https://github.com/rakshitha123/SETAR_Trees.

The remainder of this paper is organized as follows: Section 2 reviews the
relevant prior work. Section 3 explains the theoretical concepts of SETAR
models, and the proposed tree and forest algorithms in detail. Section 4
explains the experimental framework, including the datasets, error metrics,
benchmarks and statistical testing. An analysis of the results is presented in
Section 5. Section 6 concludes the paper and discusses possible future
research.

##  2 Related work

In the following, we discuss the related prior work in the areas of TAR
models, GFMs and the state-of-the-art tree-based forecasting models.

###  2.1 Threshold autoregressive models

TAR models (Tong, 1978, 1993) are piece-wise linear models that model the
state space of a given prediction problem using multiple AR models. TAR models
are also known as AR regime switching models which use separate (same or
different order) AR functions to model different regimes.

TAR models have variants such as SETAR and STAR (Terasvirta, 1994) which are
different in terms of the transition function and the method of defining
regimes. The simplest version of the TAR model is known as the SETAR model
(Tong, 1993) which defines two or more regimes based on a particular lagged
value of the time series itself. Thus, SETAR models are required to estimate
the optimal lag, and the optimal threshold corresponding with that chosen
optimal lag to define regimes. The most common way of defining these optimal
lag and threshold values is using a grid search. The optimal thresholds can be
also identified by analysing the places of time series showing significant
increasing or decreasing rates (Ghosh et al., 2006) and by using memetic
algorithms (Bergmeir et al., 2012). The transition function used by the SETAR
model is generally a step function. Thus, during prediction, SETAR models use
only one AR model. Versions of the SETAR model exist that use statistical
tests such as a Lagrange Multiplier (LM) test (Breusch and Pagan, 1980) to
determine the linearity captured by the AR models (Terasvirta, 1994). In STAR
models (Terasvirta, 1994), the transition can happen either using a past value
of a series or an external variable. Unlike the SETAR models, STAR models use
a transition function such as exponential or logistic, resulting in ESTAR and
LSTAR models, respectively, and thus, the transition happens smoothly. During
prediction, STAR models use multiple AR models to determine the final output
where the corresponding weight of each contributing AR model is determined
based on the transition function.

Aznarte and Benítez (2010) and Aznarte et al. (2010) establish the equivalence
between TAR models and certain types of fuzzy rule-based systems. Inspired by
this work, we exploit equivalences between TAR models and decision trees, and
extend TAR models to global hierarchical models which train across series.

TAR models have been used to address many real-world forecasting problems such
as stock market forecasting (Narayan, 2006), exchange rate forecasting
(Pippenger and Goering, 1998) and electricity price forecasting (Rambharat et
al., 2005). To the best of our knowledge, none of these TAR models have been
used as global models across time series yet.

###  2.2 Global time series forecasting

Global time series forecasting (Januschowski et al., 2020) is a relatively
recent trend in forecasting that became popular after contributing to the
winning methods of the M4 (Makridakis et al., 2018) and M5 (Makridakis et al.,
2022) forecasting competitions. In particular, many top solutions of the M5
competition used global (tree-based) models. GFMs build a single forecasting
model across many series sharing the same parameters across all series. Thus,
the model has the capability to explore cross-series information with a fewer
amount of parameters compared to traditional univariate forecasting models.

Many works use Recurrent Neural Networks (RNN, Hewamalage et al., 2021b) as
the baseline model when developing global forecasting frameworks. The winning
method of the M4 competition (Smyl, 2020) uses an ensembling approach,
ensemble of specialists (Smyl, 2020), that combines the predictions of a set
of globally trained RNNs. Smyl and Kuber (2016) use globally trained RNNs for
forecasting. One of the state-of-the-art probabilistic forecasting algorithms,
DeepAR (Salinas et al., 2020), uses a global AR RNN. Godahewa et al. (2022)
propose a weekly forecasting framework that uses global RNNs.

In particular, the winning method of the M5 forecasting competition uses a
globally trained LightGBM model (Ke et al., 2017). GFMs have also contributed
to the winning methods of many other recently held Kaggle competitions (Bojer
and Meldgaard, 2020). Nowadays, GFMs are used with many real-world
applications such as energy optimisation (Godahewa et al., 2022), sales demand
forecasting (Bandara et al., 2019) and emergency medical services demand
forecasting (Bandara et al., 2020).

Godahewa et al. (2021) and Bandara et al. (2020) show that a GFM should be
trained with a set of related time series. Those authors further claim that if
a collection of time series has heterogeneous series, then the series should
be first grouped based on their similarity and then, separate GFMs should be
trained per each series group. In line with this concept, our proposed SETAR-
Tree model trains multiple GFMs, one per leaf node, because the instances
belonging to a particular leaf node are similar in terms of their values
corresponding with the optimal lag identified using the SETAR methodology. As
the baseline GFM, we use a linear AR model trained across a set of time series
that is known as a PR model (Gelman and Hill, 2006).

###  2.3 Tree-based forecasting models

There are many popular tree-based algorithms used for forecasting, including
regression trees (Loh, 2011), gradient boosted trees (Friedman, 2001), XGBoost
(Chen and Guestrin, 2016), CatBoost (Prokhorenkova et al., 2018), and LightGBM
(Ke et al., 2017).

As forecasting is usually a regression task with numerical values, regression
trees (Loh, 2011) are used for forecasting, as opposed to decision trees for
classification. Spiliotis (2022) provides an overview of how decision trees
can be used for time series forecasting. In particular, the author highlights
that in the forecasting domain, a regression tree provides a piece-wise
approximation to a single continuous function in standard regression.
Regression trees use every possible binary split on every input attribute to
determine the node splits. Regression trees have been successfully applied to
address real-world time series forecasting applications. Torgo and Oliveira
(2014) use bagged regression trees to handle the diverse dynamic regimes and
non-stationarities observed in real-world time series. Cerqueira et al. (2019)
use a regression tree as a base model of a pool of forecasting experts where
the forecasts of experts are dynamically combined using arbitrating to obtain
the final forecasts. Gradient boosted trees (Friedman, 2001) use a collection
of regression trees as a sequential ensemble model. The first regression tree
is trained on the actual observations of the training instances where the
remaining trees are trained to predict the residuals of the previous tree in
the sequence. XGBoost (Chen and Guestrin, 2016) is an accurate, efficient and
scalable version of gradient boosted trees. To define the best node splits,
XGBoost uses similarity scores in a way that the gain is maximised.
Furthermore, to reduce overfitting, XGBoost uses regularisation parameters
when calculating similarity scores and tree pruning. The second winning method
of the M4 forecasting competition, Feature-based Forecast Model Averaging
(FFORMA, Montero-Manso et al., 2020) uses XGBoost as a meta-learner to
determine the weights that should be used to combine the predictions of a set
of base forecasting models. CatBoost (Prokhorenkova et al., 2018) considers
the order of data points during modelling and thus, it is more suitable to
address time series forecasting problems. CatBoost more effectively deals with
categorical variables compared to other tree-based algorithms.

As stated earlier, tree-based forecasting models recently obtained popularity
in the forecasting field after contributing to the winning approach and many
other top submissions in the M5 forecasting competition (Makridakis et al.,
2022). In particular, many top submissions of the competition incorporate a
highly efficient gradient boosted tree-based model, LightGBM (Ke et al.,
2017). The LightGBM trees use the leaf-wise tree growth approach instead of
the level-wise tree growth used in traditional tree-based algorithms. It
identifies the best node splits by using Gradient-based One-Side Sampling
(GOSS). LightGBM can handle extremely large datasets as well and nowadays, it
is used with many real-world forecasting applications such as cryptocurrency
price trend forecasting (Sun et al., 2020), sales demand forecasting (Weng et
al., 2020) and wind power forecasting (Ju et al., 2019). Januschowski et al.
(2021) highlight that tree-based algorithms are blackbox learners which are
highly efficient and more robust than the existing deep learning techniques.
In particular, the LightGBM algorithm offers a large number of loss functions
and missing value handling methods where the users can choose the
corresponding options based on their application.

Even though these tree algorithms are efficient, they are not tailored to the
forecasting problem. For example, they consider the average of the training
outputs in a particular leaf node as the forecast of the test instances that
reach that node during forecasting. In contrast, our proposed SETAR-Tree
algorithm trains a global linear AR model in each leaf node, allowing each
leaf node to learn the cross-series information.

In the literature, there already exist tree models which train linear models
in leaf nodes, known as linear regression model trees (Quinlan, 1992). They
analyse every possible split to determine the node splits where the split
which maximises the difference between the standard deviation of the actual
target values of the parent node and the mean of the standard deviations of
the actual target values of child nodes, is considered as the optimal split.
The prediction for a given test instance is obtained by its corresponding leaf
node where the prediction is then adjusted based on the nodes from the root
node to its corresponding leaf node. Cubist (Kuhn and Johnson, 2013) is such a
linear model tree that performs rule-based modelling. It trains linear models
in all tree nodes where the predictions are obtained from the leaf nodes and
are smoothed afterwards. Lefakis et al. (2019) propose a piecewise linear
regression tree which trains regularised linear models in leaf nodes. The
nodes are split in a way that the least square errors of linear predictors
after splitting are minimised. Loh (2002) proposes a regression tree which
uses the significance of the chi-square test to identify the variables that
should be used during node splitting. Dutang and Guibert (2022) propose a
splitting procedure for regression trees where the optimal splits are chosen
based on the explicit likelihood. Zeileis et al. (2008) propose a regression
model which uses recursive partitioning where the nodes are only further split
if there exists an overall parameter instability. Hothorn et al. (2006)
propose CTree, a regression model that uses recursive partitioning with
conditional inference procedures. There are hardly any works in the literature
that propose tree-based models particularly for forecasting that we are aware
of, and the few works available address local forecasting, not global
forecasting models. In particular, the STR-Tree proposed by da Rosa et al.
(2008), though proposed for regression, is deemed by its authors to be
trivially extendable to forecasting, and Epprecht and Veiga (2012) propose a
tree method to predict stock market returns where in the leaf nodes AR models
with extra inputs are trained. These works use the concept of STAR models
during node splitting, which is similar to our approach but limits them
effectively to smaller sample sizes, as the smooth transitions make their
model training computationally considerably more costly. Thus, these methods
are inherently more suitable for local forecasting and using them in current
typical global forecasting scenarios brings scalability issues.

Our proposed SETAR-Tree algorithm is in particular tailored to the time series
forecasting task by using the concept of SETAR models in defining splits where
the optimal lag and the threshold that should be used for node splitting is
determined using a grid search approach, without considering every possible
split as used by the above explained linear regression model trees. Our grid
search approach selects the lag and the threshold that is corresponding with
the split which minimises the Sum of Squared Errors (SSE) at the child nodes.
Furthermore, to determine the validity of a split, we use a statistical
linearity test and the error reduction percentage that can be gained from the
split, where a node is only further split if making that split is worth
enough. In leaf nodes, we train global AR models allowing the models to learn
the cross-series information and they do not use any regularisation during
model fitting. Also, the predictions for the test instances are obtained only
from their corresponding leaf nodes where no further adjustments are applied
to the predictions as in linear regression model trees (for details, see
Section 3.3).

A collection of parallel executed tree models is known as a forest. The RF
(Breiman, 2001) algorithm is the most common type of forest algorithms which
averages the predictions provided by a set of diverse regression trees to
obtain the final predictions. The trees can be diverse in terms of the data
and features used during training. This ensembling technique is also known as
bagging in the machine learning literature and it is expected to improve the
accuracy of the final predictions compared to the predictions of the
individual trees by properly addressing data, model and parameter
uncertainties (Petropoulos et al., 2018; Bergmeir et al., 2016). Also,
Coulombe (2021) argues that the RF does an optimal pruning as the splits that
overfit are averaged out through bagging. In the context of linear model
trees, the Macroeconomic RF (MRF) proposed by Coulombe (2020) consists of a
collection of linear model trees that use the concept of TAR models during
node splitting. However, MRF is not designed for global time series
forecasting. Unlike our SETAR-Tree, the individual trees in MRF use
regularisation during parameter estimation, and use moving average factors
calculated using the lagged values of time series as regressors during
training. The trees also do not use any stopping criterion based on
statistical tests or error reduction percentages. In particular, the
individual trees in MRF do not use early stopping where they are allowed to
grow into higher depths ensuring the trees are sufficiently diversified.
Furthermore, Athey et al. (2019) propose generalized RFs, a method based on
RFs that can be used for non-parametric statistical estimation. Consequently,
we propose a forest model that uses a collection of SETAR-Trees where the
individual trees are diverse in terms of the significance of the linearity
test and the error threshold used to make the node splits (for details, see
Sections 3.3.2, 3.3.3 and 3.4).

##  3 Methodology

This section first gives a brief overview of the theory of SETAR models. Then,
we explain the terminology of GFMs. Later, the proposed SETAR-Tree and forest
algorithms are explained in detail.

###  3.1 Self exciting threshold autoregressive models

For our study, we consider a 2-regime SETAR model (Tong, 1993) which is
defined in Equation 1. Here, $y_{t}$ is the value of the series at time $t$,
$n$ is the number of past time series lags, $T$ is the threshold value, and
$\beta^{i}$ and $\epsilon^{i}$ are respectively the AR parameters and error
terms of the $i^{th}$ regime.

$y_{t}=\begin{cases}\sum_{l=1}^{n}\beta^{(1)}_{l}{y_{t-l}}+\epsilon^{(1)}_{t},&\text{if
}y_{t-l}<T\\\ \\\
\sum_{l=1}^{n}\beta^{(2)}_{l}{y_{t-l}}+\epsilon^{(2)}_{t},&\text{otherwise}\end{cases}$
(1)

As shown in Equation 1, SETAR models define regimes based on a particular
lagged value of the time series itself. One regime is defined for the training
instances where the corresponding lagged value is less than the threshold
value $T$, and the other regime is defined for the remaining instances where
the corresponding lagged value is greater than or equal to $T$. Thus, SETAR
models have as hyperparameters the optimal lag, $l$, and the optimal threshold
value, $T$ corresponding with $l$ to define the regimes.

###  3.2 Terminology of global forecasting models

For our study, we consider univariate time series forecasting which predicts
the future values of a particular time series using its own past values and
possibly, some external covariates. We use the usual out-of-sample evaluation
in forecasting that reserves a block of data from the very end of the time
series for evaluation (Bergmeir and Benítez, 2012). For that, the time series
are first split into training and test parts where the training parts are used
for model training and the test parts are the actual values corresponding with
the expected forecast horizon.

Let $TS_{1}$ and $TS_{2}$ be two time series. To train purely autoregressive
GFMs, as shown in Figure 1, an embedded matrix is constructed using $TS_{1}$
and $TS_{2}$ where each row contains a window of a time series which is
composed of a set of lagged values and the corresponding next value. The rows
of the embedded matrix that are corresponding with the training parts of the
series are used to train GFMs. The forecasts corresponding with the expected
forecast horizon are obtained for each series using the trained GFMs. During
the evaluation, the forecasts are compared with their corresponding actual
values that are in the test set and the forecasting errors are measured using
the error metrics explained in Section 4.2.

Figure 1: Visualisation of an embedded matrix created using two time series:
$TS_{1}$ and $TS_{2}$. Each row contains a window of a time series which is
composed with a set of lagged values ($L1$ to $L5$) and the corresponding next
value/true output ($y$). The windows corresponding with the training and test
sets are respectively highlighted in yellow and blue.

###  3.3 SETAR-Tree model

Figure 2 shows an overview of an example SETAR-Tree model. Level 0 refers to
the root node of the tree. The root node is an embedded matrix constructed
using the training parts of $TS_{1}$ and $TS_{2}$ as explained in Section 3.2.
Each row of the embedded matrix is used as a training instance. We use the
concept of SETAR models to split the training instances into child nodes. As
we only consider 2-regime SETAR models, a given tree node is split into two
child nodes. The optimal lag, $l$ and the optimal threshold, $T$ that are used
to split each node are determined by using a grid search approach as explained
in Section 3.3.1. The instances where the corresponding value of the Lag $l$
is less than $T$, and greater than or equal to $T$, are separately grouped
into child nodes. However, a node is only split if making that split is worth
enough. In particular, a node is only split if there exists remaining non-
linearity in the training instances of that node (Section 3.3.2) and/or it can
gain a certain training error reduction by splitting (Section 3.3.3). Thus,
our model uses the leaf-wise tree growth approach where the branches of the
tree may have different lengths. For an example, at Level 1, only one node is
further split into child nodes. The other node is not split and it has been
treated as a leaf node because further splitting that node does not fulfill
the specified criteria.

Figure 2: Overview of an example SETAR-Tree constructed using $TS_{1}$ and
$TS_{2}$. The matrices containing the past time series lags and their
corresponding true outputs are considered as nodes. Each row of a matrix is
considered as a training instance. The training instances of each node are
optimally split into child nodes based on the concept of SETAR models. Here
$l$ and $T$ mentioned with the split nodes refer to the optimal lag and the
optimal threshold that are used to split at a particular node. A node is only
split if there exists remaining non-linearity in the training instances of
that node and/or it can gain a certain training error reduction by splitting.
This example tree model has three levels. This model uses a leaf-wise tree
growth approach and thus, the branches of the tree may have different lengths

The two methods we use to determine whether a node should be split or not, a
statistical linearity test and the error reduction gained by node splitting,
are explained in detail in Sections 3.3.2 and 3.3.3. They control the growth
of the tree and thus, the model can itself determine the maximum tree depth
that it should allow. In particular, these two methods act as the stopping
criteria of the tree growth where constructing the tree is completed when none
of the nodes should be further split. The leaf nodes can be there at different
tree levels as shown in Figure 2 where one leaf node is in Level 1 and two
leaf nodes are in Level 2. In each leaf node, a global PR model is trained
using the instances corresponding with that node.

During testing, the model identifies the leaf node corresponding with a given
test instance by following the same optimal lags and thresholds it previously
used during node splitting starting from the root node. The prediction of the
test instance is then obtained using the trained PR model corresponding with
its leaf node.

The SETAR-Tree requires minimal hyperparameter tuning. A few parameters are
used for the statistical linearity test and for measuring error reduction in
node spitting. However, the model provides competitive results with their
default parameters that we provide and thus, setting externally tuned
hyperparameters is not mandatory for our model.

####  3.3.1 Finding optimal lags and thresholds

To find the lag and threshold that will be used to split a given node, we
utilise a standard grid-search approach. Given a node, we consider each lag
and a set (grid) of thresholds with which to split the node into two children.
For each lag and threshold pair in the grid, separate linear models are fitted
to the two subsets of data formed by partitioning the chosen lag at the
specified threshold. The lag and the threshold that are results in the split
with the minimum SSE at the child nodes is selected as the “optimal” lag and
threshold for splitting.

To increase the efficiency of the grid-search, we exploit the properties of
the least-squares estimates during fitting. Recall that the leaf models are
linear AR models of the form

$y=X\beta+\varepsilon$

where $X$ is the lag-matrix, $y$ are the outputs, $\beta$ the AR coefficients
and $\varepsilon$ a vector of errors. The least-squares estimates of $\beta$
are well known to be

$\hat{\beta}=(X^{T}X)^{-1}X^{T}y$ (2)

where T and -1 respectively refer to the transpose and inverse of the matrix.
The SSE of the least-squares fit is given by

${\rm SSE}(\hat{\beta})=y^{T}y-\hat{\beta}^{T}X^{T}X\hat{\beta}.$ (3)

An important property of the least-squares estimates are that they depend on
$X$ and $y$ only through their inner products. We can exploit this to fit the
models for all potential thresholds in $O(n\log n)$ time. Let $X_{l}$ be the
lag variable we will be splitting on, and let $s=(s_{1},\ldots,s_{q})$ be the
$q$ candidate thresholds for splitting. Define the sets

$S_{i}=\left\\{j:x_{j,l}\leq s_{i}\right\\},(i=1,\ldots,q)\;{\rm
and}\;S_{q+1}=\left\\{j:x_{j,l}>s_{q}\right\\}.$

These sets indicate which rows fall into each of the $q+1$ different
partitions of the data defined by the set of thresholds. Let $\bar{x}_{i}$
denote the $i$-th row of $X$; we can then define the quantities

$B_{(k)}=\sum_{i=1}^{k}\sum_{j\in
S_{i}}\bar{x}_{j}\bar{x}_{j}^{T},\;\;c_{(k)}=\sum_{i=1}^{k}\sum_{j\in
S_{i}}\bar{x}_{j}y_{j}\;\;{\rm and}\;\;d_{(k)}=\sum_{i=1}^{k}\sum_{j\in
S_{i}}y_{j}^{2}.$ (4)

These are the relevant inner products for the “left” child node model formed
by splitting the data at each of the thresholds $s_{k}$. From Equation 2, if
we choose to split at threshold $s_{k}$, the least squares estimates for the
left and right child node models can be written as

$\hat{\beta}{(L_{k})}=B_{(k)}^{-1}c_{(k)}\;\;{\rm
and}\;\;\hat{\beta}{(R_{k})}=\left(X^{T}X-B_{(k)}\right)^{-1}\left(X^{T}y-c_{(k)}\right)$
(5)

respectively. Using Equation 3, the SSE in the left and right child node
models can be then written as

${\rm
SSE}{(L_{k})}=d_{(k)}-\left(\hat{\beta}{(L_{k})}\right)^{T}B_{(k)}\left(\hat{\beta}{(L_{k})}\right)$
(6)

and

${\rm
SSE}{(R_{k})}=\left(y^{T}y-d_{(k)}\right)-\left(\hat{\beta}{(L_{k})}\right)^{T}\left(X^{T}X-B_{(k)}\right)\left(\hat{\beta}{(L_{k})}\right)$
(7)

respectively. When considering a lag $x_{l}$ to split on, we can first sort
the rows of $X$ by $x_{l}$. In practice, these sorted values can then be used
to find a set of suitable thresholds; typically we choose fifteen equispaced
quantiles of the vector $x_{l}$. We can then use Equation 4 to compute the
relevant inner products for the first threshold, i.e., $B_{(1)}$, $c_{(1)}$
and $d_{(1)}$. We can then iteratively try each potential threshold
$s_{1},\ldots,s_{q}$ in turn using Equations 5, 6 and 7, exploiting the fact
that we can update the inner products using

$B_{(k+1)}=B_{(k)}+\sum_{j\in S_{k}}\bar{x}_{j}\bar{x}_{j}^{T},$

as we move from threshold $s_{k}$ to $s_{k+1}$, with similar expressions for
$c_{(k)}$ and $d_{(k)}$ to efficiently compute the running totals of these
quantities.

####  3.3.2 Statistical linearity test

The usage of linearity tests in determining the amount of regimes and the
goodness of TAR models is well-established in the econometrics and statistics
literature (Terasvirta, 1994). However, to the best of our knowledge, the
linearity tests have not been used to determine the goodness of a SETAR model
trained across series, in particular to determine the depth of a hierarchical
SETAR model as we propose.

We use the general linear F-test (Box, 1953) to determine whether there exists
a remaining non-linearity of a set of training instances in a particular tree
node. The node is only further split if there exists a significant remaining
non-linearity.

The null hypothesis of the linear F-test, $H_{0}$ states that there exists no
significant remaining non-linearity in the training instances and thus, the
parent node is enough to model the instances. The alternative hypothesis,
$H_{1}$ states that there exists a significant remaining non-linearity in the
training instances and therefore, a node should be further split into child
nodes.

In particular, the linear F-test calculates a statistic, namely the
F-statistic, which compares a complex model and a reduced model, and
determines whether the reduced model should be rejected in favour of the
complex model. In our node splitting scenario, the reduced and complex models
respectively refer to the global PR models that are trained on parent and
child nodes corresponding with a particular split. The suitability of the
split is determined based on the value of the F-statistic calculated using the
SSE corresponding with PR models built on parent and child nodes.

The SSE of a parent node and the corresponding child nodes are respectively
defined in Equations 8 and 9. Here, $N$ is the number of training instances,
$y_{k}$ are the observed values, $yp_{k}$ are the fitted values of the parent
node model and $yc_{k}$ are the fitted values of the child node models.

$SSE(P)=\sum_{k=1}^{N}({y_{k}-yp_{k}})^{2}$ (8)

$SSE(C)=\sum_{k=1}^{N}({y_{k}-yc_{k}})^{2}$ (9)

The F-statistic corresponding with a node split is then calculated using
Equation 10 where $df_{P}$ and $df_{C}$ are the degrees of freedom of parent
and child node models, respectively.

$F^{*}=\left(\frac{SSE(P)-SSE(C)}{df_{P}-df_{C}}\right){\left(\frac{SSE(C)}{df_{C}}\right)}^{-1}$
(10)

For simplicity, we consider the same number of past lag values when training
the PR models in all nodes. Hence, $df_{P}$ and $df_{C}$ are defined as in
Equations 11 and 12 where $L$ is the number of past lags used for modelling.

$df_{P}=N-L-1$ (11)

$df_{C}=N-2L-2$ (12)

Substituting $df_{P}$ and $df_{C}$ reduces Equation 10 as follows.

$F^{*}=\left(\frac{SSE(P)-SSE(C)}{L+1}\right){\left(\frac{SSE(C)}{N-2L-2}\right)}^{-1}$
(13)

The $p$-value of the F-test is determined by comparing the value of $F^{*}$ to
the $F$ distribution. If the $p$-value is less than the considered
significance level, $\alpha$, then $H_{1}$ is accepted and the node is further
split into child nodes. Otherwise, the node is kept as a leaf node of the
tree.

As a form of multiple testing correction (Noble, 2009), we gradually decrease
the value of $\alpha$ when the tree grows level by level as defined in
Equation 14, where $\alpha_{d}$ and $\alpha_{d+1}$ are the significance levels
corresponding with the tree levels $d$ and $d+1$, respectively.

$\alpha_{d+1}=\frac{\alpha_{d}}{\mathit{significance\\_divider}}$ (14)

This reduces unnecessary tree growth and acts as a way of preventing model
overfitting. The values of $\alpha_{0}$ and $\mathit{significance\\_divider}$
are user defined, with default values of 0.05 and 2, repectively. As our trees
perform binary splits, using a $\mathit{significance\\_divider}$ of 2
corresponds to a Bonferroni correction, where the family-wise error is kept
constant across the levels of the tree. Larger values lead accordingly to more
conservative significance levels.

####  3.3.3 Error reduction in node splitting

We also consider the error reduction percentage to determine whether a node
should be further split or not. The child nodes models together use more
parameters during model fitting and thus, they always give lower training
error compared to the parent node model. Here, the purpose is to check whether
the error reduction percentage that can be gained by splitting a node is
considerably high or not. A node is only split if the error reduction
percentage is greater than or equal to a particular error threshold ($e_{t}$)
and otherwise, the node is considered as a leaf node of the tree. The value of
$e_{t}$ is user defined and based on our preliminary experiments, it is set to
3%, by default.

Algorithm 1 shows the training and testing phases of the SETAR-Tree model. The
algorithm takes the following parameters as inputs.

Algorithm 1 SETAR-Tree Algorithm

1:procedure setar_tree($training\\_set$, $forecast\\_horizon$, $lag$,
$stopping\\_criteria$, $max\\_depth=1000$ $alpha0=0.05$,
$significance\\_divider=2$, $error\\_threshold=0.03$)

2: $alpha$ $\leftarrow$ $alpha0$

3: $tree$ $\leftarrow$ $list()$

4: $opt\\_lags$ $\leftarrow$ $list()$

5: $opt\\_thresholds$ $\leftarrow$ $list()$

6: $input\\_matrix$ $\leftarrow$ $create\\_input\\_matrix(training\\_set,lag)$

7: $node\\_data$ $\leftarrow$ $input\\_matrix$

8: for $depth$ in $max\\_depth$ do

9: $level\\_significant\\_node\\_count$ $\leftarrow$ 0

10: for $node$ in $node\\_data$ do

11: $opt\\_params$ $\leftarrow$ $get\\_opt\\_params(node)$

12: $opt\\_lag$ $\leftarrow$ $opt\\_params[1]$

13: $opt\\_threshold$ $\leftarrow$ $opt\\_params[2]$

14: $child\\_nodes$ $\leftarrow$
$split\\_node(node,opt\\_lag,opt\\_threshold)$

15: if $stopping\\_criteria$ $=$ $"linearity\\_test"$ then

16: $good\\_split$ $\leftarrow$ $check\\_linearity(node,child\\_nodes,alpha)$

17: else if $stopping\\_criteria$ $=$ $"error\\_reduction"$ then

18: $good\\_split$ $\leftarrow$
$check\\_error\\_reduction(node,child\\_nodes,error\\_threshold)$

19: else

20: $good\\_split$ $\leftarrow$ $check\\_linearity(node,child\\_nodes,alpha)$
and

21:$check\\_error\\_reduction(node,child\\_nodes,error\\_threshold)$

22: end if

23: if $good\\_split$ then

24: $tree[depth]$ $\leftarrow$ $add(tree[depth],child\\_nodes)$

25: $opt\\_lags[depth]$ $\leftarrow$ $add(opt\\_lags[depth],opt\\_lag)$

26: $opt\\_thresholds[depth]$ $\leftarrow$
$add(opt\\_thresholds[depth],opt\\_threshold)$

27: $level\\_significant\\_node\\_count$ $\leftarrow$
$level\\_significant\\_node\\_count+1$

28: end if

29: end for

30: if $level\\_significant\\_node\\_count>0$ then

31: $node\\_data$ $\leftarrow$ $tree[depth]$

32: $alpha$ $\leftarrow$ $alpha/significance\\_divider$

33: else

34: $break$

35: end if

36: end for

37: $leaf\\_nodes$ $\leftarrow$ $get\\_leaf\\_nodes(tree)$

38: $leaf\\_models$ $\leftarrow$ $list()$

39: for $leaf$ in $leaf\\_nodes$ do

40: $model$ $\leftarrow$ $train\\_pr\\_model(leaf)$

41: $leaf\\_models$ $\leftarrow$ $add(leaf\\_models,model)$

42: end for

43: $predictions$ $\leftarrow$ $matrix()$

44: $test\\_set$ $\leftarrow$ $create\\_test\\_set(training\\_set,lag)$

45: for $f$ in $forecast\\_horizon$ do

46: for $s$ in $nrow(test\\_set)$ do

47: $leaf\\_node\\_index$ $\leftarrow$
$find\\_leaf\\_node\\_index(test\\_set[s,],tree,opt\\_lags,opt\\_thresholds)$

48: $predictions[s,f]$ $\leftarrow$
$predict(leaf\\_models[leaf\\_node\\_index],test\\_set[s,])$

49: end for

50: $test\\_set$ $\leftarrow$
$update\\_test\\_set(test\\_set,predictions[,f])$

51: end forreturn $predictions$

52:end procedure

  * •

training_set: contains the training parts of each series of the dataset

  * •

forecast_horizon: the length of the prediction period

  * •

lag: the number of past lags used to forecast the next series value

  * •

stopping_criteria: the criteria that are used to determine whether a node
should be split

  * •

max_depth: maximum possible tree depth (optional)

  * •

alpha0: the initial significance level to be used by the F-test (optional)

  * •

significance_divider: the significance level used by the F-test is
sequentially divided at each tree level by this value (optional)

  * •

error_threshold: minimum error reduction that should be gained to split a node
(optional)

The algorithm first creates an embedded matrix using the training parts of
each series which is considered as the root node of the tree. The function
create_input_matrix creates this embedded matrix taking training_set and lag
as inputs (line 6). It then finds the optimal lag and the threshold that
should be used to split the node into child nodes using the function
get_opt_params (line 11). This function internally conducts a grid search over
the lags of the node and selects the lag and the threshold that are
corresponding with the split which provides the minimum SSE at the child
nodes. Then, the function split_node (line 14) splits the parent node into
child nodes using the optimal lag (opt_lag) and the threshold (opt_threshold)
returned by get_opt_params. Based on the stopping_criteria, the split is
further assessed to determine whether the child nodes should be added to the
tree or not. The functions check_linearity (line 16) and check_error_reduction
(line 18) respectively check whether there exists remaining non-linearity in
the training instances of the parent node and whether the training error
reduction gain by node splitting is greater than the error_threshold. If the
split is worth to make, then the child nodes are added to the tree (line 24).
The tree building procedure is stopped if none of the nodes at the current
tree level is split (line 34). If at least one node at the current tree level
is split, then the same node splitting procedure is conducted for the next
tree level. At each tree level, the significance level used by the F-test,
alpha, is decreased by dividing it by the significance_divider (line 32).
After completing the tree construction process, our method trains a global PR
model at each leaf node using the function train_pr_model (line 40). During
the testing phase, it first creates a set of test instances, test_set, using
the function create_test_set which returns the last set of lags in each
training series (line 44). The forecasts for the required forecast horizon are
then obtained iteratively. In each iteration, the function
find_leaf_node_index (line 47) identifies the corresponding leaf node for each
test instance by following the same optimal lags (opt_lags) and thresholds
(opt_thresholds) it previously used during node splitting. The predictions for
the test instances are then obtained using the corresponding trained models at
their leaf nodes (line 48). In each iteration, the test_set is also updated
using the function update_test_set, by appending the forecasts of the previous
iteration (line 50).

####  3.3.4 Training with covariates

The proposed SETAR-Tree model can also be trained with external numerical and
categorical covariates in addition to the past lags of time series. The
numerical covariates are treated in the same way as the time series lagged
values. The categorical covariates are converted into numerical format
beforehand by applying one-hot encoding which then allows to treat them in the
same way as numerical attributes. When determining the optimal attribute and
the threshold to split a node, the past lags of series, numerical covariates
and categorical covariates are all considered together.

###  3.4 SETAR-Forest model

A forest is an ensemble model that aggregates the predictions provided by a
set of diverse tree models to obtain a prediction, often through bootstrap
aggregation (“bagging”) (Breiman, 2001).

We introduce a new forest algorithm, SETAR-Forest, which consists of a
collection of diverse SETAR-Trees. The trees are made diverse by varying the
initial significance level ($\alpha_{0}$), the significance divider used to
calculate the sequence significance levels and the error reduction percentage
threshold ($e_{t}$) used during node splitting. In a SETAR-Tree, these three
parameters are set by default based on our preliminary experiments and in
contrast to that, in the forest algorithm, each tree selects the values of
these parameters randomly. For our experiments, we consider ten SETAR-Trees
for a SETAR-Forest. An individual SETAR-Tree does not use bagging whereas the
SETAR-Forest does. A set of randomly chosen 80% of instances and all
attributes of them are used to train each SETAR-Tree in a SETAR-Forest. Thus,
the values of bagging frequency, bagging fraction and feature fraction of our
SETAR-Forest are 10, 0.8, and 1, respectively, by default. However, our
implementation takes these values as external parameters and thus, users can
modify them based on the application.

During testing, the forecasts of the test instances are obtained by averaging
the corresponding forecasts provided by all SETAR-Trees in the forest.

##  4 Experimental framework

In this section, we discuss the experimental datasets, error metrics, and
benchmarks used in our experiments.

###  4.1 Datasets

Table 1: Datasets information

Dataset Name | No. of | Forecast | Frequency | Minimum | Maximum
---|---|---|---|---|---
| Time Series | Horizon |  | Length | Length
Rossmann | 1115 | 48 | Daily | 894 | 894
Wikipedia Web Traffic | 1000 | 59 | Daily | 744 | 744
Favorita | 1000 | 16 | Daily | 1668 | 1668
M5 | 1490 | 28 | Daily | 1941 | 1941
Tourism Monthly | 366 | 24 | Monthly | 67 | 309
Tourism Quarterly | 427 | 8 | Quarterly | 22 | 122
Chaotic Logistic | 100 | 8 | - | 592 | 592
Mackey-Glass | 100 | 8 | - | 592 | 592

We use eight experimental datasets111The experimental datasets are available
at https://github.com/rakshitha123/SETAR_Trees/tree/master/datasets., six
publicly available datasets and two simulated datasets, to evaluate the
performance of our proposed SETAR-Tree and forest algorithms. The datasets are
briefly explained in the following and Table 1 provides their summary
statistics.

  * •

Rossmann Sales Dataset (Kaggle, 2015): The dataset from the Rossmann Sales
forecasting competition that shows the daily sales of a set of Rossmann
stores. The missing observations of this dataset are replaced by carrying
forward the corresponding last observations (LOCF method). Several covariates
are also available with this dataset such as the daily number of customers to
each store, status of the stores (open or closed), promotion details and
state/school holidays information.

  * •

Kaggle Wikipedia Web Traffic Dataset (Google, 2017): The first 1000 time
series from the Kaggle Wikipedia Web Traffic forecasting competition that
shows the number of daily hits for a given set of Wikipedia web pages.
Following the procedure suggested by Hewamalage et al. (2021b), the missing
observations of this dataset are replaced by zeros. Covariates are also
available with this dataset such as day of the week. The dataset was extracted
from Godahewa et al. (2021).

  * •

Favorita Dataset (Kaggle, 2018): The first 1000 time series from the
Corporación Favorita Grocery Sales forecasting competition that shows daily
unit sales of a set of items sold at different Favorita stores. The missing
observations of this dataset are replaced by zeros. Covariates are also
available with this dataset such as day of the week.

  * •

M5 Dataset (Makridakis et al., 2022): A subset of time series from the M5
forecasting competition that shows daily unit sales of a set of items sold at
different Walmart stores. To make the amount of time series comparable with
the other datasets, we select the items that belong to the $HOBBIES\\_2$
department which is one of the ten departments used in the competition.

  * •

Tourism Monthly Dataset (Athanasopoulos et al., 2011): Monthly dataset from
the Tourism forecasting competition. The dataset was extracted from Godahewa
et al. (2021).

  * •

Tourism Quarterly Dataset (Athanasopoulos et al., 2011): Quarterly dataset
from the Tourism forecasting competition. The dataset was extracted from
Godahewa et al. (2021).

  * •

Chaotic Logistic Dataset: A simulated dataset used in Hewamalage et al.
(2021a) constructed by using the Chaotic Logistic Map (May, 1976) Data
Generation Process (DGP)

  * •

Mackey-Glass Dataset: A simulated dataset used in Hewamalage et al. (2021a)
constructed by using the Mackey-Glass Equation (Mackey and Glass, 1977) DGP.

For our study, we have mostly used datasets where globally trained tree models
are known to provide good forecasts. Note that four datasets: Rossmann, Kaggle
Web Traffic, Favorita and M5 are from Kaggle competitions where tree-based
models have won the first or second positions (Bojer and Meldgaard, 2020). The
two simulated datasets also contain non-linear time series where tree-based
models have provided accurate forecasts compared to PR and neural network
models based on the results published by Hewamalage et al. (2021a). To add
some diversity into the pool of datasets, we consider the Tourism Monthly and
Quarterly datasets from the Tourism forecasting competition where the
traditional univariate forecasting models have provided better forecasts.
There are other popular benchmarking datasets in the global forecasting field,
most notably the M3 (Makridakis and Hibon, 2000) and M4 (Makridakis et al.,
2018) datasets. However, these datasets are very heterogeneous as they contain
time series from many different use cases and require special preprocessing
and/or ensembling techniques such as clustering and ensemble of specialists
(Godahewa et al., 2021; Bandara et al., 2020; Smyl, 2020) to make global
models competitive. As we deem such additional techniques to be out of the
scope of our current paper, we do not consider these datasets for this study.

###  4.2 Error metrics

We measure the performance of our models using the modified version of
symmetric Mean Absolute Percentage Error (msMAPE, Suilin, 2017) and Mean
Absolute Scaled Error (MASE, Hyndman and Koehler, 2006) which are commonly
used error metrics in the time series forecasting field. Equations 15 and 16
respectively define the msMAPE and MASE error metrics. Here, $F_{k}$ are the
forecasts, $Y_{k}$ are the actual values for the required forecast horizon,
$M$ is the number of instances in the training set, $N$ is the number of data
points in the test set and $S$ is the length of the seasonal cycle. In
Equation 15, the $\epsilon$ is set to its proposed default value of 0.1.

$msMAPE=\frac{100\%}{N}\sum_{k=1}^{N}\frac{\lvert
F_{k}-Y_{k}\rvert}{(max(\lvert Y_{k}\rvert+\lvert
F_{k}\rvert+\epsilon,0.5+\epsilon))/2}$ (15)

$MASE=\frac{\sum_{k=1}^{N}\lvert
F_{k}-Y_{k}\rvert}{\frac{N}{M-S}\sum_{k=S+1}^{M}\lvert Y_{k}-Y_{k-S}\rvert}$
(16)

For the two simulated datasets, the length of the seasonal cycle is considered
as one. Since all these error measures are defined for each time series, we
calculate the mean and median values of them across a dataset to measure the
model performance. Therefore, four error metrics are used to evaluate each
model: mean msMAPE, median msMAPE, mean MASE and median MASE.

###  4.3 Benchmarks and variants

We use seven global models: a PR model (Gelman and Hill, 2006), Feed-Forward
Neural Network (FFNN, Goodfellow et al., 2016) and five tree-based models: a
regression tree (Loh, 2011), Cubist (Kuhn and Johnson, 2013), LightGBM (Ke et
al., 2017), XGBoost (Chen and Guestrin, 2016), and CatBoost (Prokhorenkova et
al., 2018), as the main benchmarks of our study. These tree-based models are
popular well-performing models in the forecasting domain. Regression trees
have provided successful insights and foundations for machine learning
approaches in the forecasting domain (Spiliotis, 2022). Cubist is a well-known
regression model tree that performs rule-based modelling (Quinlan, 1992). Many
top solutions of the M5 forecasting competition (Makridakis et al., 2022)
including the winning method use globally trained LightGBM models.
Furthermore, the winning methods of many recently held Kaggle forecasting
competitions use tree-based models (Bojer and Meldgaard, 2020). The FFORMA
method (Montero-Manso et al., 2020) which placed second in the M4 forecasting
competition (Makridakis et al., 2018) uses XGBoost as a meta-learner. CatBoost
considers the order of data points during modelling which makes it more
suitable for time series forecasting (Prokhorenkova et al., 2018).

The R packages `glmnet` (Friedman et al., 2010), `nnet` (Venables and Ripley,
2002), `rpart` (Therneau and Atkinson, 2019), `Cubist` (Kuhn and Quinlan,
2022), `lightgbm` (Ke et al., 2020), `xgboost` (Chen et al., 2020) and
`catboost` (Prokhorenkova et al., 2018) are respectively used to implement the
PR, FFNN, regression tree, Cubist, LightGBM, XGBoost and CatBoost models. The
hyperparameters of FFNN, LightGBM and XGBoost are tuned using a grid search
approach. For the FFNN, the number of nodes in the hidden layers is varied
from 10 to 60 in steps of 10, and the learning rate decay is varied from 0.01
to 0.1 in steps of 0.01. For LightGBM, the minimum number of instances in a
leaf node is varied from 50 to 200 in steps of 50, and the learning rate is
varied from 0.01 to 0.1 in steps of 0.01. For XGBoost, the maximum tree depth
is varied from 3 to 10 in steps of 1, and $\eta$ is varied from 0.1 to 0.5 in
steps of 0.1. The regression tree, Cubist and CatBoost models are used with
their default parameters that are assumed to provide good results. The PR
model does not require any hyperparameter tuning.

We also consider two traditional univariate forecasting models: ETS (Hyndman
et al., 2008) and ARIMA (Box et al., 2015) that are commonly used in the
forecasting space, as the benchmarks of our study. They are implemented using
the R package, `forecast` (Hyndman and Khandakar, 2008) under the default
configurations.

The SETAR (Tong, 1993) and STAR (Terasvirta, 1994) models that are trained per
series are also considered as benchmarks. They are implemented using the R
package, `tsDyn` (Narzo et al., 2022) under the default configurations.

We further consider RF (Breiman, 2001) as a benchmark to evaluate the
performance of our SETAR-Forest algorithm. For that, the implementation of RF
in the R package `lightgbm` (Ke et al., 2020) is used where the
hyperparameters are tuned using a grid search approach. The minimum number of
instances in a leaf node is varied from 50 to 200 in steps of 50, and the
learning rate is varied from 0.01 to 0.1 in steps of 0.01. The bagging
frequency, bagging fraction and feature fraction are respectively set to the
same default values used in our SETAR-Forest algorithm: 10, 0.8 and 1 to
facilitate the comparison between RF and the proposed SETAR-Forest model.

Three variants of the proposed SETAR-Tree model (Section 3.3) are used as
follows.

Tree.Lin.Test Uses the significance of the linearity test as the stopping
criterion of node splitting.

Tree.Error.Red Uses the error reduction percentage that can be gained by node
splitting as the stopping criterion.

Tree.Lin.Test.Error.Red Uses both the significance of the linearity test and
the error reduction percentage gained by node splitting as the stopping
criteria.

Three variants of the proposed SETAR-Forest model (Section 3.4) are used as
follows.

Forest.Significance The trees are randomised in terms of the significance
level and the sequence significance divider considered in the linearity test.

Forest.Error.Red The trees are randomised in terms of the error threshold used
to measure the error reduction percentage of node splitting.

Forest.Significance.Error.Red The trees are randomised in terms of the
significance level and the sequence significance divider considered in the
linearity test as well as the error threshold used to measure the error
reduction percentage of node splitting.

The number of lagged values used in all GFMs including our SETAR-Tree and
forest variants are determined using the heuristics suggested by Hewamalage et
al. (2021b) where the number of lags used with a dataset depends on its
seasonal cycle length or forecast horizon. In particular, for daily and
monthly datasets, we consider the heuristic, $seasonality\times 1.25$ to
determine the number of lags. For daily datasets: Rossmann, Kaggle Web
Traffic, Favorita and M5, this heuristic suggests 8.75 lags, where this value
is rounded (to steps of 5) and 10 lags are finally considered with those
datasets. For the Tourism Monthly dataset, 15 lags are considered as suggested
by the above heuristic. For Chaotic Logistic, Mackey-Glass and Tourism
Quarterly datasets, the seasonality-based heuristic suggests small values, so
that for these datasets we use the horizon-based heuristic, and consider 10
lags which is equal to $horizon\times 1.25$ of those datasets. For the three
datasets with external covariates: Rossmann, Kaggle Web Traffic and Favorita,
the GFMs are separately executed with and without covariates where only the
past lags are used when training without covariates, and the past lags and
covariates are together used when training with covariates.

###  4.4 Statistical testing of the results

The non-parametric Friedman rank-sum test is used to assess the statistical
significance of the results provided by different forecasting models across
time series considering a significance level of $\alpha$ = 0.05 (García et
al., 2010). Based on the corresponding msMAPE errors, the methods are ranked
on every series of the eight primary experimental datasets. In line with
traditional forecasting modelling and evaluation (Koning et al., 2005), we
treat each series as a separate entity for statistical testing, so that the
ranks per series are used as inputs for the statistical tests. This seems
reasonable as the methods are univariate and – once trained – forecast each
series in isolation. In this way, the statistical test has much more
information instead of averaging the results over the datasets before applying
the test. All primary experimental datasets have a comparable amount of series
and thus, none of the datasets will dominate the evaluation and we obtain a
representative result across all primary experimental datasets. The best
method according to the average rank is chosen as the control method. To
further characterise the statistical differences, Hochberg’s post-hoc
procedure is used (García et al., 2010).

##  5 Results

This section discusses the results in terms of main accuracy, statistical
significance and computational performance.

###  5.1 Main accuracy results

This section first presents the results and a comparison of our SETAR-Tree and
SETAR-Forest variants. The best SETAR-Tree and SETAR-Forest variants are then
compared with the benchmarks.

####  5.1.1 Comparison of SETAR-Tree and SETAR-Forest variants

Table 2: Results of SETAR-Tree and SETAR-Forest variants across all
experimental datasets. The best performing variants in each group are
italicized and the overall best performing variants are highlighted in
boldface

| Without Covariates | With Covariates
---|---|---
|  Ross- | Kag- | Favo- | M5 | Tour | Tour | Cha- | Mackey- | Ross- | Kag- | Favo-
|  mann | gle | rita |  | (M) | (Q) | otic | Glass | mann | gle | rita
Mean msMAPE |  |
Tree.Lin.Test | 39.05 | 71.04 | 83.01 | 45.28 | 23.94 | 17.04 | 41.83 | 0.00372 | 11.53 | 66.64 | 94.77
Tree.Error.Red | 54.93 | 44.74 | 85.04 | 53.92 | 22.62 | 19.04 | 49.36 | 0.00661 | 15.20 | 45.01 | 97.01
Tree.Lin.Test.Error.Red | 41.90 | 44.74 | 85.04 | 53.92 | 21.52 | 15.59 | 41.98 | 0.00372 | 12.09 | 44.88 | 97.01
Forest.Significance | 41.65 | 48.13 | 85.06 | 53.91 | 21.17 | 15.59 | 41.30 | 0.00296 | 12.06 | 46.85 | 96.67
Forest.Error.Red | 43.03 | 43.97 | 82.44 | 54.13 | 25.61 | 16.57 | 41.55 | 0.00307 | 12.07 | 47.94 | 95.32
Forest.Significance.Error.Red | 40.73 | 43.80 | 82.36 | 54.13 | 22.16 | 15.97 | 41.14 | 0.00296 | 11.93 | 47.83 | 95.28
Median msMAPE |  |
Tree.Lin.Test | 36.86 | 69.72 | 81.88 | 41.48 | 19.61 | 13.44 | 36.37 | 0.00082 | 10.02 | 63.78 | 90.09
Tree.Error.Red | 42.90 | 41.52 | 87.60 | 49.96 | 17.43 | 13.23 | 41.00 | 0.00118 | 11.64 | 40.95 | 91.33
Tree.Lin.Test.Error.Red | 41.07 | 41.52 | 87.60 | 49.96 | 17.26 | 12.59 | 36.30 | 0.00082 | 10.52 | 40.90 | 91.33
Forest.Significance | 40.72 | 45.70 | 87.66 | 49.96 | 16.96 | 12.62 | 37.41 | 0.00124 | 10.59 | 41.64 | 91.12
Forest.Error.Red | 39.27 | 40.83 | 82.80 | 50.08 | 18.40 | 12.70 | 36.07 | 0.00082 | 10.31 | 43.20 | 90.85
Forest.Significance.Error.Red | 39.11 | 40.81 | 82.80 | 50.08 | 17.71 | 12.57 | 36.38 | 0.00124 | 10.38 | 43.17 | 90.81
Mean MASE |  |
Tree.Lin.Test | 0.546 | 2.605 | 0.835 | 1.330 | 1.726 | 1.698 | 0.682 | 0.00641 | 0.434 | 3.058 | 3.940
Tree.Error.Red | 0.676 | 0.960 | 0.784 | 1.442 | 1.572 | 1.667 | 0.810 | 0.01150 | 0.571 | 1.428 | 1.015
Tree.Lin.Test.Error.Red | 0.522 | 0.960 | 0.784 | 1.442 | 1.575 | 1.559 | 0.687 | 0.00641 | 0.459 | 1.204 | 1.015
Forest.Significance | 0.511 | 0.950 | 0.784 | 1.442 | 1.546 | 1.550 | 0.670 | 0.00512 | 0.453 | 0.946 | 1.003
Forest.Error.Red | 0.506 | 0.895 | 0.757 | 1.445 | 1.627 | 1.589 | 0.676 | 0.00531 | 0.447 | 0.976 | 0.907
Forest.Significance.Error.Red | 0.456 | 0.872 | 0.777 | 1.445 | 1.601 | 1.593 | 0.667 | 0.00512 | 0.446 | 0.966 | 0.906
Median MASE |  |
Tree.Lin.Test | 0.443 | 1.872 | 0.737 | 1.037 | 1.532 | 1.428 | 0.650 | 0.00135 | 0.351 | 1.590 | 0.773
Tree.Error.Red | 0.518 | 0.776 | 0.714 | 1.138 | 1.409 | 1.365 | 0.742 | 0.00185 | 0.416 | 0.787 | 0.839
Tree.Lin.Test.Error.Red | 0.460 | 0.776 | 0.714 | 1.138 | 1.408 | 1.320 | 0.653 | 0.00135 | 0.372 | 0.787 | 0.839
Forest.Significance | 0.443 | 0.820 | 0.713 | 1.138 | 1.401 | 1.291 | 0.657 | 0.00206 | 0.379 | 0.840 | 0.828
Forest.Error.Red | 0.389 | 0.755 | 0.688 | 1.139 | 1.451 | 1.341 | 0.641 | 0.00128 | 0.370 | 0.871 | 0.777
Forest.Significance.Error.Red | 0.383 | 0.755 | 0.688 | 1.139 | 1.416 | 1.323 | 0.646 | 0.00206 | 0.371 | 0.869 | 0.777

Table 2 shows the results of the considered SETAR-Tree and SETAR-Forest
variants across all experimental datasets for mean msMAPE, median msMAPE, mean
MASE and median MASE. The tree variants and forest variants are separately
grouped in Table 2. The results of the best performing variants in each group
are italicized, and the overall best performing variants across the datasets
are highlighted in boldface. The preliminary experiments are also conducted
with other stopping criteria such as using Akaike Information Criterion (AIC),
however, the performance of the SETAR-Tree was not improved with them
confirming the claims of prior similar studies (da Rosa et al., 2008). Thus,
the corresponding results are not included here.

Across the tree model variants, on mean msMAPE, the Tree.Lin.Test model
variant shows the best performance across five primary datasets where the
Tree.Lin.Test.Error.Red model variant shows the best performance across four
primary datasets. However, on all other error metrics, Tree.Lin.Test.Error.Red
shows the overall best performance across the majority of the experimental
datasets from the tree model variants. Compared with the Tree.Lin.Test and
Tree.Lin.Test.Error.Red model variants, the Tree.Error.Red model variant shows
worse performance on all error metrics except across the Kaggle Web Traffic
dataset. This shows that using the significance of the statistical linearity
test individually or together with the error reduction percentage gained by
node splitting are better options rather than using the error reduction
percentage on its own as the stopping criterion of the SETAR-Tree to determine
its maximum depth.

In the SETAR-Forest model variants, as the stopping criteria of each
individual SETAR-Tree, we consider both the significance of the statistical
linearity test and the error reduction percentage gained by node splitting
which is the stopping criterion of our best SETAR-Tree variant,
Tree.Lin.Test.Error.Red. Out of the three SETAR-Forest variants,
Forest.Significance.Error.Red shows an overall better performance compared to
the other two forest variants across all experimental datasets on all error
metrics. This shows that when the individual SETAR-Trees of the forest are
more diversified, the resultant forecasts are more accurate. We see that
overall, the SETAR-Forest variants show the best performance across the
majority of datasets compared to the SETAR-Tree variants on all error metrics.
The forests properly address the data, model, and parameter uncertainties
through bagging compared to the individual trees and that can be considered as
the major reason for this phenomenon.

Our SETAR-Tree and forest model variants also provide interesting results when
trained with external covariates. As discussed in Section 3.3.4, the
categorical covariates of the Rossmann (status of the store, promotions, state
holiday and school holiday flags), Kaggle Web Traffic (day of the week) and
Favorita (day of the week) datasets are converted into the numerical format by
applying a one-hot encoding mechanism before model training. Furthermore, with
the Rossmann dataset, the daily number of customers to each store is
considered as an external numerical covariate. In particular, the overall
result of the Rossmann dataset has considerably improved across our proposed
SETAR-Tree and forest variants after using external covariates on all error
metrics. However, the results of the Kaggle Web Traffic and Favorita datasets
do not show such improvements after using covariates across our proposed model
variants.

####  5.1.2 Comparison of proposed models with benchmarks

Table 3 shows the results of the benchmarks and our best SETAR-Tree and SETAR-
Forest variants, Tree.Lin.Test.Error.Red and Forest.Significance.Error.Red,
across all experimental datasets for mean msMAPE, median msMAPE, mean MASE and
median MASE.

Table 3: Results across all experimental datasets. The best performing models
in each group are italicized and the overall best performing models are
highlighted in boldface

| Without Covariates | With Covariates
---|---|---
|  Ross- | Kag- | Favo- | M5 | Tour | Tour | Cha- | Mackey- | Ross- | Kag- | Favo-
|  mann | gle | rita |  | (M) | (Q) | otic | Glass | mann | gle | rita
Mean msMAPE |  |
ETS | 43.98 | 46.24 | 87.67 | 78.22 | 19.02 | 15.07 | 50.33 | 1.02983 | - | - | -
ARIMA | 45.34 | 47.96 | 87.82 | 77.81 | 19.73 | 16.58 | 48.71 | 11.12100 | - | - | -
SETAR | 62.20 | 46.75 | 94.56 | 58.18 | 31.30 | 36.14 | 52.93 | 0.04079 | - | - | -
STAR | 72.89 | 46.82 | 96.30 | 95.01 | 32.58 | 34.08 | 44.82 | 0.02094 | - | - | -
PR | 64.45 | 111.48 | 85.04 | 53.92 | 21.56 | 17.07 | 52.27 | 0.01949 | 43.02 | 68.78 | 99.22
Cubist | 38.77 | 55.69 | 85.75 | 146.12 | 19.96 | 16.02 | 43.03 | 0.26995 | 13.07 | 55.67 | 85.63
FFNN | 197.35 | 164.74 | 119.40 | 94.97 | 199.47 | 199.77 | 42.78 | 60.53347 | 197.35 | 164.74 | 115.14
Regression Tree | 55.48 | 61.88 | 101.21 | 65.94 | 64.34 | 115.02 | 44.72 | 3.42950 | 46.58 | 61.88 | 101.21
CatBoost | 49.39 | 49.66 | 90.73 | 57.33 | 23.75 | 25.37 | 42.09 | 0.65735 | 39.91 | 47.97 | 90.90
LightGBM | 56.16 | 55.63 | 96.83 | 32.60 | 22.18 | 19.72 | 42.53 | 0.56777 | 42.73 | 59.25 | 98.06
XGBoost | 48.29 | 69.73 | 86.93 | 54.99 | 23.48 | 18.84 | 44.40 | 0.45676 | 48.41 | 65.20 | 89.41
RF | 61.95 | 49.63 | 103.02 | 104.49 | 32.55 | 27.13 | 42.62 | 2.85191 | 46.53 | 49.90 | 101.62
Tree.Lin.Test.Error.Red | 41.90 | 44.74 | 85.04 | 53.92 | 21.52 | 15.59 | 41.98 | 0.00372 | 12.09 | 44.88 | 97.01
Forest.Significance.Error.Red | 40.73 | 43.80 | 82.36 | 54.13 | 22.16 | 15.97 | 41.14 | 0.00296 | 11.93 | 47.83 | 95.28
Median msMAPE |  |
ETS | 44.25 | 41.72 | 86.59 | 72.62 | 17.16 | 12.89 | 38.37 | 0.79343 | - | - | -
ARIMA | 45.89 | 45.03 | 84.84 | 73.08 | 18.00 | 13.13 | 37.26 | 9.78094 | - | - | -
SETAR | 57.45 | 41.55 | 88.75 | 46.39 | 25.27 | 23.36 | 51.64 | 0.03111 | - | - | -
STAR | 61.04 | 41.07 | 89.65 | 76.06 | 23.66 | 17.34 | 39.53 | 0.01078 | - | - | -
PR | 64.88 | 110.64 | 87.60 | 49.96 | 18.60 | 13.42 | 43.81 | 0.01255 | 41.20 | 64.64 | 91.88
Cubist | 42.91 | 53.62 | 87.79 | 148.76 | 17.24 | 13.21 | 39.93 | 0.21272 | 11.18 | 53.62 | 87.69
FFNN | 197.29 | 168.20 | 114.39 | 91.12 | 199.84 | 199.94 | 37.18 | 62.12558 | 197.29 | 168.20 | 114.22
Regression Tree | 52.32 | 53.57 | 92.15 | 62.22 | 57.22 | 120.69 | 37.46 | 2.97342 | 44.85 | 53.57 | 92.15
CatBoost | 48.38 | 45.94 | 89.55 | 47.81 | 20.38 | 19.31 | 36.70 | 0.48910 | 38.74 | 44.06 | 89.77
LightGBM | 57.16 | 51.71 | 89.37 | 26.46 | 19.39 | 14.88 | 38.09 | 0.41013 | 40.92 | 55.78 | 89.23
XGBoost | 46.90 | 67.29 | 83.10 | 47.26 | 21.05 | 16.19 | 37.39 | 0.34459 | 47.43 | 61.00 | 84.82
RF | 60.71 | 45.08 | 95.50 | 102.75 | 27.16 | 20.74 | 37.13 | 2.16739 | 45.11 | 45.38 | 92.89
Tree.Lin.Test.Error.Red | 41.07 | 41.52 | 87.60 | 49.96 | 17.26 | 12.59 | 36.30 | 0.00082 | 10.52 | 40.90 | 91.33
Forest.Significance.Error.Red | 39.11 | 40.81 | 82.80 | 50.08 | 17.71 | 12.57 | 36.38 | 0.00124 | 10.38 | 43.17 | 90.81
Mean MASE |  |
ETS | 0.584 | 2.582 | 0.786 | 1.796 | 1.526 | 1.592 | 0.770 | 1.61780 | - | - | -
ARIMA | 0.596 | 1.407 | 0.782 | 1.807 | 1.589 | 1.782 | 0.743 | 17.59567 | - | - | -
SETAR | 1.260 | 0.897 | 0.848 | 1.387 | 2.225 | 2.591 | 0.851 | 0.06783 | - | - | -
STAR | 1.380 | 0.902 | 0.857 | 1.540 | 2.110 | 1.910 | 0.788 | 0.03498 | - | - | -
PR | 1.333 | 1.483 | 0.784 | 1.442 | 1.681 | 1.654 | 0.802 | 0.03165 | 0.638 | 4.040 | 1.106
Cubist | 0.602 | 1.031 | 0.780 | 4.407 | 1.593 | 1.693 | 0.696 | 0.43350 | 0.483 | 1.036 | 0.779
FFNN | 3.167 | 1.901 | 0.977 | 2.259 | 9.747 | 12.568 | 0.699 | 77.63547 | 3.167 | 1.901 | 0.954
Regression Tree | 0.960 | 3.530 | 1.156 | 1.789 | 8.080 | 82.250 | 0.713 | 5.44867 | 0.667 | 3.530 | 1.156
CatBoost | 0.858 | 1.517 | 0.866 | 1.562 | 1.877 | 2.851 | 0.693 | 1.06160 | 0.449 | 1.584 | 0.860
LightGBM | 1.030 | 1.977 | 0.976 | 1.128 | 1.769 | 2.246 | 0.699 | 0.90972 | 0.539 | 2.660 | 1.060
XGBoost | 0.796 | 2.719 | 0.896 | 1.545 | 1.924 | 2.094 | 0.722 | 0.75285 | 0.710 | 2.753 | 0.850
RF | 1.167 | 1.901 | 1.325 | 2.522 | 2.596 | 3.431 | 0.705 | 4.46711 | 0.663 | 1.913 | 1.175
Tree.Lin.Test.Error.Red | 0.522 | 0.960 | 0.784 | 1.442 | 1.575 | 1.559 | 0.687 | 0.00641 | 0.459 | 1.204 | 1.015
Forest.Significance.Error.Red | 0.456 | 0.872 | 0.777 | 1.445 | 1.601 | 1.593 | 0.667 | 0.00512 | 0.446 | 0.966 | 0.906
Median MASE |  |
ETS | 0.552 | 0.800 | 0.683 | 1.454 | 1.276 | 1.275 | 0.697 | 1.32630 | - | - | -
ARIMA | 0.585 | 0.808 | 0.684 | 1.403 | 1.337 | 1.388 | 0.686 | 16.15742 | - | - | -
SETAR | 1.144 | 0.779 | 0.765 | 1.087 | 1.840 | 2.316 | 0.844 | 0.04834 | - | - | -
STAR | 1.250 | 0.761 | 0.771 | 1.482 | 1.689 | 1.664 | 0.706 | 0.01861 | - | - | -
PR | 1.291 | 1.368 | 0.714 | 1.138 | 1.502 | 1.442 | 0.778 | 0.02055 | 0.480 | 1.531 | 0.886
Cubist | 0.562 | 0.893 | 0.697 | 3.399 | 1.406 | 1.432 | 0.705 | 0.35068 | 0.391 | 0.894 | 0.705
FFNN | 2.981 | 1.705 | 0.808 | 1.848 | 9.105 | 11.614 | 0.688 | 78.37429 | 2.981 | 1.705 | 0.798
Regression Tree | 0.807 | 1.108 | 0.873 | 1.484 | 4.440 | 19.586 | 0.685 | 4.71614 | 0.557 | 1.108 | 0.873
CatBoost | 0.822 | 0.896 | 0.742 | 1.235 | 1.668 | 2.000 | 0.664 | 0.75808 | 0.371 | 0.854 | 0.737
LightGBM | 1.017 | 1.055 | 0.793 | 0.842 | 1.569 | 1.711 | 0.674 | 0.70062 | 0.431 | 1.238 | 0.833
XGBoost | 0.735 | 1.760 | 0.761 | 1.220 | 1.777 | 1.737 | 0.683 | 0.56261 | 0.659 | 1.479 | 0.734
RF | 1.095 | 0.865 | 0.956 | 2.035 | 2.080 | 2.174 | 0.682 | 3.86985 | 0.554 | 0.865 | 0.889
Tree.Lin.Test.Error.Red | 0.460 | 0.776 | 0.714 | 1.138 | 1.408 | 1.320 | 0.653 | 0.00135 | 0.372 | 0.787 | 0.839
Forest.Significance.Error.Red | 0.383 | 0.755 | 0.688 | 1.139 | 1.416 | 1.323 | 0.646 | 0.00206 | 0.371 | 0.869 | 0.777

The models in Table 3 are grouped based on the sub-experiments. The results of
the best performing models in each group are italicized, and the overall best
performing models across the datasets are highlighted in boldface. The first
group contains the traditional univariate forecasting models, ETS, ARIMA,
SETAR and STAR, where these models are only executed with our eight primary
experimental datasets without using external covariates. The second group
contains benchmark GFMs including PR, FFNN, state-of-the-art tree-based models
and RF. The last group contains our best SETAR-Tree and SETAR-Forest variants
that are discussed in Section 5.1.1. The GFMs are executed with and without
covariates for three datasets: Rossmann, Kaggle Web Traffic and Favorita,
separately.

In the first group, ETS shows an overall better performance compared to ARIMA,
SETAR and STAR across all error metrics. ETS shows the best performance across
the Tourism Monthly dataset on all error metrics, Tourism Quarterly dataset on
mean msMAPE and median MASE, and the Favorita dataset on median MASE. However,
GFMs including our proposed models show a better performance than these
traditional univariate forecasting models across all experimental datasets on
all error metrics except for these seven cases.

In the second group, overall, PR and Cubist show a better performance compared
to the other benchmark GFMs across all error metrics. The Cubist and CatBoost
models show a better performance when training with external covariates.

The third group shows the results of our best SETAR-Tree variant,
Tree.Lin.Test.Error.Red, and the best SETAR-Forest variant,
Forest.Significance.Error.Red. Our proposed models show an overall better
performance than other benchmark GFMs listed in the second group. In
particular, our proposed models outperform the state-of-the-art tree-based and
forest-based forecasting algorithms such as the regression tree, Cubist,
LightGBM, XGBoost, CatBoost and RF across all experimental datasets on all
error metrics except for twelve cases: Rossmann on mean msMAPE (Cubist), M5 on
all error metrics (LightGBM), M5 on median msMAPE (CatBoost and XGBoost),
Tourism Monthly on mean msMAPE, median msMAPE and median MASE (Cubist),
Favorita with covariates on median msMAPE (LightGBM) and on all error metrics
(Cubist, Catboost and XGBoost). This further shows training a global AR model
in the leaf nodes leads to more accurate forecasts rather than simply
considering the average of the training instances in the leaf nodes. We also
see the results of PR and Tree.Lin.Test.Error.Red are the same across the
Favorita dataset without covariates and the M5 dataset on all error metrics.
This means the corresponding SETAR-Tree have only one node which is the parent
node where the parent node has not been further split based on the considered
stopping criterion. The SETAR-Tree and the PR model have overall provided
better results across the above two datasets compared to most of the other
benchmark GFMs. It shows limiting the tree to a single node is a worthwhile
decision in these cases. We see that overall, Forest.Significance.Error.Red
shows the best performance across the majority of datasets compared to all
other considered benchmarks listed in the first and second groups as well as
Tree.Lin.Test.Error.Red.

Overall, our proposed models show better performance when training with
covariates compared to the benchmarks including PR, FFNN and the state-of-the-
art tree-based algorithms. As explained in Section 5.1.1, the overall result
of the Rossmann dataset has considerably improved across our proposed models
after using external covariates on all error metrics while considerably
outperforming all benchmarks. The results of the Kaggle Web Traffic and
Favorita datasets do not show such improvements after using covariates across
our proposed models as well as all other considered benchmarks. However,
across the Kaggle Web Traffic dataset, Tree.Lin.Test.Error.Red shows the
overall best performance compared to the other benchmarks when training with
covariates.

####  5.1.3 Analysis of SETAR-Forest size

Table 4: Results of Forest.Significance.Error.Red for Chaotic Logistic dataset
across four different sizes. The best performing models are highlighted in
boldface

No: of Trees | Mean msMAPE | Median msMAPE | Mean MASE | Median MASE
---|---|---|---|---
5 | 41.65 | 36.57 | 0.676 | 0.644
10 | 41.14 | 36.38 | 0.667 | 0.646
20 | 41.34 | 35.65 | 0.670 | 0.644
50 | 41.25 | 35.60 | 0.668 | 0.647

We analyse the effect of the size of the SETAR-Forest, i.e., the number of
SETAR-Trees in the forest, for the final forecasting accuracy. Table 4 shows
the results of our best SETAR-Forest variant, Forest.Significance.Error.Red,
for the Chaotic Logistic dataset across four different sizes: 5, 10, 20 and 50
for mean msMAPE, median msMAPE, mean MASE and median MASE.

As expected, based on the results, overall, the performance of the SETAR-
Forest gets better when the number of SETAR-Trees is increased. The SETAR-
Forest containing 10 SETAR-Trees shows the best performance on mean msMAPE and
mean MASE and the second best performance on median MASE. Adding more trees
does not lead to considerable accuracy gains. Thus, we deem ten SETAR-Trees
sufficient for a SETAR-Forest to obtain accurate forecasts in our experiments.

###  5.2 Statistical testing results

Table 5: Results of statistical testing Model | $p_{Hoch}$
---|---
Forest.Significance.Error.Red | -
Tree.Lin.Test.Error.Red | $0.003$
LightGBM | $<10^{-30}$
ARIMA | $<10^{-30}$
ETS | $<10^{-30}$
CatBoost | $<10^{-30}$
SETAR | $<10^{-30}$
XGBoost | $<10^{-30}$
STAR | $<10^{-30}$
Cubist | $<10^{-30}$
PR | $<10^{-30}$
Regression Tree | $<10^{-30}$
RF | $<10^{-30}$
FFNN | $<10^{-30}$

Table 5 shows the results of the statistical testing evaluation, namely the
adjusted $p$-values calculated from the Friedman test with Hochberg’s post-hoc
procedure considering a significance level of $\alpha=0.05$ (García et al.,
2010). We consider the msMAPE of each series provided by the benchmarks and
our best SETAR-Tree and SETAR-Forest variants, Tree.Lin.Test.Error.Red and
Forest.Significance.Error.Red, for the eight primary experimental datasets,
namely the Rossmann, Kaggle Web Traffic, Favorita, M5, Tourism Monthly,
Tourism Quarterly, Chaotic Logistic, and Mackey-Glass datasets. The datasets
are not considered with the external covariates during statistical testing as
ETS, ARIMA, SETAR and STAR models are only executed without covariates.

The overall $p$-value of the Friedman rank sum test is less than $10^{-30}$
which is highly significant. Forest.Significance.Error.Red performs the best
on ranking over msMAPE per each series of the eight primary experimental
datasets and thus, it is used as the control method as mentioned in the first
row. A horizontal line is used to separate the models that perform
significantly worse than Forest.Significance.Error.Red. All benchmarks and
Tree.Lin.Test.Error.Red are significantly worse than the control method as
they report $p_{Hoch}$ values less than $\alpha$.

###  5.3 Computational performance

Table 6: Computational times (in minutes) of all models across all
experimental datasets

| Without Covariates | With Covariates
---|---|---
|  Ross- | Kag- | Favo- | M5 | Tour | Tour | Cha- | Mackey- | Ross- | Kag- | Favo-
|  mann | gle | rita |  | (M) | (Q) | otic | Glass | mann | gle | rita
ETS | 7.52 | 5.59 | 8.08 | 14.43 | 5.00 | 1.06 | 0.07 | 0.22 | - | - | -
ARIMA | 163.80 | 11.16 | 123.00 | 25.41 | 47.00 | 6.78 | 0.18 | 0.30 | - | - | -
SETAR | 7.55 | 0.90 | 2.01 | 5.17 | 0.59 | 0.27 | 0.38 | 1.29 | - | - | -
STAR | 99.68 | 70.33 | 133.87 | 140.29 | 10.04 | 1.86 | 4.81 | 5.11 | - | - | -
PR | 0.55 | 0.63 | 0.65 | 1.58 | 0.03 | 0.02 | 0.02 | 0.02 | 0.90 | 0.25 | 0.55
Cubist | 1.51 | 7.40 | 2.55 | 7.44 | 0.52 | 0.27 | 0.10 | 0.06 | 1.86 | 2.16 | 2.73
FFNN | 7.57 | 18.88 | 223.80 | 24.01 | 0.74 | 0.25 | 25.00 | 4.61 | 10.16 | 24.22 | 352.20
Regression Tree | 0.98 | 0.48 | 1.31 | 1.84 | 0.02 | 0.02 | 0.03 | 0.02 | 1.22 | 0.37 | 1.60
CatBoost | 2.13 | 2.17 | 1.12 | 2.07 | 0.24 | 0.13 | 0.15 | 0.43 | 6.07 | 1.20 | 2.99
LightGBM | 3.24 | 12.73 | 5.25 | 8.91 | 0.48 | 5.03 | 3.84 | 6.27 | 48.02 | 5.44 | 12.22
XGBoost | 81.60 | 45.00 | 231.60 | 70.26 | 5.23 | 3.75 | 7.37 | 4.72 | 341.40 | 21.93 | 59.53
RF | 1.49 | 58.30 | 2.76 | 4.71 | 6.89 | 0.10 | 0.28 | 0.12 | 10.43 | 6.44 | 14.29
Tree.Lin.Test.Error.Red | 13.85 | 5.25 | 0.24 | 0.53 | 0.34 | 0.12 | 0.09 | 19.25 | 38.78 | 14.52 | 3.22
Forest.Significance.Error.Red | 124.16 | 68.22 | 35.56 | 5.53 | 7.14 | 1.79 | 0.77 | 103.82 | 310.60 | 175.88 | 300.15

Table 7: Times (in seconds) required by the proposed models to predict a data
point (given the fully trained models). The prediction time of
Forest.Significance.Error.Red is ten times higher than the prediction time of
Tree.Lin.Test.Error.Red as we consider ten SETAR-Trees for a SETAR-Forest.

| Without Covariates | With Covariates
---|---|---
|  Ross- | Kag- | Favo- | M5 | Tour | Tour | Cha- | Mackey- | Ross- | Kag- | Favo-
|  mann | gle | rita |  | (M) | (Q) | otic | Glass | mann | gle | rita
Tree.Lin.Test.Error.Red | 1.20 | 0.65 | 0.01 | 0.01 | 0.29 | 0.32 | 0.06 | 0.16 | 1.52 | 0.96 | 0.87
Forest.Significance.Error.Red | 12.00 | 6.50 | 0.10 | 0.10 | 2.90 | 3.20 | 0.60 | 1.60 | 15.20 | 9.60 | 8.70

To compare the computational performance of the benchmarks and our proposed
models, all experiments are executed in a controlled environment, namely an
Intel(R) Core(TM) i7 processor (2.6GHz) and 32GB of main memory.

Table 6 shows the computational times (in minutes) of all benchmarks and our
proposed tree and forest models across all experimental datasets. The reported
computational times include the model training, hyperparameter tuning and
forecast computation times.

From Table 6, we can see that the simple models such as PR and regression tree
show the lowest computational times. Cubist, CatBoost and LightGBM models
generally show lower computational times compared to XGBoost. Our proposed
SETAR-Tree model, Tree.Lin.Test.Error.Red, overall shows lower computational
times than Cubist, LightGBM and XGBoost across all primary datasets except
Rossmann and Mackey-Glass. The most accurate model which is our proposed
SETAR-Forest model, Forest.Significance.Error.Red, generally shows higher
computational times than the benchmarks as it executes multiple SETAR-Trees
and combines their forecasts to produce the final forecasts where a single
tree internally does a lot of automatic parameter tuning. However, our SETAR-
Forest shows lower computational times across the M5, Tourism Quarterly and
Chaotic Logistic datasets than the LightGBM and XGBoost models. RF shows lower
computational times than the SETAR-Forest model, however, as our SETAR-Forest
is more accurate than the RF, the higher computational times seem justified.

Table 7 shows the times (in seconds) required by the proposed models, once
trained, to predict a data point. Across the majority of datasets, the
proposed SETAR-Tree model takes less than a second for prediction. The
proposed SETAR-Forest model takes a few seconds for the same prediction, given
the trained models. Thus, we deem our models computationally feasible.

Based on our results, we recommend to use our best SETAR-Tree variant,
Tree.Lin.Test.Error.Red, and the best SETAR-Forest variant,
Forest.Significance.Error.Red, over the state-of-the-art tree-based
algorithms.

##  6 Conclusions and future research

Globally trained tree-based models have recently become popular in the
forecasting field due to their contribution to most top-performing methods in
the M5 forecasting competition. However, the models used are general-purpose
models not specific to the forecasting task, and they only consider the
average of the training instances at a leaf node as their prediction.

In this paper, we have proposed a new forecasting-specific globally trained
tree model and a forest model, the SETAR-Tree and SETAR-Forest, which use
time-series-specific splitting and stopping procedures. The SETAR-Tree uses
the underlying concept of SETAR models during node splitting. In contrast to
the state-of-the-art tree-based algorithms, our proposed tree model trains a
global PR model in each leaf node allowing the models to learn the cross-
series information. The SETAR-Tree internally controls the tree depth by
conducting a statistical linearity test and measuring the error reduction
percentage at each node split while making the requirement of external
hyperparameter tuning to a minimal state. The SETAR-Forest model uses a
collection of SETAR-Trees during forecasting where the trees are made diverse
by randomising the significance of the statistical linearity test and the
threshold of the error reduction percentage used during node splitting. The
proposed tree and forest models can also be trained with external covariates.
Across eight experimental datasets, we have shown that our proposed tree and
forest models can significantly outperform a set of state-of-the-art tree-
based GFMs and traditional univariate forecasting models. As our proposed
models show a better forecasting accuracy, with less parameter tuning with a
competitive computational performance compared to the benchmarks, we recommend
SETAR-Tree and SETAR-Forest as strong tree-based models for global time series
forecasting.

From our experiments, we conclude that training a global model in the leaf
nodes often leads to more accurate results in tree-based algorithms rather
than considering the average of the training instances in the leaf nodes in
the forecasting context. We further conclude that considering the significance
of the statistical linearity test and the error reduction percentage gained by
node splitting together is a good option in automatically determining the
maximum tree depth. The SETAR-Forest shows the best performance across the
majority of the datasets as it adequately addresses the data, model, and
parameter uncertainties compared to the individual tree models. The results
get further improved when the individual trees in the forest are more
diversified.

The success of this approach encourages as future work to extend the proposed
tree models to work with multiple regimes, to work with boosting techniques
and to provide probabilistic forecasts. It will be also interesting to study
whether incorporating tree pruning can further increase the forecasting
accuracy of SETAR-Trees. Analysing the effects of training a GFM in the leaf
nodes of the state-of-the-art tree-based algorithms such as LightGBM and
CatBoost is also worth to study.

## Declarations

  * •

Funding - This work was supported by the Australian Research Council under
grant DE190100045, a Facebook Statistics for Improving Insights and Decisions
research award and Monash University Graduate Research funding.

  * •

Competing interests - The authors have no relevant financial or non-financial
interests to disclose.

  * •

Ethics approval - Not applicable

  * •

Consent to participate - Not applicable

  * •

Consent for publication - The paper is the authors’ original work and has not
been published nor has been submitted simultaneously elsewhere. All the
authors have checked the manuscript and have agreed to the submission.

  * •

Availability of data and materials - All experimental datasets are publicly
available at https://github.com/rakshitha123/SETAR_Trees/tree/master/datasets.
All experimental results are mentioned and explained in the manuscript.

  * •

Code availability - All implementations of our work are publicly available at
https://github.com/rakshitha123/SETAR_Trees.

  * •

Authors’ contributions - Conceptualisation: Rakshitha Godahewa, Geoffrey I.
Webb, Christoph Bergmeir and Daniel Schmidt; Methodology: Rakshitha Godahewa
and Christoph Bergmeir; Software: Rakshitha Godahewa and Daniel Schmidt;
Formal analysis and investigation: Rakshitha Godahewa and Christoph Bergmeir;
Writing - original draft preparation: Rakshitha Godahewa; Writing - review and
editing: Rakshitha Godahewa, Geoffrey I. Webb, Christoph Bergmeir and Daniel
Schmidt; Supervision: Geoffrey I. Webb, Christoph Bergmeir and Daniel Schmidt.

## References

  *   * Athanasopoulos et al. (2011) Athanasopoulos, G., Hyndman, R. J., Song, H., & Wu, D. C. (2011).  The tourism forecasting competition.  International Journal of Forecasting, 27(3), 822–844.
  * Athey et al. (2019) Athey, S., Tibshirani, J., & Wager, S. (2019).  Generalized random forests.  The Annals of Statistics, 47(2), 1148–1178.
  * Aznarte and Benítez (2010) Aznarte, J. L. & Benítez, J. M. (2010).  Equivalences between neural-autoregressive time series models and fuzzy systems.  IEEE Transactions on Neural Networks, 21(9), 1434–1444.
  * Aznarte et al. (2007) Aznarte, J. L., Benítez, J. M., & Castro, J. L. (2007).  Smooth transition autoregressive models and fuzzy rule-based systems: Functional equivalence and consequences.  Fuzzy Sets and Systems, 158(24), 2734–2745.
  * Aznarte et al. (2010) Aznarte, J. L., Medeiros, M. C., & Benítez, J. M. (2010).  Linearity testing for fuzzy rule-based models.  Fuzzy Sets and Systems, 161(13), 1836–1851.
  * Bandara et al. (2020) Bandara, K., Bergmeir, C., Campbell, S., Scott, D., & Lubman, D. (2020).  Towards accurate predictions and causal ’what-if’ analyses for planning and policy-making: A case study in emergency medical services demand.  In International Joint Conference on Neural Networks.
  * Bandara et al. (2020) Bandara, K., Bergmeir, C., & Smyl, S. (2020).  Forecasting across time series databases using recurrent neural networks on groups of similar series: A clustering approach.  Expert Systems with Applications, 140, 112896.
  * Bandara et al. (2019) Bandara, K., Shi, P., Bergmeir, C., Hewamalage, H., Tran, Q., & Seaman, B. (2019).  Sales demand forecast in e-commerce using a long short-term memory neural network methodology.  In 26th International Conference on Neural Information Processing, pp. 462–474.
  * Bergmeir and Benítez (2012) Bergmeir, C. & Benítez, J. M. (2012).  On the use of cross-validation for time series predictor evaluation.  Information Sciences, 191, 192–213.
  * Bergmeir et al. (2016) Bergmeir, C., Hyndman, R. J., & Benítez, J. M. (2016).  Bagging exponential smoothing methods using STL decomposition and Box–Cox transformation.  International Journal of Forecasting, 32(2), 303 – 312\.
  * Bergmeir et al. (2012) Bergmeir, C., Triguero, I., Molina, D., Aznarte, J. L., & Benítez, J. M. (2012).  Time series modeling and forecasting using memetic algorithms for regime-switching models.  IEEE Transactions on Neural Networks and Learning Systems, 23(11), 1841–1847.
  * Bojer and Meldgaard (2020) Bojer, C. S. & Meldgaard, J. P. (2020).  Kaggle forecasting competitions: An overlooked learning opportunity.  International Journal of Forecasting, 37(2), 587–603.
  * Box (1953) Box, G. E. P. (1953).  Non-normality and tests on variances.  Biometrika, 40(3-4), 318–335.
  * Box et al. (2015) Box, G. E. P., Jenkins, G. M., Reinsel, G. C., & Ljung, G. M. (2015).  Time series analysis: Forecasting and control.  John Wiley and Sons.
  * Breiman (2001) Breiman, L. (2001).  Random forests.  Machine Learning, 45, 5–32.
  * Breusch and Pagan (1980) Breusch, T. S. & Pagan, A. R. (1980).  The lagrange multiplier test and its applications to model specification in econometrics.  The Review of Economic Studies, 47(1), 239–253.
  * Cerqueira et al. (2019) Cerqueira, V., Torgo, L., Pinto, F., & Soares, C. (2019).  Arbitrage of forecasting experts.  Machine Learning, 108(6), 913–944.
  * Chen and Guestrin (2016) Chen, T. & Guestrin, C. (2016).  XGBoost: a scalable tree boosting system.  In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD ’16, New York, NY, USA, pp. 785–794. Association for Computing Machinery.
  * Chen et al. (2020) Chen, T., He, T., Benesty, M., Khotilovich, V., Tang, Y., Cho, H., Chen, K., Mitchell, R., Cano, I., Zhou, T., Li, M., Xie, J., Lin, M., Geng, Y., & Li, Y. (2020).  xgboost: extreme gradient boosting.  R package version 1.2.0.1.
  * Coulombe (2020) Coulombe, P. G. (2020).  The macroeconomy as a random forest.  CoRR, abs/2006.12724.
  * Coulombe (2021) Coulombe, P. G. (2021).  To bag is to prune.  CoRR, abs/2008.07063.
  * da Rosa et al. (2008) da Rosa, J. C., Veiga, A., & Medeiros, M. C. (2008).  Tree-structured smooth transition regression models.  Computational Statistics and Data Analysis, 52(5), 2469–2488.
  * Dutang and Guibert (2022) Dutang, C. & Guibert, Q. (2022).  An explicit split point procedure in model-based trees allowing for a quick fitting of GLM trees and GLM forests.  Statistics and Computing, 32(6).
  * Epprecht and Veiga (2012) Epprecht, C. & Veiga, A. (2012).  Evaluating the predictability of stock market returns via STARX-Tree models.  China-USA Business Review, 11, 1–21.
  * Friedman et al. (2010) Friedman, J., Hastie, T., & Tibshirani, R. (2010).  Regularization paths for generalized linear models via coordinate descent.  Journal of Statistical Software, 33(1), 1–22.
  * Friedman (2001) Friedman, J. H. (2001).  Greedy function approximation: A gradient boosting machine.  The Annals of Statistics, 29(5), 1189–1232.
  * García et al. (2010) García, S., Fernández, A., Luengo, J., & Herrera, F. (2010).  Advanced non-parametric tests for multiple comparisons in the design of experiments in computational intelligence and data mining: Experimental analysis of power.  Information Sciences, 180(10), 2044–2064.
  * Gelman and Hill (2006) Gelman, A. & Hill, J. (2006).  Data analysis using regression and multilevel/hierarchical models.  Analytical Methods for Social Research. Cambridge University Press.
  * Ghosh et al. (2006) Ghosh, H., Sunilkumar, G., & Prajneshu. (2006).  Self exciting threshold autoregressive models for describing cyclical data.  Calcutta Statistical Association Bulletin, 58(1-2), 115–132.
  * Godahewa et al. (2021) Godahewa, R., Bandara, K., Webb, G. I., Smyl, S., & Bergmeir, C. (2021).  Ensembles of localised models for time series forecasting.  Knowledge-Based Systems, 233, 107518.
  * Godahewa et al. (2021) Godahewa, R., Bergmeir, C., Webb, G. I., Hyndman, R. J., & Montero-Manso, P. (2021).  Monash time series forecasting archive.  In Neural Information Processing Systems Track on Datasets and Benchmarks.
  * Godahewa et al. (2022) Godahewa, R., Bergmeir, C., Webb, G. I., & Montero-Manso, P. (2022).  An accurate and fully-automated ensemble model for weekly time series forecasting.  International Journal of Forecasting.
  * Godahewa et al. (2022) Godahewa, R., Deng, C., Prouzeau, A., & Bergmeir, C. (2022).  A generative deep learning framework across time series to optimise the energy consumption of air conditioning systems.  IEEE Access, 10, 6842–6855.
  * Goodfellow et al. (2016) Goodfellow, I., Bengio, Y., & Courville, A. (2016).  Deep learning.  MIT Press.
  * Google (2017) Google. (2017).  Web traffic time series forecasting.  https://www.kaggle.com/c/web-traffic-time-series-forecasting.
  * Hewamalage et al. (2021a) Hewamalage, H., Bergmeir, C., & Bandara, K. (2021a).  Global models for time series forecasting: A simulation study.  Pattern Recognition, 108441.
  * Hewamalage et al. (2021b) Hewamalage, H., Bergmeir, C., & Bandara, K. (2021b).  Recurrent neural networks for time series forecasting: Current status and future directions.  International Journal of Forecasting, 37(1), 388–427.
  * Hothorn et al. (2006) Hothorn, T., Hornik, K., & Zeileis, A. (2006).  Unbiased recursive partitioning: A conditional inference framework.  Journal of Computational and Graphical Statistics, 15(3), 651–674.
  * Hyndman and Khandakar (2008) Hyndman, R. J. & Khandakar, Y. (2008).  Automatic time series forecasting: The forecast package for R.  Journal of Statistical Software, 27(3), 1–22.
  * Hyndman and Koehler (2006) Hyndman, R. J. & Koehler, A. B. (2006).  Another look at measures of forecast accuracy.  International Journal of Forecasting, 22(4), 679–688.
  * Hyndman et al. (2008) Hyndman, R. J., Koehler, A. B., Ord, J. K., & Snyder, R. D. (2008).  Forecasting with exponential smoothing: The state space approach.  Springer.
  * Januschowski et al. (2020) Januschowski, T., Gasthaus, J., Wang, Y., Salinas, D., Flunkert, V., Bohlke-Schneider, M., & Callot, L. (2020).  Criteria for classifying forecasting methods.  International Journal of Forecasting, 36(1), 167–177.
  * Januschowski et al. (2021) Januschowski, T., Wang, Y., Torkkola, K., Erkkilä, T., Hasson, H., & Gasthaus, J. (2021).  Forecasting with trees.  International Journal of Forecasting.
  * Ju et al. (2019) Ju, Y., Sun, G., Chen, Q., Zhang, M., Zhu, H., & Rehman, M. (2019).  A model combining convolutional neural network and LightGBM algorithm for ultra-short-term wind power forecasting.  IEEE Access, 7, 28309–28318.
  * Kaggle (2015) Kaggle. (2015).  Rossmann store sales.  http://www.kaggle.com/c/rossmann-store-sales.
  * Kaggle (2018) Kaggle. (2018).  Corporación favorita grocery sales forecasting.  https://www.kaggle.com/c/favorita-grocery-sales-forecasting.
  * Ke et al. (2017) Ke, G., Meng, Q., Finley, T., Wang, T., Chen, W., Ma, W., Ye, Q., & Liu, T. (2017).  LightGBM: A highly efficient gradient boosting decision tree.  In Proceedings of the 31st International Conference on Neural Information Processing Systems, NIPS’17, Red Hook, NY, USA, pp. 3149–3157. Curran Associates Inc.
  * Ke et al. (2020) Ke, G., Soukhavong, D., Lamb, J., Meng, Q., Finley, T., Wang, T., Chen, W., Ma, W., Ye, Q., & Liu, T. Y. (2020).  lightgbm: light gradient boosting machine.  R package version 3.1.1.
  * Koning et al. (2005) Koning, A. J., Franses, P. H., Hibon, M., & Stekler, H. O. (2005).  The M3 competition: Statistical tests of the results.  International Journal of Forecasting, 21(3), 397–409.
  * Kuhn and Johnson (2013) Kuhn, M. & Johnson, K. (2013).  Applied predictive modeling.  New York, NY, USA: Springer.
  * Kuhn and Quinlan (2022) Kuhn, M. & Quinlan, R. (2022).  Cubist: rule- and instance-based regression modeling.  R package version 0.4.0.
  * Lefakis et al. (2019) Lefakis, L., Zadorozhnyi, O., & Blanchard, G. (2019).  Efficient regularized piecewise-linear regression trees.  CoRR, abs/1907.00275.
  * Loh (2002) Loh, W. Y. (2002).  Regression trees with unbiased variable selection and interaction detection.  Statistica Sinica, 12(2), 361–386.
  * Loh (2011) Loh, W. Y. (2011).  Classification and regression trees.  WIREs Data Mining and Knowledge Discovery, 1(1), 14–23.
  * Mackey and Glass (1977) Mackey, M. C. & Glass, L. (1977).  Oscillation and chaos in physiological control systems.  Science, 197, 287–289.
  * Makridakis and Hibon (2000) Makridakis, S. & Hibon, M. (2000).  The M3-competition: Results, conclusions and implications.  International Journal of Forecasting, 16(4), 451–476.
  * Makridakis et al. (2018) Makridakis, S., Spiliotis, E., & Assimakopoulos, V. (2018).  The M4 competition: Results, findings, conclusion and way forward.  International Journal of Forecasting, 34(4), 802–808.
  * Makridakis et al. (2022) Makridakis, S., Spiliotis, E., & Assimakopoulos, V. (2022).  The M5 accuracy competition: Results, findings and conclusions.  International Journal of Forecasting, 38(4), 1346–1364.
  * May (1976) May, R. M. (1976).  Simple mathematical models with very complicated dynamics.  Nature, 261, 459–467.
  * Medeiros and Veiga (2000) Medeiros, M. C. & Veiga, A. (2000).  A hybrid linear-neural model for time series forecasting.  IEEE Transactions on Neural Networks and Learning Systems, 11(6), 1402–1412.
  * Montero-Manso et al. (2020) Montero-Manso, P., Athanasopoulos, G., Hyndman, R. J., & Talagala, T. S. (2020).  FFORMA: Feature-based forecast model averaging.  International Journal of Forecasting, 36(1), 86–92.
  * Montero-Manso and Hyndman (2021) Montero-Manso, P. & Hyndman, R. J. (2021).  Principles and algorithms for forecasting groups of time series: Locality and globality.  International Journal of Forecasting, 37(4), 1632–1653.
  * Narayan (2006) Narayan, P. K. (2006).  The behaviour of US stock prices: Evidence from a threshold autoregressive model.  Mathematics and Computers in Simulation, 71(2), 103–108.
  * Narzo et al. (2022) Narzo, A. F. D., Aznarte, J. L., & Stigler, M. (2022).  tsDyn: nonlinear time series models with regime switching.  R package version 10-1.2.
  * Noble (2009) Noble, W. (2009).  How does multiple testing correction work?  Nature Biotechnology, 27, 1135–1137.
  * Oreshkin et al. (2020) Oreshkin, B. N., Carpov, D., Chapados, N., & Bengio, Y. (2020).  N-BEATS: Neural basis expansion analysis for interpretable time series forecasting.  In 8th International Conference on Learning Representations (ICLR).
  * Petropoulos et al. (2018) Petropoulos, F., Hyndman, R. J., & Bergmeir, C. (2018).  Exploring the sources of uncertainty: Why does bagging for time series forecasting work?  European Journal of Operational Research, 268(2), 545 – 554.
  * Pippenger and Goering (1998) Pippenger, M. K. & Goering, G. E. (1998).  Exchange rate forecasting: Results from a threshold autoregressive model.  Open Economies Review, 9, 157–170.
  * Prokhorenkova et al. (2018) Prokhorenkova, L., Gusev, G., Vorobev, A., Dorogush, A. V., & Gulin, A. (2018).  CatBoost: Unbiased boosting with categorical features.  In S. Bengio, H. Wallach, H. Larochelle, K. Grauman, N. Cesa-Bianchi, and R. Garnett (Eds.), Advances in Neural Information Processing Systems, Volume 31. Curran Associates, Inc.
  * Quinlan (1992) Quinlan, J. R. (1992).  Learning with continuous classes.  In 5th Australian Joint Conference on Artificial Intelligence, pp. 343–348. World Scientific.
  * Rambharat et al. (2005) Rambharat, B. R., Brockwell, A. E., & Seppi, D. J. (2005).  A threshold autoregressive model for wholesale electricity prices.  Journal of the Royal Statistical Society: Series C (Applied Statistics), 54(2), 287–299.
  * Salinas et al. (2020) Salinas, D., Flunkert, V., Gasthaus, J., & Januschowski, T. (2020).  DeepAR: Probabilistic forecasting with autoregressive recurrent networks.  International Journal of Forecasting, 36(3), 1181–1191.
  * Smyl (2020) Smyl, S. (2020).  A hybrid method of exponential smoothing and recurrent neural networks for time series forecasting.  International Journal of Forecasting, 36(1), 75–85.
  * Smyl and Kuber (2016) Smyl, S. & Kuber, K. (2016).  Data preprocessing and augmentation for multiple short time series forecasting with recurrent neural networks.  In 36th International Symposium on Forecasting.
  * Spiliotis (2022) Spiliotis, E. (2022).  Decision trees for time-series forecasting.  Foresight: The International Journal of Applied Forecasting, 64, 30–44.
  * Suilin (2017) Suilin, A. (2017).  kaggle-web-traffic.  https://github.com/Arturus/kaggle-web-traffic.
  * Sun et al. (2020) Sun, X., Liu, M., & Sima, Z. (2020).  A novel cryptocurrency price trend forecasting model based on LightGBM.  Finance Research Letters, 32, 101084.
  * Terasvirta (1994) Terasvirta, T. (1994).  Specification, estimation and evaluation of smooth transition autoregresive models.  Journal of the American Statistical Association, 89(425), 208–218.
  * Therneau and Atkinson (2019) Therneau, T. & Atkinson, B. (2019).  rpart: recursive partitioning and regression trees.  R package version 4.1-15.
  * Tong (1978) Tong, H. (1978).  On a threshold model.  Pattern Recognition and Signal Processing, 575–586.
  * Tong (1993) Tong, H. (1993).  Non-linear time series: A dynamical system approach.  Oxford: Clarendon Press.
  * Torgo and Oliveira (2014) Torgo, L. & Oliveira, M. (2014).  Ensembles for time series forecasting.  In Asian Conference on Machine Learning, pp. 360–370.
  * Venables and Ripley (2002) Venables, W. N. & Ripley, B. D. (2002).  Modern applied statistics with S (Fourth ed.).  Springer.
  * Weng et al. (2020) Weng, T., Liu, W., & Xiao, J. (2020).  Supply chain sales forecasting based on lightGBM and LSTM combination model.  Industrial Management & Data Systems, 120(2), 265–279.
  * Zeileis et al. (2008) Zeileis, A., Hothorn, T., & Hornik, K. (2008).  Model-based recursive partitioning.  Journal of Computational and Graphical Statistics, 17(2), 492–514.