§ INTRODUCTION

Forecasting is the procedure of creating predictions based on past and current data. Subsequently, these predictions can be compared against what happened. For example, one might estimate future infections and then compare them against the actual outcomes. Forecasting might refer to specific formal statistical methods employing time-series, cross-sectional or longitudinal data, less traditional judgmental methods or the process of prediction and resolution itself [1].

Although challenging, forecasting time-series is an essential task to which substantial research effort has been applied. Makridakis et al. [6] emphasise two facts about the field: first, no one has prophetic powers to predict the future accurately, and second, all predictions are subject to uncertainty, especially within social contexts.

Forecasting has applications in many fields where estimates of future conditions are helpful. Depending on the area, accuracy varies significantly. If the factors that relate to the forecast are known and well understood and a significant amount of data can be used, the final value will likely be close to the estimates. If this is not the case or if the actual outcome is affected by the forecast, the dependability of the predictions can be significantly lower. Some forecasting applications include Climate change, preventative maintenance, anomalies, stock returns, epidemics, economic trends, and the development of conflict situations [2, 3, 5, 4].

§.§ Model Fusion and Meta-learning

In this section we refine the concepts of meta-learning and model fusion and present a brief framework to better contextualise the later literature review as it relates to the presented forecasting techniques. First, we give a definition of model fusion as it relates to forecasting and how it differs from traditional data fusion. Next we describe what meta-learning is within the context of model fusion for time-series forecasting.

Data fusion is the multi-modal, multi-resolution, multi-temporal process of integrating data sources to produce more consistent, accurate, and useful information than the sources individually. Traditional multi-modal data fusion approaches are grouped into classes, based on the processing level at which the fusion occurs. Early fusion or data-level fusion, late fusion or decision level fusion and intermediate fusion which combines late and early approaches [61].

Analogously, we define model fusion as integrating base learners to produce a lower biased and variance; and more robust, consistent, and accurate model than the learners individually. These base learners can either be homogeneous from the same hypothesis class (e.g. decisions trees in a random forest) or be heterogeneous from different hypothesis classes (e.g. neural network with a support vector machine for the classification). Drawing further from multi-modal data fusion which sub-categorises based on the whether or not the fusion process is early or late we present definitions for early model fusion, late model fusion and incremental model fusion. We define these terms as follows:

Early model fusion integrates the base learners before training. The combined model is then trained as a single fused model.

Late model fusion first trains the base learners individually. The now pre-trained base learners are then integrated without further modification.

Incremental model fusion performs model integration while training the base learners incrementally. Each combined base learner's parameters remains fixed once trained.

An important aspect of model fusion is that of meta-learning. Meta-learning is typically classified into metric-based, model-based, and optimisation-based meta-learning. Of interest in our discussion is model-based meta-learning. Model-based meta-learning uses metadata/meta-features about the problem together with a meta-model to improve overall performance. For instance, stacked generalisation works by using linear-regression to weight heterogeneous base learners [21]. We are now in a position to define the model integration process:

Meta-model fusion uses model-based meta-learning to perform the model integration process.

Aggregation fusion or just aggregation uses a simple aggregation scheme, like weighted averaging, to perform the model integration.

Combing the above discussion with the commonly used appellative for the classes of techniques allows us to arrive at the following incomplete taxonomy shown in Figure <ref>. In the following section we review relevant literature for forecasting model-fusion in the context of the presented taxonomy.

An incomplete taxonomy of model fusion. Showing Ensemble Learning, Boosting, Parallel Ensembles, Stacking, Deep Learning, and Hybrid-models [24, 29].

§.§ Related Literature

Despite the diversity of the forecasting models, the no free lunch theorem holds that no single algorithm universally outperforms any other for all problems [20]. As a result, an essential question is raised about choosing the best algorithm for a specific time-series. These many techniques are underpinned by different principles to perform forecasting. However, numerous forecasting studies have shown that a fusion of models generally outperforms the individuals [12, 13, 14].

Arinze [21] first introduced meta-learning as stacking technique to select one forecasting model among six others based on the learnt performance for six time-series meta-features (e.g., autocorrelation, trend, coefficients). Raftery et al. [22] point out that single model selection is open to issues of model uncertainty and proposed using Bayesian model averaging for stacking instead. Bayesian model averaging makes weighted average predictions as the posterior probability of each model being correct given the same dataset. Stacking can also learn the weights from the predictions of heterogeneous base learners using cross-validation instead of using posterior probabilities [23]. Clarke [23] reported cross-validation based stacking as a more robust method than Bayesian model averaging for computations that involve sensitive changes over their variables.

Stacking methods' performance might also be improved by utilising additional time-series statistics as meta-features in the learning process. Cawood and van Zyl [1], for example, proposed a feature-weighted stacking approach that improves the regression from base learner predictions by supplementing them with time-series statistics extracted from the input series. Recent work has established many valuable statistics for meta-learning, including, but not limited to, the linearity, curvature, autocorrelation, stability, and entropy [25, 26, 27].

Ribeiro and Coelho [24] studied contemporary ensemble methods in forecasting time-series within agribusiness. In their work, they compare approaches of bagging with RF, boosting with gradient boosting machines, XGBoost and a stacked generalisation of machine learning models. Ribeiro and Coelhoa's findings suggest that gradient boosting methods generally produce the lowest prediction errors. The top two submissions of the M4 forecasting competition both implement ensemble learning to produce state of the art results.

The runner-up in the M4 Competition, FFORMA, implements stacking using extreme gradient boosting to learn the weightings for the base learners based on several meta-features extracted from the input time-series [27]. The weightings learned by the trained meta-model are then used to combine the models' predictions as a feature-weighted average.

Hybrid approaches combine linear and nonlinear models since time-series are rarely pure linear or nonlinear in reality [29]. Hybrid forecasting models were first proposed by Zhang [29], who showed that the hybridisation of the ARIMA and MLP models yields improved accuracy since it takes advantage of combining the strength of both linear- and nonlinear models. The authors note that the main advantage of neural networks is their ability to do nonlinear modelling. In their implementation, the MLP learns from the residuals of the linear ARIMA method to make the final integrated predictions. The ES-RNN is a hybrid forecasting model and winning submission of the M4 Competition that merges a modified Holt-Winters and dilated LSTM stacks [28].

Hybrid methods have gained traction in recent time-series forecasting research [30, 31, 32]. The most common approach is to combine deep learning with a linear method, similar to the initial work by Zhang. Numerous hybrids implement ANN or combine them with traditional models [7, 8, 9]. Zhang et al. [10] identify this attraction to ANN as their ability to self-adapt and model underlying relationships within the data, even when the relationships are unknown. However, some have proposed hybridising fuzzy logic techniques with Machine Learning models as well [33, 34]. Further, others have also reported success in integrating Machine Learning models with genetic algorithms like particle swarm, and ant colony optimisation [35, 37].

Hybrid machine Learning techniques are currently the preferred family for research in forecasting time-series, owing to their predictive accuracy [36]. As it stands, N-BEATS achieves state of the art forecasting accuracy by 3% over the ES-RNN [38]. N-BEATS is a hybrid forecasting method that integrates a fully connected neural network with traditional time-series decomposition. N-BEATS proposes a novel doubly residual topology that stacks blocks of different architectures to model the various time-series components. The N-BEATS stacking topology resembles the DenseNet architecture proposed by Huang et al. [40]; however, it has two residual branches that run over both foreword and backward prediction branches of each layer. The neural network thus produces expansion coefficients that are projected to the block basis functions in each layer to produce partial forecasts and backcasts that assist downstream blocks in ignoring components of their inputs that are not useful.

