This is a preprint copy that has been accepted for publication in IEEE
Transactions on Cybernetics.

Please cite this article as:
Yukun Bao, Tao Xiong, Zhongyi Hu, “PSO-MISMO Modeling Strategy
for Multi-Step-Ahead Time Series Prediction”. IEEE Transactions on
Cybernetics. 2013, accepted. doi:10.1109/TCYB.2013. 2265084.

Copyright Claim: It should be noted that the copyright is reserved by
IEEE. This preprint copy is only for personal use.

Copyright © 2013 IEEE. All rights reserved.

> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) <

1

PSO-MISMO Modeling Strategy for
Multi-Step-Ahead Time Series Prediction
Yukun Bao, Member, IEEE, Tao Xiong, Zhongyi Hu

Abstract—Multi-step-ahead time series prediction is one of the
most challenging research topics in the field of time series
modeling and prediction, and is continually under research.
Recently, the multiple-input several multiple-outputs (MISMO)
modeling strategy has been proposed as a promising alternative
for multi-step-ahead time series prediction, exhibiting advantages
compared with the two currently dominating strategies, the
iterated and the direct strategies. Built on the established MISMO
strategy, this study proposes a particle swarm optimization
(PSO)-based MISMO modeling strategy, which is capable of
determining the number of sub-models in a self-adaptive mode,
with varying prediction horizons. Rather than deriving crisp
divides with equal-size s prediction horizons from the established
MISMO, the proposed PSO-MISMO strategy, implemented with
neural networks, employs a heuristic to create flexible divides
with varying sizes of prediction horizons and to generate
corresponding sub-models, providing considerable flexibility in
model construction, which has been validated with simulated and
real datasets.
Index Terms—Multi-step-ahead time series prediction,
multiple-output models, particle swarm optimization, genetic
algorithm.

I. INTRODUCTION

D

espite intensive research efforts in time series modeling

and prediction, most models continue to be limited to
one-step-ahead prediction rather than multi-step-ahead
prediction, primarily because of the challenges arising from
increased uncertainty from longer prediction horizons. The
most common modeling strategies for multi-step-ahead time
series prediction rely either on iterated or direct strategies [1-3].
An iterated strategy first constructs a one-step-ahead prediction
model, and then uses the predicted values as known data to
Manuscript received September 12, 2012. This work was supported in part
by the Natural Science Foundation of China under Grant No. 70771042, the
Fundamental Research Funds for the Central Universities (2012QN208-HUST),
the MOE (Ministry of Education in China) Project of Humanities and Social
Science (Project No.13YJA630002), and a grant from the Modern Information
Management Research Center at Huazhong University of Science and
Technology.
Yukun Bao (corresponding author), is with school of management,
Huazhong University of Science and Technology, Wuhan, P.R.China (Tel:
86-27-87558579; fax: 86-27-87556437; e-mail: yukunbao@hust.edu.cn).
Tao Xiong, is with school of management, Huazhong University of Science
and Technology, Wuhan, P.R.China (e-mail: taoxiong@hust.edu.cn).
Zhongyi Hu, is with school of management, Huazhong University of
Science and Technology, Wuhan, P.R.China (e-mail: huzhyi21@hust.edu.cn).

predict the next ones. In the direct strategy, a set of H
prediction models are solved to predict values H steps ahead.
The accumulation of errors in the iterated case deteriorates the
accuracy of the prediction, and the direct strategy requires a
sharply increased computational expense, especially in the case
of a longer prediction horizon. Based upon these two
mainstream modeling strategies, there has been a revival of
interest in developing novel multi-step-ahead prediction
strategies such as the DirRec strategy [4], the multiple-input
multiple-output (MIMO) strategy [5] and multiple-input
several multiple-outputs (MISMO) strategy [6, 7] (DIRMO in
[8]). A review and comparative study on strategies for
multi-step-ahead time series has indicated that the
multiple-output strategies, including MIMO and MISMO
(DIRMO), are the best performing approaches, which has been
validated with NN5 forecasting competition data [8].
This study focuses on a MISMO strategy, improving it in
terms of modeling flexibility and more accurate prediction.
MISMO is achieved by transforming the original H
-step-ahead prediction task into n = H s subtasks, yielding n
sub-models with a fixed number of multiple outputs, s , for
each, where 1 ≤ s ≤ H . The value of s is chosen through
cross-validation, analyzing the performance of the MISMO
model for different values of s on the learning set, and then
using the best value to estimate the outputs [6, 7]. In the special
case where s = H , MISMO can be regarded as MIMO, and
when s = 1 , MISMO is the same as a direct strategy. The
superior out-of-sample performance of the MISMO strategy
presented in [7] depends on the trade-off between heterogeneity
and computational complexity. The MISMO strategy exploits
the heterogeneity across subtasks efficiently, with less
computational complexity than the direct prediction technique.
An important feature of the MISMO strategy is that the original
prediction modeling task is divided into several subtasks with
an equal number of outputs (prediction horizons). Handling
indivisibility in the problem is an explicit difficulty with the
method. In addition, an underlying problem arises from the
crisp divides on the prediction horizons with equal size s that
induce a latent separation of dependencies among the steps in
the different sub models. This problem may prevent the
MISMO strategy from considering complex dependencies
between prediction steps within different sub models and
consequently reduce the prediction accuracy. Thus a strategy
with more flexibility and a mechanism to self-adaptively
determine the prediction horizon divides while modeling is
appealing. This study proposes an improved PSO-MISMO

> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) <
modeling strategy for multi-step-ahead prediction that
incorporates a heuristic based on particle swarm optimization
(PSO) into the MISMO modeling process to self-adaptively
determine the number of sub-models with varying prediction
horizons.
For the purpose of justification, this study compares the rank
of the proposed modeling strategy with the four
well-established strategies (i.e., iterated, direct, MIMO and
MISMO) with neural networks (NNs). In addition, as a
straightforward alternative, binary genetic algorithm-based
MISMO modeling strategy (GA-MISMO) is selected as
counterpart against the proposed PSO-MISMO and compared
as well. Both simulated (i.e., Logistic and Mackey-Glass time
series) and real (i.e., NN3 competition data) datasets are used
for the comparisons.
The paper is structured as follows. Section II provides a brief
review of the multiple-output modeling strategies, including
MIMO and MISMO, with limited details. Section III describes
the basic concepts of FNN for multi-step-ahead prediction and
the implementation of the FNN-based prediction model in the
current study. The concept of the binary PSO is explained in
Section IV. In Section V, details of the proposed PSO-MISMO
strategy are presented. Section VI details the experimental
setting on data, accuracy measure, genetic algorithm as
alternative, and experimental procedure. Experimental results
are discussed in Section VII, and concluding remarks are
provided in Section VIII.
II. MIMO AND MISMO STRATEGIES
Multi-step-ahead prediction can be described as an estimation
of ϕ N + h , ( h = 1, 2, , H ) , where H is an integer greater than
one, given the current and previous observation
ϕt , (t = 1, 2, , N ) . In this section, two competing modeling
strategies, MIMO and MISMO, with multiple-input
multiple-output structures for multi-step-ahead prediction are
described, to enable understanding of the proposed
PSO-MISMO strategy in Section V.
A. MIMO Strategy
MIMO was first proposed in [5] and characterized as a
multiple-input multiple-output approach, where the predicted
value is not a scalar quantity but a vector of future values
(ϕ N +1 ,ϕ N + 2 , , ϕ N + H ) of the time series ϕt , (t = 1, 2, , N ) .
Compared with the direct strategy, which estimates
ϕ N + h , ( h = 1, 2, , H ) using H models, MIMO employs only
one multiple-output model, preserving the temporal stochastic
dependency hidden in the predicted time series.
The MIMO modeling strategy learns one multiple-output
forecasting model as follows:

=
(ϕi +1 , , ϕi + H ) f (ϕi , , ϕi − d +1 ) + ε

(1)

where d is the maximum embedding order and ε is a vector
noise term of zero mean and nondiagonal covariance.
After the learning process, the estimations of the H next
values are returned according to the following equation:
(2)
(ϕˆˆi +1 , , ϕi + H ) = fˆ (ϕi , , ϕi − d +1 )

2

The limitation of MIMO lies in forcing all the prediction
horizons to be predicted with the same model structure (for
instance, the same inputs (ϕi , , ϕi − d +1 ) ), resulting in a

reduction of the modeling effort, but likely to bias the returned
model results.
B. MISMO Strategy
A solution to the shortcomings of MIMO was proposed in [6,
7], where the constraint of the fixed modeling structure is
relaxed by tuning an integer parameter s , which calibrates the
dimensionality of the outputs on the basis of a validation
criterion. For a given s , MISMO generates n = H s
sub-models, which may have different inputs from each other
but a fixed size of outputs at s , resulting in slightly more
modeling flexibility.
The MISMO modeling strategy learns n multiple-output
sub-models as follows:
 , ϕi + ks −1 ) f k (ϕi , , ϕi − d +1 ) + ε, k ∈ {1, , n} . (3)
