text:
\section{Introduction}\label{sec_introduction}
In August 2015, the United States Department of Transportation (DoT) released a test plan \cite{dot_testplan} for assessing the coexistence of Dedicated Short Range Communications (DSRC) and unlicensed devices in the 5.850-5.925 GHz (5.9 GHz) band. As suggested by the Congress in September 2015, the Federal Communications Commission (FCC), in its latest public notice \cite{fcc1668a1}, now seeks to refresh the record of its pending 5.9 GHz rulemaking to provide potential sharing solutions between Wi-Fi and DSRC at 5.9 GHz.

Wi-Fi provides short-range, high-speed, unlicensed wireless connections. DSRC uses short-range wireless communication links to facilitate information transfer for vehicle to infrastructure (V2I) and vehicle to vehicle (V2V), based on IEEE 802.11p for its PHY and MAC layers. The current focus of the FCC's solicitation in \cite{fcc1668a1} is two-fold: (i) prototype of interference-avoiding devices for testing; (ii) test plans to evaluate electromagnetic compatibility between unlicensed devices and DSRC.

Recent work discusses coexistence between DSRC system and Wi-Fi \cite{infocom_12}-\cite{gaurang17}. Especially the results of experiments and simulations in \cite{gaurang17} suggest a method of allocating channels for DSRC and Wi-Fi. However, despite the thorough study, the implications on the coexistence problem are unreliable due to discrepancy between the experimental and simulation results.

This leads to a conclusion that we need a method to accurately describe the performance of a DSRC network to accomplish the coexistence test plans solicited in \cite{fcc1668a1}. Another recent study \cite{priyabrata_twc16} provides an extensive analysis for a DSRC network, yet it lacks consideration of DSRC-Wi-Fi coexistence. Moreover, we found that using the current metrics makes it difficult to precisely evaluate the performance of a DSRC network. The main objective of DSRC networking is to support safety-critical applications that utilize a functionality called basic safety message (BSM). BSMs are periodically broadcast from a vehicle, which necessitates a unique metric to compute the network-wise performance. Many of the prior studies including \cite{infocom_6}-\cite{priyabrata_tits16} rely on the typical metrics that only myopically capture the broadcast nature of DSRC networks. Examples of these metrics are packet delivery rate (PDR) and packet delay/latency.

An advanced metric that addresses this problem is proposed in \cite{vanet06}, namely inter-reception time (IRT). The typical packet latency is defined as the time spent by a successful packet to travel from its source to final destination. But this metric is not suitable to capture the performance of broadcast-based safety applications since latency is measured only for successful packets. in other words, such a latency does not capture the impact of packet losses and collisions on the latency perceived by applications that are based on periodic broadcast of packets. Defined as the time elapsed between two successive successful receptions of packets broadcast by a specific transmitter (TX) vehicle, IRT can more accurately display the performance of a BSM-based safety-critical application on a DSRC network.

\begin{figure}[t]
\vspace{-0.3in}
\centering
\includegraphics[width = 0.5\linewidth]{fig_geometry_topology.pdf}
\caption{Geometric analysis of a broadcast-based vehicular network}
\label{fig_geometry_topology}
\vspace{-0.1in}
\end{figure}

However, IRT still shows shortcomings. Limitations of IRT are highlighted when taking a geometric view to a broadcast-based vehicular network, as depicted in Fig. \ref{fig_geometry_topology}. For instance, when a packet is transmitted from a TX vehicle, all the vehicles located in the transmission range of it become potential receivers (RXs). When a packet collides with a packet from another vehicle that is at the ``left'' side of the TX vehicle, the vehicles that are located at the ``right'' side are still able to receive the packet without corruption. This geometric insight cannot be captured using the aforementioned classical metrics.

A one-dimensional geometric analysis is provided in \cite{elsevier14}. However, this paper extends the framework to a two-dimensional space, and proposes a metric that can capture partial receptions when a packet is collided. More importantly, this paper has a clear contribution to addressing the 5.9 GHz coexistence problem. We apply the proposed metric for suggesting proper values of MAC parameters when DSRC coexists with Wi-Fi--namely inter-broadcast interval (IBI) and contention window (CW).

In this paper, a new metric is therefore proposed that provides a geometric point of view to measure the performance of a broadcast-based application on a DSRC network. Based on this new method, we characterize DSRC-Wi-Fi coexistence. The contributions of this paper can be formulated as:
\begin{enumerate}
\item{This paper proposes a novel metric that evaluates the performance of a broadcast-based DSRC network. The main advantage of our method is that it can illustrate a wider perspective than the typical metrics that provide only myopic observation at a RX. Considering the broadcast nature, some RXs are still able to receive a packet from a TX even if other RXs are not. The new metric can capture such partial possibility of packet reception, which leads to a more rigorous analysis}
\item{Considering that a significant challenge that DSRC faces is co-channel coexistence with Wi-Fi, we apply the novel performance metric for characterization of DSRC-Wi-Fi coexistence. Our network simulation results suggest adequate IBI and CW values to improve the performance of a broadcast-based safety-critical application.}
\end{enumerate}

\section{System Model}\label{sec_model}
\begin{figure}[t]
\vspace{-0.3in}
\centering
\includegraphics[width = 0.7\linewidth]{fig_simulation_topology_5page_wifi.pdf}
\caption{Geometry of the simulation environment}
\label{fig_simulation_topology_5page_wifi}
\vspace{-0.1in}
\end{figure}

Fig. \ref{fig_simulation_topology_5page_wifi} illustrates the environment that comprises a 6-lane, 1000-meter (m) road with a dividing strip of 5 m between the two sides (the upper side for right-to-left and the lower side for left-to-right direction). There is a Wi-Fi system that is composed of an access point (AP) and a station (STA) right above the road, where the AP acts as the TX.

A snapshot of DSRC-Wi-Fi coexistence geometry is illustrated in Fig. \ref{fig_simulation_topology_5page_wifi} in a Cartesian coordinate system. Position of the $i$th vehicle is formally written as $\mathtt{x}_{i}=\left(x_{i},y_{i}\right) \in \mathbf{R}^2_{\boldsymbol{\mathsf{D}}}$. Since the region $\mathbf{R}_{\boldsymbol{\mathsf{D}}}^2$ is composed of 6 lanes, there are 6 values of $y_{i}$. For each lane, X axis ranges $x_{i} = \left[0, 1000\right]$, in which vehicles are distributed according to an independent and homogeneous one-dimensional Poisson Point Process (PPP) of intensity $\lambda_{dsrc}$ vehicles per 1000 m per lane. Note that for ease of analysis, we assume homogeneous traffic across all 6 lanes so that each lane has the same intensity of traffic. Thus, the number of vehicles distributed in $\mathbf{R}_{\boldsymbol{\mathsf{D}}}^2$ that is provided in our results in Section \ref{sec_evaluation} are given by $6\lambda_{dsrc}$. According to the uniformity property of a homogeneous point process \cite{daley}, we can assume that the DSRC vehicles are uniformly randomly scattered in each lane.

The DSRC network shares the 5.9 GHz spectrum with a Wi-Fi system in a co-channel manner. Position of the AP is defined as $\mathtt{x}_{ap}=\left(x_{ap},y_{ap}\right) \in \mathbf{R}^2_{\boldsymbol{\mathsf{W}}}$. Then, a Wi-Fi network is formed around the AP, which is formulated as $\left|\mathtt{x}_{ap}-\mathtt{x}_{sta}\right| \le r_{\boldsymbol{\mathsf{W}},tx}$ where $\mathtt{x}_{sta}=\left(x_{sta},y_{sta}\right)$ is position of a Wi-Fi STA, and $r_{\boldsymbol{\mathsf{W}},tx}$ is the transmission range of a Wi-Fi AP. Distribution of STAs is represented as another PPP in $\mathbf{R}_{\boldsymbol{\mathsf{W}}}^2$ with $\lambda_{sta}\left(>0\right)$.

This paper focuses on evaluating the performance of a DSRC network; the only role of Wi-Fi transmissions is external interference to the DSRC. We assume the Wi-Fi to be in a saturated situation, where the AP transmits every slot. Under this assumption, the value of $\lambda_{sta}$ does not matter as long as the AP continues transmission; thus we assume $\lambda_{sta} = 1$. Also, the Wi-Fi AP adopts an omni-directional antenna; once a position of AP is set, the only factor that determines the level of the interference at a DSRC vehicle is the vehicle's distance from the Wi-Fi AP.

There are four possible results of a packet transmission including \textit{successful delivery (DLVY)}; \textit{expiration (EXP)}; \textit{synchronized transmission (SYNC)}; and \textit{hidden node problem (HN)} \cite{elsevier14}. A DLVY does not undergo either contention nor collision. Also, with a DLVY, we assume that every RX vehicle in a TX vehicle's transmission range successfully receives the packet. On the other hand, an EXP defines a situation where a packet is not even transmitted due to contention with other vehicles or a Wi-Fi TX.

There are two main types of packet \textit{collision}: A SYNC refers to a situation where more than one TXs start transmission at the same time due to the same value of backoff in CSMA/CA. On the other hand, a HN occurs in relation to carrier-sense threshold.

\section{Proposed Performance Evaluation Method}\label{sec_rgb}
The three types of packet failure--EXP, SYNC, and HN--do not necessarily incur a lost message for all the RX vehicles in a TX vehicle's transmission range, $r_{tx}$. Based on that main idea, this section introduces a geometric method that precisely evaluates the performance of broadcast-based DSRC network.

We focus on the geometric analysis that can reflect a wider perspective of a network as opposed to gauging the performance of a certain TX-RX pair. We propose a metric that more accurately displays reception of the broadcast packets and thus represents the performance of a vehicular network more accurately. This metric, \textit{Reception Geometry for Broadcast} ($\boldsymbol{\mathsf{RGB}}$) indicator, is defined at an arbitrary TX vehicle as
\begin{align}\label{eq_rgb}
\boldsymbol{\mathsf{RGB}} &= \frac{\text{Area where a packet can be received}}{\text{Area of the transmission range}}\nonumber\\
&=\frac{\mathsf{A}_{rx}}{\pi r_{tx}^{2}}.
\end{align}
As implied in (\ref{eq_rgb}), an $\boldsymbol{\mathsf{RGB}}$ is a normalized quantity with a range of [0,1]. For discussing how to compute $\boldsymbol{\mathsf{RGB}}$, we start with DLVY and EXP since they are straightforward. As a DLVY does not undergo either contention nor collision, $\boldsymbol{\mathsf{RGB}}$ is given by
\begin{align}\label{eq_rgb_suc}
\mathsf{A}_{rx} = \pi r_{tx}^{2} \Rightarrow \boldsymbol{\mathsf{RGB}}_{dlvy}=1.
\end{align}
On the other hand, $\boldsymbol{\mathsf{RGB}}$ for an EXP is given by
\begin{align}\label{eq_rgb_exp}
\mathsf{A}_{rx} = 0 \Rightarrow \boldsymbol{\mathsf{RGB}}_{exp}=0.
\end{align}
The reason for this is that a TX vehicle cannot even start transmission once it gets backed off by a contention with other TXs in the carrier-sense range. As a result, there is no RX vehicle in the network that can receive the packet.

\begin{figure}
\vspace{-0.3in}
\subfigure[Minimum $\mathsf{A}_{rx}$]
{
\centering
\includegraphics[height = 1.5in]{fig_geometry_sync_min.pdf}
\label{fig_geometry_sync_min}
}
\centering
\subfigure[Maximum $\mathsf{A}_{rx}$]
{
\centering
\includegraphics[height = 1.2in]{fig_geometry_sync_max.pdf}
\label{fig_geometry_sync_max}
}
\caption{Geometry of intersection for SYNC}
\label{fig_geometry_sync}
\vspace{-0.1in}
\end{figure}