§.§ Makridakis Competitions (M-Competitions)

The Makridakis Competitions (or M-Competitions) are undoubtedly the most influential initiative driving continuous research in improving the aforementioned forecasting methods. The M-competitions are a standard benchmark for establishing an innovatory of forecasting methods [6, 11]. Makridakis et al. [6] point out that machine learning models have delivered poor performances for the forecasting of social sciences when considering the M-competitions prior to the recent M5 Competition.

In the M4 Competition, Makridakis et al. [19] presented a large scale comparative study of $100,000$ time-series and 61 forecasting methods. They compare the accuracy of statistical methods e.g., ARIMA and Holt-Winters versus that of Machine Learning methods e.g., MLP and RNN. Similar to previous studies, they find that utilising hybrid-model fusion of heterogeneous base learners produces the best results [12, 13, 14]. For a comprehensive study of statistical and Machine Learning forecast models, refer to the extensive comparative paper by Makridakis et al. [19].

More recently, the M5 Competition was the first in the series of M-competitions with hierarchical time-series data [11]. Subsequently, the M5 Competition led research efforts to confirm the dominance of meta-model fusion in forecasting homogeneous data. The top entries of the M5 Competition, methods that adopt parallel ensemble learning and gradient boosting produce state of the art forecasting results [15]. Ensembles are known to outperform their constituent models for many machine learning problems [16]. Likewise, integrating models with different architectures has been reported to deliver superior results [17, 18].

For the above reasons, this study returns focused to the M4 Competition to investigate if the results from the M5 Competition relating to fusion strategies hold over the forecasting of univariate time-series.

§.§ Contributions of this Paper

A review of time-series forecasting literature highlights the need to determine the improvement one might expect to achieve by implementing late data fusion with the state of the art hybrid time-series forecasting models. Moreover, research is required to compare the performance of different ensembles using the forecasts from the same pool of forecasting models. There is also a requirement in forecasting research that necessitates one to validate the utility of ensemble learning with an advanced hybrid model. Our work is novel in that it is the first study that adopts late meta-model fusion to integrate traditional base learners with a hybrid-model. We contribute to the forecasting literature by:

* presenting a novel taxonomy for organising the current literature around forecasting model fusion;

* studying the potential improvement of the predictive power of any state of the art forecasting model;

* contrasting the performance of multiple ensembling techniques from different architectures; and

* delivering an equitable comparison of techniques by providing validation results of the ensembles over five runs of ten-fold cross-validation.

§ MATERIALS AND METHODS

In this section we first give a brief overview of the M4 Competition followed by each base learner utilised. The base learners include four statistical ones and the ES-RNN, a hybrid deep learning model. We reused the original base learner forecasts from the M4 submissions as the input forecasts for the ensembles. The later part of this section gives an overview of the different ensemble techniques that we implemented for our experiments and how they differ in their complexity and meta-learning strategies.

§.§.§ The M4 Forecasting Competition