(ϕi + ( k −1) s ,=
where d is the maximum embedding order and ε is a zero-mean
vector of size s .
After the learning process, the estimations of the H next
values are returned by the following:
(4)
fˆk (ϕi , , ϕi − d +1 ), k ∈ {1, , n} .
(ϕˆˆ=
i + ( k −1) s ,  , ϕi + ks −1 )
The MISMO approach highlights the trade-off between the
property of preserving the stochastic dependency among the
predicted values and the flexibility of the modeling procedure
[8]. However, the increased flexibility of the MISMO strategy
requires an additional parameter s, which gives rise to several
issues in itself. An explicit issue is handling indivisibility when
computing n = H s , where all the parameters should be
integers. Furthermore, the number of outputs of every
sub-model derived from MISMO is fixed at s . This induces a
latent separation of dependencies among the steps that may
prevent the MISMO strategy from considering complex
dependencies between prediction steps within different
sub-models and consequently reduce the prediction accuracy.
Further study has been solicited to determine the divides of
prediction horizons adaptively while modeling rather than
fixing them at s .
III. FEED-FORWARD NEURAL NETWORKS
Although many types of neural networks (NNs) have been
proposed, the most popular one for time series prediction is the
feed-forward NNs (FNNs). Because the focus of this study is
the evaluation of the proposed modeling strategy relative to the
other competitors for multi-step-ahead prediction, a simple
three-layer FNN is implemented to model the selected
strategies, capable of creating a common and reliable
benchmark. In addition, the proposed modeling strategy and the
four well-established strategies can all be implemented in the
FNN because of the flexible structure of the FNN. Thus, the
standard three-layer FNN with nodes in adjacent layers fully
connected is used in this study. In this section, a short
introduction to feed-forward neural networks for
multi-step-ahead prediction is provided, and then the

> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) <
implementation of the FNN-based prediction model in the
current study is described.

H

MSE= ∑ (ϕˆi +h -ϕi +h )

2

H

3

(7)

h =1

A. FNN for Multi-Step-Ahead Prediction
Fig. 1 shows an example of a typical three-layer FNN model
with four input nodes, four hidden nodes (one hidden layer),
and two output neurons used for multi-step-ahead prediction.
The input nodes are the previous lagged observations while the
outputs provide the forecast for the future values in a
multi-step-ahead fashion. Hidden nodes with appropriate
nonlinear activation functions are used to process the
information received by the input nodes. Bias (or intercept)
terms are used both in the hidden and the output layers. To
facilitate the understanding of the FNN for multi-step-ahead
time series prediction, the FNN with the MIMO modeling
strategy will be used as an example in this section.

ϕt +1

ϕt +2

Output Layer

Hidden Layer

Input Layer

ϕt

ϕt -1

ϕt -3

ϕt -2

Fig. 1. A three-layer feed-forward neural network

=
D
Given the training sample

{( x , y ) ∈ (ℜ
i

i

d

}

× ℜH )

N

i =d

, the

FNN model for H -step-ahead prediction can be written as
follows:
 k

 n

(5)
y=  ∑ w0jhσ  ∑ wijr xi +b rj  +bh0  , h =1, …, H
 i =1

 j =1

where x is the input vector (ϕi , , ϕi − d +1 ) , y is the output
r
0
r
0
vector (ϕi +1 , , ϕi + H ) , w , w , b , and b are the weights

and biases of hidden and output layers, n and H are the sizes
of input and output vectors, k is the number of hidden units,
and σ ( ⋅) is a nonlinear activation function.

Functionally, equation (5) can also be written as follows:

(6)
(ϕi +1 , , ϕi + H ) =θf (ϕi , , ϕi − d +1 , ) +ε
where f ( ⋅) represents the FNN model and θ is a vector that
contains all the parameters in equation (5). The determination
of parameter set θ is based on the local error between the
measured (ϕi +1 , , ϕi + H ) and predicted values (ϕˆˆi +1 , , ϕi + H ) .

In the current study, the mean square error (MSE) is used as the
index of multi-step-ahead prediction performance to be
minimized, defined as follows:

B. Implementation of the FNN-based Prediction Model
Following the procedure from [9], the design tasks for the
FNN-based prediction model can be roughly divided into four
parts: data preprocessing, FNN designing, FNN
implementation and validation. The implementation of the
FNN-based model for multi-step-ahead prediction is presented
step-by-step below:
 Data preprocessing. Normalization is a standard
requirement for time series modeling and prediction. Thus,
the data sets are first scaled by linear transference to map
onto a range of [0, 1]. After the linear transference,
deseasonalization and detrending are performed when
necessary. Deseasonalization is performed with the most
recent X-12-ARIMA seasonal adjustment procedure.
Detrending is performed by fitting a polynomial time trend
to the data, and then subtracting the estimated trend from
the series when trends are detected by the Mann-Kendall
test.
 FNN designing. Selecting an appropriate architecture for
the FNN is normally the first step when designing an
FNN-based prediction system. The standard three-layer
FNN with fully connected nodes in adjacent layers is used
in this study. The number of input nodes, hidden nodes,
output nodes, and the type of activation functions are
defined in the present study as follows.
1) Selecting the number of input nodes: The number of
input nodes is determined by the input selection. The filter
method, which selects a set of inputs by optimizing a
criterion over different combinations of inputs by means of
a search algorithm, is employed for input selection in the
current study. The filter method requires the setting of two
elements: the criterion, i.e., the statistic which estimates the
quality of the selected variables, and the search algorithm,
which describes the policy used to explore the input space.
With respect to the criterion, the partial mutual information
[10] is used for the models using the iterated or direct
strategies, and an extension of the Delta test [6] is used for
the models with MIMO, MISMO, GA-MISMO, or
PSO-MISMO strategies. With respect to the search
algorithm, a forward-backward selection method that offers
the flexibility to reconsider input variables previously
discarded and to discard input variables previously selected
is used. The maximum embedding order, d , is set to 15.
2) Selecting the number of hidden nodes: The number of the
hidden nodes cannot be determined in advance. Thus
empirical experimentations are needed to determine this
number [11]. Because of the small sample size of many of
the series, the experimentation is limited to five possible
values of hidden nodes, i.e., 2, 4, 6, 8, or 10. The best
number of hidden nodes is determined by using the original
Akaike’s information criterion (AIC).
3) Selecting the number of output nodes: The number of
output nodes is determined by the modeling strategy. Most
of the modeling strategies found in the literature can be
classified into two categories, single-output structure and
multiple-output structure. In the single-output structure

> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) <
strategies, such as the iterated and direct strategies, multiple
inputs map to a single output. The MIMO, MISMO,
GA-MISMO, and PSO-MISMO strategies are all based on
the multiple-output structure that maps multiple inputs to
multiple outputs. Therefore, for the iterated and direct
strategies, only one output node is used for the FNN model.
When employing the MIMO strategy, the number of output
nodes is equal to the number of prediction horizons. In the
MISMO strategy, the original H -step-ahead prediction
tasks are separated into n = H s subtasks, resulting in n
sub-models with fixed multiple outputs, s , for each, where
1 ≤ s ≤ H . Unlike the MISMO strategy, which uses a fixed
number of outputs s for each sub-model, the GA-MISMO
and PSO-MISMO strategies self-adaptively determines the
number of sub-models with varying number of outputs, as
described in Section V.
4) Selecting the type of activation functions: A sigmoid
transfer function is used at each hidden node and a linear
transfer function is used at each output node.
 FNN implementation. The Levenberg and Marquardt
algorithm (LMA) provided by the MATLAB (Version
2009b) NN toolkit is used for the training. Node biases are
applied at the hidden and output layers. For the stopping
criterion, the number of learning epochs is chosen as 1000,
having no prior knowledge of the appropriate value. To
determine the optimal parameters for the FNN, the common
practice of fivefold cross-validation is used.
 Validation: For each series of the simulated and real
dataset, the last 18 observations are separated for ex-ante
performance assessment (out-of-sample), statistically
independent from the parameter estimation process. The
remainder of the data (in-sample) is used for parameter
estimation and the FNN modeling.
IV. BINARY PARTICLE SWARM OPTIMIZATION
Particle swarm optimization (PSO) [12] is a
population-based self-adaptive search algorithm that exploits a
population of individuals to probe promising regions of the
search space. In the PSO, the population is referred to as a
swarm that consists of a number of particles. Each particle
represents a potential solution of the optimization task and has a
position represented by a vector. Each particle moves to a new
position according to both local information (the particle
memory) and global information (the knowledge of all the
particles). Thus, the PSO has the ability to converge to local
and/or global optimal solutions in a small number of
generations. In this section, a brief introduction to a variant of
the original PSO, i.e., binary PSO, is presented, to understand
the mechanism of the proposed PSO-MISMO strategy applied
in this study.
The PSO was initially designed to solve continuous
optimization problems. Since the introduction of PSO,
successful applications to several optimization problems [13-17]
have demonstrated its potential. However, the major obstacle of
applying a PSO successfully is its continuous nature. The first
variant of the PSO algorithm for solving problems with
binary-value solution elements (binary PSO) [18] was also
developed by the creator of the original PSO. There are two key

4

