<figure><img src="image_1.png" /></figure>

HAL
open science

# Pattern matching for multivariate time series forecasting

Noé Lebreton, Julien Ah-Pine, Julien Jacques, Matthieu Neveu

► To cite this version:

Noé Lebreton, Julien Ah-Pine, Julien Jacques, Matthieu Neveu. Pattern matching for multivariate time series forecasting. 2025. hal-05106642v1

HAL Id: hal-05106642

https://hal.science/hal-05106642v1

Preprint submitted on 11 Jun 2025 (v1), last revised 4 Dec 2025 (v2)

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.
# Pattern matching for multivariate time series forecasting

Noé Lebreton<fnref n="1" />,<fnref n="2" />, Julien Ah-Pine<fnref n="3" />,<fnref n="4" />, Julien Jacques<fnref n="1" />, Matthieu Neveu<fnref n="2" />

## Abstract

This article presents a new approach to multivariate time series forecasting. While most existing techniques in the literature focus on forecasting a single time series, forecasting multiple time series is a common goal in many applications. To deal with this, we introduce a new method, Weighted Nearest Neighbours for multivariate time series (*WNN<sub>multi</sub>*). This method forecasts the future of a given series by identifying similar patterns not only in its own past but also in the past of related time series. Once the *k* nearest neighbours are identified, forecasts are made by averaging their future values. We evaluate the proposed approach on several real-world datasets and compare its performance against state-of-the-art forecasting techniques. The results demonstrate that our method achieves comparable or significantly improved performance, showcasing its effectiveness.

**Keywords** Multivariate Time Series, Multi-step forecasting, Pattern matching, *k* nearest neighbours.

## 1 Introduction

Time series are sets of data collected and structured along a temporal dimension, typically recorded at regular intervals. They play a vital role in many sectors, such as economics, healthcare, energy, and commerce, due to their ability to capture the evolution of phenomena over time. Each point in a time series represents an observation recorded at a specific moment, providing a chronological view of a system's or phenomenon's progression. Time series analysis often aims to extract key information such as trends, recurrent patterns, or seasonal behaviours. These elements enhance understanding of the underlying processes and support various applications. Forecasting time series is a central task [1, 2, 3], involving the projection of future values based on past observations. This field faces several challenges, including the presence of noise in the data, uncertainty over long horizons, and heterogeneity among time series.

Time series forecasting is traditionally approached from a univariate perspective, where the goal is to forecast the future of a single series using only its own past values. However, in many real-world scenarios, multiple time series are observed simultaneously, each measuring different

<footnote n="1">¹Université Lumière Lyon 2, Universite Claude Bernard Lyon 1, ERIC, 69007, Lyon, France</footnote>

<footnote n="2">²Data-Prisme, Saint Paul, 97460, France</footnote>

<footnote n="3">³ Université Clermont Auvergne, Clermont Auvergne INP, Mines Saint-Etienne, CNRS, LIMOS, Clermont-Ferrand, 63000, France</footnote>

<footnote n="4">⁴ Université Clermont Auvergne, IRD, CNRS, CERDI, Clermont-Ferrand, 63000, France</footnote>
<figure><img src="image_2.png" /><figcaption>Figure 1: Example of multivariate time series. From top to bottom, the graphs show the price of gold, silver and oil over 20 years. The graph on the right shows the 3 nearest neighbours of a section of the silver serie in black.</figcaption></figure>

but potentially related phenomena. A common example is meteorology, where temperature, rainfall, and wind speed are recorded at the same weather station. If the objective is to forecast temperature, the other variables, such as rainfall and wind speed, may provide valuable additional information.

Traditional methods like AutoRegressive Integrated Moving Average (ARIMA) [4] or eXtreme Gradient Boosting (XGBoost) are often designed for univariate, single-horizon forecasting tasks and struggle to capture complex relationships between multiple series. These models typically aim to learn a function $f(\cdot)$ that maps a sequence of past observations $(x_{T-w}^i, \dots, x_T^i)$ to future values $(\hat{x}_{T+1}^i, \dots, \hat{x}_{T+h}^i)$, where $T$ denotes the current time, $w$ the size of the input window and $h$ the forecast horizon. However, when the dependencies between past and future are unclear or highly non-linear, such models may become less effective.