Interestingly, an $\boldsymbol{\mathsf{RGB}}$ is evaluated differently between SYNC and HN. Explanation of the rationale refers to the stochastic geometry that was discussed in Section \ref{sec_model}. The geometry is illustrated in Fig. \ref{fig_geometry_sync}. A SYNC occurs when a colliding TX is in the carrier-sense range of a TX but they are assigned the same value of backoff. More specifically, $\mathsf{A}_{rx}$ of an $r_{tx}$ is determined as the ``rest'' of the intersection of two circles of transmission ranges of the TX and colliding TXs as in Fig. \ref{fig_geometry_sync}, which is given by
\begin{align}\label{eq_arxsync}
&\mathsf{A}_{rx}\left(\mathsf{d}\left(\mathtt{x}_{T},\mathtt{x}_{C}\right)\right)\nonumber\\
&= \mathsf{A}_{tx} - \mathsf{A}_{col}\left(\mathsf{d}\left(\mathtt{x}_{T},\mathtt{x}_{C}\right)\right)\nonumber\\
&= \pi r_{tx}^2 - \Biggl( 2r_{tx}^2 \cos^{-1} \left(\frac{\mathsf{d}\left(\mathtt{x}_{T},\mathtt{x}_{C}\right)}{2r_{tx}}\right)\nonumber\\
&{\rm{~~~~~~~~~~~~~~}}- \frac{\mathsf{d}\left(\mathtt{x}_{T},\mathtt{x}_{C}\right)}{2}\sqrt{4r_{tx}^2 - \mathsf{d}^2\left(\mathtt{x}_{T},\mathtt{x}_{C}\right)} \Biggr).
\end{align}
We denote by $\mathsf{A}_{col}$ the collision area and is given by $\mathsf{A}_{col} = \mathsf{A}_{tx} - \mathsf{A}_{rx}$. Also, $\mathtt{x}_{T}$ and $\mathtt{x}_{C}$ are the positions of the TX and colliding TXs, respectively. $r_{tx}$ denotes the transmission range of a vehicle. The distance between the TX vehicle and the colliding TX is defined as a function of positions $\mathtt{x}_{T}$ and $\mathtt{x}_{C}$ as
\begin{align}\label{eq_rgb_sync_inequality}
d_{avg} \le \mathsf{d}\left(\mathtt{x}_{T},\mathtt{x}_{C}\right) \le r_{cs}.
\end{align}
As Fig. \ref{fig_geometry_sync_min} shows, a $\min \mathsf{A}_{rx}$ occurs with $\max \mathsf{d}\left(\mathtt{x}_{T},\mathtt{x}_{C}\right)$; a colliding TX is placed at the border of the TX vehicle's carrier-sense range. In contrast, as in Fig. \ref{fig_geometry_sync_max}, a $\max \mathsf{A}_{rx}$ occurs with $\min \mathsf{d}\left(\mathtt{x}_{T},\mathtt{x}_{C}\right)$. Since infinitely many cases are possible for $\min d$, we consider $d_{avg}$, the average value of inter-vehicle distance, to be the minimum separation between two arbitrary vehicles. This bounds the $\boldsymbol{\mathsf{RGB}}$ for a SYNC in $\mathbf{R}_{\boldsymbol{\mathsf{D}}}^2$ as
\begin{align}\label{eq_rgb_sync}
\frac{\mathsf{A}_{rx}\left(d_{avg}\right)}{\pi r_{tx}^2} \le \boldsymbol{\mathsf{RGB}}_{sync}\left(\mathtt{x}_{T},\mathtt{x}_{C}\right) \le \frac{\mathsf{A}_{rx}\left(r_{cs}\right)}{\pi r_{tx}^2}.
\end{align}

Now, the $\boldsymbol{\mathsf{RGB}}$ for a HN can be formulated in a similar logic to (\ref{eq_arxsync})-(\ref{eq_rgb_sync}). While the way to obtain $\mathsf{A}_{rx}$ refers to (\ref{eq_arxsync}), $\mathsf{d}\left(\mathtt{x}_{T},\mathtt{x}_{C}\right)$ is computed as
\begin{align}\label{eq_rgb_hn_inequality}
r_{cs} < \mathsf{d}\left(\mathtt{x}_{T},\mathtt{x}_{C}\right) < 2r_{cs}.
\end{align}
A $\max \mathsf{d}\left(\mathtt{x}_{T},\mathtt{x}_{C}\right)$ causing $\min \mathsf{A}_{rx}$ occurs when the $r_{tx}$'s of the TX and the colliding TX contact at the border. A $\min \mathsf{d}\left(\mathtt{x}_{T},\mathtt{x}_{C}\right)$ yielding $\max \mathsf{A}_{rx}$ occurs when the colliding TX is right outside of $r_{tx}$ of the TX vehicle. This yields
\begin{align}\label{eq_rgb_hn}
\frac{\mathsf{A}_{rx}\left(r_{cs}\right)}{\pi r_{tx}^2} \le \boldsymbol{\mathsf{RGB}}_{hn}\left(\mathtt{x}_{T},\mathtt{x}_{C}\right) \le \frac{\mathsf{A}_{rx}\left(2r_{cs}\right)}{\pi r_{tx}^2}.
\end{align}

Finally, one can obtain a \textit{mean} $\boldsymbol{\mathsf{RGB}}$ that is averaged over (i) all the packets transmitted by a TX vehicle within a unit time, (ii) all the possible positions of the TX vehicle, and (iii) all the possible positions of the RX vehicles:
\begin{align}\label{eq_rgb_mean}
\overline{\boldsymbol{\mathsf{RGB}}} &= \frac{1}{\left|\mathbf{R}_{\boldsymbol{\mathsf{D}}}^2\right|^2} \int_{\mathtt{x}_{T} \in \mathbf{R}_{\boldsymbol{\mathsf{D}}}^2} \int_{\mathtt{x}_{C} \in \mathbf{R}_{\boldsymbol{\mathsf{D}}}^2} \nonumber\\
&{\rm{~~~~~~}}\sum_{k=1}^{N_{pkt}}\sum_{j \in \mathcal{S}_{pb}}{\boldsymbol{\mathsf{RGB}}_{k,j}\left(\mathtt{x}_{T},\mathtt{x}_{C}\right)}{\mathsf{P}_{k,j}} d\mathtt{x}_{T}d\mathtt{x}_{C}.
\end{align}
Denote the number of packets that are generated in a unit time by $N_{pkt}$. Alos, $\mathcal{S}_{pb}$ is the set of the four possible results of a packet transmission, $\mathcal{S}_{pb} = \left\{\text{DLVY, EXP, SYNC, HN}\right\}$. $\boldsymbol{\mathsf{RGB}}_{k,j}$ denotes each of (\ref{eq_rgb_suc}), (\ref{eq_rgb_exp}), (\ref{eq_rgb_sync}), and (\ref{eq_rgb_hn}) for the $k$th packet, and $\mathsf{P}_{k,j}$ means the probability that each of $\mathcal{S}_{pb}$ occurs.

\begin{table}[t]
\vspace{-0.3in}
\caption{Parameters}
\centering
\scriptsize
\begin{tabular}{| l | l |}
\hline
\multicolumn{1}{|c|}{\textbf{Parameter}} & \multicolumn{1}{|c|}{\textbf{Value}} \\ \hline \hline
\multicolumn{2}{|c|}{\cellcolor{gray!20}Road}\\ \hline
Length of the road & 1000 m \\ \hline
Number of lanes & 6 \\ \hline
Lane width & 4 m \\ \hline
Dividing strip width & 5 m \\ \hline
\multicolumn{2}{|c|}{\cellcolor{gray!20}DSRC}\\ \hline
Number of vehicles in $\mathbf{R}_{\boldsymbol{\mathsf{D}}}^2$ (= $6\lambda_{dsrc}$) & \{30,120,210,270\} \\ \hline
Inter-vehicle distance & \{200, 50, 28.57, 22.22\} m \\ \hline
Velocity & \{42.07, 17.49, 11.93, 9.72\} m/sec \\ \hline
Receiver sensitivity & -91 dBm\\ \hline
Transmission antenna gain & 0 dBm \\ \hline
Packet length & 500 bytes \\ \hline
Data rate & 6 Mbps \\ \hline
Slot time & 13 $\mu$s \\ \hline
\multicolumn{2}{|c|}{\cellcolor{gray!20}Wi-Fi}\\ \hline
Receiver sensitivity & -82 dBm\\ \hline
Transmission power & 30 dBm \\ \hline
CW Size & [15,1023] \\ \hline
Packet length & 10,800 bytes \\ \hline
Data rate & 54 Mbps \\ \hline
Slot time & 9 $\mu$s \\ \hline
\end{tabular}
\label{table_simulation_parameters}
\vspace{-0.1in}
\end{table}

\begin{figure*}[t]
\vspace{-0.3in}
\centering
\subfigure[With CW = 63]
{
\centering
\includegraphics[width=0.44\linewidth]{fig_proportion_cw63}
\label{fig_proportion_63}
}
\hspace{0.1in}
\subfigure[With CW = 255]
{
\centering
\includegraphics[width=0.44\linewidth]{fig_proportion_cw255}
\label{fig_proportion_255}
}
\caption{$\mathsf{P}_{j,k}$ according to vehicle density}
\label{fig_proportion}
\vspace{-0.1in}
\end{figure*}

\section{Performance Evaluation}\label{sec_evaluation}
We evaluate the performance of a DSRC-based vehicular network through the \textit{network simulator 2 (ns-2)}. The simulation topology refers to Fig. \ref{fig_simulation_topology_5page_wifi}. Each result is an average from 180 runs of 20-second experiments on \textit{ns-2}, based on the settings described in Table \ref{table_simulation_parameters}. The 180 runs are composed of 30 runs for each of the 6 lanes, in order to remove bias that could occur if not all lanes were considered.

\subsection{Setting}
As mentioned in Section \ref{sec_model}, the DSRC vehicles are uniformly distributed on the road. Their mobility setting is summarized in Table \ref{table_simulation_parameters}. The total number of vehicles situated in region $\mathbf{R}_{\boldsymbol{\mathsf{D}}}^2$ is the number of lanes multiplied by the number of vehicles per lane, i.e., $6\lambda_{dsrc}$. Inter-vehicle distance is obtained by $\lambda_{dsrc} / \text{(Road length)} = \lambda_{dsrc} / 1000$. The inter-vehicle distance accordingly determines the vehicle velocity. A set of empirical data \cite{distance_velocity} matches the velocity with the minimum inter-vehicle distance required for avoiding a collision. Since the data set was a discrete matrix that omitted some values that we needed, we fitted the data using MATLAB based on a fitting technique called \textit{shape-preserving fitter}. The resulting values of velocities were used in the \textit{ns-2} experiments.

To display the three types of packet failure--EXP, SYNC, and HN, we modified the C++ stack of the \textit{ns-2}. More specifically, we modified the packet trace function inside the stack in such a way that a colliding packet is declared to be HN/SYNC if the TX of the colliding packet is outside/within the TX range of a vehicle, respectively.

\begin{figure*}[t]
\vspace{-0.3in}
\centering
\subfigure[Without interference from Wi-Fi]
{
\centering
\includegraphics[width=0.44\linewidth]{fig_rgb_ibi_wo80211ac.eps}
}
\hspace{0.1in}
\subfigure[With interference from Wi-Fi]
{
\centering
\includegraphics[width=0.44\linewidth]{fig_rgb_ibi_w80211ac.eps}
}
\caption{$\boldsymbol{\mathsf{RGB}}$ according to IBI}
\label{fig_rgb_ibi}
\end{figure*}

\begin{figure*}[t]
\centering
\subfigure[Without interference from Wi-Fi]
{
\centering
\includegraphics[width=0.44\linewidth]{fig_rgb_cw_wo80211ac.eps}
\label{fig_rgb_cw_wo80211ac}
}
\hspace{0.1in}
\subfigure[With interference from Wi-Fi]
{
\centering
\includegraphics[width=0.44\linewidth]{fig_rgb_cw_w80211ac.eps}
\label{fig_rgb_cw_w80211ac}
}
\caption{$\boldsymbol{\mathsf{RGB}}$ according to CW}
\label{fig_rgb_cw}
\vspace{-0.1in}
\end{figure*}

\subsection{Results}
The results provide the values of $\mathsf{P}_{j,k}$ and $\boldsymbol{\mathsf{RGB}}_{j,k}$ that are given in (\ref{eq_rgb_mean}).

\subsubsection{$\mathsf{P}_{j,k}$ Analysis}
Fig. \ref{fig_proportion} demonstrates $\mathsf{P}_{j,k}$ for DLVY, EXP, SYNC, and HN, according to CW and presence of Wi-Fi TX. There are two bars at each value of the X-axis. The left and right bars indicate the results with and without interference from Wi-Fi. Note that comparison of Fig. \ref{fig_proportion} to Figs. 3 through 5 in \cite{elsevier14} validates our results.

Note that two representative values of CW are selected: 63 and 255. First, CW = 63 is chosen since it is a proper value to highlight impact of presence of Wi-Fi interference in a balanced manner among SYNC and HN. We found that smaller values of CW did not highlight the impact on HN due to too many SYNCs. Also, selection of CW = 255 is to compare the impact of a larger CW. In fact, Fig. \ref{fig_proportion_255} shows dramatic increase in EXP and decrease in SYNC.

Two important observations need to be made. First, EXP increases with a larger value of CW, higher traffic density, and presence of Wi-Fi TX. Second, SYNC increases with higher traffic density and presence of Wi-Fi TX, but decreases as CW increases. The reason is that now there is a wider choice of CW values for DSRC TXs and as a result $\mathsf{P}_{sync}$ gets smaller. Third, HN increases as traffic density increases and as a Wi-Fi TX is introduced. It implies that the Wi-Fi TX acts as a hidden node.