differences between the binary PSO and the original PSO. The
first difference is the representation of the particle. In the binary
PSO, every particle is characterized by a binary solution
representation (Note that the representation is primarily related
to a coding issue in the proposed PSO-MISMO strategy). The
other difference is that the velocity of a particle in the binary
PSO is a probability vector, where each probability element
determines the likelihood of that binary variable having the
value of one. Fig. 2 lists the pseudocode algorithm for the basic
binary PSO.
Algorithm 1 Binary particle swarm optimization
Initialize a population of particles with random positions and
velocities, throughout the input space.
While not sufficiently good performance or maximum number
of iterations do
for each particle i in {1, , I } do
if f ( Pi ) < f ( pbesti ) then

pbesti = Pi
end if
end
pbestmin =

=
f ( pbest )
{ pbest
i*

i*

}

min ( f ( pbesti ) , i ∈ {1, , I } )

if f ( pbestmin ) < f ( gbest ) then

gbest = pbestmin
end if
for each particle i in {1, , I } do
Update velocity Vi according to equation (8).
Update position Pi according to equation (11).
end
end while
Fig. 2. Pseudocode for the binary PSO algorithm

The binary PSO updates the velocity according to equation
(8) as follows:
Vi t = w ⋅ Vi t −1 + c1 ⋅ rand1 ⋅ ( pbestit −1 − Pi t −1 )
(8)
+ c2 ⋅ rand 2 ⋅ ( gbestit −1 − Pi t −1 )
where c1 and c2 are the cognitive and interaction coefficients,
rand1 and rand 2 are random real numbers uniformly
distributed between 0 and 1. w , the inertia weight, a
user-specified parameter controls the momentum of the particle.
A larger inertia weight pushes towards global exploration while
a smaller inertia weight helps in fine-tuning the current search
area. The following weighting function is usually utilized:
T −t
(9)
+ wmin
w = ( wmax − wmin ) ×
t
where wmax and wmin are initial weight and final weight,
respectively, T is the maximum number of allowable iterations
and t is the current iteration.
The velocity is constrained to the interval [0, 1] by using the
following sigmoid transformation:
1
(10)
S (Vi t ) =
1 − exp ( −Vi t )

> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) <
where

, is the probability that bit

avoid

approaching 0 or 1, a constant

4 such that

is equal to 1. To
is often set at

Each particle, at each time step, changes its current position
according to equation (10), based on equation (12), as follows:
(12)

, then the following relationships are

applied:

5

where
is a random real number uniformly distributed
between 0 and 1.
(11)

V. PROPOSED PSO-MISMO MODELING STRATEGY

Fig. 3. Flowchart of the proposed PSO-MISMO modeling strategy

To address the limitations of MISMO mentioned first in
Section I and illustrated further in Section II.B, an improved
binary PSO-based MISMO is developed that self-adaptively
determines the number of sub-models with varying prediction
horizons. The proposed strategy is abbreviated PSO-MISMO.
Fig. 3 depicts the flowchart of the proposed PSO-MISMO
modeling strategy. As shown in Fig. 3, the detailed
PSO-MISMO algorithm consists of three major operations:
initialization, evaluation, and update. Before describing details
of these operations, coding issues will be addressed. The
overall learning process in Fig. 3 is elaborated step by step
below.
composed of
 Coding. Consider a time series
observations for which the next
observations are to be
predicted. To avoid the restriction associated with adopting a
fixed in the original MISMO method, a flexible vector
is used, where

is the number of prediction

tasks, and

represents the number of consecutive outputs to

be predicted at a time by the
th

model (Note that

model refers to the j sub model with
of

and

outputs). The values

are determined by the binary vector
that represents the segmentation point (see

Algorithm 2). Each particle

, representing a

possible solution to the optimization problem, is determined
by the binary PSO algorithm. The values of
components of the binary vector
are either
0 or 1; if the value of a variable is 0, then the original task will
not be divided by this segmentation point. If the value of a
variable is 1, then the prediction task will be divided by this
segmentation point. Fig. 4 illustrates the solution
representation.

> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) <
Algorithm 2. Computing the value of S = ( s1 , , sJ )

Prediction horizon

1

2

p2

p1



3

p3



6

H −1

pH − 2

H

pH −1

Input: H , prediction horizon.
Input: P = ( p1 , , pH −1 ) , binary vector.
Output: S = ( s1 , , sJ ) , s j represents the number of
consecutive outputs to be predicted at a time by the FNN j
model
H −1

J = ∑ pi + 1

Segmentation point

i =1

Fig. 4. Solution representation

H −1

1,=
i 1, 2, , H − 1} , the strategy
In the special case that { p=
i

pi 0,=
i 1, 2, , H − 1} ,
equals the direct strategy, and when { =
it equals the MIMO strategy.
For further illustration, assume a 10-steps-ahead prediction
task ( H =10 ). Given the binary vector P = ( 0,0,1,0,0,0,0,1,0 ) ,

N = ∑ pi
i =1

SP is the vector of size N
n =1

for i in {1, , H − 1} do
if pi == 1 then

then S = ( 3,5, 2 ) and J = 3 (see Algorithm 2). Thus, the
original task is separated into three prediction subtasks as
follows:
=
(ϕi +1 , ϕi + 2 , ϕi +3 ) f1 (ϕi , ϕi −1 , , ϕi − d +1 ) + ε1

ϕi + 5 , ϕi + 6 , ϕi + 7 , ϕi +8 ) f 2 (ϕi , ϕi −1 , , ϕi − d +1 ) + ε 2
(ϕi + 4 ,=
=
(ϕi + 9 , ϕi +10 ) f3 (ϕi , ϕi −1 , , ϕi − d +1 ) + ε3
where {ε i , i = 1, 2,3} is a zero-mean vector of size 3, 5, 2.
After the learning process, the estimation of the ten next
values is returned by the three steps as follows:
ˆ
(ϕˆˆˆ
i +1 , ϕi + 2 , ϕi + 3 ) = f1 (ϕi , ϕi −1 ,  , ϕi − d +1 )
ˆ
(ϕˆˆˆˆ,̂
i + 4 ϕi + 5 , ϕi + 6 , ϕi + 7 , ϕi + 8 ) = f 2 (ϕi , ϕi −1 ,  , ϕi − d +1 )
(ϕˆˆi + 9 , ϕi +10 ) = fˆ3 (ϕi , ϕi −1 , , ϕi − d +1 )

Prediction horizon H =10

1

2

0

3

0

4

1

5

0

6

0

7

0

8

0

9

1

n =n + 1
end if
end

s1 =SP [1]

sJ =H − SP [ N ]
for j in {2, , J − 1} do

s j =SP [ j ] − SP [ j − 1]
end
Return S = ( s1 , , sJ )
Fig. 6. Pseudocode for computing the value of S = ( s1 , , sJ )

The prediction horizons of these prediction sub-models (or
the number of outputs of the three FNN models) are 3, 5, and 2,
respectively. This difference, or heterogeneity, is the highlight
of the proposed strategy. For clarity, the flow of this example is
illustrated in Fig. 5.

Original task

SP [ n ] =i

10

0

 Initialization. All the particles are randomly generated with 0
or 1 using a 50% probability. All the velocity components are
assigned the initial value of 0. Concerning the selection of
parameters (i.e., inertial weight w , coefficients c1 and c2 ,
swarm size I , and number of iterations T ) in binary PSO, it is
yet another challenging model selection task [19]. Fortunately,
several empirical and theoretical studies have been performed
about the parameters of PSO from which valuable information
can be obtained [20-23]. In this study, the parameters are
determined according to the recommendations in these studies
and selected based on the prediction performance and
computational time in a trial-error fashion. Table I summaries
the final parameters.
TABLE I
PARAMETER SELECTION OF BINARY PSO
Parameters
Values
20
Swarm size I
100
Number of iterations T
2.0
Cognitive coefficient c

Segmentation point

1

Prediction horizon H = 5 Prediction horizon H = 2

Prediction horizon H = 3

Subtasks

1

2

Subtask І

3

4

5

6

Subtask Π

Fig. 5. An example of solution representation

7

8

9

10

Subtask Ш

Interaction coefficient c2

2.0

Initial weight wmax

0.9

Final weight wmin

0.4

()
 Evaluation. In the current generation t , for each particle pi ,
t

the

position

P = ( p1 , , pH −1 )

is

transformed

to

> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) <

S = ( s1 , , sJ ) (see Algorithm 2). Then, each s j ∈ {s1 , , sJ }
()
derived from the particle pi , is trained using the training
t

FNN ij model and the dataset
=
Dj
where: x j ⊂ {ϕ n -d +1 , , ϕ n } ,

{( x , y ) ∈ (ℜ
j

j

m

×ℜ

sj

)}

N

,

s0 =0 .
After the training process, the estimations of the next values
of s j are returned. Once all H next values for particle

{

pi( ) ∈ p1( ) , , pI( )
t

t

t

pi( )
t

particle

} are obtained, the fitness value of each
is

calculated

MSE i( ) = ∑ (ϕˆ Ni +h -ϕ Ni +h )
t

H

2

according

to

H , where ϕˆ Ni +h is the h- step-ahead

h =1

()
forecast for particle pi , ϕ Ni +h is the true value for particle
t

pi(t ) , and H is the prediction horizon.
 Update. In the current generation, once the fitness values of