In certain applications (as in [5]), time series may exhibit correlations and evolve in similar ways. However, this similarity does not necessarily span the entire observation period; it may be limited to specific time intervals. In the literature, the term pattern is often used to describe characteristic shapes or behaviours that time series can present. Some series, while differing in their overall trends, may still share common structures over short periods. Identifying these local patterns allows for the detection of meaningful relationships between time series, even when their global behaviour diverges. Figure 1 illustrates this concept using data from the Swiss National Bank (https://data.snb.ch), showing the evolution of gold and silver prices (in USD per ounce) and oil prices (in USD per barrel) from September 2004 to September 2024. The gold and silver series display a strong correlation, suggesting potential benefits from joint modeling. In contrast, the oil price series appears more volatile and less clearly correlated with the others. However, by focusing on local patterns, we observe short-term similarities across all three series. The graph on the right of Figure 1 highlights in black the last six values of the silver price and identifies the three most similar historical patterns across the dataset. The closest match appears in the gold series, followed by silver itself, and then oil, demonstrating that even seemingly unrelated series can exhibit shared local dynamics.

From a mathematical point of view, we have a set of $d$ time series called $\mathbf{x}^i$ with $i \in \{1, \dots, d\}$:

$$ (\mathbf{x}^1, \mathbf{x}^2, \dots, \mathbf{x}^d)^\mathrm{T} \qquad (1) $$
where $x^i = (x_{t_1}^i, \dots, x_T^i)$. In our context, time series can start at different times, so $t_1^i$ can change between series. On the other hand, the last time step $T$ is common to all the series. In this paper, the objective is to forecast the $d$ time series on a horizon of size $h$:

$$ \hat{x}_h^i = (\hat{x}_{T+1}^i, \dots, \hat{x}_{T+h}^i) \text{ for } i = 1, \dots, d. \qquad (2) $$

The general idea of our approach is to explore the past of all time series in order to identify patterns (i.e., sequences of values) similar to the recent sequence $(x_{T-w+1}^i, \dots, x_T^i)$ for each series $i$, where $w$ denotes the size of the historical window preceding the forecasting point. In [6], the authors proposed the Weighted Nearest Neighbours (WNN) method, originally designed for a univariate framework. We extend this approach to the multivariate context, enabling the capture of interactions between multiple time series. The WNN<sub>multi</sub> method provides a natural and intuitive solution for multi-horizon forecasting ($h > 1$), as it directly yields a forecast vector. Moreover, it avoids the limitations of traditional approaches by bypassing the need for explicit modeling of the relationship between past and future.

The rest of the paper is structured as follows. In Section 2, we present the WNN method for univariate time series and discuss its limitations. Section 3 present WNN<sub>multi</sub>. Several aspects of this approach will be discussed. Section 4 illustrates the methods presented in Section 2 using several real-world datasets. These datasets have various characteristics, such as different forecast horizons, hyperparameter settings, and numbers of time series. Finally, Section 5 concludes the paper.

# 2 Related work

## 2.1 Methods for multivariate time series forecasting

Multivariate time series involve interrelated variables observed over time, requiring analysis that captures both temporal changes and interactions between variables. Statistical models such as Vector AutoRegressive (VAR) and AutoRegressive Integrated Moving Average with exogenous inputs (ARIMAX) offer structured approaches for handling these relationships. The VAR model [7] models $x_{T+1}^i$ as a linear combination of the past $w-1$ observations from all the series, denoted $(x_{T-j}^i)_{1 \le i \le d, 0 \le j < w-1}$. However, its complexity and the risk of overfitting increase as the number of variables grows. ARIMAX, on the other hand, extends the ARIMA model by incorporating exogenous variables into the forecast of a time series. It models $x_{T+1}^i$ as a function of its own past values and errors, as well as the current and past values of one or more external series. However, ARIMAX also assumes linearity, which limits its capacity to model non-linear relationships often found in real-world systems.

Deep learning methods represent a significant advancement in time series forecasting, particularly for large, non-linear datasets or those with long-term dependencies. Long Short-Term Memory (LSTM) networks [8] effectively capture long-term relationships due to their memory mechanism. Transformers [9], with their attention mechanism, handle long-range dependencies, making them ideal for processing long and complex time series. However, they often require large volumes of data to be effective. Additionally, their "black box" nature limits the ability to understand the interactions between variables. For these reasons, we have chosen not to use neural network models.
Pattern matching is a technique used in time series analysis to identify recurrent patterns in historical data for forecasting future values. The underlying principle is simple: if a specific sequence of values has occurred in the past, it is likely that a similar sequence will follow under comparable conditions. One commonly used algorithm for implementation is k Nearest Neighbours (kNN), which identifies historical patterns closest to the target series and uses their future trajectories to forecast future values [10]. The MV-kWNN method [5] adapts the kNN framework to multivariate time series by defining neighbours as a set of synchronous observations across all series at a given time index. This makes it possible to forecast all the variables simultaneously and obtain a forecast matrix of dimension $d \times h$. This approach is computationally efficient as it avoids forecasting future values step by step, where each prediction depends on the previous one. In cases of strong heterogeneity, the inclusion of irrelevant patterns can degrade the distance metrics, leading to reduced forecast accuracy. Another notable approach is Pattern Sequence-based Forecasting (PSF) [11], which operates in two stages: clustering sub-series to assign cluster labels and then performing sequence-based forecasting by identifying similar historical sequences of cluster labels. Extensions to PSF, such as those in [12], adapt the method to multivariate time series.

Traditional methods represent time series as discrete sequences $(x_1^i, \dots, x_T^i)$, which may overlook temporal continuity and smooth variations. Functional Data Analysis (FDA), by modeling series as continuous functions $x(t)$, provides a more accurate representation of trends and correlations. Gaussian processes (GPs) [13] are one way to incorporate the functional nature of time series. This approach offers several advantages, such as the ability to quantify uncertainty through confidence intervals. To address the challenges associated with multiple-horizon forecasting using Gaussian Processes while incorporating multiple time series, the authors of [14] proposed an extension of classical GPs, named MAGMA. Their approach introduces a shared mean process across series, improving forecasts by leveraging collective information. However, this approach faces challenges with heterogeneous time series. MagmaClust [15] tackles this by incorporating clustering, but this introduces additional computational complexity.

## 2.2 Forecasting with Weighted Nearest Neighbours

The methodology proposed in this paper extends to the multivariate case the univariate WNN approach [6], which we give an overview in this section. Consequently, since we are considering univariate time series for the moment, we omit the exponent $i$ of $x^i$. Let start by defining a target pattern:

$$x_{T-w+1,T} = (x_{T-w+1}, \dots, x_T) \quad (3)$$

which corresponds to the $w$ last values of the time series at time $T$. Our goal is to forecast the future of this pattern. The next step is to segment the whole history of the time series into two sets of sub-series. The first set corresponds to past sub-series used to identify similar patterns (boxes with solid lines in Figure 2), defined as:

$$\mathcal{E}_p = \{x_{j-w+1,j}, \text{ with } j = w, \dots, T-h\}, \quad (4)$$

where each $x_{j-w+1,j}$ correspond to:

$$x_{j-w+1,j} = (x_{j-w+1}, \dots, x_j). \quad (5)$$
The second set $\mathcal{E}_f$ is composed of future sub-series:

$$ \mathcal{E}_f = \{\mathbf{x}_{j+1,j+h}, \text{ with } j = w, \dots, T-h\}, \qquad (6) $$

where each $\mathbf{x}_{j+1,j+h}$ is the future pattern of $\mathbf{x}_{j-w+1,j}$ (boxes with dashed lines in Figure 2):

$$ \mathbf{x}_{j+1,j+h} = (x_{j+1}, \dots, x_{j+h}). \qquad (7) $$

Let's now define $\mathcal{D}$ the set of distances $d_j$ between the target pattern $\mathbf{x}_{T-w+1,T}$ and each of the

<figure><img src="image_3.png" /><figcaption>Figure 2: Segmentation and nearest neighbour selection process for WNN.</figcaption></figure>

possible past pattern of $\mathcal{E}_p$:

$$ d_j = \|\mathbf{x}_{T-w+1,T} - \mathbf{x}_{j-w+1,j}\|_2 \qquad (8) $$

with $j = w, \dots, T-h$ and where $\| \cdot \|_2$ refers to the Euclidean distance. To identify the most similar patterns, we select the indexes of the $k$ smallest distances from $\mathcal{D}$, corresponding to the $k$ nearest neighbors of $\mathbf{x}_{T-w+1,T}$:

$$ \mathcal{J}_k = \{j_1, j_2, \dots, j_k\} \qquad (9) $$

where $j_i$ denotes the index of the $j^{th}$ nearest neighbor, denoted by $\mathbf{x}_{j_i-w+1,j_i}$. The forecast can then be obtained as the average of the future of the $k$ nearest neighbors weighted according to their distance:

$$ \hat{\mathbf{x}}_h = \frac{1}{\sum_{j \in \mathcal{J}_k} \alpha_{\mathbf{x}_{j+1,j+h}}} \sum_{j \in \mathcal{J}_k} \alpha_{\mathbf{x}_{j+1,j+h}} \mathbf{x}_{j+1,j+h} \qquad (10) $$

with

$$ \alpha_{\mathbf{x}_{j+1,j+h}} = \frac{1}{(d_j)^2}. \qquad (11) $$

The choice of weights for forecasting is crucial and can be decisive depending on the application. The authors of [16] explore several weighting strategies and show that the weight defined by equation (11) is particularly effective in giving strong importance to the sub-series closest to the target pattern.
### 3 A new method for multivariate time series forecasting:

#### WNN<sub>multi</sub>

Forecasting methods for multivariate time series are still limited and often require specific adaptations, particularly due to the large number of time series. One of the main challenges of multivariate approaches lies in the simultaneous management of multiple series, which also increases the complexity of analyzing interactions between variables. To overcome some of these limitations, we propose a novel non-parametric approach, called WNN<sub>multi</sub>. The goal of this method is to forecast the future of a time series over multiple horizons by leveraging information from several time series, if they are relevant. The core idea of this approach is based on an extension of the classic WNN method, which is based on pattern matching. To achieve this, the search domain for nearest neighbours is significantly expanded. This generalisation makes it possible to use information from multiple time series independently, without imposing synchronization or explicitly modeling the interactions between them. Such an approach can be particularly useful in applications where there are many time series or where the series have a relatively short history.

Multivariate time series can be challenging to compare, as they may have different scales or units. To address this issue, we normalize each time series by dividing it by its maximum value (as in [5]), ensuring that all values range from 0 to 1.

With WNN<sub>multi</sub> (Figure 3), our objective is to forecast the future values of each series x<sup>i</sup>. To achieve this, we identify the k nearest neighbours for each series from the complete set of patterns across all series. Let x<sup>i</sup><sub>T-w+1,T</sub> be the target pattern of series i. We define E<sup>all</sup><sub>p</sub> (12) and E<sup>all</sup><sub>f</sub> (13) as the two sets, past and future.

$$E_{p}^{all} = \{x_{j-w+1,j}^{i}, \text{ with } j = w, \dots, T-h \text{ and } i = 1, \dots, d\}, \quad (12)$$

$$E_{f}^{all} = \{x_{j+1,j+h}^{i}, \text{ with } j = w, \dots, T-h \text{ and } i = 1, \dots, d\}. \quad (13)$$

Let's now define D<sup>all</sup> the set of distances d<sub>j</sub> between the target pattern x<sup>i</sup><sub>T-w+1,T</sub> and each of

<figure><img src="image_4.png" /><figcaption>Figure 3: Segmentation and nearest neighbour selection process for WNN<sub>multi</sub>.</figcaption></figure>
the possible past pattern of $\mathcal{E}_p^{all}$:

$$d_j = \|\mathbf{x}_{T-w+1,T}^i - \mathbf{x}_{j-w+1,j}^{i'}\|_2 \quad (14)$$

where $i' = 1, \dots, d$ and with $j = w, \dots, T-h$. To find the most similar patterns, we extract the indexes of the $k$ smallest distances from $\mathcal{D}^{all}$, corresponding to the $k$ nearest neighbours of $\mathbf{x}_{T-w+1,T}^i$:

$$\mathcal{J}_k^{all} = \{(i'_1, j_1), (i'_2, j_2), \dots, (i'_k, j_k)\} \quad (15)$$

where $j_k$ is the index of the nearest $k^{th}$ neighbour in time, and $i'_k$ is the time series from which this neighbour comes, so that the corresponding pattern is given by $\mathbf{x}_{j_k-w+1,j_k}^{i'_k}$. The forecast is then obtained as the average of the futures of the $k$ nearest neighbours. The futures are weighted according to their distance:

$$\hat{\mathbf{x}}_h^i = \frac{1}{\sum_{(i',j) \in \mathcal{J}_k^{all}} \alpha_{\mathbf{x}_{j+1,j+h}^{i'}} \sum_{(i',j) \in \mathcal{J}_k^{all}} \alpha_{\mathbf{x}_{j+1,j+h}^{i'}} \mathbf{x}_{j+1,j+h}^{i'}} \quad (16)$$

In this formula, the $\alpha_{\mathbf{x}_{j+1,j+h}^{i'}}$ weights are calculated using the same rule as in the univariate case (11) since it gives a significant weight to the nearest neighbours [16].

# 4 Numerical experiments

In this section, we describe the experimental setup, the baseline methods used for comparison, the chosen hyperparameters, the computation times for each approach, and the results obtained on the two datasets.

## 4.1 Experimental setup

In order to compare $WNN_{multi}$ to alternative forecasting methods, temporal cross validation is performed, with 60% of the data is used for training, 20% for hyperparameter tuning, and the remaining 20% for testing. This setup, illustrated in Figure 4, preserves the chronological order of the data. Additionally, we maintain a fixed training window throughout the evaluation to ensure consistent and comparable results over time. We consider the two following evaluation criteria: Root Mean Squared Error (RMSE) and Mean Absolute Error (MAE) (17). These metrics offer distinct perspectives on the accuracy of forecasts. MAE measures the mean absolute error between forecasted and actual values, providing a straightforward interpretation of average accuracy. In contrast, RMSE gives greater weight to larger errors by squaring them, making it more sensitive to large deviations.

$$\text{RMSE} = \sqrt{\frac{1}{h} \|\mathbf{x}_h^i - \hat{\mathbf{x}}_h^i\|_2^2} \quad \text{and} \quad \text{MAE} = \frac{1}{h} \|\mathbf{x}_h^i - \hat{\mathbf{x}}_h^i\|_1 \quad (17)$$

where $\| \cdot \|_1$ correspond to the sum of the absolute values.

To illustrate our method, we use a first dataset that includes a limited number of individuals, representing electricity consumption and price in Spain (used in [5]). This dataset consists of two
<figure><img src="image_5.png" /><figcaption>Figure 4: Example of splitting between training, validation, and test.</figcaption></figure>

primary variables: consumption (measured in MW) and price (expressed in Euro/MWh). The goal is to forecast both electricity consumption and price over a full day. Data is recorded hourly, setting the forecast horizon at $h = 24$ hours. To ensure reliable results, we selected a sufficiently large sample to have an adequate number of days for both the validation and test sets. The dataset covers the period from January 1, 2018, to April 6, 2024. Of these, 1,372 days are used for model training, 457 for hyperparameter validation, and 457 for model evaluation on the test set.

The second dataset used in this study contains daily exchange rates (used in [17]) for eight countries: Australia, the UK, Canada, Switzerland, China, Japan, New Zealand, and Singapore, covering the period from 1990 to 2016. Each country has a total of 7588 observations, providing a substantial historical record for analysis. As this study also focuses on long-term forecasting, the forecast horizon is set to $h = 14$ days. Regarding sample sizes, the training set consists of 4550 days, while 1512 days are allocated for both hyperparameter validation and model evaluation.

## 4.2 Competitors

We evaluate several methods to determine whether our model is competitive with standard approaches commonly used for similar tasks.

For machine learning methods, there are several ways to make predictions. In this article we propose two approaches, one univariate and one multivariate. For the univariate approach, we decided to stay consistent with our pattern matching framework, defining features as fixed-length segments of the time series using a sliding window of size $w$. Each training example is defined as a pair consisting of a segment and its next value:

$$ \{(\mathbf{x}_{1,w}, x_{w+1}), \dots, (\mathbf{x}_{T-w,T-1}, x_T)\} \qquad (18) $$

To produce forecasts of length $h$, we adopt a recursive forecasting strategy [18, 19], in which each predicted value is iteratively used as an input for forecasting the next step, since the models are designed to predict one step at a time. The use of XGBoost and Random Forest in the univariate framework will be noted $XGB_{uni}$ and $RF_{uni}$.

The multivariate approach generalizes the univariate method by incorporating multiple related time series as covariates. The objective of this model is to predict the next value, denoted by
$x_{T+1}^i$, of the time series $i$ from $\{\mathbf{x}_{T-w+1,T}^i\}$ for $j = 1, \dots, d$. This step is performed for each time series, $i = 1, \dots, d$. The training sample for the time series $i$ therefore corresponds to:

$$ \{ (\{\mathbf{x}_{1,w}^1, \dots, \mathbf{x}_{1,w}^d\}, x_{w+1}^i), \dots, (\{\mathbf{x}_{T-w,T-1}^1, \dots, \mathbf{x}_{T-w,T-1}^d\}, x_T^i) \} \quad (19) $$

This formulation yields a total of $w \times d$ explanatory variables, where $d$ is the number of time series. This richer representation enables the multivariate models to use patterns from multiple time series in an asynchronous and independent manner. However, it also increases the dimensionality, which can make model training more complex. As with the univariate approach, forecasts are obtained using a recursive strategy. The use of XGBoost and Random Forest with these covariates is designed by $XGB_{multi}$ and $RF_{multi}$. For consistency and speed, we have decided to set $w$ for machine learning approaches according to the optimal value identified for pattern matching methods. For the XGBoost model, three main hyperparameters are adjusted: the number of estimators, which controls the maximum number of trees, the maximum depth of the trees and the learning rate, which reduces the impact of each weight update to limit overfitting. For Random Forest, we adjusted the number of trees and their maximum depth. The number of $w$ features used in machine learning models is aligned with the window sizes optimized in the pattern matching methods.

For pattern matching approaches, $WNN_{multi}$ is described in Section 3 and MV-kWNN in [5]. Two hyperparameters need to be adjusted: the pattern size ($w$) and the number of neighbours taken into account ($k$).

As MagmaClust does not scale to large datasets, we propose the $NN_{MAGMA}$ method, which selects a small, homogeneous subset of nearest neighbours and then apply MAGMA on this subset of neighbors. This allows MAGMA to be used more efficiently, even on large datasets. Unlike previous approaches that rely on access to future values of other series, $NN_{MAGMA}$ uses pattern matching to clearly distinguish between the past and the future of each series. This separation enables simultaneous forecasting of all series without future information. Therefore, we can forecast large and heterogeneous sets of time series. For $NN_{MAGMA}$, we have opted to use the default hyperparameters (a zero mean function and a Squared Exponential Kernel). The aim of this decision is to concentrate our efforts on the hyperparameters specific to pattern matching ($w$ and $k$) in order to simplify the analysis of the results.

## 4.3 Hyperparameters grid

As mentioned previously, hyperparameters play a key role in model performance, but their number can be substantial, particularly due to the diversity of the methods tested. The ranges of values considered for these hyperparameters were defined by taking into account the specific characteristics of the datasets, such as the forecast horizon, the number of series, and the length of the time series. For electricity data, we test three values for $w$: $W = [12, 24, 48]$, to capture relevant patterns, along with five values for $k$: $K = [5, 10, 20, 50, 100]$. For exchange rate data, the values tested are $W = [7, 14, 28, 42]$. For $k$, we explore values similar to those tested previously. For machine learning methods, we will test these different values: $n_{estimators} \in [50, 100, 150]$, $\max depth \in [5, 10, 15]$) and learning rate $\in [0.01, 0.1, 0.3]$)
## 4.4 Computation time between methods

To evaluate performance in terms of computation time, we selected the most time-consuming configurations and applied them to each exchange rate time series. This approach simulates the scenario with the longest computation time. Table 1 illustrates significant differences between the methods. The fastest approach is MV-kWNN and our approach is very close with a difference of one second. Multi-task GP ($NN_{MAGMA}$) give the worst performance in terms of computation time. On the other hand, concurrent methods systematically take more time.

<table><caption>Table 1: Performance for 5 iterations per model</caption><thead><tr><td>Model</td><td>Seconds</td></tr></thead><tbody><tr><td>$NN_{multi}$</td><td>4</td></tr><tr><td>MV-kWNN</td><td>2.9</td></tr><tr><td>$NN_{MAGMA}$</td><td>10080</td></tr><tr><td>XGBoost (uni)</td><td>20.04</td></tr><tr><td>XGBoost (multi)</td><td>101.64</td></tr><tr><td>RF (uni)</td><td>64.75</td></tr><tr><td>RF (multi)</td><td>761.88</td></tr></tbody></table>

## 4.5 Spanish electricity market data

The optimal hyperparameter combinations are reported in Table 6. We evaluated all the methods on a new test dataset (Table 2) of the same size as the validation set. The previously observed trends are confirmed: for both price and demand forecasting, $NN_{multi}$ achieves the best average and median performance. The selected machine learning models still perform significantly worse, and while the competing pattern-matching method is similar to ours, it delivers slightly lower results. When forecasting time series, it is useful to complement traditional global performance metrics with a visual analysis of the results. Since the proposed methods are based on similarity, we chose to directly examine the differences between the actual time series and the forecasts obtained with $NN_{multi}$. The graphs illustrating price trends are presented in Figure 5, while those related to demand are shown in Figure 6. Each figure includes two plots: one corresponding to the best RMSE and another showing the worst RMSE. For both variables, the best forecasts closely follow the actual series, demonstrating the model's ability to capture the underlying dynamics. Finally, the forecasts with the worst RMSEs show large differences from the real values, revealing the method's limitations in some cases.

## 4.6 Exchange rate data

With this dataset, given that $d = 8$, presenting all the results becomes complex, especially given the large number of methods tested. The optimal hyperparameter combinations are reported in Table 7.

In this section, the results are presented in Table 3 and Table 4, summarising forecast performance across all time series during the test period, using each method with its optimal hyperparameters. Performance varies depending on the time series; for instance, exchange rates in
<table><caption>Table 2: Means/medians(standard deviation) with best hyperparameters for test data</caption><thead><tr><th>Data</th><th>Model</th><th>RMSE</th><th>MAE</th></tr></thead><tbody><tr><td rowspan="7">Prices</td><td>MV-kWNN</td><td>20.42/17.08 (13.13)</td><td>17.16/14.08 (11.72)</td></tr><tr><td>WNN<sub>multi</sub></td><td><strong>20.23/16.96</strong> (13.28)</td><td><strong>16.9/13.85</strong> (11.92)</td></tr><tr><td>NN<sub>MAGMA</sub></td><td>21.9/18.003 (15.18)</td><td>18.37/14.53 (13.71)</td></tr><tr><td>XGBoost (uni)</td><td>26.18/22.07 (14.41)</td><td>21.82/18.52 (12.33)</td></tr><tr><td>XGBoost (multi)</td><td>22.15/17.79 (14.19)</td><td>18.27/14.83 (12.24)</td></tr><tr><td>RF (uni)</td><td>26.85/22.68 (14.26)</td><td>22.27/18.67 (12.18)</td></tr><tr><td>RF (multi)</td><td>25.73/22.15 (14.86)</td><td>21.25/18.09 (12.8)</td></tr><tr><td rowspan="7">Demand</td><td>MV-kWNN</td><td>1694.92/1508.66 (880.82)</td><td>1455.49/1246.21 (832.02)</td></tr><tr><td>WNN<sub>multi</sub></td><td><strong>1636.83/1356.54</strong> (960.4)</td><td><strong>1406.47/1136.81</strong> (909.86)</td></tr><tr><td>NN<sub>MAGMA</sub></td><td>1738.94/1462.031 (1045.88)</td><td>1478.31/1215.98 (968.91)</td></tr><tr><td>XGBoost (uni)</td><td>2324.69/2087.62 (1113.98)</td><td>1964.4/1757.04 (987.12)</td></tr><tr><td>XGBoost (multi)</td><td>2209.73/1910.8 (1264.43)</td><td>1855.04/1581.38 (1111.68)</td></tr><tr><td>RF (uni)</td><td>2371.23/2121.42 (1143.05)</td><td>2006.7/1825.38 (997.47)</td></tr><tr><td>RF (multi)</td><td>2370.23/2133.69 (1295.55)</td><td>1991.66/1757.04 (1126.22)</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>Figure 5: Comparison between the WNN<sub>multi</sub> forecast and the actual electricity price series. The graph on the left shows the day with the best RMSE and the graph on the right shows the day with the worst RMSE.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: Comparison between the WNN<sub>multi</sub> forecast and the actual electricity demand series. The graph on the left shows the day with the best RMSE and the graph on the right shows the day with the worst RMSE.</figcaption></figure>
Japan are forecasted much more accurately than those in Australia. However, when comparing methods within the same time series, the performance differences are less pronounced than those observed with electricity data. The only clear exception is the MV-kWNN method, which gives errors up to twice as large as the other methods for all time series. Since numerical results alone make it difficult to identify the best method, we assess whether there are significant differences in the average performance of each method using an ANOVA on the RMSEs. This statistical test generalises the t-test for comparing multiple means. The results of the ANOVA, shown in Table 5, indicate that for all time series, the p-value exceeds 0.05. In conclusion, no method clearly performs better than the others, except MV-kWNN, which consistently gives worse results.

<table><caption>Table 3: Means/medians (standard deviation) for NN methods</caption><thead><tr><th>Data</th><th>Model</th><th>RMSE</th><th>MAE</th></tr></thead><tbody><tr><td rowspan="3">Australia</td><td>MV-kWNN</td><td>0.016/0.012 (0.012)</td><td>0.014/0.011 (0.011)</td></tr><tr><td>WNN<sub>multi</sub></td><td>0.01/0.0085 (0.006)</td><td>0.0084/0.0071 (0.005)</td></tr><tr><td>NN<sub>MAGMA</sub></td><td>0.01/0.0093 (0.005)</td><td>0.0084/0.008 (0.0038)</td></tr><tr><td rowspan="3">British</td><td>MV-kWNN</td><td>0.029/0.024 (0.019)</td><td>0.027/0.022 (0.019)</td></tr><tr><td>WNN<sub>multi</sub></td><td>0.016/0.013 (0.01)</td><td>0.014/0.0111 (0.01)</td></tr><tr><td>NN<sub>MAGMA</sub></td><td>0.017/0.0093 (0.005)</td><td>0.014/0.012 (0.0078)</td></tr><tr><td rowspan="3">Canada</td><td>MV-kWNN</td><td>0.015/0.011 (0.01)</td><td>0.014/0.011 (0.01)</td></tr><tr><td>WNN<sub>multi</sub></td><td>0.0085/0.0068 (0.005)</td><td>0.0074/0.0068 (0.005)</td></tr><tr><td>NN<sub>MAGMA</sub></td><td>0.0095/0.0088 (0.0042)</td><td>0.0077/0.0071 (0.0034)</td></tr><tr><td rowspan="3">Switzerland</td><td>MV-kWNN</td><td>0.013/0.011 (0.01)</td><td>0.014/0.011 (0.01)</td></tr><tr><td>WNN<sub>multi</sub></td><td>0.0098/0.0076 (0.008)</td><td>0.0085/0.0064 (0.007)</td></tr><tr><td>NN<sub>MAGMA</sub></td><td>0.01/0.009 (0.0065)</td><td>0.0086/0.0075 (0.0051)</td></tr><tr><td rowspan="3">China</td><td>MV-kWNN</td><td>0.0013/0.0009 (0.002)</td><td>0.0011/0.0008 (0.001)</td></tr><tr><td>WNN<sub>multi</sub></td><td>0.00078/0.0004 (0.002)</td><td>0.00058/0.00034 (0.001)</td></tr><tr><td>NN<sub>MAGMA</sub></td><td>0.00076/0.0004 (0.002)</td><td>0.00054/0.00032 (0.0011)</td></tr><tr><td rowspan="3">Japan</td><td>MV-kWNN</td><td>0.00023/0.00016 (0.0002)</td><td>0.00022/0.00014 (0.0002)</td></tr><tr><td>WNN<sub>multi</sub></td><td>0.00012/0.000098 (0.0001)</td><td>0.00011/0.000081 (0.0001)</td></tr><tr><td>NN<sub>MAGMA</sub></td><td>0.00013/0.00011 (0.00008)</td><td>0.00011/0.000092 (0.000068)</td></tr><tr><td rowspan="3">New Zealand</td><td>MV-kWNN</td><td>0.015/0.013 (0.01)</td><td>0.014/0.011 (0.01)</td></tr><tr><td>WNN<sub>multi</sub></td><td>0.01/0.0093 (0.005)</td><td>0.0089/0.0074 (0.005)</td></tr><tr><td>NN<sub>MAGMA</sub></td><td>0.011/0.01 (0.0045)</td><td>0.0092/0.0081 (0.004)</td></tr><tr><td rowspan="3">Singapore</td><td>MV-kWNN</td><td>0.0075/0.0062 (0.005)</td><td>0.007/0.0057 (0.005)</td></tr><tr><td>WNN<sub>multi</sub></td><td>0.0049/0.0041 (0.003)</td><td>0.0042/0.0036 (0.002)</td></tr><tr><td>NN<sub>MAGMA</sub></td><td>0.0049/0.0044 (0.002)</td><td>0.004/0.0037 (0.0016)</td></tr></tbody></table>

# 5 Conclusion

The approach presented in this paper is an efficient method for forecasting several time series at the same time, which has the particularity of exploiting the links that may exist between the series. This capability is particularly useful when complementary sources of information, even from other individuals, contribute to a better understanding of the phenomenon to be forecast. Its flexibility
<table><caption>Table 4: Means/medians (standard deviation) for XGBoost and Random Forest</caption><thead><tr><th>Data</th><th>Model</th><th>RMSE</th><th>MAE</th></tr></thead><tbody><tr><td rowspan="4">Australia</td><td>XGBoost (uni)</td><td>0.012/0.01 (0.0074)</td><td>0.01/0.0086 (0.0063)</td></tr><tr><td>XGBoost (multi)</td><td>0.012/0.01 (0.0083)</td><td>0.011/0.0088(0.0072)</td></tr><tr><td>RF (uni)</td><td>0.012/0.0098 (0.008)</td><td>0.01/0.0086(0.0071)</td></tr><tr><td>RF (multi)</td><td>0.011/0.0098 (0.0066)</td><td>0.0097/0.0084 (0.0056)</td></tr><tr><td rowspan="4">British</td><td>XGBoost (uni)</td><td>0.017/0.013 (0.011)</td><td>0.015/0.011(0.01)</td></tr><tr><td>XGBoost (multi)</td><td>0.016/0.013 (0.011)</td><td>0.014/0.011(0.0099)</td></tr><tr><td>RF (uni)</td><td>0.017/0.013 (0.01)</td><td>0.014/0.011 (0.009)</td></tr><tr><td>RF (multi)</td><td>0.017/0.014 (0.012)</td><td>0.015/0.012 (0.012)</td></tr><tr><td rowspan="4">Canada</td><td>XGBoost (uni)</td><td>0.0087/0.0068 (0.006)</td><td>0.0075/0.0059 (0.0055)</td></tr><tr><td>XGBoost (multi)</td><td>0.0083/0.0071 (0.0048)</td><td>0.0072/0.0058 (0.0043)</td></tr><tr><td>RF (uni)</td><td>0.0084/0.0068 (0.0052)</td><td>0.0072/0.0058(0.0047)</td></tr><tr><td>RF (multi)</td><td>0.0081/0.0069(0.0049)</td><td>0.0069/0.0058(0.0044)</td></tr><tr><td rowspan="4">Switzerland</td><td>XGBoost (uni)</td><td>0.014/0.0068 (0.0099)</td><td>0.012/0.0096(0.0085)</td></tr><tr><td>XGBoost (multi)</td><td>0.015/0.012 (0.014)</td><td>0.013/0.0096 (0.011)</td></tr><tr><td>RF (uni)</td><td>0.013/0.011 (0.0096)</td><td>0.011/0.0089(0.0076)</td></tr><tr><td>RF (multi)</td><td>0.013/0.011 (0.0097)</td><td>0.011/0.0094 (0.0078)</td></tr><tr><td rowspan="4">China</td><td>XGBoost (uni)</td><td>0.00096/0.00041 (0.0026)</td><td>0.00059/0.00033(0.00088)</td></tr><tr><td>XGBoost (multi)</td><td>0.00087/0.00052 (0.002)</td><td>0.00064/0.00046 (0.00071)</td></tr><tr><td>RF (uni)</td><td>0.00086/0.00042 (0.0021)</td><td>0.00056/0.00035(0.00072)</td></tr><tr><td>RF (multi)</td><td>0.00075/0.0004 (0.002)</td><td>0.00054/0.00035 (0.00082)</td></tr><tr><td rowspan="4">Japan</td><td>XGBoost (uni)</td><td>0.00013/0.0001 (0.0001)</td><td>0.00011/0.000084(0.000096)</td></tr><tr><td>XGBoost (multi)</td><td>0.00014/0.00011 (0.0001)</td><td>0.00012/0.000093 (0.000096)</td></tr><tr><td>RF (uni)</td><td>0.00013/0.0001 (0.00011)</td><td>0.00012/0.000084 (0.000098)</td></tr><tr><td>RF (multi)</td><td>0.00013/0.000096 (0.00011)</td><td>0.00011/0.000081 (0.000098)</td></tr><tr><td rowspan="4">New Zealand</td><td>XGBoost (uni)</td><td>0.011/0.0095 (0.0058)</td><td>0.0096/0.0083(0.003)</td></tr><tr><td>XGBoost (multi)</td><td>0.011/0.0097 (0.0062)</td><td>0.0097/0.008 (0.0056)</td></tr><tr><td>RF (uni)</td><td>0.011/0.0094 (0.006)</td><td>0.0096/0.0082 (0.0056)</td></tr><tr><td>RF (multi)</td><td>0.011/0.0093 (0.0062)</td><td>0.0092/0.0078 (0.0057)</td></tr><tr><td rowspan="4">Singapore</td><td>XGBoost (uni)</td><td>0.0052/0.004 (0.0033)</td><td>0.0045/0.0035(0.003)</td></tr><tr><td>XGBoost (multi)</td><td>0.0056/0.0045 (0.0034)</td><td>0.0049/0.004 (0.0032)</td></tr><tr><td>RF (uni)</td><td>0.0053/0.0044 (0.0033)</td><td>0.0046/0.0036 (0.003)</td></tr><tr><td>RF (multi)</td><td>0.0049/0.0038 (0.0032)</td><td>0.0042/0.0032(0.0029)</td></tr></tbody></table>
<table><caption>Table 5: ANOVA test for RMSE</caption><thead><tr><td>Data</td><td>F</td><td>p-value</td></tr></thead><tbody><tr><td>Australia</td><td>0.8</td><td>0.57</td></tr><tr><td>British</td><td>1.15</td><td>0.33</td></tr><tr><td>Canada</td><td>0.65</td><td>0.69</td></tr><tr><td>Switzerland</td><td>1.42</td><td>0.2</td></tr><tr><td>China</td><td>0.67</td><td>0.68</td></tr><tr><td>Japan</td><td>0.13</td><td>0.993</td></tr><tr><td>New Zealand</td><td>0.11</td><td>0.996</td></tr><tr><td>Singapore</td><td>1.3</td><td>0.25</td></tr></tbody></table>

is a major asset, allowing it to adapt to various contexts. This adaptability is coupled with notable simplicity: the method relies entirely on similarity-based principles and does not require complex modeling of the data. Another key advantage lies in its ability to make multi-horizon forecasting more natural. This stands in contrast to traditional machine learning or time series approaches, which often require multiple forecasts. Finally, the results confirm its effectiveness: it outperforms other methods on electricity consumption data and achieves comparable performance on exchange rate data, while remaining simpler and less computationally demanding. This combination of performance and flexibility makes it a promising approach for time series forecasting.

The flexibility of the kNN framework opens the door to several extensions. Firstly, the approach is not limited to Euclidean distance; other similarity measures can be used. For example, Dynamic Time Warping (DTW) may be a useful alternative in contexts where time shifts between patterns are frequent. Secondly, once the nearest neighbours have been identified, various forecasting methods can be employed beyond the weighted average. These include machine learning algorithms trained specifically on the selected neighbours. This flexibility provides a rich exploration space and could improve the adaptability of the model in different datasets.
# A Hyperparameters

## A.1 Spanish electricity market data

<table><caption>Table 6: Best hyperparameters for validation data. Format: w_best / k_best / n_estimators / max_depth / lr, with - for non-applicable values.</caption><thead><tr><th>Data</th><th>Model</th><th>Hyperparameters</th></tr></thead><tbody><tr><td rowspan="7">Prices</td><td>MV-kWNN</td><td>12 / 10 / - / - / -</td></tr><tr><td>WNN<sub>multi</sub></td><td>12 / 10 / - / - / -</td></tr><tr><td>NN<sub>MAGMA</sub></td><td>24 / 5 / - / - / -</td></tr><tr><td>XGBoost (uni)</td><td>(12) / - / 100 / 5 / 0.1</td></tr><tr><td>XGBoost (multi)</td><td>(12) / - / 150 / 5 / 0.1</td></tr><tr><td>RF (uni)</td><td>(12) / - / 50 / 0.1 / -</td></tr><tr><td>RF (multi)</td><td>(12) / - / 150 / 15 / -</td></tr><tr><td rowspan="7">Demand</td><td>MV-kWNN</td><td>12 / 10 / - / - / -</td></tr><tr><td>WNN<sub>multi</sub></td><td>24 / 20 / - / - / -</td></tr><tr><td>NN<sub>MAGMA</sub></td><td>24 / 5 / - / - / -</td></tr><tr><td>XGBoost (uni)</td><td>(12) / - / 150 / 10 / 0.1</td></tr><tr><td>XGBoost (multi)</td><td>(12) / - / 150 / 10 / 0.1</td></tr><tr><td>RF (uni)</td><td>(12) / - / 100 / 15 / -</td></tr><tr><td>RF (multi)</td><td>(12) / - / 50 / 15 / -</td></tr></tbody></table>

## A.2 Exchange rate data
<table><caption>Table 7: Best hyperparameters per country and model. Format: `w_best/k_best/n_estimators/max_depth/lr`, with - for non-applicable values.</caption><thead><tr><th>Model</th><th>Australia</th><th>British</th><th>Canada</th><th>Switzerland</th></tr></thead><tbody><tr><td>MV-kWNN</td><td>7/10/-/-/-</td><td>7/10/-/-/-</td><td>7/10/-/-/-</td><td>7/10/-/-/-</td></tr><tr><td>WNN<sub>multi</sub></td><td>7/100/-/-/-</td><td>7/50/-/-/-</td><td>7/50/-/-/-</td><td>7/10/-/-/-</td></tr><tr><td>NN<sub>MAGMA</sub></td><td>7/20/-/-/-</td><td>7/50/-/-/-</td><td>7/20/-/-/-</td><td>7/50/-/-/-</td></tr><tr><td>XGBoost (uni)</td><td>(7)/-50/10/0.3</td><td>(7)/-50/5/0.3</td><td>(7)/-50/5/0.3</td><td>(7)/-100/10/0.1</td></tr><tr><td>XGBoost (multi)</td><td>(7)/-150/5/0.1</td><td>(7)/-50/5/0.1</td><td>(7)/-100/5/0.1</td><td>(7)/-100/5/0.1</td></tr><tr><td>RF (uni)</td><td>(7)/-50/10/-</td><td>(7)/-100/15/-</td><td>(7)/-100/10/-</td><td>(7)/-50/15/-</td></tr><tr><td>RF (multi)</td><td>(7)/-50/15/-</td><td>(7)/-150/5/-</td><td>(7)/-50/15/-</td><td>(7)/-50/10/-</td></tr></tbody></table><table><thead><tr><th>Model</th><th>China</th><th>Japan</th><th>New Zealand</th><th>Singapore</th></tr></thead><tbody><tr><td>MV-kWNN</td><td>28/10/-/-/-</td><td>7/10/-/-/-</td><td>7/10/-/-/-</td><td>7/10/-/-/-</td></tr><tr><td>WNN<sub>multi</sub></td><td>7/100/-/-/-</td><td>7/100/-/-/-</td><td>7/10/-/-/-</td><td>7/10/-/-/-</td></tr><tr><td>NN<sub>MAGMA</sub></td><td>7/50/-/-/-</td><td>7/50/-/-/-</td><td>7/10/-/-/-</td><td>7/50/-/-/-</td></tr><tr><td>XGBoost (uni)</td><td>(7)/-100/10/0.1</td><td>(7)/-100/5/0.1</td><td>(7)/-100/5/0.1</td><td>(7)/-100/5/0.1</td></tr><tr><td>XGBoost (multi)</td><td>(7)/-50/10/0.3</td><td>(7)/-150/5/0.1</td><td>(7)/-50/5/0.1</td><td>(7)/-100/5/0.1</td></tr><tr><td>RF (uni)</td><td>(7)/-150/15/-</td><td>(7)/-50/10/-</td><td>(7)/-150/10/-</td><td>(7)/-50/10/-</td></tr><tr><td>RF (multi)</td><td>(7)/-50/15/-</td><td>(7)/-150/15/-</td><td>(7)/-50/10/-</td><td>(7)/-150/15/-</td></tr></tbody></table>

## References

[1] R. Hyndman and G. Athanasopoulos. *Forecasting: Principles and Practice*. OTexts, Australia, 3rd edition, 2021.

[2] S. Makridakis, S. Wheelwright, and R. Hyndman. *Forecasting: Methods and Applications*, volume 35. 1984.

[3] J. Armstrong. *Principles of Forecasting: A Handbook for Researchers and Practitioners*, volume 69. 2001.

[4] G.E.P. Box, G.M. Jenkins, G.C. Reinsel, and G.M. Ljung. *Time series analysis: forecasting and control*. John Wiley & Sons, 2015.

[5] R. Talavera-Llames, R. Pérez-Chacón, A. Troncoso, and F. Martínez-Álvarez. Mv-kwnn: A novel multivariate and multi-output weighted nearest neighbours algorithm for big data time series forecasting. *Neurocomputing*, 353:56-73, 2019.

[6] R. Talavera-Llames, R. Pérez-Chacón, M. Martínez Ballesteros, A. Troncoso, and F. Martínez-Álvarez. A nearest neighbours-based algorithm for big time series data forecasting. volume 9648, pages 174-185, 2016.

[7] C. A. Sims. Macroeconomics and reality. *Econometrica*, 48(1):1-48, 1980.

[8] S. Hochreiter and J. Schmidhuber. Long Short-Term Memory. *Neural Computation*, 9(8):1735-1780, 1997.
[9] A. Vaswani, N. Shazeer, N. Parmar, J. Uszkoreit, L. Jones, A. N Gomez, Ł. Kaiser, and I. Polosukhin. Attention is all you need. *Advances in neural information processing systems*, 30, 2017.

[10] A. Troncoso, J. M. R. Santos, A.G. Exposito, J.L.M. Ramos, and J. C. Riquelme. Electricity market price forecasting based on weighted nearest neighbors techniques. *IEEE Transactions on Power Systems*, 22(3):1294-1301, 2007.

[11] F. Martinez Alvarez, A. Troncoso, J.C. Riquelme, and J.S. Aguilar Ruiz. Energy time series forecasting based on pattern sequence similarity. *IEEE Transactions on Knowledge and Data Engineering*, 23(8):1230-1243, 2011.

[12] Y. Fujimoto and Y. Hayashi. Pattern sequence-based energy demand forecast using photo-voltaic energy records. In *2012 International Conference on Renewable Energy Research and Applications (ICRERA)*, pages 1-6, 2012.

[13] C. E. Rasmussen and C. K. I. Williams. Gaussian processes for machine learning. Adaptive computation and machine learning. MIT Press, 2006.

[14] A. Leroy, P. Latouche, B. Guedj, and S. Gey. Magma: inference and prediction using multi-task gaussian processes with common mean. *Machine Learning*, 111(5):1821-1849, 2022.

[15] A. Leroy, P. Latouche, B. Guedj, and S. Gey. Cluster-specific predictions with multi-task gaussian processes. *Journal of Machine Learning Research*, 24(5):1-49, 2023.

[16] R. Talavera-Llames, R. Pérez-Chacón, A. Troncoso, and F. Martínez-Álvarez. Big data time series forecasting based on nearest neighbours distributed computing with spark. *Knowledge-Based Systems*, 161:12-25, 2018.

[17] Guokun Lai, Wei-Cheng Chang, Yiming Yang, and Hanxiao Liu. Modeling long- and short-term temporal patterns with deep neural networks. *The 41st International ACM SIGIR Conference on Research & Development in Information Retrieval*, 2017.

[18] S. Ben Taieb and G. Bontempi. Recursive multi-step time series forecasting by perturbing data. In *2011 IEEE 11th International Conference on Data Mining*, pages 695-704, 2011.

[19] L.J. Herrera, H. Pomares, I. Rojas, A. Guillén, A. Prieto, and O. Valenzuela. Recursive prediction for long term time series forecasting using advanced models. *Neurocomputing*, 70(16):2870-2880, 2007.