\subsubsection{$\boldsymbol{\mathsf{RGB}}$ Analysis}
Fig. \ref{fig_rgb_ibi} illustrates $\overline{\boldsymbol{\mathsf{RGB}}}$ that is given in (\ref{eq_rgb_mean}), according to (i) traffic density, (ii) IBI, and (iii) presence of interfering Wi-Fi TX. There are two noteworthy interpretations. First, the $\boldsymbol{\mathsf{RGB}}$ decreases as IBI increases. This implies that traffic congestion is a predominating factor that determines the performance of a DSRC network. When each vehicle increases IBI, the total amount of packets over the air decreases. Second, with interference from Wi-Fi, the $\boldsymbol{\mathsf{RGB}}$ is decreased for all traffic density cases.

Fig. \ref{fig_rgb_cw} shows $\overline{\boldsymbol{\mathsf{RGB}}}$ as a function of (i) traffic density, (ii) CW, and (ii) presence of Wi-Fi interference. Note the following two observations. First, a smaller value of CW is recommended as traffic density increases. This is osbserved for higher vehicle densities--210 and 270 vehicles in $\mathbf{R}_{\boldsymbol{\mathsf{D}}}^2$. Comparison of boxes (A) and (B) in Figs. \ref{fig_rgb_cw_wo80211ac} and \ref{fig_rgb_cw_w80211ac} highlights the tendency. This can be explained follows: the $\boldsymbol{\mathsf{RGB}}$ is most dependent on $\mathsf{P}_{dlvy}$, since the $\boldsymbol{\mathsf{RGB}}$ value for a DLVY is the largest (Recall that $\boldsymbol{\mathsf{RGB}}_{dlvy} = 1 > \boldsymbol{\mathsf{RGB}}_{sync}, \boldsymbol{\mathsf{RGB}}_{hn}, \boldsymbol{\mathsf{RGB}}_{exp}$). Because $\mathsf{P}_{dlvy}$ depends on both expired and collided messages, a balance needs to be found between the two quantities. After a certain point, the number of expired packets begins to prevail and the reception rate starts a slower, but steady, decrease. When $\mathsf{P}_{exp}$ goes beyond its peak (and the number of collisions starts increasing again), $\mathsf{P}_{dlvy}$ goes through another phase, where its value remains almost constant. Finally, when $\mathsf{P}_{exp}$ becomes low enough, the effect of collisions becomes predominant again and the number of total receptions starts decreasing once more. Second, Wi-Fi interference can be interpreted as additional traffic. This can be seen from the fact that the peak moves to lower CW values between Figs. \ref{fig_rgb_cw_wo80211ac} and \ref{fig_rgb_cw_w80211ac}.

\section{Conclusion}\label{sec_conclusion}
This paper proposes a new metric that accurately measures the performance of a DSRC-based vehicular network that broadcasts BSMs for safety-critical applications. Compared to traditional metrics, it more finely evaluates reception of a broadcast packet that is collided, by taking into account the RX vehicles that still are able to receive the packet. Based on the metric, it analyzes impacts of the external interference from IEEE 802.11ac-based Wi-Fi. The results suggest on the coexistence that (i) higher IBI and (ii) smaller CW are recommended in cases of higher vehicle density.


text:
\section{Introduction}

Time series forecasting is a growing field of interest playing an important role in nearly all fields of science and engineering, such as economics, finance, meteorology and telecommunication~\citep{compTimeSeries2005}. Unlike one-step ahead forecasting, multi-step ahead forecasting tasks are more difficult~\citep{tiaotsay}, since they have to deal with various additional complications, like accumulation of errors, reduced accuracy, and increased uncertainty~\citep{santafe,longterm}.

The forecasting domain has been influenced, for a long time, by linear statistical methods such as ARIMA models. However, in the late 1970s and early 1980s, it became increasingly clear that linear models are not adapted to many real applications~\citep{GH06}. In the same period, several useful nonlinear time series models were proposed such as the bilinear model~\citep{PT86}, the threshold autoregressive model~\citep{TongRS1980,Tong83,Tong00} and the autoregressive conditional heteroscedastic~(ARCH) model~\citep{Engle1982}~(see \citep{GH06}  and \citep{DeGooijer1992135} for a review). Nowadays, Monte Carlo simulation or bootstrapping methods are used to compute nonlinear forecasts. Since no assumptions are made about the distribution of the error process, the latter approach is preferred~\citep{CFS04,GH06}. However, the study of nonlinear time series analysis and forecasting is still in its infancy compared to the development of linear time series~\citep{GH06}.

In the last two decades, machine learning models have drawn attention and have established themselves as serious contenders to classical statistical models in the forecasting community~\citep{empiricalMLTimeSeries,compTimeSeries2005,ZPH98}. These models, also called black-box or data-driven models~\citep{mlmitchel}, are examples of nonparametric nonlinear models which use only historical data to learn the stochastic dependency between the past and the future. For instance, Werbos found that Artificial Neural Networks~(ANNs) outperforms the classical statistical methods such as linear regression and Box-Jenkins approaches~\citep{Werbos74,Werbos1988}. A similar study has been conducted by Lapedes and Farber~\citep{Lapedes87} who conclude that ANNs can be successfully used for modeling and forecasting nonlinear time series. Later, others models appeared such as decision trees, support vector machines and nearest neighbor regression~\citep{hastie_09_elements-of.statistical-learning,alpaydin-04-ml}. Moreover, the empirical accuracy of several machine learning models has been explored in a number of forecasting competitions under different data conditions~(e.g. the NN3, NN5, and the annual ESTSP competitions~\citep{nn3,nn5,estsp07,estsp08}) creating interesting scientific debates in the area of data mining and forecasting~\citep{debateMiningForecasting,debateMiningForecastingCommentPrice,debateMiningForecastingCommentCrone}.

In the forecasting community, researchers have paid attention to several aspects of the forecasting procedure such as model selection~\citep{Aha97,msejor06,msnn99,Chapelle00modelselection}, effect of deseasonalization~\citep{modellingSeas,fma,nelson99,Zhang2005501}, forecasts combination~\citep{Bates69thecombination,Clemen89,forecastcomb} and many other critical topics~\citep{GH06}. However, approaches for generating multi-step ahead forecasts for machine learning models did not receive as much attention, as pointed out by Kline: ``\emph{One issue that has had limited investigation is how to generate multiple-step-ahead forecasts}''~\citep{klineMethods}.

To the best of our knowledge, five alternatives~(or strategies) have been proposed in the literature to tackle an $H$-step ahead forecasting task. The \emph{Recursive} strategy~\citep{santafe, longterm,multistep,tiaotsay,klineMethods,Hamzacebi20093839} iterates, $H$ times, a one-step ahead forecasting model to obtain the $H$ forecasts. After estimating the future series value, it is fed back as an input for the following forecast.

In contrast to the previous strategy which use a single model, the \emph{Direct} strategy~\citep{santafe, longterm,multistep,tiaotsay,klineMethods,Hamzacebi20093839} estimates a set of $H$ forecasting models, each returning a forecast for the $i$-th value~($i \in \{1,\dots,H\} $).

A combination of the two previous strategies, called \emph{DirRec} strategy has been proposed in~\citep{sorjamaa_esann_06}. The idea behind this strategy is to combine aspects from both, the Direct and the Recursive strategies. In other words, a different model is used at each step but the approximations from previous steps are introduced into the input set.

In order to preserve, between the predicted values, the stochastic dependency characterizing the time series, the Multi-Input Multi-Output~(\emph{MIMO}) strategy has been introduced and analyzed in ~\citep{mimo,mimoijf,klineMethods}. Unlike the previous strategies where the models return a scalar value, the MIMO strategy returns a vector of future values in a single step.

The last strategy, called \emph{DIRMO}~\citep{mismo}, aims to preserve the most appealing aspects of both the DIRect and miMO strategies. This strategy aims to find a trade-off between the property of preserving the stochastic dependency between the forecasted values and the flexibility of the modeling procedure.

In the literature, these five forecasting strategies have been presented separately, sometimes, using different terminologies. The \emph{first contribution} of this paper is to present a thorough unified review as well as a theoretical comparative analysis of the existing strategies for multi-step ahead forecasting.

Despite the fact that many studies have compared between the different multi-step ahead approaches, the collective outcome of these studies regarding forecasting performance has been inconclusive.  So the modeler is still left with little guidance as to which strategy to use. For example, research from \citep{Bontempi99BirattariBersini-ICML99,Weigend92HubermanRumelhart} provide experimental evidence in favor of Recursive strategy against Direct strategy. However, results from \citep{Zhang94Hutchinson,longterm,Hamzacebi20093839} support the fact that the Direct strategy is better than the Recursive strategy.  The work by \citep{sorjamaa_esann_06} shows that the DirRec strategy gives better performance than Direct and Recursive strategies. The Direct and Recursive strategies have been theoretically and empirically compared in \citep{Atiya99acomparison}. In this study the authors obtained theoretical and experimental evidence in favor of Direct strategy. Concerning the MIMO strategy, Kline~\citep{klineMethods} and Cheng et al~\citep{multistep} support the idea that MIMO strategy provides worse forecasting performance than Recursive and Direct strategies. However, in \citep{mimo,mimoijf}, the comparison between MIMO, Recursive, and Direct was in favor of MIMO. Finally, ~\citep{mismo,BenTaieb20101950} show that the DIRMO strategy gives better forecasting results than Direct and MIMO strategies when the parameter controlling the degree of dependency between forecasts is correctly identified.
These previous comparisons have been performed with different datasets in different configurations using different forecasting methods, such as Multiple Linear Regression, Artificial Neural Networks, Hidden Markov Models and Nearest Neighbors.

All the contradictory findings of these studies make it all the more necessary to investigate further to find the truth concerning the relative performance of these strategies.
The \emph{second contribution} of this paper is an experimental comparison of the different multi-step ahead forecasting strategies on the $111$ time series of the NN5 international forecasting competition benchmark. These time series pose some of the realistic problems that one usually encounters in a typical multi-step ahead forecasting task, for example the existence of several times series of possibly related dynamics, outliers, missing values, and multiple overlying seasonalities. This experimental comparison is performed for a variety of different configurations (regarding seasonality, input selection and combination), in order to have the comparison as encompassing as can be. In addition, the methodology used for this experimental comparison is based on the guidelines and recommendations advocated in some of the methodological papers \citep{Demsar_statisticalcomparisons,extensionDemsar}.

In other words, the aim of this paper is not to make a comparison of machine learning algorithms for forecasting~(which was already conducted in \citep{empiricalMLTimeSeries}) but rather to show for a given learning algorithm, how the choice of the forecasting strategy can sensibly influence the performance of the multi-step ahead forecasts. In this work, we adopted the Lazy Learning algorithm~\citep{Aha97}, a particular instance of local learning, which has been successfully applied to many real-world forecasting tasks~\citep{Sauer94,Bontempi98BirattariBersini_Leuven,McNames98}.

Last but not least, the paper proposes also a Lazy Learning entry to the NN5 forecasting competition~\citep{nn5}. The goal is to assess how this model fares compared to the other computational intelligence models that were proposed for the competition~\citep{mimoijf}. This will give us an idea about the potential of this approach.

The paper is organized as follows. The next section presents a review of the different forecasting strategies. Section \ref{sectionLazy} describes the Lazy Learning model and the associated algorithms for the different forecasting strategies. Section \ref{sectionExperiments} gives a detailed presentation of the datasets and the methodology applied for the experimental comparison. Section \ref{sectionResDiscussion} presents the results and discusses them. Finally, Section \ref{sectionConclusion} gives a summary and concludes the work.

\section{Strategies for Multi-Step-Ahead Time Series Forecasting}

\label{sectionStrategies}

A multi-step ahead~(also called long-term) time series forecasting task consists of predicting the next $H$ values $[y_{N+1},\dots,y_{N+H}]$ of a historical time series $[y_1,\dots,y_N]$ composed of $N$ observations, where $H>1$ denotes the forecasting horizon.

This section will first give a presentation of the five forecasting strategies and next, a subsection will be devoted to a comparative analysis of these strategies in terms of number and types of models to learn as well as forecasting properties.

We will use a common notation where $f$ and $F$ denote the functional dependency between past and future observations, $d$ refers to the embedding dimension~\citep{Casdagli91EubankFarmerGibson} of the time series, that is the number of past values used to predict future values and $w$ represents the term that includes modeling error, disturbances and/or noise.

\subsection{Recursive strategy}
\label{RecSection}

The oldest and most intuitive forecasting strategy is the \emph{Recursive}~(also called \emph{Iterated} or \emph{Multi-Stage}) strategy~\citep{santafe, longterm,multistep,tiaotsay,klineMethods,Hamzacebi20093839}. In this strategy, a single model $f$ is trained to perform a {\it one-step ahead} forecast, i.e.

\begin{equation} y_{t+1}=f(y_t,\dots,y_{t-d+1})+w \textnormal{,} \label{stochRec} \\ \end{equation}