all the particles in the swarm are calculated, the pbest of each
particle, and the gbest of the swarm are obtained. Then, all the
particles’ velocities and positions are updated according to
equations (8) and (12), respectively. If the termination
conditions are met, then output gbest (that is, P* ) and the
corresponding H - step-ahead forecast

{ϕˆˆ

*
N +1

, ,ϕ N* +H } for

particle P* ; otherwise, go back to the previous step.
VI. EXPERIMENTAL SETTING
A. Datasets Description
To evaluate the performances of the proposed PSO-MISMO
and the counterparts in terms of the forecast accuracy, two
simulated time series, i.e., Logistic and Mackey-Glass time
series, and a real world dataset, i.e., NN3 competition dataset,
are used in this present study.
Logistic and Mackey-Glass time series are recognized as
benchmark time series that have been commonly used and
reported by a number of studies related to time series modeling
and forecasting [24-26].
The Logistic map is a demographic model that was
popularized by May [27] as an example of a simple non-linear
system that exhibits complex, chaotic behavior. It is drawn
from equation (13).

=
ϕi 4ϕi −1 (1 − ϕi −1 )

Systems Toolbox 1 from the MATLAB software.
TABLE II
INITIALIZATION AND SAMPLE SIZE OF THE SIMULATED TIME
SERIES
No.
DGP
Sample size

i =d





y j = ϕ
, , ϕ j  ,
j
s( k -1)
n + ∑ sk
 n +1+ ∑

k =1
k =1

(13)

The Mackey-Glass time series is approximated from the
differential equation (14) (see [28]).
dϕ i
0.2ϕi −τ
(14)
=
− 0.1ϕi
di 1 + ϕi −τ 10
For each data-generating process (DGP), i.e., Logistic and
Mackey-Glass process, we simulate twenty time series with
different initialization and sample size, as was shown in Table
II. The data for these time series are generated by the Chaotic

7

( ϕ1 )

Logistic
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20

0.100
0.125
0.150
0.175
0.200
0.225
0.260
0.275
0.300
0.325
0.350
0.375
0.400
0.425
0.450
0.475
0.510
0.525
0.550
0.575

Mackey-Glass ( ϕ1 , τ )
1.000
1.200
1.400
1.600
1.800
2.000
1.000
1.200
1.400
1.600
1.800
2.000
1.000
1.200
1.400
1.600
1.800
2.000
1.000
1.200

15
15
15
15
15
15
16
16
16
16
16
16
17
17
17
17
17
17
18
18

485
496
523
548
674
692
726
758
779
791
821
843
869
889
912
926
946
964
987
1002

The NN3 competition was organized in 2007, targeting
computational-intelligence forecasting approaches. The
competition dataset of 111 monthly time series drawn from
homogeneous population of real business time series is used for
evaluation 2.
As such, three datasets of 20 Logistic time series, 20
Mackey-Glass time series, and 111 NN3 time series are used
for evaluating the performances of the proposed PSO-MISMO
and the counterparts in this study. Each series is split into an
estimation sample and a hold-out sample. The last 18
observations are saved for evaluating and comparing the
out-of-sample forecast performances of the various
multi-step-ahead prediction strategies. All performance
comparisons are based on these 18 × 20 out-of-sample points
for Logistic and Mackey-Glass datasets and 18 × 111
out-of-sample points for NN3 datasets.
B. Accuracy Measure
For each prediction horizon h , three alternative forecast
accuracy measures are considered: the mean absolute
percentage error (MAPE), the symmetric mean absolute
percentage error (SMAPE), and the mean absolute scaled error
(MASE). The definitions of each are as follows:
1 M ϕ m -ϕˆ m
(15)
MAPE h = ∑ N +h m N +h ∗100
M m =1 ϕ N +h
SMAPE h =

1
2

1
M

ϕ Nm+h -ϕˆ Nm+h
∗100
∑
m
ˆ Nm+h
m =1 ϕ N +h +ϕ
M

(16)

Available at http://www.mathworks.com/matlabcentral/fileexchange/1597
Available at http://www.neural-forecasting-competition.com/NN3/datasets.htm

> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) <

8

PSO-MISMO, as was presented in Section V.
MASE h =

1
M

M

∑

m =1

ϕ

m
N +h
N

-ϕˆ

m
N +h

1
∑ ϕim -ϕim-1
N -1 i = 2

Start

(17)

where ϕˆ Nm+h is the h- step-ahead forecast for time series m ,

ϕ Nm+h is the true time series value for series m , H is the
prediction horizon (in this case, H =18 ) and M is the number
of time series in the datasets (in this case, M = 20 for Logistic
and Mackey-Glass datasets and M =111 for NN3 datasets). For
preprocessed data, including normalization, deseasonalization
and detrending, we convert the outputs from the neural
networks back to their original scales. Performance prediction
can then be compared directly.
C. Genetic Algorithm as Alternative
To further assess the performance of the proposed
PSO-MISMO modeling strategy for multi-step-ahead
prediction, we compare the experimental results of the
PSO-MISMO with those produced by other discrete
evolutionary algorithms, such as genetic algorithm (GA) (for
detailed introduction to GA, please refer to [29-31]). This
subsection briefly presents the implementation of GA-MISMO,
a straightforward alternative against PSO-MISMO.
Fig. 7 depicts the flowchart of the GA-MISMO modeling
strategy. As shown in Fig. 7, the detailed GA-MISMO
algorithm consists of three major parts: initialization,
evaluation, and operatation. Before describing details of these
parts, coding issues will be addressed. The overall learning
process in Fig. 7 is elaborated step by step below.
 Coding. The chromosome design in GA-MISMO is the same
with the particle design in PSO-MISMO, as was presented in

Section V. Actually, the individual P = ( p1 , , pH −1 ) may be
regarded as chromosome in terms of GA-MISMO, or as
particle in terms of PSO-MISMO. The values of H − 1
components of the binary vector P = ( p1 , , pH −1 ) are either 0

or 1: if the value of a variable is 0, then the original task will not
be divided by this segmentation point; if the value of a variable
is 1, then the prediction task will be divided by this
segmentation point. The coding issue in GA-MISMO is the
same as PSO-MISMO and will not be presented in details here
to save space.
TABLE III
PARAMETER SELECTION OF GA
Parameters
Values
Population size
20
Number of iterations
100
Crossover probability
0.90
Mutation probability
0.02

 Initialization. All the chromosomes are randomly generated
with 0 or 1 using a 50% probability. In setting the parameters of
GA, it is required to do extensive simulations to find suitable
values for various parameters. In this study, GA’s parameters
are determined through preliminary simulation and selected
according to the recommendations in [32-34] in a trial-error
fashion. Table III summarizes the final parameter of GA.
 Evaluation. The evaluation in GA-MISMO is the same as

Step 1：
Initialization

Set population size I, maximum number of
generations T, and t=1
Distribute chromosomes randomly across the
design space
The termination
conditions are met?

Yes

No

Step 2：
Evaluation

Step 3：
Operatation

For each chromosome Pi (t ) in the population,
calculate the fitness value
Apply GA operators:
selection, crossover, and mutation
t=t+1

(

Output P* and ϕˆˆN* +1 , ,ϕ N* +H

)

End
Fig. 7. Flowchart of the GA-MISMO modeling strategy

 Operatation. Designing the operators for selection, crossover,
and mutation is a major issue in a GA implementation and is
always done ad hoc. In each generation, after the fitness values
of all the chromosomes in the same population are calculated,
the selection operator that chooses chromosomes from the
current generation’s population for inclusion in the next
generation’s population is executed. There are many ways for
choosing parents [35]. In this study, the roulette method is used
for Logistic time series and top percent method is adopted for
Mackey-Glass and NN3 time series through preliminary
simulation in a trial-error fashion. After selection of parents, the
crossover operator should be applied to produce two offspring
by exchanging some genetic information between the parents.
Crossover occurs during evolution according to the crossover
probability. The crossover probability is how often a crossover
will be performed. In this study, the probability of crossover is
set to 0.9, as shown in Table III; two-point operator is used.
This process of selection and reproduction is repeated until the
number of offspring becomes equal to the number of eliminated
chromosomes; by adding the offspring to population, its size
becomes equal to the initial size [34]. The mutation is the third
operator, which is applied on the population (parents and their
offspring) excluding the best chromosome. A simple mutation
operator that selects a few percent of the genes of the
population randomly and flips their values from zero to one and
vice versa is used in this study [36]. In this study, the
probability of mutation is set to 0.02, as shown in Table III. The
above process is repeated until a termination conditions are
reached. Once the termination conditions are met, the best
chromosome P* and the corresponding H - step-ahead forecast

> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) <

{ϕˆˆ

*
N +1

, ,ϕ N* +H } for chromosome P* are obtained; otherwise,

go back to the previous step.

Start
m=1
Input mth time series
Split series into estimation sample and
hold-out sample

Conduct input selection
Perform Levenburg-Marquardt training algorithm and
five-fold cross validation for model selection of FNN
Conduct multi-step-ahead prediction using

Iterated
strategy

MIMO
strategy

Direct
strategy

No

MISMO
strategy

GA-MISMO
strategy

PSO-MISMO
strategy

m>M
Yes