One instance of the M-Competitions is the M4 forecasting competition. The M4 dataset contains $100,000$ time-series of different frequencies, i.e., yearly, quarterly, monthly, weekly, daily and hourly. The minimum number of observations varies for the different subsets, e.g., 13 for the yearly series and 16 for quarterly. The data is freely available on Github [https://github.com/Mcompetitions/M4-methods/tree/master/Dataset] and Table <ref> provides a summary of the number of series per frequency and domain. Economic, Finance, Demographics, and Industries are among the domains, with data from Tourism, Trade, Labor and Wage, Real Estate, Transportation, Natural Resources, and the Environment also included.

Number of series per data frequency and domain [19].

Data subset
Yearly        6,538      3,716      3,903      6,519      1,088     1,236     23,000

Quarterly     6,020      4,637      5,315      5,305      1,858     865       24,000

Monthly       10,975     10,017     10,016     10,987     5,728     277       48,000

Weekly        112        6          41         164        24        12        359

Daily         1,476      422        127        1,559      10        633       4,227

Hourly        0          0          0          0          0         414       414

Total         25,121     18,798     19,402     24,534     8,708     3,437     100,000

The forecast horizon lengths for each frequency are six steps ahead for yearly data, eight steps ahead for quarterly data, 18 steps ahead for monthly data, 13 steps ahead for weekly data, 14 steps ahead for daily data, and 48 forecasts for hourly data. The M4 Competition comprises one essential considerations relating to OWA as a metric for comparison detailed further next.

The OWA computes the average of the two most popular accuracy measures to evaluate the performance of the PFs relative to the Naïve model 2. The OWA, therefore, combines the two metrics of the sMAPE [41] and mean absolute scaled error MASE [42], each calculated as follows:

\begin{equation}
\begin{aligned}
\label{formula:smape}
sMAPE = \frac{1}{h} \sum^{h}_{t=1} \frac{2|Y_t - \hat{Y_t}|}{|Y_t| + |\hat{Y_t}|} \times 100\%,
\end{aligned}
\end{equation}

\begin{equation}
\begin{aligned}
\label{formula:mase}
MASE = \frac{1}{h} \frac{\sum^{n+h}_{t=n+1}|Y_t - \hat{Y_t}|}{\frac{1}{n-m}\sum^{n}_{t=m+1}|Y_t - Y_{t-m}|},
\end{aligned}
\end{equation}

where $Y_t$ is the actual time-series value at time step $t$, $\hat{Y}_t$ is the predicted value of $Y_t$ and $h$ the length of the forecasting horizon. The denominator and scaling factor of the MASE formula are the in-sample mean absolute error from one-step-ahead predictions of the Naïve model 2, and $n$ is the number of data points. The term $m$ defines the time interval between each successive observation, i.e., $12$ for time-series that have a monthly frequency, four for those with a quarterly frequency, $24$ for hourly series and one for the other frequencies that are nonseasonal series. The OWA error is then computed as the average of the MASE and sMAPE errors relative to the Naïve model 2 predictions as follows:

\begin{equation}
\begin{aligned}
\label{formula:owa}
OWA = \frac{1}{2} \ \left( \frac{MASE }{MASE_{Naive_2}} + \frac{sMAPE }{sMAPE_{Naive_2}} \right)
\end{aligned}
\end{equation}

§.§ Statistical Base Learners

This section gives a brief overview of the statistical base learners used in our experimentation.

Auto-ARIMA a standard method for comparing forecast methods' performances. We use the forecasts from an Auto-ARIMA method that uses maximum-likelihood estimation to approximate the parameters [43].

Theta the best method of the M3 competition [44]. Theta is a simple forecasting method that averages the extrapolated Theta-lines, computed from two given Theta-coefficients, applied to the second differences of the time-series.

Damped  Holt's exponential smoothing with a trend component modified with a damping parameter $\phi$ imposed on the trend component [45, 46].

Comb (or COMB S-H-D) the arithmetic average of the three exponential smoothing methods: Single, Holt-Winters and Damped exponential smoothing [47]. Comb was the winning approach for the M2 competition, and was used as a benchmark in the M4 Competition.

§.§ ES-RNN Base Learner

The ES-RNN method uses late fusing of ES models with LSTM networks to produce more accurate forecasts than either approach. Smyl [28] describes the three main elements of the ES-RNN as deseasonalisation plus adaptive normalisation, generation of forecasts and ensembling.

The first element of Smyl's approach normalises then deseasonalises the series on the fly by ES decomposition. ES decomposition also computes the level, seasonality, and second seasonality components to integrate them with the RNN forecasts.

Second, Smyl's methodology allocates the predictions of the time-series trends to the RNN method since it can model nonlinear relationships. Hybrid forecasting models often exploit the benefits of linear and nonlinear methods by integrating them [29, 48].

The final ensembling of the forecasts is from multiple instances of the hybrid model to mitigate parameter uncertainty [49] and take advantage of averaging over model combinations [50]. In Smyl's methodology, the hybrid instances are trained with the same architecture for independent randomly initialised runs and different subsets of the time-series if computationally feasible. The final forecast for a given series is the average of the forecasts produced by the top-N best models.

§.§.§ Preprocessing

In Smyl's methodology, deseasonalisation is achieved on the fly using ES models for time-series that are nonseasonal (yearly and daily frequencies), single-seasonal (monthly, quarterly and weekly frequencies) and double-seasonal (hourly frequency) [51]. In addition, the formulas are updated to remove the linear trend component, which is modelled using the RNN instead. The updated formulas are as follows. For the Nonseasonal models:

\begin{equation}
\begin{aligned}
\label{formula:es_nonseasonal}
l_t = \alpha y_t + (1 - \alpha) l_{t-1},
\end{aligned}
\end{equation}

for the Single-seasonality models:

\begin{equation}
\begin{aligned}
\label{formula:es_singleseasonal}
l_t &= \alpha y_t / s_t + (1 - \alpha) l_{t-1},\\
s_{t + K} &= \beta y_t / l_t + (1 - \beta) s_t,
\end{aligned}
\end{equation}

and the Double-seasonality models:

\begin{equation}
\begin{aligned}
\label{formula:es_doubleseasonal}
l_t &= \alpha y_t / (s_t u_t) + (1 - \alpha) l_{t-1},\\
s_{t + K} &= \beta y_t / (l_t u_t) + (1 - \beta) s_t,\\
u_{t + L} &= \gamma y_t / (l_t s_t) + (1 - \gamma) u_t,
\end{aligned}
\end{equation}

where $y_t$ is the value of the series at time step $t$; $l_t$, $s_t$ and $u_t$ are the level, seasonality and second-seasonality components, respectively; $K$ denotes the number of seasonal observations (i.e., four for quarterly, 12 for monthly and 52 for weekly) and $L$ is the number of double-seasonal observations ($168$ for the hourly frequency data.)

Smyl adopts normalisation using the typical approach of constant size, rolling input and output windows to normalise the level and seasonality components produced by the ES methods in Eqs. <ref> - <ref>. The input window size is defined after experimentation, and the output window size is equal to the length of the forecasting horizon. The values of the input and output windows are then divided by the last value of the level of the input window and, if the series is seasonal, additionally divided by the seasonality component. Lastly, the log function is applied to counter the effects of outliers on the learning process. Furthermore, the time-series domain (e.g., micro, macro and finance) are one-hot encoded and presented as the only meta-features to the RNN model.

§.§.§ Forecasts by the RNN

The RNN receives the normalised, deseasononalised and squashed values of level and seasonality together with the meta-features as inputs. The outputs of the RNN are then integrated to complete the ES-RNN hybridisation in the following way. For the Nonseasonal models:

\begin{equation}
\begin{aligned}
\label{formula:rnn_nonseasonal}
\hat{y}_{t + 1, \dots, t + h} = \exp{(\operatorname{RNN}(\textbf{x}))} \times l_t,
\end{aligned}
\end{equation}

for the Single-seasonal models:

\begin{equation}
\begin{aligned}
\label{formula:rnn_singleseasonal}
\hat{y}_{t + 1, \dots, t + h} = \exp{(\operatorname{RNN}(\textbf{x}))} \times l_t \times s_{t + 1  , \dots  , t + h},
\end{aligned}
\end{equation}

and Double-seasonal models:

\begin{equation}
\begin{aligned}
\label{formula:rnn_doubleseasonal}
\hat{y}_{t + 1, \dots, t + h} = \exp{(\operatorname{RNN}(\textbf{x}))} \times l_t \times s_{t + 1  , \dots  , t + h} \times u_{t + 1 , \dots  , t + h},
\end{aligned}
\end{equation}

where $\operatorname{RNN}(\textbf{x})$ models the linear trend component from the preprocessed input vector $\textbf{x}$, and $h$ is the forecasting horizon. The $l$, $s$ and $u$ components are from the outputs of the ES models during the preprocessing step.

§.§.§ The Architecture

The RNNs of the ES-RNN are constructed from dilated LSTM [52] stacks, and in some cases, followed by a nonlinear layer and always a final linear layer. Smyl refers to the linear layer as the "adapter" layer since it adapts the size of the last layer to the size of the forecasting horizon, or twice the size of the forecasting horizon for prediction interval (PI) models. The dilated LSTMs improve performance using significantly fewer parameters [52]. Smyl also extends the dilated LSTM with an attention mechanism that exposes the hidden states to the weights of the previous states- a horizon equal to the dilation. To achieve this, Smyl embeds a linear two-layer network into the LSTM.

§.§.§ Loss Function and Optimiser

The ES-RNN implements a pinball loss function to fit the models using Stochastic Gradient Descent. The loss is defined as follows:

\begin{equation}
\begin{split}
\label{formula:pinball_loss}
L_t &= ( y_t - \hat{y}_t )\tau  \text{, \ if} \ y_t \geq \hat{y}_t \\
& = (\hat{y}_t  - y_t  )(1 - \tau)  \text{, \ if} \ \hat{y}_t > y_t,
\end{split}
\end{equation}

where $\tau$ is configured typically between $0.45$ and $0.49$. Smyl notes that the pinball function is asymmetric, and that it penalises values outside a quantile range differently to deal with any biasing, and minimising it produces quantile regression

A level variability penalty is implemented as a regulariser to smooth the level values in the loss function. Smyl notes that this drastically improves the performance of the hybrid method as it can concentrate on modelling the trend instead of over-fitting seasonality-related patterns.
The ES-RNN's level variability penalty is found by: 1) compute the log change, i.e., $d_t = \log(y_{t+1}/y_t)$; 2) compute the difference of the changes: $e_t = d_{t+1} - d_t$; 3) square and average the differences; and 4) lastly, the level variability penalty is multiplied by a constant parameter in the range of $50$ - $100$ before adding it to the loss functions.

§.§ AVG

As a simple benchmark, this study implements a model averaging technique that averages the forecasts of the weak learners [53]. The combined forecast for a set of forecast models $M$, might thus be expressed by the following notation:

\begin{align}
\hat{y_t} = \frac{1}{n} \cdot \sum_{m=1}^{n} y_{m_t},
\end{align}

where $n$ is the number of models in $M$ and $y_{m_t}$ is the forecast of a model at time step $t$.

§.§ FFORMA

The FFORMA framework adopts a feature-weighted model averaging strategy [27]. A meta-learner learns how effectively a pool of forecasting models are at their task for different regions of the meta-feature space and then combines their predictions based on the learned model weightings. In addition, the FFORMA implements the gradient tree boosting model of XGBoost, which the authors note is computationally efficient and has shown promising results for other problems [54].

Montero-Manso et al.'s [27] meta-learning methodology use inputs of nine models and features extracted from the time-series, which measures the characteristics of a time-series: including, but not limited to, features of lag, correlation, the strength of seasonality and spectral entropy.

Montero-Manso et al. [27] implement a custom objective function for the XGBoost to minimise. XGBoost requires both a gradient and hessian of the objective function to fit the model. The functions for their model is derived as follows. The term $p_m(f_n)$ is firstly defined as the output of the meta-learner for model $m$. A softmax transformation is applied to the numeric values to compute the model weights as the probability that each model is the best as:

\begin{equation}
\begin{aligned}
\label{formula:fforma_softmax}
w_m(f_n) = \frac{\exp{(p_m(f_n))}}{\sum_m \exp{(p_m(f_n))}},
\end{aligned}
\end{equation}

where $w_m(f_n)$ is the weight produced by the GB meta-learner for base learner $m \in M$. For each $n$ time-series, the contribution of each method for the OWA error is denoted as $L_{nm}$. The weighted average loss function is computed as:

\begin{equation}
\begin{aligned}
\label{formula:fforma_average_loss}
\Bar{L}_n = \sum^{M}_{m=1} w_m(f_n) L_{nm},
\end{aligned}
\end{equation}

where $\Bar{L}_n$ is the weighted average loss over the set of base learners $M$. The gradient of $\Bar{L}_n$ is then computed as follows:

\begin{equation}
\begin{aligned}
\label{formula:fforma_gradient}
G_{nm} = \frac{\partial \Bar{L}_n}{\partial p_m(f_n)} =
w_{nm}(L_{nm} - \Bar{L}_n).
\end{aligned}
\end{equation}

The hessian of the objective function is then finally derived as follows:

\begin{equation}
\begin{aligned}
\label{formula:fforma_hessian}
H_{nm} = \frac{\partial G_n}{\partial p_m(f_n)} \approx w_{nm}(L_{nm}(1 - w_{nm}) - G_{nm}).
\end{aligned}
\end{equation}

In order to minimise the objective function $\Bar{L}$, the functions $G$ and $H$ are passed to XGBoost, and the model's hyperparameters are found using Bayesian optimisation on a limited search space that is determined based on preliminary results and rules-of-thumb.

Montero-Manso et al.'s [27] methodology combines the predictions of their pool of forecasting models using Algorithm <ref>.
The forecasts are produced from the meta-learner's estimated model weightings, given the meta-features of a new time-series. Thus, the fusion is achieved using the vector of weights $w(f_x) \in \mathbb{R}^M$ with a set of $M$ forecasting models for each time-series $x$. Figure <ref> depicts this experiment's FFORMA forecasting pipeline.

The FFORMA forecasting pipeline.

OFFLINE PHASE: TRAINING

$\{x_1, x_2, \ldots, x_N\}:$ N observed time-series from the reference set.

$\{f_1, f_2, \ldots, f_F\}:$ set of F functions to calculate the meta-features.

$\{m_1, m_2, \ldots, m_M\}:$ set of M forecasting models.

Meta-learner model

Prepare the metadata

$n \gets 1$ to $N$
1. Split $x_n$ into training and test series.

2. Calculate  the meta-features $f_n \in F$.

3. Fit each base forecasting method $m \in M$ and generate forecasts.

4. Compute forecast losses $L_{nm}$ over test period.

Train the meta-learner model using the meta-data and base learner forecasting errors, by minimising:

\[ \underset{w}{\argmin{}}\sum^N_{n=1}\sum^M_{m=1} w(f_n)_m L_{nm}.\]

ONLINE PHASE: FORECASTING

The meta-learner model from the offline phase.

$\{x_1, x_2, \ldots, x_N\}:$ N observed time-series from validation set.

$\{m_1, m_2, \ldots, m_M\}:$ set of M forecasting models.

$\{y_1, y_2, \ldots, y_N\}:$  Forecast for each series in test set.

$n \gets 1$ to $N$
1. Calculate meta-features $f_n \in F$.

2. Use meta-learner to produce $w(f_n)$, an M-vector of weights.

3. Generate forecasts for each $m \in M$.

4. Combine the forecasts using $w$.

3-4. FFORMS-G: Select from $M$ the model with the highest allocated $w$ to produce  final forecast.

FFORMA's forecast combination [27].

*FFORMS-G modification of FFORMA steps 3 & 4.

The meta-feature for each time-series is computed by a function $f$. The FFORMA extracts $43$ meta-features from each time-series using the R package called tsfeatures [55]. For nonseasonal time-series, features that only apply to seasonal time-series are set to zero. Unlike the ES-RNN, the domain-specific features supplied with the M4 dataset are not utilised in the FFORMA.

§.§ FFORMS-R

FFORMS-R, the precursor to FFORMA, learns to select a single model from a pool of forecasting models according to their varying performance observed over some meta-data feature space [56]. To this end, FFORMS-R uses a random forest ensemble learner [57] to classify a single model as the most relevant for some time-series features extracted from the reference series. Following the original FFORMS paper, we utilise the Gini impurity to determine the quality of a split. The same meta-features are used as per the FFORMA approach.

§.§ FFORMS-G

The FFORMS-R framework demonstrated the benefit of selecting a single model from a pool of forecasting models [56]. We additionally propose replacing the weighted-averaging of the FFORMA with the selection of the model with the highest allocated weighting. This replacement allows us to compare how effective model selection against weighted model averaging. This change is highlighted in Algorithm <ref>.

§.§ NN-STACK

Cawood and van Zyl [1] proposed model stacking to do forecasting of nonseasonal time-series. Their implementation performs regression over a feature space consisting of both the ensemble's model forecasts and a set of statistics extracted from the time-series. This paper builds on their tudy by experimenting with more meta-features and a more extensive dataset of different domains and seasonality.

The forecasts of a single timestep is combined using regression over both the ensemble's base learner forecasts and the meta-features extracted from the input series. The MLP model performs regression over inputs of the model forecasts and the extracted meta-features, with the target variables as the predicted series' actual values. A neural network of a basic architecture is adopted, and each layer is transformed using the ReLU [60] activation function. The model is fit using mini-batches, and the Adam [58] algorithm with a mean absolute error loss function.

This research extends the original NN-STACK implementation by including a more extensive set of meta-features for use by the FFORMA framework. Feature selection is achieved using Spearman's rank correlation coefficient ($\rho$) [59], of the meta-features' correlation with the change in each model performance. This procedure helps reduce problems with overfitting, as only the most correlating features are included in the ensemble.

§.§ FFORMA-N

We propose an additional MLP meta-learner that takes the time-series statistics as inputs and a one-hot encoded vector of the best performing model as the model's targets. The FFORMA-N is a deep-learning approach to the FFORMA's feature-weighted model averaging methodology. FFORMA-N adopts a softmax activation function in the MLP's output layer to normalise the network's output to a probability distribution, i.e., the probability of each base learner's adequacy to model a time-series. The same meta-features and feature selection procedure are used as per the NN-STACK approach.

The neural network is of a deep architecture and adopts a ReLU activation function at each layer except for the softmax activated output layer. The model is fit using mini-batches and the Adam stochastic gradient descent algorithm with a categorical cross-entropy loss function.

The predictions of the forecasting models are combined, similar to the FFORMA's fusion technique. The forecasts from the base learners are summed after they are weighted according to the meta-learner's estimated probability distribution.

§.§ N-BEATS

Oreshkin et al. [38] reported a 3% accuracy improvement over the forecasts of the ES-RNN. We treat the N-BEATS method as a state of the art benchmark to compare the performance of the ensembles. Therefore, the N-BEATS method is excluded from the ensembles' pool of base learners. N-BEATS is a pure deep learning methodology that does not rely on feature engineering or input-scaling[38].

§ EXPERIMENTAL METHOD AND RESULTS

All ensembles are evaluated using ten-fold cross-validation. This validation process is repeated five times, and the scores are averaged over all fifty validation sets to rule out chance from random initialisation. A pseudorandom number generator [39] is configured to produce the indices for splitting the data into the training and validation sets for each of the five runs. The seed numbers used were one, two, three, four and five.

All algorithms were implemented in Python and run on a 2.60GHz Intel Core i7 PC with 16 16GB RAM and 1,365 MHz Nvidia RTX 2060 GPU with 6GB GDDR6 memory. The source code of our experiments is available on GitHub [https://github.com/Pieter-Cawood/FFORMA-ESRNN].

§.§ Hyper-parameter Tuning

The FFORMA hyperparameters

Hyperparameter     H     D     W     M     Y     Q

n-estimators          2000      2000     2000     1200     1200     2000

min data in leaf      63        200       50      100      100      50

number of leaves      135       94       19       110      110      94

eta                   0.61      0.90     0.46     0.20     0.10     0.75

max depth             61        9       17        28       28       43

subsample             0.49      0.52     0.49     0.50     0.50     0.81

colsample bytree      0.90      0.49     0.90     0.50     0.50     0.49

Subset names H, D, W, M, Y and Q of Tables <ref> - <ref> correspond to the Hourly, Daily, Weekly, Monthly, Yearly and Quarterly data subsets. The ensembles' architectures and hyperparameters were found using the training set of the first fold of the first cross-validation run. The two ensembles that use neural networks (NN-STACK and FFORMA-N) were tuned using backtesting, and the gradient boosting methods (FFORMA and FFORMS-G) were tuned using Bayesian optimisation. The parameter search space for the Bayesian optimisation was limited based on some initial results, and the Gaussian process method was configured to estimate the parameters that minimise the OWA error over $300$ runs of parameter observations.

The FFORMA and FFORMS-G architecture and hyper-parameters are identical and were automatically tuned using Bayesian optimisation over $300$ runs of a limited parameter search space. For the model's training, early stopping with patience ten was employed using the validation set loss. The validation set is a quarter of the training dataset. Table <ref> presents the hyperparameter settings used to train the FFORMA ensembles for the M4 Competition's dataset.

For the FFORMS-R, we used the same hyperparameters to model all subsets of data. The forests were limited to $100$ trees, with each tree having a maximum of sixteen nodes.

The NN-STACK's architectures and hyperparameters were determined using backtesting on the training dataset of the first fold of the first run of cross-validation. The MLP architectures for the Hourly and Weekly subsets were deep with $11$ hidden layers with neuron  in each layer of [$100$, $100$, $100$, $100$, $100$, $50$, $50$, $50$, $50$, $20$ and $20$] sequentially. A more shallow network was used to model the Daily, Monthly, Yearly and Quarterly subsets with three hidden layers with ten neurons each.

Early stopping with patience $15$ and max epochs $300$ was used, except the Daily subset, which had max epochs of $600$. All batch sizes were $225$ except the Daily subset, set at $1200$. Similarly, a learning rate of $0.0003$ was used everywhere except the Daily subset, which used $0.0001$.

Similar to the stacking method, the FFORMA-N's architecture and hyperparameters were determined using backtesting on the training dataset of the first fold of the first run of cross-validation. The MLP architecture was deep, with $11$ hidden layers with number of neurons [$100$, $100$, $100$, $100$, $100$, $50$, $50$, $50$, $50$, $20$ and $20$] sequentially.

Early stopping was configured to terminate the training loop after a configured patience interval of three except for the Daily subset, set at $15$. A maximum number of epochs of $300$ was used except for the Daily subset set to $600$. All batch sizes were $52$ except the Weekly subset at $225$ and the Daily subset, set at $1200$. A learning rate of $0.0001$ was used throughout.

§.§ Detailed Results

The OWA error distributions on the M4 test set. Where FFORMS represents the similar distributions of FFORMS-R and FFORMS-G.

Average OWA errors on the M4 test set. Where lower values are better.

& Q

8cBase Learners

ARIMA    0.577   1.047   0.925     0.903   0.892   0.898    0.899

Comb     1.556   0.981   0.944     0.920   0.867   0.890    0.899

Damped   1.141   0.999   0.992     0.924   0.890   0.893    0.907

ES-RNN   0.440   1.046   0.864     0.836   0.778   0.847    0.825

Theta    1.006   0.998   0.965     0.907   0.872   0.917    0.901

FFORMA                0.415   0.983      0.725   0.800   0.732   0.816    0.788

FFORMS-R              0.423   0.981      0.740      0.817      0.752      0.830       0.805

FFORMS-G$^\ddagger$   0.427      0.984      0.740      0.810   0.745   0.826       0.798

AVG                   0.847      0.985      0.860      0.863      0.804      0.856       0.847

FFORMA-N$^\ddagger$   0.428      0.979   0.718   0.813      0.746      0.828       0.801

NN-STACK              1.427      0.927   0.810      0.833      0.771      0.838       0.819

8cState of the Art

N-BEATS$^\dagger$   -   -   -   0.819   0.758   0.800    0.799

8l$\dagger$ reproduced for comparison [38].

8l$\ddagger$ proposed methods.

Median OWA errors on the M4 test set. Where lower values are better.

Rank [63]

8cBase Learners

ARIMA    0.332    0.745    0.688    0.670    0.612    0.633     9

Comb     1.121    0.718    0.623    0.684    0.595    0.648     8

Damped   0.940    0.727    0.637    0.691    0.602    0.644     9

ES-RNN   0.370    0.764    0.546    0.637    0.550    0.610     5

Theta    0.817    0.723    0.673    0.692    0.617    0.686     11

FFORMA                 0.318    0.723       0.529    0.602    0.491    0.580     1

FFORMS-R               0.311    0.711       0.552       0.615       0.511       0.589     2

FFORMS-G$^\ddagger$    0.328       0.722       0.538    0.610    0.506    0.596        2

AVG                    0.738       0.656    0.632       0.652       0.557       0.619        7

FFORMA-N$^\ddagger$    0.326       0.720       0.539       0.614       0.508       0.594        2

NN-STACK               0.685       0.646    0.606       0.637       0.535       0.595        5

8l$\ddagger$ proposed methods.

Tables <ref>-<ref> presents the average and median OWA error performance for each forecasting method for each of the six data subsets of the M4 forecasting competition. Column "M, Y & Q" presents the mean over the three larger subsets. The inclusion of this mean is to provide a comparison to N-BEATS. The number of time-series of each data subset is provided in parenthesis, the best results are given in bold, and the top two are highlighted in grey.

Egrioglu and Fildes [62] showed that the M4 competition's ranking was flawed and that the average OWA error measurement is calculated from non-symmetric error distributions and suggest that the errors should be ranked based on their median values instead. Egrioglu and Fildes [62] also argue that evaluating the combined results is of little value as different methods are ranked best across different subsets. In order to present these more robust results, we present the median OWA errors in Table <ref>. In addition, Figures <ref> - <ref> depict the OWA error distributions and visualises the data quartiles and extreme values. The visualisations of the error distributions are used to analyse the consistency of each forecast method's accuracy. The violin plots' upper extreme values were clipped to $3.5$, a high OWA measurement representing model failure.

§.§.§ The Hourly Subset

For the Hourly subset, the ES-RNN produced the only stable forecasts amongst the pool of base learners (see Figure <ref>). Subsequently, the NN-STACK failed as an ensembling method and produced the experiment's worst average ensemble result. This failure was partly due to the small dataset ($414$ series) and the most extensive forecasting horizon requirement of $48$ points.

The FFORMA outperformed all other forecasting methods and produced at least twice as good as the AVG benchmark. It was noted here that the ES-RNN's median OWA error might slightly ($0.059$) be improved by utilising the hybrid model with the FFORMA ensemble technique.

§.§.§ The Daily Subset

The distributions of base learner errors (see Figure <ref>) are exceptionally similar. Moreover, this was the only subset where the ES-RNN failed to outperform other base learners. Subsequently, the gradient boosting methods failed to gain the advantage of the ES-RNN's forecasts as per the other subsets, and the NN-STACK method was the most successful ensemble and the best approach.

The results suggest that the NN-STACK method is a more suitable ensemble for situations where it is difficult for the ensemble to learn the weightings of superiority amongst the pool of similar performing base learners.

§.§.§ The Weekly Subset

Considering the OWA distributions depicted in Figure <ref>, the ES-RNN produced the lowest OWA errors amongst the other base learners, with considerably fewer forecasts reaching the extreme value of $3.0$. This exceptional performance of at least one base learner allowed all ensemble methods to outperform all base learners (see Table <ref>).

The FFORMA ensemble improved the ES-RNN's median OWA error by a small ($0.017$) margin. The FFORMA performed the best of all methods, whereas the FFORMS-G was the second-best.

§.§.§ The Monthly Subset

The Monthly subset consisted of $48,000$ time-series and was the largest- and most domain-balanced dataset of all the experiments. Figure <ref> depicts the distribution of the OWA errors, and the ES-RNN performed best amongst the base learners, with a notably smaller number of model failures. These conditions are ideal for analysing the general performance of the ensemble methods, and the results show similar performance in median OWA errors and a slight ($0.037$) improvement when the FFORMA learns from the predictions of the ES-RNN. Nevertheless, considering the median OWA error, the FFORMA remained the superior method and outperformed all other approaches.

§.§.§ The Yearly Subset

The Yearly subset was one of the experiments' larger and nonseasonal datasets. It was noted that all ensembles, except for the AVG benchmark, slightly outperformed the ES-RNN. Similar OWA distributions were observed between the gradient boosting and neural network ensembles (see Figure <ref>.) However, the FFORMA showed the greatest ($0.056$) improvement, and it is the preferred ensemble for time-series with a Yearly seasonality.

§.§.§ The Quarterly Subset

The Quarterly subset is similar to the Yearly one in size, balance and forecast horizon. The only exception is that the Quarterly time-series contain a higher seasonality. Considering the ensembles' OWA distributions visualised in Figure <ref>, similar ensembling performance was observed to that of the ensembling performance on the Yearly subset. Likewise, no base learner was distinguishable in terms of the OWA distributions. However, the more extensive dataset allowed the ensemble methods to avoid the case of the minor Daily dataset (with similar performing base learners), and the FFORMA ensemble produced a $0.03$ improvement over the median accuracy of the ES-RNN.

The N-BEATS approach showed the best average OWA result, and it was able to outperform the FFORMA for this subset only by a small margin of $0.058$.

§.§ Overall Results

N-BEATS outperformed all other ensembles only for the Quarterly subset. Therefore, the utility of N-BEATS might have notably boosted the ensembles' accuracy for the Quarterly subset. The median OWA error gap between the AVG and other ensemble methods (e.g., between the FFORMA: $0.101$) indicates that Machine Learning is a robust solution to late data fusion.

Although all the ensembles outperform both the ES-RNN's and N-Beats average error when considering the OWA error distributions (Figures <ref> - <ref>) of the three larger subsets (Monthly, Yearly and Quarterly), there is still generally a large amount of uncertainty when considering the (Hourly, Daily and Weekly) data sets. Nevertheless, the ensembles were still able to lower the upper quantiles of the OWA error distributions by a notable margin. A more significant improvement was observed in the OWA error distributions of the three smaller subsets (Hourly, Daily and Weekly) (see Figures <ref> - <ref>), where the weak learners showed more remarkable performance (considering their OWA errors.)

The outstanding ensemble learning performance by both our presented methods and previously by N-BEATS provides evidence that:

* the performance of ensemble learning is dependent on the performances of its weak learners;

* ensemble learning's performance is dependent on the diversity of the data (i.e., the similarity of the time-series from the different domains.);

* for smaller subsets, the traditional methods do better but are still outperformed by the ensemble methods and even more so for larger data sets, where more cross-learning can be exploited;

* ensembles of hybrids can still lead to improved performance; and

* we reaffirm that there is no free lunch concerning the selection of ensembling methods and that validation will still be required for model selection.

When considering the ensemble learning methods, we note that the gradient boosting ensembles generally outperformed those of neural networks. Furthermore, both gradient boosting ensembles outperformed all other ensemble methods on all subsets except for the Daily subset. Table <ref> statistically shows that for all datasets, it is always possible to improve a standalone time-series forecast model by adding it to an ensemble with other statistical forecast models. Lastly, since the FFORMA (the second-place M4 submission) utilises the ES-RNN as a base learner, the FFORMA could outperform the standalone ES-RNN (the M4 winner) for all six subsets of data.

§ CONCLUSIONS

Our primary objective was to empirically compare ensemble methods against the state of the art forecasting methods, the ES-RNN and N-BEATS, to determine whether we might 1) boost the accuracy of contemporary hybrid models and 2) single out a method as the state of the art ensembling technique for future research efforts.