\noindent with $t \in \{d,\dots,N-1\}$.

When forecasting $H$ steps ahead, we first forecast the first step
by applying the model. Subsequently, we use the value just forecasted
as part of the input variables for forecasting the next step (using the
same one-step ahead model).
We continue in this manner until we have forecasted
the entire horizon.

Let the trained one-step ahead model be $\hat f$. Then the forecasts
are given by:

\begin{equation}
\label{equationRec}
\hat y_{N+h} =
\begin{cases}
\hat{f}(y_{N},\dots,y_{N-d+1}) & \text{if } h=1 \\
\hat{f}(\hat y_{N+h-1}, \dots,\hat y_{N+1},y_{N},\dots,y_{N-d+h}) & \text{if } h \in \{2,\dots,d\} \\
\hat{f}(\hat y_{N+h-1}, \dots,\hat y_{N+h-d}) & \text{if } h \in \{d+1,\dots,H\} \\
\end{cases}
\end{equation}

Depending on the noise present in the time series and the forecasting horizon, the recursive strategy may suffer from low performance in multi-step ahead forecasting tasks. Indeed, this is especially true if the forecasting horizon $h$ exceeds the embedding dimension $d$, as at some point all the inputs are forecasted values instead of actual observations~(Equation \ref{equationRec}). The reason for the potential inaccuracy is that the Recursive strategy is sensitive to the accumulation of errors with the forecasting horizon.
Errors present in intermediate forecasts will propagate forward as these forecasts
are used to determine subsequent forecasts.

In spite of these limitations, the Recursive strategy has been successfully used to forecast many real-world time series by using different machine learning models, like recurrent neural networks~\citep{rnnTimeSeries} and nearest-neighbors~\citep{McNames98,Bontempi99BirattariBersini-ICML99}.

\subsection{Direct strategy}

The \emph{Direct}~(also called Independent) strategy~\citep{santafe, longterm,multistep,tiaotsay,klineMethods,Hamzacebi20093839} consists of forecasting each horizon independently from the others. In other terms, $H$ models $f_h$ are learned~(one for each horizon) from the time series $[y_1,\dots,y_N]$ where

\begin{equation} y_{t+h}=f_h(y_t,\dots,y_{t-d+1})+w \label{stochDirect} \textnormal{,} \\ \end{equation}

with $t \in \{d,\dots,N-H\}$ and $h \in \{1,\dots,H\}$.

The forecasts are obtained by using the $H$ learned models $\hat{f}_h$ as follows:

\begin{equation}
\label{equationDir}
\hat{y}_{N+h}=\hat{f}_h(y_{N},\dots,y_{N-d+1}) \textnormal{.}
\end{equation}

This implies that the Direct strategy does not use any approximated values to compute the forecasts~(Equation \ref{equationDir}), being then immune to the accumulation of errors. However, the $H$ models are learned independently inducing a conditional independence of the $H$ forecasts. This affects the forecasting accuracy as it prevents the strategy from considering complex dependencies between the variables $\hat{y}_{N+h}$\citep{mimo,mimoijf,klineMethods}. For example consider a case where the best forecast is a linear or mildly nonlinear trend. The direct method could yield a broken curve because of the ``uncooperative" way the $H$ forecasts are generated. Also, this strategy demands a large computational time since there are as many models to learn as the size of the horizon.

Different machine learning models have been used to implement the Direct strategy for multi-step ahead forecasting tasks, for instance neural networks~\citep{klineMethods}, nearest neighbors~\citep{longterm} and decision trees~\citep{directTrees}.

\subsection{DirRec strategy}

The \emph{DirRec} strategy~\citep{sorjamaa_esann_06} combines the architectures and the principles underlying the Direct and the Recursive strategies. DirRec computes the forecasts with different models for every horizon~(like the Direct strategy) and, at each time step, it enlarges the set of inputs by adding variables corresponding to the forecasts of the previous step~(like the Recursive strategy). However, note that unlike the two previous strategies, the embedding size $d$ is not the same for all the horizons. In other terms, the DirRec strategy learns $H$ models $f_h$ from the time series $[y_1,\dots,y_N]$ where

\begin{equation} y_{t+h}=f_h(y_{t+h-1},\dots,y_{t-d+1})+w \textnormal{,} \\  \label{stochDirRec} \end{equation}

\noindent with $t \in \{d,\dots,N-H\}$ and $h \in \{1,\dots,H\}$.

To obtain the forecasts, the $H$ learned models are used as follows:

\begin{equation}
\label{equationDirRec}
\hat y_{N+h} =
\begin{cases}
\hat{f_h}(y_{N},\dots,y_{N-d+1}) & \text{if } h=1 \\
\hat{f_h}(\hat y_{N+h-1}, \dots,\hat y_{N+1},y_{N},\dots,y_{N-d+1}) & \text{if } h \in \{2,\dots,H\} \\
\end{cases}
\end{equation}

This strategy outperformed the Direct and the Recursive strategies on two real-world time series: Santa Fe and Poland Electricity Load data sets~\citep{sorjamaa_esann_06}. Few research has been done regarding this strategy, so there is a need for further evaluation.
\subsection{MIMO strategy}

The three previous strategies~(Recursive, Direct and DirRec) may be considered as Single-Output strategies~\citep{BenTaieb20101950} since they model the data as a (multiple-input) single-output function~(see Equations \ref{equationRec}, \ref{equationDir} and \ref{equationDirRec}).

The introduction of the \emph{Multi-Input Multi-Output}~(MIMO) strategy~\citep{mimo,mimoijf}~(also called Joint strategy~\citep{klineMethods}) has been motivated by the need to avoid the modeling of single-output mapping, which neglects the existence of stochastic dependencies between future values and consequently affects the forecast accuracy~\citep{mimo,mimoijf}.

The MIMO strategy learns one multiple-output model $F$ from the time series $[y_1,\dots,y_N]$ where

\begin{equation}  [y_{t+H},\dots,y_{t+1}]=F(y_{t},\dots,y_{t-d+1})+ \mathbf{w} \textnormal{,} \\  \label{stochMIMO} \end{equation}

\noindent with $t \in \{d,\dots,N-H\}$, $F: \mathbb{R}^d \rightarrow \mathbb{R}^H$ is a vector-valued function~\citep{vector-valued}, and $\mathbf{w} \in \mathbb{R}^H$ is a noise vector with a covariance that is not necessarily diagonal~\citep{matias} .

The forecasts are returned in one step by a multiple-output model $\hat F$ where

\begin{equation} \label{equationMIMO} [\hat{y}_{t+H},\dots,\hat{y}_{t+1}]=\hat{F}(y_{N},\dots,y_{N-d+1}) \textnormal{.} \\ \end{equation}

The rationale of the MIMO strategy is to preserve, between the predicted values, the stochastic dependency characterizing the time series. This strategy avoids the conditional independence assumption made by the Direct strategy as well as the accumulation of errors  from which plagues the Recursive strategy. So far, this strategy has been successfully applied to several real-world multi-step ahead time series forecasting tasks~\citep{mimo,mimoijf,mismo,BenTaieb20101950}.

However, the need to preserve the stochastic dependencies by using one model has a drawback as it constrains all the horizons to be forecasted with the same model structure. This constraint could reduce the flexibility of the forecasting approach~\citep{mismo}. This was the motivation for the introduction of a new multiple-output strategy: DIRMO~\citep{mismo,BenTaieb20101950}, presented next.

\subsection{DIRMO strategy}

The DIRMO strategy~\citep{mismo} aims to preserve the most appealing aspects of both the DIRect and miMO strategies. Taking a middle approach, DIRMO forecasts the horizon $H$ in blocks, where each block is forecasted in a MIMO fashion. Thus, the $H$-step-ahead forecasting task is decomposed into $n$ multiple-output forecasting tasks~($n=\frac{H}{s}$), each with an output of size $s$~($s \in \{1,\dots,H\}$).

When the value of the parameter $s$ is $1$, the number of forecasting tasks $n$ is equal to $H$ which correspond to the Direct strategy. When the value of the parameter $s$ is $H$, the number of forecasting tasks $n$ is equal to $1$ which correspond to the MIMO strategy. There are intermediate configurations between these two extremes depending on the value of a parameter $s$.

The tuning of the parameter $s$ allows us to improve the flexibility of the MIMO strategy by calibrating the dimensionality of the outputs~(no dependency in the case $s=1$ and maximal dependency for $s=H$). This provides a beneficial trade off between the preserving a larger degree of the stochastic dependency between future values and having a greater flexibility of the predictor.

The DIRMO strategy, previously called MISMO strategy~\citep{mismo}~(renamed for clarity reason), learns $n$ models $F_p$ from the time series $[y_1,\dots,y_N]$ where

\begin{equation}  [y_{t+p*s},\dots,y_{t+(p-1)*s+1}]=F_p(y_{t},\dots,y_{t-d+1})+ \mathbf{w} \textnormal{,} \\ \label{stochDIRMO} \end{equation}

\noindent with $t \in \{d,\dots,N-H\}$, $p \in \{1,\dots,n\}$ and $F_p: \mathbb{R}^d \rightarrow \mathbb{R}^s$ is a vector-valued function if $s>1$.

The $H$ forecasts are returned by the $n$ learned models as follows:

\begin{equation} \label{equationDIRMO}  [\hat{y}_{N+p*s},\dots,\hat{y}_{N+(p-1)*s+1}]=\hat{F}_p(y_{N},\dots,y_{N-d+1}) \textnormal{.} \\ \end{equation}

The DIRMO strategy has been successfully applied to two forecasting competitions: ESTSP'07~\citep{mismo} and NN3~\citep{BenTaieb20101950}.

\subsection{Comparative Analysis}

To summarize, there are five possible forecasting strategies that perform a multi-step ahead forecasting task: \emph{Recursive}, \emph{Direct}, \emph{DirRec}, \emph{MIMO} and \emph{DIRMO} strategies. Figure \ref{strategiesGraph} shows the different forecasting strategies with links indicating their relationships.

\begin{figure}[h]
\centering
\includegraphics[scale=0.8]{strategiesGraph.pdf}
\caption{\label{strategiesGraph} The different forecasting strategies with the links showing their relationship.}
\end{figure}

As we see, the DirRec strategy is a combination of the Direct and the Recursive strategy, while the DIRMO strategy is a combination of the Direct and the MIMO strategy.

Contingent on the selected strategy, a different number and type of models will be required. Before presenting the general comparison of the multi-step ahead forecasting strategies, let us highlight using an example the differences between the forecasting strategies.

Consider a multi-step ahead forecasting task for the time series $[y_1,\dots,y_N]$ where the forecasting horizon $H$ is $4$. Table \ref{exampleStrategies} shows, for each strategy, the different input sets and forecasting models involved in the calculation of the four forecasts $[\hat{y}_{N+1},\dots,\hat{y}_{N+4}]$.

\begin{table}[h]
\begin{tabular}{|l|l|l|l|l|}
\hline
&   $   \hat{y}_{N+1}$ & $\hat{y}_{N+2}$ & $\hat{y}_{N+3}$ & $\hat{y}_{N+4}   $ \\
  \hline
  \hline
Recursive &      $ {\scriptstyle   \hat{f}(y_{N},\dots,y_{N-d+1})   } $  & ${\scriptstyle    \hat{f}(\hat{y}_{N+1},y_{N},\dots,y_{N-d+2})   }$   & ${\scriptstyle   \hat{f}(\hat{y}_{N+2},\hat{y}_{N+1},\dots,y_{N-d+3})   }$ & ${\scriptstyle   \hat{f}(\hat{y}_{N+3},\hat{y}_{N+2},\dots,y_{N-d+4})  }$  \\
    \hline
Direct &      ${\scriptstyle  \hat{f}_1(y_{N},\dots,y_{N-d+1})  }$ & ${\scriptstyle   \hat{f}_2(y_{N},\dots,y_{N-d+1})  }$ & ${\scriptstyle  \hat{f}_3(y_{N},\dots,y_{N-d+1})  }$ & ${\scriptstyle  \hat{f}_4(y_{N},\dots,y_{N-d+1}) }$ \\
    \hline
DirRec &  $ {\scriptstyle \hat{f_1}(y_{N},\dots,y_{ N-d+1}) }$ & ${\scriptstyle \hat{f}_2(\hat{y}_{N+1},y_{N},\dots,y_{N-d+1})  }$ & ${\scriptstyle  \hat{f}_3(\hat{y}_{N+2},\hat{y}_{N+1},\dots,y_{N-d+1})  }$ & ${\scriptstyle  \hat{f}_4(\hat{y}_{N+3},\hat{y}_{N+2},\dots,y_{N-d+1})  }$  \\
  \hline
MIMO &  \multicolumn{4}{|c|}{${\scriptstyle  \hat{F}(y_{N},\dots,y_{N-d+1})  }$} \\
  \hline