Calculate MAPE h , SMAPE h, and MASE h for each prediction horizon and dataset

No

identify the significantly different prediction models in
multiple pair-wise comparisons.
VII. RESULTS AND DISCUSSION

D. Experimental Procedure

m=m+1

9

Repeated 10
times?

Yes
Calculate the mean of each accuracy measure

Compare the quality of obtained results
Perform ANOVA and HSD test
End

Fig. 8. Experiment procedure

Fig. 8 shows the experimental procedure using the simulated
and real time series. Each series is first split into the estimation
sample and the hold-out sample. Then, the input selection and
model selection for each series are determined using the filter
method, the Levenberg and Marquardt algorithm, and fivefold
cross-validation with iterated, direct, MIMO, MISMO,
GA-MISMO, and PSO-MISMO strategies. Finally the models
are tested with the hold-out samples and the MAPE h ,

SMAPE h and MASE h are computed for each prediction
horizon h (in this case h =1,2, … ,18 ) over datasets (i.e.,
Logistic, Mackey-Glass, and NN3 datasets). The modeling
process for each series is repeated ten times. Upon the
termination of this loop, performance of the examined models
with selected strategies at each prediction horizon and dataset is
judged in terms of the mean, averaged by ten, of the MAPE h ,

SMAPE h and MASE h . Analysis of variance (ANOVA) test
procedures are used to determine if the means of performance
measures are significantly different among the three models for
each prediction horizon and datasets. If so, Tukey’s honesty
significant difference (HSD) tests [37] are then used to further

A. PSO-MISMO vs. the Others in Prediction Accuracy
The prediction performances of all six modeling strategies
examined in terms of the three accuracy measures (MAPE,
SMAPE, and MASE) and average rank for three datasets are
shown in Table IV. The columns labeled as ‘Average 1- h ’
show the average accuracy measures over the prediction
horizon 1 to h . The last column shows the average ranking for
each model over all prediction horizons of the out-of-sample
prediction performance. For each column of Table IV, the entry
with the smallest value is set in boldface and marked with an
asterisk, and the entry with second smallest value is set in
boldface type.
The results in Table IV lead to the following conclusions:
 The rankings from best to worst are: PSO-MISMO and
GA-MISMO are almost a tie, then MISMO, MIMO, direct, and
iterated strategies, regardless of the accuracy measures
considered. Thus our findings are robust to the choice of
accuracy measures.
 The modeling strategies based on the multiple-output
structure (i.e., PSO-MISMO, GA-MISMO, MISMO, and
MIMO strategies) outperform those based on the single-output
structure (i.e., iterated and direct strategies), which is in
agreement with [7].
 Despite the popularity of the iterated strategy in
multi-step-ahead prediction literature, it is consistently the
worst performing strategy for prediction.
 Comparing the iterated strategy with the direct strategy, the
direct strategy is better regardless of the accuracy measures
considered, which demonstrates that the accumulation of errors
in the case of iterated strategy drastically deteriorates the
accuracy of the prediction.
 The MIMO strategy consistently achieves more accurate
forecasts than either the iterated or the direct strategies for all
prediction horizons. It is conceivable that the reason for the
superiority of the MIMO strategy is that it preserves, among the
predicted values, the stochastic dependency characterizing the
time series.
 The MISMO strategy seems to produce forecasts which are
more accurate than those of the MIMO strategy (though only
marginally). It is conceivable that the reason for the superiority
of the MISMO strategy is that it trades off the property of
preserving the stochastic dependency between future values
with a greater flexibility of the predictor.
 When comparing the results of two heuristic-based
modeling strategies (i.e., PSO-MISMO and GA-MISMO) with
those of MISMO strategy, the heuristic-based strategies are
generally better, which may indicate that the fixed number of
outputs of every sub-model derived from the MISMO strategy
tends to prevent the MISMO strategy from considering
complex dependencies between prediction steps within
different sub models and consequently reduces the prediction
accuracy.

> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) <

10

TABLE IV
PREDICTION ACCURACY MEASURES FOR HOLD-OUT SAMPLE
Dataset

Strategy

Average Rank
Prediction horizon( h )
Average 1- h
1
2
3
6
12
18
1-3
1-12
1-18
MAPE
ITER
39.94
39.41
43.28
35.21
51.41
42.41
43.18
46.18
5.888
47.89
DIR
47.93
37.18
36.27
41.71
36.14
49.22
40.46
42.95
44.87
5.000
MIMO
47.96
38.42
37.81
40.18
34.81
47.81
41.39
42.54
42.18
4.000
MISMO
47.96
37.64
35.54
41.21
34.15
45.12
40.38
41.54
41.84
3.111
GA-MISMO
47.88*
37.04*
35.28
39.84
31.01*
43.01
40.06*
39.24
39.14
1.555
PSO-MISMO
47.95
37.17
35.17*
38.28*
32.73
42.18*
40.09
38.91*
38.51*
1.444*
SMAPE
ITER
31.10
33.42
36.74
39.74
45.12
39.21
33.75
39.13
40.14
5.888
DIR
32.02
37.18
36.18
45.85
37.44
33.39
38.67
38.79
5.055
30.97
Logistic
MIMO
31.18
32.59
36.88
34.79
44.14
35.42
33.55
38.04
38.17
3.944
MISMO
36.31
34.82
42.74
34.18
33.05
37.87
37.49
3.166
30.87*
31.98
GA-MISMO
31.09
31.86*
36.14
32.97
40.81*
32.95
33.03
35.41*
35.82
1.555
PSO-MISMO
31.11
32.03
35.74*
33.58*
41.58
32.07*
32.96*
36.13
35.18*
1.500*
MASE
ITER
0.706
0.521
0.911
0.914
1.184
1.054
0.713
0.941
1.124
5.944
DIR
0.705
0.894
0.881
0.914
0.867
0.700
0.887
0.914
4.944
0.501
MIMO
0.696
0.539
0.814
0.784
0.846
0.814
0.683
0.761
0.824
4.111
MISMO
0.541
0.795
0.718
0.819
0.748
0.675
0.745
0.794
3.000
0.691*
GA-MISMO
0.705
0.505
0.731*
0.624*
0.718*
0.642
0.647
0.687*
0.678*
1.444*
PSO-MISMO
0.645*
0.695
0.499*
0.742
0.648
0.743
0.618*
0.702
0.704
1.555
MAPE
ITER
5.15
5.26
6.34
7.15
10.18
9.46
5.58
9.01
9.22
5.944
DIR
4.93
5.87
6.92
8.55
8.49
5.30
7.83
8.13
5.055
5.11*
MIMO
5.16
5.02
5.19
6.31
6.24
6.87
5.12
6.37
6.38
3.777
MISMO
5.23
5.87
5.97
6.05
5.08
5.89
5.64
3.222
5.13
4.88
GA-MISMO
5.16
4.87*
4.61
4.87*
4.28*
4.22*
4.88
4.78*
4.71*
1.444*
PSO-MISMO
5.15
4.91
4.50*
5.06
4.58
4.31
4.85*
4.96
4.87
1.555
SMAPE
ITER
5.31
5.48
8.10
10.18
11.18
5.35
8.91
10.52
5.555
5.28*
DIR
5.30
5.11
7.29
8.46
9.85
5.05
7.84
9.86
5.000
4.75
5.25
5.62
6.57
6.87
7.46
5.38
6.05
8.05
4.055
Mackey-Glass MIMO
5.29
MISMO
5.31
5.32
6.18
5.94
6.28
5.22
5.87
6.91
3.111
5.05*
GA-MISMO
5.10
5.29
4.75
5.40
4.91
4.31*
5.04
5.16
5.19
1.611*
PSO-MISMO
5.31
5.07
4.64*
5.23*
4.72*
4.45
5.01*
5.12*
5.03*
1.777
MASE
ITER
1.785
1.683
1.628
1.935
2.148
2.184
1.698
1.924
2.015
5.722
DIR
1.783
1.531
1.910
1.924
1.924
1.662
1.875
1.987
4.833
1.672
MIMO
1.694
1.654
1.875
1.816
1.794
1.705
1.846
1.914
4.000
1.769*
MISMO
1.681
1.615
1.849
1.748
1.649
1.690
1.789
1.841
3.055
1.775
GA-MISMO
1.781
1.663*
1.522
1.746*
1.601
1.478
1.655*
1.710
1.647
1.666
PSO-MISMO
1.777
1.675
1.520*
1.759
1.588*
1.462*
1.657
1.703*
1.638*
1.611*
MAPE
ITER
14.93
13.69
17.05
19.76
18.48
29.49
15.10
22.38
23.27
5.944
DIR
12.23
17.04
18.87
17.08
23.84
14.71
17.29
18.74
5.000
14.87*
MIMO
11.68
17.09
18.01
16.07
21.54
14.55
16.41
17.57
4.055
14.88
MISMO
17.27
17.89
15.96
21.31
14.56
15.95
16.12
2.944
14.88
11.62
GA-MISMO
11.63
14.88
16.67*
17.20
15.04*
20.34
14.39*
15.68
16.03
1.611
PSO-MISMO
14.89
11.54*
16.84
17.15*
15.15
20.22*
14.42
15.61*
15.91*
1.444*
SMAPE
ITER
15.16
13.31
14.97
20.57
17.84
20.68
14.48
19.97
21.59
5.944
DIR
11.84
14.34
18.37
15.67
19.97
13.76
17.75
19.84
5.000
15.10
NN3
MIMO
15.12
11.17
14.22
16.37
14.83
18.34
13.50
16.58
18.34
4.055
MISMO
10.97
15.97
14.56
18.16
13.14
16.13
17.10
2.944
15.10
13.37
GA-MISMO
15.09*
10.88
13.21*
15.86
13.43*
17.24*
13.06*
15.61*
16.47*
1.500*
PSO-MISMO
15.11
13.46
10.68*
15.91*
13.61
17.38
13.08
15.72
16.51
1.555
MASE
ITER
0.835
0.845
1.104
1.297
0.184
1.418
0.928
1.164
1.187
5.944
DIR
0.859
1.197
1.168
0.108
1.387
0.962
0.109
1.046
5.055
0.829*
MIMO
0.746
1.019
1.035
0.981
1.197
0.866
0.943
0.981
3.888
0.834
MISMO
0.996
1.007
0.975
1.184
0.857
0.923
0.934
3.055
0.834
0.741
GA-MISMO
0.835
0.731*
0.987
0.912
0.921*
1.083
0.856
0.910
0.917
1.555
PSO-MISMO
0.742
0.834
0.974*
0.907*
0.926
1.061*
0.850*
0.903*
0.911*
1.500*
Note: For each column of table, the entry with the smallest value is set in boldface and marked with an asterisk, and the entry with second smallest value is set in
boldface type.

> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) <

dataset

Measure

MAPE h

Logistic

SMAPE h

MASE h

MAPE h

Mackeyglass

SMAPE h

MASE h

MAPE h

SMAPE h
NN3

MASE h

Note:

TABLE V
MULTIPLE COMPARISON RESULTS WITH RANKED STRATEGIES FOR HOLD-OUT SAMPLE
Prediction
Rank of strategies
1
2
3
4
horizon ( h )
4, 5, 8, 13, 16, 17
PSO-MISMO
GA-MISMO
MISMO
MIMO
<
<
<*
<
<
7, 9-11, 15
GA-MISMO
PSO-MISMO
MISMO
MIMO
<*
<
18
PSO-MISMO
GA-MISMO
MISMO
MIMO
<*
<*
<
<
12, 14
GA-MISMO
PSO-MISMO
MISMO
MIMO
<*
<
<
<
6
PSO-MISMO
GA-MISMO
MIMO
MISMO
<
<
<
2
GA-MISMO
PSO-MISMO
MISMO
DIR
<
<
<
3
PSO-MISMO
GA-MISMO
MISMO
DIR
<
<
6, 11, 14, 15, 18
PSO-MISMO
GA-MISMO
MISMO
MIMO
<*
<
<
5, 8, 13, 16
GA-MISMO
PSO-MISMO
MISMO
MIMO
<*
<
7
PSO-MISMO
GA-MISMO
MISMO
MIMO
<*
<*
<
<
<
10
PSO-MISMO
GA-MISMO
MIMO
MISMO
<
12, 17
GA-MISMO
PSO-MISMO
MISMO
MIMO
<*
<*
<
<
<
4, 9
PSO-MISMO
GA-MISMO
MISMO
DIR
<
<
5, 10, 14, 15
PSO-MISMO
GA-MISMO
MISMO
MIMO
<*
<
<
7, 8, 11, 16, 17
GA-MISMO
PSO-MISMO
MISMO
MIMO
<*
<
9, 13
PSO-MISMO
GA-MISMO
MISMO
DIR
<*
<*
<
<
<
4, 12
GA-MISMO
PSO-MISMO
MISMO
MIMO
<
<
18
PSO-MISMO
GA-MISMO
MISMO
MIMO
<*
<
<
<
3
GA-MISMO
PSO-MISMO
MISMO
MIMO
<
<
6
PSO-MISMO
GA-MISMO
MISMO
MIMO
<*
<
<
5, 7, 10, 15
PSO-MISMO
GA-MISMO
MISMO
MIMO
<*
<
<
8, 9, 13, 16
GA-MISMO
PSO-MISMO
MISMO
MIMO
<*
<
4, 11, 14
PSO-MISMO
GA-MISMO
MISMO
MIMO
<*
<*
<
<
<
12, 17
GA-MISMO
PSO-MISMO
MISMO
MIMO
<
<
18
GA-MISMO
PSO-MISMO
MIMO
MISMO
<*
<
<
6
GA-MISMO
PSO-MISMO
MISMO
MIMO
<*
<
<
3
PSO-MISMO
GA-MISMO
MIMO
MISMO
<*
<
<
6, 9, 11, 13, 15
PSO-MISMO
GA-MISMO
MISMO
MIMO
<*
<
<
4, 7, 8, 14, 16, 17
GA-MISMO
PSO-MISMO
MISMO
MIMO
<*
<
12, 10
PSO-MISMO
GA-MISMO
MISMO
MIMO
<*
<*
<
18
GA-MISMO
PSO-MISMO
MISMO
MIMO
<*
<*
<
<
5
GA-MISMO
PSO-MISMO
MISMO
MIMO
<*
<
<
<
3
PSO-MISMO
GA-MISMO
DIR
MISMO
<
<
7, 10, 11, 14, 16
PSO-MISMO
GA-MISMO
MISMO
MIMO
<*
<
<
5, 8, 13
GA-MISMO
PSO-MISMO
MISMO
MIMO
<*
<
18, 15
PSO-MISMO
GA-MISMO
MISMO
MIMO
<*
<*
<
<
<
4, 17
GA-MISMO
PSO-MISMO
MISMO
MIMO
<
<
6, 9
GA-MISMO
PSO-MISMO
MISMO
MIMO
<*
<
3
PSO-MISMO
GA-MISMO
DIR
MISMO
<*
<*
<
<
12
PSO-MISMO
GA-MISMO
MISMO
MIMO
<*
<
<
4, 6, 13,15,18
PSO-MISMO
GA-MISMO
MISMO
MIMO
<*
<
<
5, 8, 10-12
GA-MISMO
PSO-MISMO
MISMO
MIMO
<*
<
<
7, 16, 17
PSO-MISMO
GA-MISMO
MISMO
MIMO
<*
<
<
9, 14
GA-MISMO
PSO-MISMO
MISMO
MIMO
<*
<
<
<
2,
PSO-MISMO
MISMO
GA-MISMO
MIMO
<
<
5, 7-9, 14,17
PSO-MISMO
GA-MISMO
MISMO
MIMO
<*
<
<
4, 10, 11, 12, 16
GA-MISMO
PSO-MISMO
MISMO
MIMO
<*
<
<
<
6
PSO-MISMO
GA-MISMO
MISMO
MIMO
<
<
<
13
PSO-MISMO
MISMO
GA-MISMO
MIMO
<
<
3, 15
GA-MISMO
PSO-MISMO
MISMO
MIMO
<*
<
<
<
18
GA-MISMO
PSO-MISMO
MISMO
MIMO
<
<
2
PSO-MISMO
GA-MISMO
MISMO
MIMO
<
<
<
5, 7, 8, 11
PSO-MISMO
GA-MISMO
MISMO
MIMO
<*
<
<
10, 13, 14, 16
GA-MISMO
PSO-MISMO
MISMO
MIMO
<*
<
<
9, 15, 18
PSO-MISMO
GA-MISMO
MISMO
MIMO
<*
<
<
4, 12, 17
GA-MISMO
PSO-MISMO
MISMO
MIMO
<*
<
<
<
2
GA-MISMO
MISMO
PSO-MISMO
MIMO
<
<
<
6
PSO-MISMO
GA-MISMO
MISMO
MIMO
<
<
<
3
PSO-MISMO
GA-MISMO
MISMO
MIMO

* indicates the mean difference between the two adjacent strategies is significant at the 0.05 level

<*
<*
<*
<
<
<
<
<*
<*
<*
<
<
<
<*
<*
<
<*
<
<
<*
<*
<*
<*
<*
<*
<*
<*
<*
<*
<*
<*
<*
<*
<*
<*
<*
<*
<
<
<*
<*
<*
<*
<*
<*
<*
<*
<*
<*
<
<*
<*
<*
<*
<*
<
<*
<*
<*

5
DIR
DIR
DIR
ITER
DIR
MIMO
MIMO
DIR
DIR
DIR
DIR
ITER
MIMO
DIR
DIR
MIMO
DIR
DIR
DIR
DIR
DIR
DIR
DIR
ITER
DIR
DIR
DIR
DIR
DIR
DIR
ITER
ITER
ITER
DIR
DIR
DIR
ITER
DIR
ITER
ITER
DIR
DIR
DIR
DIR
DIR
DIR
DIR
DIR
DIR
ITER
DIR
DIR
DIR
DIR
DIR
DIR
DIR
DIR
ITER

11

<*
<*
<*
<
<*
<*
<*
<*
<*
<*
<*
<
<*
<*
<*
<*
<*
<
<*
<
<*
<*
<*
<*
<*
<
<*
<*
<*
<*
<*
<
<*
<*
<*
<
<*
<
<
<
<*
<*
<
<
<*
<*
<*
<*
<*
<
<*
<*
<*
<*
<
<*
<
<*
<