The experiment included four traditional base learners, namely ARIMA, Comb, Damped and Theta and an advanced and hybrid machine learning
model, the ES-RNN.
The ensembles of the experiment adopt different ensembling strategies of model averaging, stacking, and selection. We experimented with a random forest meta-learner, namely the FFORMS-R, two gradient boosting ones, namely the FFORMA, FFORMS-G; and two neural network approaches, namely NN-STACK and FFORMA-N; and a naive arithmetic average as a benchmark. We additionally compare the results of the state of the art benchmark, namely N-BEATS using the results reported in their paper.

Regarding the obtained results, we show that the FFORMA is a state of the art ensemble technique that might be used to boost the predictive power of powerful methods like the ES-RNN and N-BEATS. However, considering the case of the Daily subset, the NN-STACK approach might be a more suitable ensemble when the pool of forecasting models has similar performance.

Further, we show that weighted model averaging is a superior ensemble approach as both the FFORMA and FFORMA-N generally outperform versions of their same architecture that do model stacking or selection.

NN-STACK performs regression over single points of the base learner forecasts. NN-STACK, therefore, ignores the temporal element of the data, and the meta-learner instead learns a temporally-blind fusion function. Consequently, we recommend that future research consider investigating the performance of a stacking ensemble that is conscious of the data's temporality. Furthermore, considering the equal-weighted averaging adopted by the top submissions of the M5 Competition, further research is needed to assess the performance of feature weighted model averaging for multivariate time-series data.