DIRMO~($s=2$) &  \multicolumn{2}{|c|}{${\scriptstyle  \hat{F}_1(y_{N},\dots,y_{N-d+1})  }$}  & \multicolumn{2}{|c|}{${\scriptstyle  \hat{F}_2(y_{N},\dots,y_{N-d+1})  }$}  \\
  \hline
\end{tabular}
\caption{ \label{exampleStrategies} The different forecasting models used by each strategy to obtain the 4 forecasts needed.}
\end{table}

Let $T_{SO}$ and $T_{MO}$ denote the amount of computational time needed to learn~(with a given learning algorithm) a Single-Output model and a Multiple-Output model, respectively. For a given $H$-step ahead forecasting task, we can see in Table \ref{modelsToLearn} for each strategy the number and type of models to learn, the size of the output for each model as well as the computational time.

\begin{table}[h!]
	\begin{center}
		\begin{tabular}{|c|c|c|c|c|}
			\hline
			 & Number of Models & Types of models & Size of output & Computational time\\
			 \hline
			Recursive & 1 & SO & 1 & $1 \times T_{SO}$  \\
			\hline
			Direct & $H$ & SO & 1 &  $ H \times T_{SO}$ \\
			\hline
			DirRec & $H$ & SO & 1 & $ H \times  (T_{SO}+\mu) $\\
			\hline
			MIMO & 1 & MO & $H$ & $1 \times T_{MO}$ \\
			\hline
			DIRMO & $\frac{H}{s}$ & MO & $s$ & $\frac{H}{s} \times T_{MO}$\\
			\hline
		\end{tabular}
		\caption{For each forecasting strategy: the number and type of models~(Single-Output or Multiple-Output) to learn  and the size of the output for each model.}
		\label{modelsToLearn}
	\end{center}
\end{table}

Suppose $T_{MO} = T_{SO}  + \delta $, which is a reasonable assumption because learning a model with a vector-valued output takes more time than learning a model with a single-valued output. This allows us to rank the forecasting strategies according to their computation time for training given in Table \ref{modelsToLearn}. Indeed, we have

 \begin{equation} \underbrace{1 \times T_{SO} }_{Recursive} < \underbrace{1 \times T_{MO}}_{MIMO} < \underbrace{\frac{H}{s} \times T_{MO} }_{DIRMO} < \underbrace{H \times T_{SO}}_{Direct} < \underbrace{H \times (T_{SO}+\mu) }_{DirRec} \text{ ,} \label{CtimeStrategies}\end{equation}

 \noindent where we suppose that the parameter $s$ of DIRMO is not equal to $1$ or $H$.

Note, in one hand, that the time needed to learn a SO model of the DirRec strategy equals $T_{SO}+\mu$ because the input size of each SO task increases at each step. On the other hand, if we need to select the value of the parameter $s$ by on some tuning, the DIRMO strategy will take more time and hence will be the slowest one.

In the following, we conclude this section by summarizing the pros and cons of the five forecasting strategies as depicted on Table \ref{summaryStrategies}.

\begin{landscape}
\begin{table}
	\begin{center}
		\begin{tabular}{|M{2cm}|M{4cm}|M{4cm}|c|}
			\hline
			 & Pros & Cons & Computational time needed  \tabularnewline
			 \hline
			Recursive & Suitable for noise-free time series (e.g. chaotic)  & Accumulation of errors & $+$  \tabularnewline
			\hline
			Direct & No accumulation of errors & Conditional independence assumption & $++++$  \tabularnewline
			\hline
			DirRec & Trade-off between Direct and Recursive & Input set grows linearly with $H$ & $+++++$ \tabularnewline
			\hline
			MIMO & No conditional independence assumption & Reduced flexibility: same model structure for all the horizons & $++$  \tabularnewline
			\hline
			DIRMO & Trade-off between total dependence and total independence of forecasts & One additional parameter to estimate & $+++$ \tabularnewline
			\hline
		\end{tabular}
		\caption{A Summary of the Pros and Cons of the Different Multi-step Forecasting Strategies}
		\label{summaryStrategies}
	\end{center}
\end{table}
\end{landscape}

\clearpage

\section{Lazy Learning for Time Series Forecasting}
\label{sectionLazy}

Each of the forecasting strategies introduced in Section \ref{sectionStrategies} demands the definition of a specific forecasting model or learning algorithm to estimate either the scalar-valued function $f$~(see Equations \ref{stochRec}, \ref{stochDirect} and \ref{stochDirRec}) or the vector-valued function $F$~(see Equations \ref{stochMIMO} and \ref{stochDIRMO}) which represent the temporal stochastic dependencies. As the goal of the paper is not to compare forecasting models~(as in \citep{empiricalMLTimeSeries}) but rather multi-step ahead forecasting strategies, the choice of a underlying forecasting model is required to setup the experiments. In this paper, we adopted the Lazy Learning algorithm, which is a particular instance of local learning models, since it has been shown to be particularly effective in time series forecasting tasks~\citep{Bontempi98BirattariBersini_Leuven,phdbontempi,mimo,mimoijf,mismo,BenTaieb20101950}.

Next section gives a general comparison of global models with local models. Section \ref{LLalgo} presents the Lazy Learning Algorithm in terms of learning properties. Section \ref{SOalgo} and \ref{MOalgo} describe two Lazy Learning algorithms for two types of learning tasks, namely the Single-Output and Multiple-Output Lazy Learning algorithms. Finally, the a discussion is presented on the model combination.

\subsection{Global vs local modeling for supervised learning}

Forecasting the future values of a time series using past observations can be reduced to a supervised learning problem or, more precisely, to a regression problem. Indeed, the time series can be seen as a dataset made of pairs where the first component, called input, is a past temporal pattern and the second, called output, is the corresponding future pattern. Being able to predict the unknown output for a given input is equivalent to forecasting the future values given the last observations of the time series.

\emph{Global modeling} is the typical approach to the supervised learning problem. Global models are parametric models that describe the relationship between the inputs and the output values as an analytical function over the whole input domain. Examples of global models are linear models~\citep{MontgomeryBook2006}, nonlinear statistical regressions~\citep{Seber89Wild} and Neural Networks~\citep{Rumelhart86HintonWilliams}.

Another approach is the \emph{divide-and-conquer} modeling which consists in relaxing the global modeling assumptions by dividing a complex problem into simpler problems, whose solutions can be combined to yield a solution to the original problem~\citep{phdbontempi}. The divide-and-conquer has evolved in two different paradigms: the \emph{modular architectures} and the \emph{local modeling} approach~\citep{phdbontempi}.

\emph{Modular techniques} replace a global model with different modules covering different parts of the input space. Examples based on this approach are Fuzzy Inference Systems~\citep{Takagi85Sugeno}, Radial Basis Functions~\citep{Moody89Darken,Poggio90Girosi}, Local Model Networks~\citep{MurraySmith94}, Trees~\citep{BreimanFriedmanOlshenStone84} and Mixture of Experts~\citep{Jordan94Jacobs}. The modular approach is in the intermediate scale between the two extremes, the global and the local approach. However, their identification is still performed on the basis of the whole dataset and requires the same procedures used for generic global models.

\emph{Local modeling} techniques are at the extreme end of divide-and-conquer methods. They are nonparametric models that combine excellent theoretical properties with a simple and flexible learning procedure. Indeed, they do not aim to return a complete description of the input/output mapping but rather to approximate the function in a neighborhood of the point to be predicted~(also called the query point). There are different examples of local models, for example nearest neighbor, weighted average, and locally weighted regression~\citep{Atkeson97locallyweighted}. Each of these models use data points near the point to be predicted for estimating the unknown output. \emph{Nearest neighbor} models simply find the closest point and uses its output value. \emph{Weighted average} models combines the closest points by averaging them with weights inversely proportional to their distance to the point to be predicted. \emph{Locally weighted regression} models fit a model to nearby points with a weighted regression where the weights are function of distances to the query point.

The effectiveness of local models is well-known in the time series and computational intelligence community. For example, the method proposed by Sauer~\citep{Sauer94} gave good performance and ranked second best forecast for the Santa Fe A dataset from a forecasting competition organized by Santa Fe institute. Moreover, the two top-ranked entries of the KULeuven competition used local learning methods~\citep{Bontempi98BirattariBersini_Leuven,McNames98}.

In this work, we will restrict to consider a particular instance of local modeling algorithms: the Lazy Learning algorithm~\citep{Aha97}.

\subsection{The Lazy Learning algorithm}
\label{LLalgo}

It is possible to encounter different degree of ``laziness'' in local learning algorithms. For instance, a $k$ Nearest Neighbor ($k$-NN) algorithm, which learns the best value of $k$ before the query is requested, is hardly a lazy approach since, after the query is presented, it requires only a reduced amount of learning procedure, only the computation of the neighbors and the average. On the contrary, a local method, which depends on the query to select the number of neighbors or other structural parameters presents a higher degree of ``laziness''.

The Lazy Learning~(LL) algorithm, extensively discussed in \citep{Birattari99BontempiBersini,bontempi}, is a \emph{query-based} local modeling technique where \emph{the whole learning procedure} is deferred until a forecast is required. When the query is requested, the learning procedure may start to select the best value of the number of neighbors~(or other structural parameters) and next, the dataset is searched for the nearest neighbors of the query point. The nearest neighbors are then used for estimating a local model, which returns a forecast. The local model is then discarded and the procedure is repeated from scratch for subsequent queries.

The LL algorithm has  a number of attractive features~\citep{Aha97}, namely, the reduced number of assumptions, the online learning capability and the capacity to model nonstationarity. LL assumes no a priori knowledge on the process underlying the data, which is particularly relevant in real datasets. These considerations motivate the adoption of the LL algorithm as a learning model in a multi-step ahead forecasting context.

Local modeling techniques require the definition of a set of model parameters namely \emph{the number $k$ of neighbors}, \emph{the kernel function}, \emph{the parametric family} and \emph{the distance metric}[REF]. In the literature, different methods have been proposed to select automatically the adequate configuration~\citep{Atkeson97locallyweighted,Birattari99BontempiBersini}. However, in this paper, we will limit the search on only the selection of the number of neighbors~(also called or equivalent to the bandwidth selection). This is essentially the most critical parameter, as it controls the bias/variance trade-off. Bandwidth selection is usually performed by rule-of-thumb techniques~\citep{Fan95Gijbels}, plug-in methods~\citep{Ruppert95SheatherWand} or cross-validation strategies~\citep{Atkeson97MooreSchaal}. Concerning the other parameters, we use the tricubic kernel~\citep{Cleveland198887} as kernel function, a constant model for the parametric family and the euclidean distance as metric.

Note that in order to apply local learning to a time series, we need to embed it into a dataset made of pairs where the first component is a temporal pattern of length $d$ and the second component is either the future value~(in the case of Single-Output Modeling) or the consecutive temporal pattern of length $H$~(in the case of Multiple-Output Modeling). In the following sections, $D$ will refer to the embedded time series with $M$ input/output pairs.

\subsection{Single-Output Lazy Learning algorithm}
\label{SOalgo}

In the case of Single-Output learning~(i.e with scalar output), the Lazy Learning procedure consists of a sequence of steps detailed in Algorithm \ref{LLSO}. The algorithm assesses the generalization performance of different local models and compares them in order to select the best one in terms of generalization capability. To perform that, the algorithm associate a Leave-One-Out~(LOO) error $e_{LOO}(k)$ to the estimation $y_{q}^{(k)}$ obtained with $k$ neighbors~(lines \ref{LOOLLSO2} and \ref{LOOLLSO}).

The LOO error can provide a reliable estimate of the generalization capability. However the disadvantage of such an approach is that it requires to repeat $k$ times the training process, which means a large computational effort. Fortunately, in the case of linear models there exists a powerful statistical procedure to compute the LOO cross-validation measure at a reduced computational cost: the PRESS~(Prediction Sum of Squares) statistic~\citep{Allen1974}.

In case of constant model, the LOO error $e_{LOO}(k)$ for the estimation $y_{q}^{(k)}$ of the query point $\mathbf{x}_{q}$ is calculated as follows~\citep{phdbontempi}:

\begin{equation} e_{LOO}(k)=\frac{1}{k}\sum_{j=1}^{k} (e_{j}(k))^2  \label{pressconst} \text{,}\end{equation}

\noindent where $e_{j}(k)$ designates the error obtained by setting aside the $j$-th neighbor of $\mathbf{x}_{q}$~($j \in \{1,\dots,k\}$). If we define the output of the $k$ closest neighbors of $\mathbf{x}_{q}$ as $\{y_{[1]},\dots,y_{[k]}\}$ then, $e_{j}(k)$ is defined as