6
ITER
ITER
ITER
DIR
ITER
ITER
ITER
ITER
ITER
ITER
ITER
DIR
ITER
ITER
ITER
ITER
ITER
ITER
ITER
ITER
ITER
ITER
ITER
DIR
ITER
ITER
ITER
ITER
ITER
ITER
DIR
DIR
MIMO
ITER
ITER
ITER
DIR
ITER
MIMO
DIR
ITER
ITER
ITER
ITER
ITER
ITER
ITER
ITER
ITER
DIR
ITER
ITER
ITER
ITER
ITER
ITER
ITER
ITER
DIR

> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) <

B. PSO-MISMO vs. GA-MISMO in Convergence and
Computational Time
Up to now, we have compared the proposed PSO-MISMO
with other five competitors (i.e., iterated, direct, MIMO,
MISMO, and GA-MISMO) in terms of prediction accuracy. In
this subsection, the convergence and computational time of the
two heuristic-based prediction strategies (i.e., PSO-MISMO
and GA-MISMO) are further examined.
It should be noted that although we have used a number of
time series for each datasets, i.e., Logistic, Mackey-Glass, and
NN3 datasets, general results do not change much within the
time series for each datasets. Therefore, to save space, we

report only the results for randomly selected three time series
from Logistic, Mackey-Glass, and NN3 datasets, respectively.
They are the No.1 series from Logistic time series with
ϕ1 = 0.100 and length=485, the No.18 series from
Mackey-Glass time series with ϕ1 = 2 , τ = 17 , and length=964,
and the No.55 series from NN3 dataset. In addition, in the
following experiment, the swarm/population size is fixed to 20
for both PSO-MISMO and GA-MISMO strategies, and
generation number is set to 100, as described in Section V. Each
strategy is run ten times over each of these time series, in an
attempt to eliminate the influence of a lucky initial solution. All
the numerical experiments are performed on a personal
computer, Inter(R) Core(TM) 2 Duo CPU 2.50 GHz, 1.87-GB
memory, and MATLAB environment (Version R2009b).
As was mentioned above, the learning process is repeated 10
times for both PSO-MISMO and GA-MISMO for each time
series. Thus, typical convergences of the best and average
fitness values as a function of generations for both
PSO-MISMO and GA-MISMO on three aforementioned time
series are shown in Fig. 9(a)-(c), respectively. As it can be seen,
in a typical PSO-MISMO optimization process the best fitness
value decreases rapidly and converges at about 48 generations
for No.1 Logistic time series, 65 generations for No.18
Mackey-Glass time series, and 59 generations for No. 55 time
series in NN3 datasets; whereas for GA-MISMO it takes about
63 generations for No.1 Logistic time series, 77 generations for
No.18 Mackey-Glass time series, and 68 generations for No. 55
time series in NN3 datasets. To show how the evolution process
is going on for both PSO-MISMO and GA-MISMO, the
convergence of the average fitness values are also shown in Fig.
9. Looking at Fig. 9, it is clear that PSO-MISMO seems to
perform better than GA-MISMO. Thus, for the present problem
the performance of the PSO-MISMO is better than
GA-MISMO from an evolutionary point of view.
0.9098
Fitness of PSO-MISMO and GA-MISMO (MSE)

 As far as the comparison between PSO-MISMO and
GA-MISMO is concerned, they are almost a tie and the results
are mixing among the prediction measures examined. In terms
of MAPE, PSO-MISMO wins in Logistic dataset and NN3
dataset but loses in Mackey-Glass dataset. In terms of SMAPE,
PSO-MISMO loses in Mackey-Glass dataset and NN3 dataset
but wins in Logistic dataset. In terms of MASE, PSO-MISMO
wins in Mackey-Glass dataset and NN3 dataset but loses in
Logistic dataset.
Following the experimental procedure presented in Fig. 8, an
ANOVA procedure is performed to determine if there exists a
statistically significant difference among the six modeling
strategies in the hold-out sample for each of the performance
measures and the prediction horizon. The results are not
included in detail to save space. All ANOVA results are
significant at the 0.05 level (with a few exceptions), suggesting
that there are significant differences among the six modeling
strategies. To further identify the significant difference between
any two strategies, the Tukey’s HSD test is used to compare all
pair-wise differences simultaneously. Tukey’s HSD test is a
post-hoc test, meaning that Tukey’s HSD test should not be
performed unless the results of the ANOVA procedure are
positive. Table V shows the results of the multiple comparison
tests for three datasets. For each accuracy measure and
prediction horizon, the strategies are rank ordered from 1 (the
best) to 6 (the worst). Several observations can be made from
Table V:
 For each accuracy measure and dataset, PSO-MISMO and
GA-MISMO significantly outperforms MISMO for the
overwhelming majority of prediction horizons.
 Considering the two heuristic-based modeling strategies, we
can see that, whatever the dataset used, whatever the accuracy
measure considered, and whatever the prediction horizon
examined, the difference in prediction performance between
PSO-MISMO and GA-MISMO is not significant at the 0.05
level.
 As far as the comparison MISMO vs. MIMO is concerned,
the difference in prediction performance is not significant at the
0.05 level, even with a few exceptions
 Generally, the strategies based on the single-output structure
perform significantly worse than those based on the
multiple-output structure.
 Concerning the current two leading strategies, the direct
strategy significantly outperforms the iterated strategy for the
majority of prediction horizons.
 The iterated strategy performs the poorest at 95% statistical
confidence level in most cases, even with some exceptions.

12

PSO-MISMO Best
PSO-MISMO Average
GA-MISMO Best
GA-MISMO Average

0.6998

0.5386

0.4141

0.3185

0.2451
0.0612
0

10

20

30

40

50
60
Generation
(a)

70

80

90

100

> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) <

13

VIII. CONCLUSIONS
Fitness of PSO-MISMO and GA-MISMO (MSE)

0.1688

PSO-MISMO Best
PSO-MISMO Average
GA-MISMO Best
GA-MISMO Average

0.1062
0.0667
0.0420
0.0264
0.0166
0.0104
0.0066
0.0041
0.0026
0

10

20

30

40

50
60
Generation

70

80

90

100

(b)

Fitness of PSO-MISMO and GA-MISMO (MSE)

2308912

PSO-MISMO Best
PSO-MISMO Average
GA-MISMO Best
GA-MISMO Average

2025361

1776633

Multi-step-ahead time series prediction has usually proved to
be an intractable task due to the growing amount of
uncertainties arising from various sources. For instance, an
accumulation of errors and lack of information make
multi-step-ahead prediction more difficult. Thus, modeling
strategies for multi-step-ahead prediction are major research
topics with significant practical implications. The contribution
of this study is an extension to a well-established MISMO
modeling strategy by incorporating a heuristic based on binary
particle swarm optimization into the MISMO modeling process
to self-adaptively determine the number of sub-models with
varying prediction horizons, and conduct a large scale
comparative study with neural networks for validation.
Quantitative and comprehensive assessments are performed
with the simulated and real time series on the basis of the
prediction accuracy, convergence, and computational time.
Experimental results and comparisons demonstrate the
superiority of the proposed PSO-MISMO modeling strategy for
multi-step-ahead time series prediction.

1558450

REFERENCES
[1]

1367061

1199177

1051909
0

[2]
10

20

30

40

50
60
Generation
(c)

70

80

90

100

[3]

Fig. 9. Convergence of best and average fitness for PSO-MISMO and
GA-MISMO strategies. (a) No.1 Logistic time series with ϕ1 = 0.100 and

[4]

length=485. (b) No.18 Mackey-Glass time series with ϕ1 = 2 , τ = 17 , and
length=964.(c) No.55 time series in NN3 dataset.

[5]

The second test performed is designed to measure the
average computational time needed per time series from all the
datasets within these two strategies. The result is shown in
Table VI, from which we can see that the PSO-MISMO needs
generally less time than GA-MISMO.

[6]

TABLE VI
COMPUTATIONAL TIME FOR PSO-MISMO AND GA-MISMO
Dataset
Strategy
CPU time (M)
PSO-MISMO
24.514
Logistic dataset
GA-MISMO
38.421
Mackey-Glass dataset

PSO-MISMO
GA-MISMO

30.824
41.725

NN3 dataset

PSO-MISMO
GA-MISMO

11.854
19.627

[7]

[8]

[9]

Generally speaking, thus, from the above experiments, we
can draw the following three main conclusions. (1) Two
heuristic-based modeling strategies (i.e., PSO-MISMO and
GA-MISMO) are capable of obtaining high-quality
multi-step-ahead time series forecasts related to the four
well-established strategies. (2) The performance of the
PSO-MISMO is better than that of GA-MISMO from an
evolutionary point of view. (3) PSO-MISMO showed to be
faster than GA-MISMO in terms of average running time.

[10]

[11]

[12]

[13]