Conceptualization, P.C. and T.VZ.; methodology, P.C.; software, P.C.; validation, P.C and T.VZ.; formal analysis, P.C.; investigation, P.C and T.VZ.; resources, P.C.; data curation, P.C.; writing—original draft preparation, P.C.; writing—review and editing, P.C. and T.VZ.; visualization, P.C. and T.VZ; supervision, T.VZ.; project administration, P.C. and T.VZ; All authors have read and agreed to the published version of the manuscript.

This research received no external funding.

Not applicable.

The M4 Competition's dataset and the base learner forecasts can be accessed via the organisers' GitHub repository at <https://github.com/Mcompetitions/M4-methods> (accessed on 23 May 2022.)

The authors declare no conflict of interest.

The following abbreviations are used in this manuscript:

ARIMA     Autoregressive integrated moving average

AVG     Model averaging

ES     Exponential smoothing

ES-RNN     Exponential-smoothing-recurrent neural network

FFORMA     Feature-based forecast model averaging

FFORMS     Feature-based forecast model selection

LSTM     Long short-term memory

MASE     Mean absolute scaled error

MLP     Multi layer perceptron

N-BEATS     Neural-basis expansion analysis

FFORMA-N     Neural network model averaging

NN-STACK     Neural network model stacking

OWA     Overall weighted average