\begin{align}
e_{j}(k) &= y_{[j]} - \frac{\sum_{i=1(i \neq j)}^{k} y_{[i]}}{k-1} \label{looNormal}\\
&= \frac{(k-1)y_{[j]} - \sum_{i=1(i \neq j)}^{k} y_{[i]}}{k-1} \\
&= \frac{(k-1)y_{[j]} +y_{[j]} - y_{[j]} - \sum_{i=1(i \neq j)}^{k} y_{[i]}}{k-1} \\
&= \frac{ky_{[j]} - \sum_{i=1}^{k} y_{[i]}}{k-1} \\
&=  \boxed{ k \frac{y_{[j]}-y_{q}^{(k)}}{k-1} }\text{.} \label{looPRESS}
\end{align}

Note that if we use Equation \ref{looNormal} to calculate the LOO error~(Equation \ref{pressconst}), the training process is repeated $k$ times since the sum in Equation \ref{looNormal} is performed for each index $j$. However, by using the PRESS statistic~(Equation \ref{looPRESS}), we avoid this large computational effort since the sum is replaced by the previously computed $y_{q}^{(k)}$ which was already calculated. This makes the PRESS statistic an efficient method to compute the LOO error.

After evaluating the performance of local models with different number of neighbors $k$~(lines \ref{LLSOFORSTART} to \ref{LLSOFOREND}), the best one which minimizes the LOO error~(having index $k^{*}$) is selected (lines \ref{bestlook} and \ref{bestlook2}). Finally, the prediction of the output of $\mathbf{x}_{q}$ is returned~(line \ref{LOOLLSOreturn}).

\IncMargin{0.1cm}
\begin{algorithm}[h]
\SetKwInOut{Input}{Input}
\SetKwInOut{Output}{Output}
\caption{Single-Output Lazy Learning \label{LLSO}}
\Input{$D=\{(\mathbf{x}_i,y_{i}) \in (\mathbb{R}^{d}\times \mathbb{R}) \text{ with } i \in \{1,\dots,M\}\}$, dataset.}
\Input{$\mathbf{x}_{q} \in \mathbb{R}^{d}$, query point.}
\Input{$Kmax$, the maximum number of neighbors.}
\Output{$\hat{y}_{q}$, the prediction of the (scalar) output of the query point $\mathbf{x}_{q}$.}
\BlankLine
Sort increasingly the set of vectors $\{\mathbf{x}_{i}\}$ with respect to the distance to $\mathbf{x}_{q}$. \\
$[j]$ designate the index of the $j$th closest neighbor of $\mathbf{x}_q$. \\
\BlankLine

 \For{$k \in \{2,\dots,Kmax\}$}{ \nllabel{LLSOFORSTART}
$y_{q}^{(k)}= \frac{1}{k}\sum_{j=1}^{k} y_{[j]}$.\nllabel{LOOLLSO2}  \\
Calculate $e_{LOO}(k)$ which is defined in Equation \ref{pressconst}. \nllabel{LOOLLSO}  \\
} \nllabel{LLSOFOREND}
 $k^{*} = \text{arg min$_{k \in \{2,\dots,Kmax\}}$ } e_{LOO}(k)$. \nllabel{bestlook} \\
 $\hat{y}_{q}=y_{q}^{(k^{*})}$. \nllabel{bestlook2}  \\
\BlankLine
\BlankLine
\Return $\hat{y}_{q}$.\nllabel{LOOLLSOreturn}  \\
\end{algorithm}
\DecMargin{0.1cm}

\subsection{Multiple-Output Lazy Learning algorithm}
\label{MOalgo}

The adoption of Multiple-Output strategies requires the design of multiple-output~(or equivalently multi-response) modeling techniques~\citep{matias,multiresponse,vector-valued} where the output is no more a scalar quantity but a vector of values. Like in the Single-Output case, we need criteria to assess and compare local models with different number of neighbors. In the following, we present two criteria: the first one is  an extension of the LOO error for the Multiple-Output case~(Algorithm \ref{LLMO})~\citep{mimo,BenTaieb20101950} and the second one is a criterion proper to the Multiple-Output modeling~(Algorithm \ref{LLDMO})~\citep{BenTaieb20101950,mimoijf}. Note that, in the two algorithms, the output is a vector of size $l$~(e.g. $l$ will equal $H$ with the MIMO strategy or $s$ in the DIRMO strategy).

\IncMargin{0.1cm}
\begin{algorithm}[h]
\SetKwInOut{Input}{Input}
\SetKwInOut{Output}{Output}

\caption{Multiple-Output Lazy Learning~(\emph{LOO criterion}): MIMO-LOO \label{LLMO}}
\Input{$D=\{(\mathbf{x}_i,\mathbf{y}_{i}) \in (\mathbb{R}^{d}\times \mathbb{R}^{l}) \text{ with } i \in \{1,\dots,M\} \}$, dataset.}
\Input{$\mathbf{x}_{q} \in \mathbb{R}^{d}$, query point.}
\Input{$Kmax$, the maximum number of neighbors.}
\Output{$\mathbf{\hat{y}}_{q}$, the prediction of the (vectorial) output of the query point $\mathbf{x}_{q}$.}
\BlankLine
Sort increasingly the set of vectors $\{\mathbf{x}_{i}\}$ with respect to the distance to $\mathbf{x}_{q}$. \\
$[j]$ will designate the index of the $j$th closest neighbor of $\mathbf{x}_q$. \\
\BlankLine

\For{$k \in \{2,\dots,Kmax\}$}{
$\mathbf{y}_{q}^{(k)}= \frac{1}{k}\sum_{j=1}^{k} \mathbf{y}_{[j]}$. \\
$ E_{LOO}(k)=\frac{1}{l} \sum_{h=1}^{l} e_{LOO}^{l}(k)$ where $e_{LOO}^{l}(k) $ is defined in Equation \ref{pressconst}. \nllabel{LOOLLMO} \\

}
$k^{*} = \text{arg min$_{k \in \{2,\dots,Kmax\}}$ } E_{LOO}(k)$. \\
$\mathbf{\hat{y}}_{q}=\mathbf{y}_{q}^{(k^{*})}$
\BlankLine
\BlankLine
\Return $\mathbf{\hat{y}}_{q}$. \\
\end{algorithm}
\DecMargin{0.1cm}

Algorithm \ref{LLMO} is an extension of the Algorithm \ref{LLSO} for vectorial outputs. We still use the LOO cross-validation measure as a criterion to estimate the generalization capability of the model but here, the LOO error is an aggregation of the errors obtained for each output~(line \ref{LOOLLMO}). Note that the same number of neighbors is selected for all the outputs~(e.g. MIMO strategy) unlike what could happen with different Single-Output tasks~(e.g. Direct strategy).
\clearpage

The second criterion uses the fact that the forecasting horizon $H$ is supposed to be large~(multi-step ahead forecasting) and hence we have enough samples to estimate some descriptive statistics. Then, instead of using the Leave-One-Out error, we can use as criterion a measure of stochastic discrepancy between the forecasted values and the training time series. The lower the discrepancy between the descriptors of the forecasts and the training time series, the better is the quality of the forecasts~\citep{mimoijf}.

Several measures of discrepancy can be defined, both linear and non-linear. For example, the autocorrelation can be used as linear statistics and the maximum likelihood as a non-linear one. In this work, we will consider only one linear measure using both the autocorrelation and the partial correlation.

The assessement of the quality of the estimation $y_{q}^{(k)}$ of the query point $\mathbf{x}_{q}$ is calculated as follows

\begin{equation} E_{\Delta}(k)= \underbrace{1-|cor[\rho(ts \cdot \mathbf{y}_{q}^{(k)} ),\rho(ts)]| }_{}+ \underbrace{ 1-|cor[\pi(ts \cdot \mathbf{y}_{q}^{(k)}),\pi(ts)]| }_{} \label{dCriterion} \text{,} \end{equation}

\noindent where the symbol ``$\cdot$'' represents the concatenation, $ts$ represent the training time series and $cor$ is the Pearson correlation. This discrepancy measure is composed of two parts where the first part uses the autocorrelation~(noted $\rho$) and the second uses the partial autocorrelation~(noted $\pi$).

For each part, we calculate the discrepancy~(estimated with the correlation, noted $cor$) between, on one hand, the autocorrelation~(or partial autocorrelation) of the concatenation of the training time series $ts$  and the forecasted sequence $\mathbf{y}_{q}^{(k)}$ and, on the other hand, the autocorrelation~(or partial autocorrelation) of the training time series $ts$~\citep{mimo,mismo}.

In Algorithm \ref{LLDMO}, after evaluating the performance of local models with different number of neighbors $k$ (lines \ref{LLSOFORSTART}
to \ref{LLSOFOREND}), the best one, which minimizes the discrepancy between the forecasting sequence and the training time series~(having index $k^*$), is selected~(lines \ref{bestd} and \ref{bestd2}). In other words, the goal is to select the best number of neighbors $k^{*}$ which preserve the stochastic properties of the time series in the forecasted sequence. Finally, the prediction of the output of $\mathbf{x}_{q}$ is returned (line \ref{bestpred}).

\IncMargin{0.1cm}
\begin{algorithm}[h]
\SetKwInOut{Input}{Input}
\SetKwInOut{Output}{Output}

\caption{Multiple-Output Lazy Learning~(\emph{discreprancy criterion}): MIMO-ACFLIN \label{LLDMO}}
\Input{$ts=[ts_1,\dots,ts_N]$, time series.}
\Input{$D=\{(\mathbf{x}_i,\mathbf{y}_{i}) \in (\mathbb{R}^{d}\times \mathbb{R}^{l}) \text{ with } i \in \{1,\dots,M\} \}$, dataset.}
\Input{$\mathbf{x}_{q} \in \mathbb{R}^{d}$, query point.}
\Input{$Kmax$, the maximum number of neighbors.}
\Output{$\mathbf{\hat{y}}_{q}$, the prediction of the (vectorial) output of the query point $\mathbf{x}_{q}$.}
\BlankLine
Sort increasingly the set of vectors $\{\mathbf{x}_{i}\}$ with respect to the distance to $\mathbf{x}_{q}$. \\
$[j]$ will designate the index of the $j$th closest neighbor of $\mathbf{x}_q$. \\
\BlankLine

\For{$k \in \{2,\dots,Kmax\}$}{ \label{DFORSTART}
$\mathbf{y}_{q}^{(k)}= \frac{1}{k}\sum_{j=1}^{k} \mathbf{y}_{[j]}$. \\
Calculate $E_{\Delta}(k)$ which is defined in Equation \ref{dCriterion}. \\
} \label{DFOREND}
$k^{*} = \text{arg min$_{k \in \{2,\dots,Kmax\}}$ } E_{\Delta}(k)$. \nllabel{bestd} \\
$\mathbf{\hat{y}}_{q}=\mathbf{y}_{q}^{(k^{*})}$ \nllabel{bestd2}
\BlankLine
\BlankLine
\Return $\mathbf{\hat{y}}_{q}$. \nllabel{bestpred} \\
\end{algorithm}
\DecMargin{0.1cm}

\clearpage

\subsection{Model selection or model averaging}
\label{modelCombination}

Considering the Algorithm \ref{LLSO}, we can see that we generate, for the query $\mathbf{x}_{q}$, a set of predictions $\{y_{q}^{(2)},y_{q}^{(3)},\dots,y_{q}^{(Kmax)}\}$, each obtained with different number of neighbors. For each of these predictions, a testing error $\{e_{LOO}(2),e_{LOO}(3),\dots,e_{LOO}(Kmax)\}$ has been calculated. Note that the next considerations are also applicable to Algorithms \ref{LLMO} and \ref{LLDMO}.

The goal of model selection is to use all this information~(set of predictions and testing errors) to estimate the final prediction $\hat{y}_{q}$ of the query point $\mathbf{x}_{q}$. There exist two main paradigms mainly the \emph{winner-take-all} and \emph{combination} approaches.

In the Algorithm \ref{LLSO}, we presented the winner-take-all approach~(noted WINNER)~\citep{racing97} which consists of comparing the set of models $y_{q}^{(k)}$ and selecting the best one in terms of testing error $e_{LOO}(k)$~(see line \ref{bestlook}).

Selecting the best model according to the testing error is intuitively the approach which should work the best. However, results in machine learning show that the performance of the final model can be improved by combining models having different structures~\citep{Raudys2006,Jacobs1991,Breiman1996,Schapire1998}.

In order to apply the model averaging, lines $7$ and $8$  of the Algorithm \ref{LLSO} can be replaced by

\begin{equation}  \hat{y}_{q}=    \frac{ p_{2} \text{ } y_{q}^{(2)} + \dots +  p_{Kmax} \text{ } y_{q}^{(Kmax)}         }  {\sum_{k=2}^{Kmax} p_{k}    }   \text{,} \end{equation}