A. F. Atiya, S. M. El-Shoura, S. I. Shaheen, and M. S. El-Sherif, "A
comparison between neural-network forecasting techniques-case
study: river flow forecasting," Neural Networks, IEEE Transactions
on, vol. 10, pp. 402-409, 1999.
G. Chevillon, "DIRECT MULTI ‐ STEP ESTIMATION AND
FORECASTING," J. Econ. Surveys, vol. 21, pp. 746-785, 2007.
A. Sorjamaa, J. Hao, N. Reyhani, Y. Ji, and A. Lendasse,
"Methodology for long-term prediction of time series,"
Neurocomputing, vol. 70, pp. 2861-2869, 2007.
A. Sorjamaa and A. Lendasse, "Time series prediction using DirRec
strategy," presented at the European Symposium on Artificial
Neural Networks Bruges, 2006.
G. Bontempi, "Long term time series prediction with multi-input
multi-output local learning," in Proceedings of the 2nd European
Symposium on Time Series Prediction (TSP), ESTSP08, Helsinki,
Finland, 2008, pp. 145-154.
S. Ben Taieb, G. Bontempi, A. Sorjamaa, and A. Lendasse,
"Long-term prediction of time series by combining direct and
MIMO strategies," in Proceedings of the 2009 IEEE International
Joint Conference on Neural Networks, Atlanta, U.S.A., 2009, pp.
3054-3061.
S. Ben Taieb, A. Sorjamaa, and G. Bontempi, "Multiple-output
modeling for multi-step-ahead time series forecasting,"
Neurocomputing, vol. 73, pp. 1950-1957, 2010.
S. Ben Taieb, G. Bontempi, A. F. Atiya, and A. Sorjamaa, "A
review and comparison of strategies for multi-step ahead time series
forecasting based on the NN5 forecasting competition," Expert Syst.
Appl., 2012.
H. S. Hippert, C. E. Pedreira, and R. C. Souza, "Neural networks for
short-term load forecasting: A review and evaluation," Power
Systems, IEEE Transactions on, vol. 16, pp. 44-55, 2001.
A. Sharma, "Seasonal to interannual rainfall probabilistic forecasts
for improved water supply management: Part 1—A strategy for
system predictor identification," J Hydrol, vol. 239, pp. 232-239,
2000.
G. P. Zhang and D. M. Kline, "Quarterly time-series forecasting
with neural networks," Neural Networks, IEEE Transactions on, vol.
18, pp. 1800-1814, 2007.
J. Kennedy and R. Eberhart, "Particle swarm optimization," in
Proceedings of the IEEE International Conference on Neural
Networks, 1995, pp. 1942-1948 vol. 4.
C. F. Juang, "A hybrid of genetic algorithm and particle swarm
optimization for recurrent network design," Systems, Man, and
Cybernetics, Part B: Cybernetics, IEEE Transactions on, vol. 34, pp.
997-1006, 2004.

> REPLACE THIS LINE WITH YOUR PAPER IDENTIFICATION NUMBER (DOUBLE-CLICK HERE TO EDIT) <
[14]

[15]

[16]

[17]
[18]

[19]

[20]

[21]
[22]

[23]

[24]

[25]

[26]

[27]
[28]
[29]
[30]

[31]
[32]

[33]

[34]

[35]
[36]

[37]

Y. P. Chen, W. C. Peng, and M. C. Jian, "Particle swarm
optimization with recombination and dynamic linkage discovery,"
Systems, Man, and Cybernetics, Part B: Cybernetics, IEEE
Transactions on, vol. 37, pp. 1460-1470, 2007.
B. Liu, L. Wang, and Y. H. Jin, "An effective PSO-based memetic
algorithm for flow shop scheduling," Systems, Man, and
Cybernetics, Part B: Cybernetics, IEEE Transactions on, vol. 37, pp.
18-27, 2007.
R. Xu, J. Xu, and D. Wunsch, "A Comparison Study of Validity
Indices on Swarm-Intelligence-Based Clustering," Systems, Man,
and Cybernetics, Part B: Cybernetics, IEEE Transactions on, vol.
42, pp. 1243-1256, 2012.
T. Nguyen, S. Yang, and C. Li, "A self-learning particle swarm
optimizer for global optimization problems," 2012.
J. Kennedy and R. C. Eberhart, "A discrete binary version of the
particle swarm algorithm," in Systems, Man, and Cybernetics, 1997.
Computational Cybernetics and Simulation., 1997 IEEE
International Conference on, Orlando, FL, 1997, pp. 4104-4108 vol.
5.
K. De Jong, "Parameter setting in EAs: a 30 year perspective," in
Parameter Setting in Evolutionary Algorithms, ed: Springer, 2007,
pp. 1-18.
I. C. Trelea, "The particle swarm optimization algorithm:
convergence analysis and parameter selection," Inform Process Lett,
vol. 85, pp. 317-325, 2003.
J. F. Kennedy, J. Kennedy, and R. C. Eberhart, Swarm intelligence:
Morgan Kaufmann Pub, 2001.
Y. Shi and R. C. Eberhart, "Empirical study of particle swarm
optimization," in Evolutionary Computation, 1999. CEC 99.
Proceedings of the 1999 Congress on, 1999.
B. Liu, L. Wang, and Y.-H. Jin, "An effective PSO-based memetic
algorithm for flow shop scheduling," Systems, Man, and
Cybernetics, Part B: Cybernetics, IEEE Transactions on, vol. 37, pp.
18-27, 2007.
G. Rubio, H. Pomares, I. Rojas, and L. J. Herrera, "A heuristic
method for parameter selection in LS-SVM: Application to time
series prediction," Int J Forecasting, vol. 27, pp. 725-739, 2011.
L.-C. Chang, P.-A. Chen, and F.-J. Chang, "Reinforced
Two-Step-Ahead Weight Adjustment Technique for Online
Training of Recurrent Neural Networks," Neural Networks and
Learning Systems, IEEE Transactions on, vol. 23, pp. 1269-1278,
2012.
Y. Liu and X. Yao, "Simultaneous training of negatively correlated
neural networks in an ensemble," Systems, Man, and Cybernetics,
Part B: Cybernetics, IEEE Transactions on, vol. 29, pp. 716-725,
1999.
R. M. May, "Simple mathematical models with very complicated
dynamics," Nature, vol. 261, pp. 459-467, 1976.
M. C. Mackey and L. Glass, "Oscillation and chaos in physiological
control systems," Science, vol. 197, pp. 287-289, 1977.
D. E. Goldberg, "Genetic algorithms in search, optimization, and
machine learning," Reading, MA: Addison-Wesley, 1989.
K. M. Sim, "BLGAN: Bayesian learning and genetic algorithm for
supporting negotiation with incomplete information," Systems, Man,
and Cybernetics, Part B: Cybernetics, IEEE Transactions on, vol.
39, pp. 198-211, 2009.
L. Davis, "Handbook of genetic algorithms," New York: Van
Nostrand Reinhold, 1991.
D. Golmohammadi, R. C. Creese, H. Valian, and J. Kolassa,
"Supplier selection based on a neural network model using genetic
algorithm," Neural Networks, IEEE Transactions on, vol. 20, pp.
1504-1519, 2009.
D. E. Goldberg and K. Deb, "A comparative analysis of selection
schemes used in genetic algorithms," Urbana, vol. 51, pp.
61801-2996, 1991.
M. V. Baghmisheh, K. Madani, and A. Navarbaf, "A discrete
shuffled frog optimization algorithm," Artif Intell Rev, vol. 36, pp.
267-284, 2011.
R. L. Haupt and S. E. Haupt, Practical genetic algorithms:
Wiley-Interscience, 2004.
E. Elbeltagi, T. Hegazy, and D. Grierson, "Comparison among five
evolutionary-based optimization algorithms," Adv Eng Inform, vol.
19, pp. 43-53, 2005.
F. Ramsay and D. Schaefer, The Statistical Sleuth: Duxbury, Boston,
Mass, 1996.

14

Yukun Bao received his B.Sc., M.Sc.,
and Ph.D in management science and
engineering from Huazhong University
of Science and Technology, P.R.China,
in 1996, 1999 and 2002 respectively.
Currently, he is an associate professor
at the department of management science
and information systems, school of
management, Huazhong University of
Science and Technology, P.R.China. He
has been the principal investigator for two research projects
funded by Natural Science Foundation of China and has served
as a referee of paper review for several IEEE journals and
international journals and a PC member for several
international academic conferences. His research interests are
time series modeling and forecasting, business intelligence and
data mining.

Tao Xiong received his B.Sc. and M.Sc.
in management science and engineering
from Huazhong University of Science
and Technology, P.R.China, in 2008 and
2010 respectively.
Currently, he is working toward the
Ph.D. degree in management science and
engineering, Huazhong University of
Science and Technology, China. His
research
interests
include
multi-step-ahead time series forecasting, interval data analysis,
computational intelligence.

Zhongyi Hu received B.Sc. degree in
information system in 2009 from
Harbin Institute of Technology at
Weihai, P.R.China. He received his
M.Sc. degree in management science
and engineering in 2011 from
Huazhong University of Science and
Technology, P.R.China.
Currently, he is working toward the
Ph.D. degree in management science and engineering,
Huazhong University of Science and Technology, P.R.China.
His research interests include support vector machines, swarm
intelligence, memetic algorithms, and time series forecasting.