RNN     Recurrent neural network

sMAPE     Symmetric mean absolute percentage error

[1]
P. Cawood and T. L. van Zyl, “Feature-weighted stacking for nonseasonal time
series forecasts: A case study of the covid-19 epidemic curves,” in
2021 8th International Conference on Soft Computing Machine
Intelligence (ISCMI), 2021, pp. 53–59.

[2]Atherfold, J. & Van Zyl, T. “A method for dissolved gas forecasting in power transformers using ls-svm.” 2020 IEEE 23rd International Conference On Information Fusion (FUSION). pp. 1-8 (2020)

[3]Mathonsi, T. & Zyl, T. “Multivariate anomaly detection based on prediction intervals constructed using deep learning.” Neural Computing And Applications. pp. 1-15 (2022)

[4]Timilehin, O. & Zyl, T. “Surrogate Parameters Optimization for Data and Model Fusion of COVID-19 Time-series Data.” 2021 IEEE 24th International Conference On Information Fusion (FUSION). pp. 1-7 (2021)

[5]Freeborough, W. & Zyl, T. “Investigating Explainability Methods in Recurrent Neural Network Architectures for Financial Time Series Data.” Applied Sciences. 12, 1427 (2022)

[6]
S. Makridakis, R. J. Hyndman, and F. Petropoulos, “Forecasting in social
settings: The state of the art,” International Journal of
Forecasting, vol. 36, no. 1, pp. 15–28, 2020.