\noindent where an average is calculated. The weights $p_{k}$ will take different values depending on the combination approach adopted. If $p_{k}$ equals $\frac{1}{Kmax}$, we are in the case of equally weighted combination and $\hat{y}_{q}$ reduces to an arithmetic mean~(noted COMB). Otherwise, if weights are assigned according to testing errors, $p_{k}$ will equal $\frac{1}{e_{LOO}(k)}$ and $\hat{y}_{q}$ reduces to a weighted mean~(noted WCOMB).

\clearpage

\section{Experimental Setup}
\label{sectionExperiments}

\subsection{Time Series Data}

In the last decade, several time series forecasting competitions~(e.g. the NN3, NN5, and the ESTSP competitions~\citep{nn3,nn5,estsp07,estsp08}) have been organized in order to compare and evaluate the performance of computational intelligence methods. Among them, the NN5 competition~\citep{nn5} is one of the most interesting one since it includes the challenges of a real-world multi-step ahead forecasting task, namely multiple time series, outliers, missing values as well as multiple overlying seasonalities, etc. Figure \ref{NN5TimeSeries} shows four time series from the NN5 dataset.

Each of the $111$ time series of this competition represents roughly two years of daily cash money withdrawal amounts~($735$ data points) at ATM machines at one of the various cities in the UK. For each time series, the competition required to forecast the values of the next $56$ days, using the given historical data points, as accurately as possible. The performance of the forecasting methods over one time series was assessed by the symmetric mean absolute percentage of error~(SMAPE) measure~\citep{nn5}, defined as

\begin{equation} \text{SMAPE} = \frac{1}{H} \sum_{h=1}^{H} \frac{|\hat{y}_h-y_h|}{(\hat{y}_h+y_h)/2} \times 100 \textnormal{,} \end{equation}

\noindent where $y_h$ is the target output and $\hat{y}_h$ is the prediction. Since this is a relative error measure, the errors can be averaged over all time series to obtain a mean SMAPE defined as

\begin{equation} \text{SMAPE}^{*} = \frac{1}{111} \sum_{i=1}^{111} \text{SMAPE}_{i} \textnormal{,} \label{smapeCriterion} \end{equation}

\noindent where SMAPE$_{i}$ denotes the SMAPE of the $i$th time series.

\clearpage

\begin{figure}[h!]
\centering
\vspace{-0.3cm}
\includegraphics[trim = 0cm 0.5cm 0cm 1.9cm, clip,scale=0.8,scale=0.6]{NN5TimeSeries1.pdf} \\
\includegraphics[trim = 0cm 0.5cm 0cm 1.9cm, clip,scale=0.8,scale=0.6]{NN5TimeSeries2.pdf} \\
\includegraphics[trim = 0cm 0.5cm 0cm 1.9cm, clip,scale=0.8,scale=0.6]{NN5TimeSeries3.pdf} \\
\includegraphics[trim = 0cm 0.5cm 0cm 1.9cm, clip,scale=0.8,scale=0.6]{NN5TimeSeries4.pdf} \\
\caption{Four time series from NN5 time series forecasting competition. \label{NN5TimeSeries}}
\end{figure}

\clearpage

\subsection{Methodology}
\label{secMethodo}

The aim of the experimental study is to compare the accuracy of the five forecasting strategies in the
context of the NN5 competition. Since the accuracy of a forecasting technique is known to be dependent on several design choices (e.g. the deseasonalization or the input selection) and we want to focus our analysis on the multi-step ahead forecasting strategies, we consider a number of different configurations in order to increase the statistical power of our comparison. Every configuration is composed of several preprocessing steps as sketched in Figure  \ref{methodoGeneral}. Since some of  these steps (e.g. deseasonalization) can be performed in alternative ways (e.g. two alternatives for the deseasonalization, two alternatives for input selection, three alternatives for  the model selection), we come up with $12$ configurations. The details about each step are given in what follows.

\begin{figure}[h!]
\centering
\includegraphics[scale=0.4]{methodoNN5.pdf}
\caption{The Different Preprocessing Steps. \label{methodoGeneral}}
\end{figure}

\noindent \textbf{Step 1: Gaps removal}

The specificity of the NN5 series requires a preprocessing step called \emph{gaps removal} where by gap we mean two types of
anomalies:  (i) zero values that indicate that no money withdrawal occurred and (ii) missing observations for which no value was recorded.
About  $2.5 \%$ of the data are corrupted by gaps. In our experiments we adopted the gap removal method proposed in~\citep{Wichard2010}: if $y_m$ is the gap sample, this method replaces the gap with the median of the set $[y_{m+365}, y_{m-365}, y_{m+7} \text{ and } y_{m-7}]$ among which are available.

\noindent \textbf{Step 2: Deseasonalization}
The adoption of deseasonalization may have a large impact on the forecasting strategies
because the NN5 time series possess a variety of periodic patterns.
For that reason we decided to consider tasks with and without deseasonalization in order to better account for
the role of the forecasting strategy.
We adopt the deseasonalization methodology discussed in \citep{amirnn5} to remove the strong day of the week seasonality as well as the moderate day of the month seasonality. Of course after we deseasonalize and apply the forecasting model
we restore back the seasonality.

\noindent \textbf{Step 3: Embedding dimension selection}

Every forecasting strategy requires the setting of the size $d$  of the embedding dimension~(see Equations \ref{stochRec} to \ref{stochDIRMO}). Several approaches have been proposed in the literature to select this value~\citep{nonlinear}. Since this aspect is not a central theme
in our paper we just applied the state-of-the-art approach reviewed in \citep{Crone:2009:ISN:1704555.1704739}, which consists of selecting the  time-lagged realizations with significant partial correlation function~(PACF). This method allows to select the value of the embedding dimension and then to identify the relevant variables within the window of past observations. We set the maximum lag of the PACF to $200$ to provide a sufficiently comprehensive pool of features. However, note that the final dimensionality of the input vectors for all the time series is on average equal to $24$.

\noindent \textbf{Step 4: Input Selection}

We considered the forecasting task with and without input variable selection step.
A variable selection procedure requires the setting of  two elements: the \emph{relevance criterion}, i.e. statistics which estimates the quality of the selected variables, and the \emph{search procedure}, which describes the policy to explore the input space.
We adopted the \emph{Delta test}(DT) as  relevance criterion. The DT has been introduced in time series forecasting domain by Pi and Peterson in \citep{Pi94} and later successfully applied to several forecasting task~\citep{mismo,BenTaieb20101950,Liitiainen07,guillen08newdelta,mateo10_neuro}. This criterion is based on applying some kind of a noise variance estimator, and then selecting the set of input variables that yield the strongest and most deterministic dependence between inputs and outputs~\citep{Mateo08}.

Concerning the search procedure, we adopted a Forward-Backward Search (FBS) procedure which is a combination of forward selection (sequentially adding input variables) and backward search (sequentially removing some input variables). This choice was motivated by the flexibility of the FBS procedure which allows a deeper exploration of the input space. Note that the search is initialized by the set of variables
defined in the previous step.

\noindent \textbf{Step 5: Model Selection}

Concerning the model selection procedure, three approaches~(see Section \ref{modelCombination}) are taken into consideration in our experiments:

\begin{description}
\item [WINNER]: This approach selects the model that gives best performance for the test set (winner-take-all approach).
\item [COMB]: This approach combines all alternative models by simple averaging.
\item [WCOMB]: This approach combines models by weighted averaging where weights are inversely proportional to the test errors.

\end{description}

\subsubsection{The Compared forecasting strategies}

Table \ref{tableStrategies} presents the eight forecasting strategies that we tested, showing also their respective acronyms.
\begin{table}[h]
\centering
\begin{tabular}{|l|M{4cm}|M{6cm}|}
\hline
\cline{1-3}
$1.$ REC  & \multicolumn{2}{l|}{The \emph{Recursive} forecasting strategy.}  \tabularnewline
\hline
\cline{1-3}
$2.$ DIR & \multicolumn{2}{l|}{The \emph{Direct} forecasting strategy.}  \tabularnewline
\hline
\cline{1-3}
$3.$ DIRREC & \multicolumn{2}{l|}{The \emph{DirRec} forecasting strategy.}  \tabularnewline
\hline
\cline{1-3}
\multicolumn{1}{|c|}{ \multirow{2}{*}{\vspace{-1cm} MIMO} }  & $4.$ MIMO-LOO & A variant of the \emph{MIMO} forecasting strategy with the LOO selection criteria. \tabularnewline
\cline{2-3}
  & $5.$ MIMO-ACFLIN & A variant of the \emph{MIMO} forecasting strategy with the autocorrelation selection criteria. \tabularnewline
\hline
\cline{1-3}
\multicolumn{1}{|c|}{ \multirow{3}{*}{\vspace{-2.4cm} DIRMO } }  & $6.$ DIRMO-SEL & The \emph{DIRMO} forecasting strategy which select the best value of the parameter $s$. \tabularnewline
  \cline{2-3}
  & $7.$ DIRMO-AVG & A variant of the \emph{DIRMO} strategy which calculates a simple average of the forecasts obtained with different values of the parameter $s$. \tabularnewline
  \cline{2-3}
  & $8.$ DIRMO-WAVG & The DIRMO-AVG with a weighted average where weights are inversely proportional to testing errors. \tabularnewline
\hline
\cline{1-3}
\end{tabular}
\caption{The five forecasting strategies with their respective variants which gives eight forecasting strategies. \label{tableStrategies}}
\end{table}

\clearpage

\subsubsection{Forecasting performance evaluation}
\label{sectionFevaluation}
This section describes the assessment procedure (Figure~\ref{PerformanceEvaluation}) of the 8 forecasting strategies.

\begin{figure}[h!]
\centering
\includegraphics[scale=0.5]{PerformanceEvaluation.pdf}
\caption{Different steps of the forecasting performance evaluation. \label{PerformanceEvaluation}}
\end{figure}

The procedure for comparing between the eight forecasting strategies is shown in Figure~\ref{PerformanceEvaluation}. The accuracy of each forecasting strategy is first measured using the SMAPE* measure calculated over the $111$ time series and defined in Equation \ref{smapeCriterion}.
To test if there are significant general differences
in performance between the different strategies,
we have to consider the problem of comparing multiple models
on multiple data sets. For such case Dem\v{s}ar ~\citep{Demsar_statisticalcomparisons,extensionDemsar} in a detailed comparative study
recommended using a two stage procedure: first to apply
Friedman's or Iman and Davenport's tests to test if the compared models
have the same mean rank. If this test rejects the null-hypothesis, then
post-hoc pairwise tests are to be performed to compare
the different models. These tests adjust the critical values higher
to ensure that there is at most a 5\% chance that one of
the pairwise differences will be erroneously found significant.

\noindent \textbf{Friedman test}

The Friedman test~\citep{friedman1937,friedman1940} is a non-parametric procedure which tests the significance of differences between multiple ranks. It ranks the algorithms for each dataset separately: the rank of 1 will be given to the best performing algorithm, the rank of 2 to the second best and so on. Note that average ranks are assigned in case of ties.

After ranking the algorithms for each dataset, the Friedman test compares the average ranks of algorithms. Let $r^{i}_{j}$ be the rank of the $j$-th of $k$ algorithms on the $i$-th of $N$ data sets, the average rank of the $j$-th algorithm is $R_{j}=\frac{1}{N}\sum_{i} r^{i}_{j}$.

The null-hypothesis states that all the algorithms are equivalent and so their ranks $R_{j}$ should be equal. Under the null-hypothesis, the Friedman statistic

\begin{equation}
Q = \frac{12N}{k(k+1)} \left[  \sum_{j} R_{j}^{2} - \frac{k(k+1)^2 }{4}  \right]
\end{equation}

\noindent is distributed according to a chi-squared with $k-1$ degrees of freedom~($\chi^{2}_{k-1}$), when $N$ and $k$ are large enough~(as a rule of a thumb, $N>10$ and $k>5$)~\citep{Demsar_statisticalcomparisons}.

Iman and Danvenport~\citep{Iman80approximationsof}, showing that Friedman's statistic is undesirably conservative, derived another
improved statistic, given by

\begin{equation}
S = \frac{(N-1) Q}{N(k-1)-Q}
\end{equation}

\noindent which is distributed, under the null-hypothesis, according to the F-distribution with $k-1$ and $(k-1)(N-1)$ degrees of freedom.

\noindent \textbf{Post-hoc test}

When the null-hypothesis is rejected, i.e. there is a significant difference between at least 2 strategies, a post-hoc test is performed to identify significant pairwise differences among all the algorithms. The test statistic for comparing the $i$-th and the $j$-th algorithm is

\begin{equation}
z=\frac{(R_{i}-R_{j})}{\sqrt{\frac{k(k+1)}{6N}}} \text{,}
\end{equation}

\noindent which is asymptotically normally distributed under the null hypothesis. After the corresponding $p$-value is calculated, it is compared with a given level of significance $\alpha$.