[7]
T. Mathonsi and T. L. van Zyl, “Prediction interval construction for
multivariate point forecasts using deep learning,” in 2020 7th
International Conference on Soft Computing & Machine Intelligence
(ISCMI).1em plus 0.5em minus 0.4emIEEE, 2020, pp. 88–95.

[8]
S. Laher, A. Paskaramoorthy, and T. L. Van Zyl, “Deep learning for financial
time series forecast fusion and optimal portfolio rebalancing,” in
2021 IEEE 24th International Conference on Information Fusion
(FUSION).1em plus 0.5em minus 0.4emIEEE, 2021, pp. 1–8.

[9]
T. Mathonsi and T. L. van Zyl, “A statistics and deep learning hybrid method
for multivariate time series forecasting and mortality modeling,”
Forecasting, vol. 4, no. 1, pp. 1–25, 2022.

[10]
G. Zhang, B. E. Patuwo, and M. Y. Hu, “Forecasting with artificial neural
networks:: The state of the art,” International journal of
forecasting, vol. 14, no. 1, pp. 35–62, 1998.

[11]
S. Makridakis, E. Spiliotis, and V. Assimakopoulos, “The m5 competition:
Background, organization, and implementation,” International Journal
of Forecasting, 2021.

[12]
N. G. Reich, L. C. Brooks, S. J. Fox, S. Kandula, C. J. McGowan, E. Moore,
D. Osthus, E. L. Ray, A. Tushar, T. K. Yamana et al., “A
collaborative multiyear, multimodel assessment of seasonal influenza
forecasting in the united states,” Proceedings of the National Academy
of Sciences, vol. 116, no. 8, pp. 3146–3154, 2019.

[13]
C. J. McGowan, M. Biggerstaff, M. Johansson, K. M. Apfeldorf, M. Ben-Nun,
L. Brooks, M. Convertino, M. Erraguntla, D. C. Farrow, J. Freeze
et al., “Collaborative efforts to forecast seasonal influenza in the
united states, 2015–2016,” Scientific reports, vol. 9, no. 1, pp.
1–13, 2019.

[14]
M. A. Johansson, K. M. Apfeldorf, S. Dobson, J. Devita, A. L. Buczak,
B. Baugher, L. J. Moniz, T. Bagley, S. M. Babin, E. Guven et al., “An
open challenge to advance probabilistic forecasting for dengue epidemics,”
Proceedings of the National Academy of Sciences, vol. 116, no. 48, pp.
24 268–24 274, 2019.

[15]
S. Makridakis, E. Spiliotis, and V. Assimakopoulos, “M5 accuracy competition:
Results, findings, and conclusions,” International Journal of
Forecasting, 2022.

[16]
L. Rokach, “Ensemble-based classifiers,” Artificial intelligence
review, vol. 33, no. 1, pp. 1–39, 2010.

[17]
Q. Wu, “A hybrid-forecasting model based on gaussian support vector machine
and chaotic particle swarm optimization,” Expert Systems with
Applications, vol. 37, no. 3, pp. 2388–2394, 2010.

[18]
M. Khashei and M. Bijari, “A new class of hybrid models for time series
forecasting,” Expert Systems with Applications, vol. 39, no. 4, pp.
4344–4357, 2012.

[19]
S. Makridakis, E. Spiliotis, and V. Assimakopoulos, “The m4 competition:
100,000 time series and 61 forecasting methods,” International Journal
of Forecasting, vol. 36, no. 1, pp. 54–74, 2020.

[20]
D. H. Wolpert and W. G. Macready, “No free lunch theorems for optimization,”
IEEE transactions on evolutionary computation, vol. 1, no. 1, pp.
67–82, 1997.

[21]
B. Arinze, “Selecting appropriate forecasting models using rule induction,”
Omega, vol. 22, no. 6, pp. 647–658, 1994.

[22]
A. E. Raftery, D. Madigan, and J. A. Hoeting, “Bayesian model averaging for
linear regression models,” Journal of the American Statistical
Association, vol. 92, no. 437, pp. 179–191, 1997.

[23]
B. Clarke, “Comparing bayes model averaging and stacking when model
approximation error cannot be ignored,” Journal of Machine Learning
Research, vol. 4, no. Oct, pp. 683–712, 2003.

[24]
M. H. D. M. Ribeiro and L. dos Santos Coelho, “Ensemble approach based on
bagging, boosting and stacking for short-term prediction in agribusiness time
series,” Applied Soft Computing, vol. 86, p. 105837, 2020.

[25]
A. C. Lorena, A. I. Maciel, P. B. de Miranda, I. G. Costa, and R. B.
Prudêncio, “Data complexity meta-features for regression problems,”
Machine Learning, vol. 107, no. 1, pp. 209–246, 2018.

[26]
S. Barak, M. Nasiri, and M. Rostamzadeh, “Time series model selection with a
meta-learning approach; evidence from a pool of forecasting algorithms,”
arXiv preprint arXiv:1908.08489, 2019.

[27]
P. Montero-Manso, G. Athanasopoulos, R. J. Hyndman, and T. S. Talagala,
“Fforma: Feature-based forecast model averaging,” International
Journal of Forecasting, vol. 36, no. 1, pp. 86–92, 2020.

[28]
S. Smyl, “A hybrid method of exponential smoothing and recurrent neural
networks for time series forecasting,” International Journal of
Forecasting, vol. 36, no. 1, pp. 75–85, 2020, m4 Competition.

[29]
G. P. Zhang, “Time series forecasting using a hybrid arima and neural network
model,” Neurocomputing, vol. 50, pp. 159–175, 2003.

[30]
N. Liu, Q. Tang, J. Zhang, W. Fan, and J. Liu, “A hybrid forecasting model
with parameter optimization for short-term load forecasting of micro-grids,”
Applied Energy, vol. 129, pp. 336–345, 2014.

[31]
J.-Z. Wang, Y. Wang, and P. Jiang, “The study and application of a novel
hybrid forecasting model–a case study of wind speed forecasting in china,”
Applied Energy, vol. 143, pp. 472–488, 2015.

[32]
Y. Qin, K. Li, Z. Liang, B. Lee, F. Zhang, Y. Gu, L. Zhang, F. Wu, and
D. Rodriguez, “Hybrid forecasting model based on long short term memory
network and deep learning neural network for wind signal,” Applied
energy, vol. 236, pp. 262–272, 2019.

[33]
A. Aksoy, N. Öztürk, and E. Sucky, “Demand forecasting for apparel
manufacturers by using neuro-fuzzy techniques,” Journal of Modelling
in Management, 2014.

[34]
W. Deng, G. Wang, and X. Zhang, “A novel hybrid water quality time series
prediction method based on cloud model and fuzzy forecasting,”
Chemometrics and Intelligent Laboratory Systems, vol. 149, pp. 39–49,

[35]
R. Rahmani, R. Yusof, M. Seyedmahmoudian, and S. Mekhilef, “Hybrid technique
of ant colony and particle swarm optimization for short term wind energy
forecasting,” Journal of Wind Engineering and Industrial
Aerodynamics, vol. 123, pp. 163–170, 2013.

[36]
G. R. Shinde, A. B. Kalamkar, P. N. Mahalle, N. Dey, J. Chaki, and A. E.
Hassanien, “Forecasting models for coronavirus disease (covid-19): a survey
of the state-of-the-art,” SN Computer Science, vol. 1, no. 4, pp.
1–15, 2020.

[37]
S. Kumar, S. K. Pal, and R. Singh, “A novel hybrid model based on particle
swarm optimisation and extreme learning machine for short-term temperature
prediction using ambient sensors,” Sustainable Cities and Society,
vol. 49, p. 101601, 2019.

[38]
B. N. Oreshkin, D. Carpov, N. Chapados, and Y. Bengio, “N-beats: Neural basis
expansion analysis for interpretable time series forecasting,” arXiv
preprint arXiv:1905.10437, 2019.

[39]
L. Blum, M. Blum, and M. Shub, “A simple unpredictable pseudo-random number
generator,” SIAM Journal on computing, vol. 15, no. 2, pp. 364–383,

[40]
G. Huang, Z. Liu, L. Van Der Maaten, and K. Q. Weinberger, “Densely connected
convolutional networks,” in Proceedings of the IEEE conference on
computer vision and pattern recognition, 2017, pp. 4700–4708.

[41]
S. Makridakis, “Accuracy measures: theoretical and practical concerns,”
International journal of forecasting, vol. 9, no. 4, pp. 527–529,

[42]
R. J. Hyndman and A. B. Koehler, “Another look at measures of forecast
accuracy,” International journal of forecasting, vol. 22, no. 4, pp.
679–688, 2006.

[43]
R. J. Hyndman and Y. Khandakar, “Automatic time series forecasting: the
forecast package for r,” Journal of statistical software, vol. 27,
no. 1, pp. 1–22, 2008.

[44]
V. Assimakopoulos and K. Nikolopoulos, “The theta model: a decomposition
approach to forecasting,” International journal of forecasting,
vol. 16, no. 4, pp. 521–530, 2000.

[45]
C. C. Holt, “Forecasting seasonals and trends by exponentially weighted moving
averages,” International journal of forecasting, vol. 20, no. 1, pp.
5–10, 2004.

[46]
E. McKenzie and E. S. Gardner Jr, “Damped trend exponential smoothing: a
modelling viewpoint,” International Journal of Forecasting, vol. 26,
no. 4, pp. 661–665, 2010.

[47]
S. Makridakis and M. Hibon, “The m3-competition: results, conclusions and
implications,” International journal of forecasting, vol. 16, no. 4,
pp. 451–476, 2000.

[48]
O. Fathi, “Time series forecasting using a hybrid arima and lstm model,”
Velvet Consulting, 2019.

[49]
F. Petropoulos, R. J. Hyndman, and C. Bergmeir, “Exploring the sources of
uncertainty: Why does bagging for time series forecasting work?”
European Journal of Operational Research, vol. 268, no. 2, pp.
545–554, 2018.

[50]
F. Chan and L. L. Pauwels, “Some theoretical results on forecast
combinations,” International Journal of Forecasting, vol. 34, no. 1,
pp. 64–74, 2018.

[51]
E. S. Gardner Jr, “Exponential smoothing: The state of the art,”
Journal of forecasting, vol. 4, no. 1, pp. 1–28, 1985.

[52]
S. Chang, Y. Zhang, W. Han, M. Yu, X. Guo, W. Tan, X. Cui, M. Witbrock,
M. Hasegawa-Johnson, and T. S. Huang, “Dilated recurrent neural networks,”
arXiv preprint arXiv:1710.02224, 2017.

[53]
G. Claeskens, N. L. Hjort et al., “Model selection and model
averaging,” Cambridge Books, 2008.

[54]
T. Chen and C. Guestrin, “Xgboost: A scalable tree boosting system,” in
Proceedings of the 22nd acm sigkdd international conference on
knowledge discovery and data mining, 2016, pp. 785–794.

[55]
R. J. Hyndman, E. Wang, and N. Laptev, “Large-scale unusual time series
detection,” in 2015 IEEE international conference on data mining
workshop (ICDMW).1em plus 0.5em minus 0.4emIEEE, 2015, pp.

[56]
T. S. Talagala, R. J. Hyndman, G. Athanasopoulos et al., “Meta-learning
how to forecast time series,” Monash Econometrics and Business
Statistics Working Papers, vol. 6, p. 18, 2018.

[57]
A. Liaw, M. Wiener et al., “Classification and regression by
randomforest,” R news, vol. 2, no. 3, pp. 18–22, 2002.

[58]
D. P. Kingma and J. Ba, “Adam: A method for stochastic optimization,”
arXiv preprint arXiv:1412.6980, 2014.

[59]
J. H. Zar, “Spearman rank correlation,” Encyclopedia of biostatistics,
vol. 7, 2005.

[60]
A. F. Agarap, “Deep learning using rectified linear units (relu),”
arXiv preprint arXiv:1803.08375, 2018.

[61]
Michelsanti, Daniel and Tan, Zheng-Hua and Zhang, Shi-Xiong and Xu, Yong and Yu, Meng and Yu, Dong and Jensen, Jesper, “An overview of deep-learning-based audio-visual speech enhancement and separation,”
IEEE/ACM Transactions on Audio, Speech, and Language Processing, col. 29, p. 1368–1396, 2021

[62]
E. EGRIOGLU and R. FİLDES, “A note on the robustness of performance of
methods and rankings for m4 competition,” Turkish Journal of
Forecasting, vol. 4, no. 2, pp. 26–32, 2020.

[63]Schulze, M. “The Schulze method of voting.” ArXiv Preprint ArXiv:1804.02973. (2018)