However, in multiple comparisons, as there are a possibly large number of pairwise comparisons, there is a relatively high chance that some pairwise test are incorrectly rejected. Several procedures exist to adjust the value of $\alpha$ to compensate for this bias, for instance Nemenyi, Holm, Shaffer as well as Bergmann and Hommel~\citep{Demsar_statisticalcomparisons}. Based on the suggestion of Garcia and Herrera~\citep{extensionDemsar} we adopted Shaffer's correction. The reason is that Garcia and Herrera~\citep{extensionDemsar} showed that Shaffer's procedure has the same complexity Holm's procedure, but with the advantage of using information about logically related hypothesis.

\subsection{Experimental phases}

In order to reproduce the same context of the NN5 forecasting competition the experimental setting is made of two phases: the \emph{pre-competition} and the \emph{competition} phase.

\subsubsection{Pre-competition phase}

The pre-competition phase is devoted to the comparison of the different forecasting strategies using the available observations of $111$ time series. The goal is to learn the different parameters and then estimate the forecasting performance and compare between the different strategies.

To estimate the forecasting performance of each strategy, we used a learning scheme with training-validation-testing sets. Each time series~(containing $735$ observations) is partitioned in three mutually exclusive sets~(A, B and C) as shown in Figure \ref{lPeriod}: training~(Day $1$ to Day $623$: $623$ values), validation~(Day $624$ to Day $679$: $56$ values) and testing~(Day $680$ to Day $735$: $56$ values).

The validation set~(B in Figure \ref{lPeriod}) is used to build and tune the models. Specifically, as we use a Lazy Learning approach, we need  to select, for each model, the range of number of neighbors($[2,\dots,Kmax]$) to use in performing the forecasting task.

The test set~(C in Figure \ref{lPeriod}) is used to measure the performances of each forecasting strategy. To make the utmost use of the available data, we adopt a multiple time origin test as suggested by Tashman in~\citep{Tashman00}, where the time origin denotes the point from which the multi-step ahead forecasts are generated.

The time origin and corresponding forecast intervals are given as:

\begin{enumerate}
\item Day $680$ to Day $735$~($56$ data points)
\item Day $687$ to Day $735$~($49$ data points)
\item Day $694$ to Day $735$~($42$ data points)
\end{enumerate}

In other words, we perform the forecast three times starting from the three different starting points, each time forecasting
a number of steps ahead till the end of the interval.
Note that we used the same test period and evaluation criterion (i.e. the SMAPE) as used by Andrawis et al in~\citep{amirnn5}. This allows us to compare our results with several other machine learning models tested in this article.

\begin{figure}[h!]
\centering
\includegraphics[trim = 0cm 0.5cm 0cm 1.9cm, clip,scale=0.8,scale=0.6]{learningPeriod.pdf}
\caption{Learning with three mutually exclusive sets for training~(A), validation~(B) and testing~(C).\label{lPeriod}}
\end{figure}

\subsubsection{Competition phase}

In the competition phase we generate the final forecasts, made up of $56$ future observations, which would have been submitted to the competition. This phase takes advantage of the lessons learned and the design choices made in the pre-competition phase. Here, we combine the training set with the test set~(A+B in Figure \ref{fPeriod}) to retrain the models of the different strategies and then generate the final forecast~(which will be submitted to the competition). The training set~(A+B on Figure \ref{fPeriod}) is now made of $679$ data points and the validation set~(C on Figure \ref{fPeriod}) is composed of the next $56$ data points, as shown in Figure \ref{fPeriod}. In other words, the $735$ values are then used to build and tune the models, which will next return the forecasted values.

\begin{figure}[h!]
\centering
\includegraphics[trim = 0cm 0.5cm 0cm 1.9cm, clip,scale=0.8,scale=0.6]{forecastingPeriod.pdf}
\caption{Forecasting with two mutually exclusive sets for training~(A+B) and validation~(C). \label{fPeriod}}
\end{figure}

\clearpage

\section{Results and discussion}
\label{sectionResDiscussion}

This section presents and discusses the prediction results of the forecasting strategies for the \emph{pre-competition} and \emph{competition} phase. For each phase, we report the results obtained in the $12$ different configurations introduced in Section~\ref{secMethodo}.

The forecasting performance of the strategies is measured using the criteria discussed in Section~\ref{sectionFevaluation} and presented by means of two tables.  The first one provides the average SMAPE as well as the ranking for each forecasting strategy. Since the null-hypothesis stating that all the algorithms are equivalent
has been rejected (using the Iman-Davenport statistic) for all the configurations, we proceeded with the post-hoc test. The second table presents the results of this post-hoc test, which partitions the set of strategies in several groups which are statistically signicantly different in terms of forecasting performance.

Note that the configurations which require the input selection do not contain the DIRMO results since combining the selection of the inputs with the selection of the parameter $s$ would have needed an excessive amount of computation time.

\subsection{Pre-competition results}

The SMAPE and ranking results for the \emph{pre-competition} phase are presented in Table~\ref{tableErrors.precompetition} while the results of the post-hoc test are summarized in Table~\ref{posthoc.precompetition}.

\input{tableErrors.precompetition}
\input{posthoc.precompetition}

\clearpage
The availability of the SMAPE* results, obtained according to the procedure used in \citep{amirnn5},  makes
possible the comparison of our pre-competition results with those of several others learning methods
available in \citep{amirnn5}. For the sake of comparison, Table \ref{amirTable} reports the forecasting errors
for some of the techniques considered in \citep{amirnn5}, notably Gaussian Process Regression~(GPR), Neural Network~(NN), Multiple Regression~(MULT-REGR), Simple Moving Average~(MOV-AVG), Holt's Exponential Smoothing and a combination (Combined) of such techniques.
\begin{table}[!htp]
\centering
\begin{tabular}{|c|c|}
\hline
Model &SMAPE*\\
\hline
\hline
GPR-ITER & 19.90\\
\hline
GPR-DIR & 21.22\\
\hline
GPR-LEV & 20.19\\
\hline
NN-ITER & 21.11\\
\hline
NN-LEV & 19.83\\
\hline
MULT-REGR1 & 19.11\\
\hline
MULTI-REGR2 & 18.96\\
\hline
MULT-REGR3 & 18.94\\
\hline
MOV-AVG & 19.55\\
\hline
Holt's Exp Sm & 23.77\\
\hline
Combined & 18.95\\
\hline
\end{tabular}
\caption{Forecasting errors for the different forecasting models. \label{amirTable}}
\end{table}
The comparison shows that the best configuration of Table \ref{tableErrors.precompetition.WithoutSeas.InputSel}, that is the MIMO-ACFLIN strategy,  is competitive with all these models with a SMAPE* amounting to $18.81\%$.

\clearpage
\subsection{Competition results}

The SMAPE and ranking results for the \emph{competition} phase are presented in Table~\ref{tableErrors.competition} while the results of the post-hoc test are summarized in Table~\ref{posthoc.competition}.

\input{tableErrors.competition}
\input{posthoc.competition}

\clearpage

The pre-competition results presented in the previous section suggest us to use the MIMO-ACFLIN strategy with the \verb+Comb+ model selection approach by removing the seasonality and applying the input selection procedure, since this configuration obtains the smallest forecasting error~($18.81 \%$).

By using the MIMO-ACFLIN strategy and the corresponding configuration in the competition phase, we would generate forecasts with a SMAPE* equals to $20.28 \%$ which is quite good compared to the best computational intelligence entries of the competition as shown in Table \ref{NN5results}. Figure \ref{NN5TimeSeriesForecasts} shows the forecasts of the MIMO-ACFLIN strategy versus the actual values for four NN5 time series to illustrate the forecasting capability of this strategy.

\begin{table}[!htp]
\centering
\begin{tabular}{|c|c|}
\hline
Model &SMAPE*\\
\hline
\hline
\textbf{MIMO-ACFLIN} & $\mathbf{20.28}$\\
\hline
Andrawis & 20.4\\
\hline
Vogel & 20.5\\
\hline
D'yakonov & 20.6\\
\hline
Rauch & 21.7\\
\hline
Luna & 21.8\\
\hline
Wichard & 22.1\\
\hline
Gao & 22.3\\
\hline
Puma-Villanueva & 23.7\\
\hline
Dang & 23.77\\
\hline
Pasero & 25.3\\
\hline
\end{tabular}
\caption{Forecasting errors for different computational intelligence forecasting models which participate to the NN5 forecasting competition. \label{NN5results}}
\end{table}

\clearpage

\begin{figure}
\centering
\vspace{-0.3cm}
\includegraphics[trim = 0cm 0.5cm 0cm 1.9cm, clip,scale=0.8,scale=0.6]{NN5TimeSeriesForecasts1.pdf} \\
\includegraphics[trim = 0cm 0.5cm 0cm 1.9cm, clip,scale=0.8,scale=0.6]{NN5TimeSeriesForecasts2.pdf} \\
\includegraphics[trim = 0cm 0.5cm 0cm 1.9cm, clip,scale=0.8,scale=0.6]{NN5TimeSeriesForecasts3.pdf} \\
\includegraphics[trim = 0cm 0.5cm 0cm 1.9cm, clip,scale=0.8,scale=0.6]{NN5TimeSeriesForecasts4.pdf} \\
\caption{ The forecast versus the actual of the MIMO-ACFLIN strategy for four NN5 time series. \label{NN5TimeSeriesForecasts}}
\end{figure}

\clearpage

\subsection{Discussion}

From all presented results one can deduce the following observations below.
These findings refer mainly to the pre-competition results. But, one can easily see
that they mostly also apply to the competition phase results.

\begin{itemize}

\item The overall best method is MIMO-ACFLIN, used with input selection, deseasonalization
and equal weight combination (COMB).

\item The Multiple-Output strategies (MIMO and DIRMO) are invariably the best strategies.
They beat the Single-Output strategies, such as DIR, REC, and DIRREC.
Both MIMO and DIRMO give comparable performancce.
For DIRMO, the selection of the parameter $s$ is critical, since it has a great impact on the performance.
Should there be an improved selection approach, this strategy would have a big potential.

\item
Both versions of MIMO are comparable. Also the versions of
DIRMO give close results, with perhaps DIRMO-WAVG a little better
than the other two versions.

\item Among the Single-Output strategies, the REC strategy has almost always a smaller SMAPE and a better
ranking than the DIR strategy. DIRREC is the worse strategy overall,
and gives especially low accuracy when no deseasonalization is performed.

\item Deseasonalization leads to consistently better results (in 38 out of 39 models).
This result is consistent with some other studies, such as \citep{Zhang2005501}.
The possible reason for this is that when no deseasonalization is performed, we are putting
a higher burden on the model to forecast the
future seasonal pattern plus the trend and the other aspects,
which apparently is hard to simultaneously satisfy.

\item Input selection is especially beneficial when we perform
a deseasonalization. Absent deseasonalization, the results are mixed
(as to whether input selection improves the results or not).
The possible explanation is that when no deseasonalization is performed,
the model needs all the previous cycle to construct the future seasonal
pattern. Performing an input selection will deprive it from essential information.

\item Concerning the model selection aspect, both combination
approaches (COMB and WCOMB) are superior to the winner-take-all (WINNER).
Both COMB and WCOMB are comparable, and the results do not
differ by much. This is consistent with much of the findings
in forecast combination literature, e.g.~\citep{amirnn5,Clemen89,forecastcomb,Andrawis2010}

\item The relative performance and ranking of the different strategies
 is persistent. Most findings that are based on the pre-competition
 results are true for the competition phase results. This is also true
 for the findings concerning the deseasonalization, input selection, and model selection.
 This persistence is reassuring, as we can have some confidence
 in relying on the test or validation results for selecting the best strategies.

 \item The best strategy based on the pre-competition data, the MIMO-ACFLIN method,
 would have topped all computational intelligence entries of the NN5 competition in the true competition
 hold-out data.

\end{itemize}

\section{Conclusion}
\label{sectionConclusion}
Forecasting a time series many steps into the future is a very
hard problem because the larger the forecast horizon, the higher is the uncertainty.
In this paper we presented a comparative review of existing strategies for multi-step ahead forecasting, together with an extensive comparison, applied on the 111 time series of the NN5 forecasting competition.
The comparison gave some interesting lessons that could help researchers channel their
experiments into the most promising approaches.
The most consistent findings are that Multiple-Output approaches are invariably better than Single-Output
approaches. Also, deseasonalization had a very considerable positive impact on the
performance. Finally, the results are clearly quite persistent.
So, selecting the best strategy based on testing performance is a very potent approach.
A possible direction for future research could therefore be
developing other new improved Multiple-Output strategies. Also, possibly tailoring
deseasonalization methods specifically for Multiple-Output strategies could also be a promising
research point.

\section*{Acknowledgments}

We would like to thanks the authors of the paper \citep{extensionDemsar} for making their methods available at \url{http://sci2s.ugr.es/keel/multipleTest.zip}.

\bibliographystyle{plainnat}

redpajama_set_name: RedPajamaArXiv