# Local Functional Coefficient Autoregressive Model for Multistep Prediction of Chaotic Time Series

Authors: Liyun Su; Chenlong Li
Journal: Discrete Dynamics in Nature and Society
 (2015)

---

## Abstract

A new methodology, which combines nonparametric method based on local functional coefficient autoregressive (LFAR) form with chaos theory and regional method, is proposed for multistep prediction of chaotic time series. The objective of this research study is to improve the performance of long-term forecasting of chaotic time series. To obtain the prediction values of chaotic time series, three steps are involved. Firstly, the original time series is reconstructed inm-dimensional phase space with a time delay τ by using chaos theory. Secondly, select the nearest neighbor points by using local method in the m-dimensional phase space. Thirdly, we use the nearest neighbor points to get a LFAR model. The proposed model’s parameters are selected by modified generalized cross validation (GCV) criterion. Both simulated data (Lorenz and Mackey-Glass systems) and real data (Sunspot time series) are used to illustrate the performance of the proposed methodology. By detailed investigation and comparing our results with published researches, we find that the LFAR model can effectively fit nonlinear characteristics of chaotic time series by using simple structure and has excellent performance for multistep forecasting.


---

## Body

## 1. Introduction

In recent decades, researchers have paid much attention to chaos motion in many fields, such as meteorology, medicine, economics, signal processing, traffic flow, power load, Sunspot prediction, and many others [1–12] and bring about lots of new models for predicting chaotic time series. In the late 1960s, researchers found it is a difficult task to forecast chaotic time series which is the evolution of a chaotic system’s observations by using traditional time series forecasting methods [1]. Then a series of theories and methods was established for understanding essence of chaos motion, such as Takes’ embedding theory [13]. Now, chaos theory has become an important part of nonlinear science and is used for forecasting chaotic time series.Up to now, modeling of chaotic systems constructed from observed data and predicting one or several future values of the time series have become an important issue [14]. There are many prediction methods that have been proposed, such as adaptive prediction [15], the support vector machine (SVM) [16–20], polynomial estimation [21–24], and neural network (NN) [25–29]. In most of the published literature, single-step prediction was considered. For multistep prediction, the direct and iterative methods are proposed as two main categories. The direct multistep prediction does not use the prediction values in the future; the iterative multistep prediction uses short-term predictor and is built through recursive prediction, which means the future values are calculated by the predictor itself. However, multistep prediction becomes a difficult task because of the limited largest Lyapunov exponent of the chaotic system. Some researchers have been focusing on multistep prediction and using NN or its extended models to improve the performance of multistep prediction [29–31]. Some researchers’ studies show that the accuracy of prediction can be improved by using hybrid technique, such as combined SVM and Neuro-Fuzzy [32] and neural network and Neuro-Fuzzy [25, 26]. Researchers’ studies also show that hybrid technique can appear to have good performance by using the prediction error, such as the combined PCA and SVM [19] and ARMA and RESN [7]. The generalized nonlinear filtering methods are investigated for 5-step prediction of chaotic time series in [33]. These methods generally prompt better results than those single models, but they are complex, affected by personal experience, and easy to overfit.In this paper, we propose to use functional coefficient autoregressive (FAR) model instead of local linear structure to approximate the local attractor in reconstructed phase space. As in [34], it is a nonparametric estimation of nonlinear dynamics. The proposed method combines chaos theory and local technique and has excellent spatial adaptation to effectively fit nonlinear characteristics of chaotic time series. Unlike RBF-AR model, the LFAR model has reasonable simple implementation and is rarely affected by personal experience. Furthermore, the LFAR model can avoid overfitting by controlling the dimension of the primary functions which are used for estimating the functional coefficients of LFAR model. In this study, an algorithm based on the dynamic least squares criterion for estimation of local functional coefficients is proposed. The effectiveness of the proposed model is demonstrated by the application to simulated data (Lorenz and Mackey-Glass systems) and real data (Sunspot time series). In these cases, we analyze and estimate the functional coefficients by using the proposed algorithm and examine the properties of iterative multistep prediction.The remainder of this paper is organized as follows. Section2 reviews the concept of the LFAR model and the optimal parameter set is established by using GCV. Section 3 uses the simulated chaotic systems and one real life time series as examples to evaluate the proposed models and discuss the properties of model’s parameters and also compares the results with published researches. Section 4 presents the conclusion of this paper.


## 2. Methodology

### 2.1. Phase Space Reconstruction for Chaotic Time Series

For a scalar chaotic time seriesx
1
,
x
2
,
…
,
x
n, the phase space can be reconstructed by Takes’ embedding theory and the reconstructed phase points are X
t
=
(
x
t
,
x
t
-
τ
,
…
,
x
t
-
(
m
-
1
)
τ
)
T, where t
=
1
+
(
m
-
1
)
τ
,
2
+
(
m
-
1
)
τ
,
…
,
n. The embedding dimension m and the time delay τ can be obtained by using Cao’s method [38]. Then, a continued vector mapping F
:
R
m
→
R
m or f
:
R
m
→
R can be described by the unknown evolution from X
t to X
t
+
1 or x
t
+
1. That is, X
t
+
1
=
F
(
X
t
) or x
t
+
1
=
f
(
X
t
).


### 2.2. The LFAR Model and Estimation Method

A chaotic time series prediction model for describing evolution fromX
t to x
t
+
1 can be written as(1)
x
t
+
1
=
f
X
t
=
f
x
t
,
x
t
-
τ
,
…
,
x
t
-
m
-
1
τ
.The continued vector mapping f is the best prediction function in the sense that f minimizes the expected prediction error:(2)
min
f
⁡
E
x
t
+
1
-
f
x
t
,
x
t
-
τ
,
…
,
x
t
-
m
-
1
τ
2
.The saturated nonparametric functionf
(
x
t
,
x
t
-
τ
,
…
,
x
t
-
(
m
-
1
)
τ
) cannot be estimated with reasonable accuracy due to the curse of dimensionality [30]. A LFAR model for chaotic reconstruction data is presented in m-dimensional reconstructed phase space. That is,(3)
x
t
+
1
=
∑
j
=
1
m
g
j
x
t
-
d
x
t
-
j
-
1
τ
,where d
∈
{
τ
,
2
τ
,
…
,
(
m
-
1
)
τ
} is the lag of the model-dependent variable x
t
-
d, m is the embedding dimension, τ is the time delay, and functional coefficients {
g
j
}
j
=
1
m are continued functions.The functional coefficients{
g
j
}
j
=
1
m are difficult to estimate because they are considered as nonparametric and do not have a conformed form. There are many nonlinear forms that can be used. Finding a good nonlinear form is hard by trying one model after another. Here the local nonparametric method is applied to obtain the estimations of unknown functional coefficients {
g
j
}
j
=
1
m. Using Taylor’s series expansion, g
j
(
u
) with r-order derivative near the point u
0 can be described as follows:(4)
g
j
u
=
∑
i
=
0
r
g
j
i
u
0
i
!
u
-
u
0
i
+
o
u
-
u
0
r
.Letf
i
(
u
)
=
(
u
-
u
0
)
i,  i
=
0,1
,
…
,
r,  u
∈
(
u
0
-
h
,
u
0
+
h
),  h
>
0 is the bandwidth for controlling the number of the nearest neighbor points, and a
j
i
=
g
j
(
i
)
(
u
0
)
/
i
!, i
=
0,1
,
…
,
r. Ignoring the higher order infinitesimal, we can approximate g
j
(
x
t
-
d
) by the r-dimensional primary functions {
f
0
,
f
1
,
…
,
f
r
} as follows:(5)
g
j
x
t
-
d
≈
∑
i
=
0
r
g
j
i
u
0
i
!
u
-
u
0
i
≡
∑
i
=
0
r
a
j
i
f
i
u
,where u is used to replace variable x
t
-
d and a
j
i represent the weight coefficients of primary functions.The LFAR model at the current state pointX
M in m-dimensional phase space can be described as follows:(6)
x
M
+
m
-
1
τ
+
1
=
∑
j
=
1
m
g
j
x
t
-
d
x
M
+
m
-
1
τ
-
j
-
1
τ
=
∑
j
=
1
m
∑
i
=
0
r
a
j
i
f
i
u
x
M
+
m
-
1
τ
-
j
-
1
τ
,where u
∈
(
x
M
+
(
m
-
1
)
τ
-
d
-
h
,
x
M
+
(
m
-
1
)
τ
-
d
+
h
), h
>
0. In order to obtain a LFAR model at the current state point X
M in the reconstructed phase space, we select q nearest neighbor points X
M
j
(
j
=
1,2
,
…
,
q
) by using the Euclidean distances d
(
i
)
=
X
i
-
X
M
(
i
=
1,2
,
…
,
M
-
1
). The estimations of the parameters θ
^
≡
(
a
^
00
,
…
,
a
^
0
r
,
…
,
a
^
m
0
,
…
,
a
^
m
r
)
(
m
+
1
)
(
r
+
1
)
×
1
T of the LFAR model at state point X
M can be obtained by solving the following weighted least squares (WLS) regression problem:(7)
θ
^
X
M
;
h
=
arg
⁡
min
θ
⁡
J
X
M
;
h
,where(8)
J
X
M
;
h
=
∑
t
=
1
q
x
M
t
+
m
-
1
τ
+
1
-
∑
j
=
1
m
∑
i
=
0
r
a
j
i
f
i
x
M
t
+
m
-
1
τ
-
d
x
M
t
+
j
-
1
τ
2
K
h
z
t
,and K
h
(
z
t
) is a kernel function and z
t
≡
d
(
t
)
-
min
⁡
d
1
,
d
2
,
…
,
d
q
(
o
r
z
t
≡
d
(
t
)
) is a nonnegative function which emphasizes neighbor observations around X
M. The parameter h is used to determine the weights of neighboring observations around X
M in estimating θ
^
(
X
M
).Letu
=
x
M
+
(
m
-
1
)
τ
-
d; we have g
^
j
(
x
t
-
d
)
=
a
^
j
0. Then we can obtain a prediction of next point:(9)
x
^
M
+
m
-
1
τ
+
1
=
∑
j
=
1
m
g
^
j
x
t
-
d
x
M
+
m
-
1
τ
-
j
-
1
τ
=
∑
j
=
1
m
a
^
j
0
x
M
+
m
-
1
τ
-
j
-
1
τ
.Set(10)
χ
=
U
1
x
M
1
U
1
x
M
1
+
τ
⋯
U
1
x
M
1
+
m
-
1
τ
U
2
x
M
2
U
2
x
M
2
+
τ
⋯
U
2
x
M
2
+
m
-
1
τ
⋮
⋮
⋱
⋮
U
q
x
M
q
U
q
x
M
q
+
τ
⋯
U
q
x
M
q
+
m
-
1
τ
q
×
m
r
+
1
,
Y
=
x
M
1
+
m
-
1
τ
+
1
,
…
,
x
M
q
+
m
-
1
τ
+
1
,
U
t
=
f
0
x
M
t
+
m
-
1
τ
-
d
,
f
1
x
M
t
+
m
-
1
τ
-
d
,
…
,
f
r
x
M
t
+
m
-
1
τ
-
d
,
t
=
1,2
,
…
,
q
.Then the WLS solution is(11)
θ
^
X
M
=
χ
T
W
χ
-
1
χ
T
W
Y
,where W is a q
×
q diagonal matrix with K
h
(
z
i
) as its ith diagonal element, which entails g
^
j
(
x
t
-
d
)
=
a
^
j
0.


### 2.3. Determination of Optimal Parameters

In the process of modeling LFAR, the parameter set(
m
,
τ
,
q
,
r
,
d
,
h
) needs to be estimated. The embedding dimension m and the time delay τ can be calculated by Cao’s method and the autocorrelation function method. For the remaining part of the parameter set, it is clear to see that the accuracy of prediction based on LFAR model is sensitive to the kernel function and the parameter h. Here we have the form of kernel function as follows:(12)
K
h
z
t
=
exp
⁡
-
z
t
a
h
b
.Conventionally, we have a
=
0.5 and b
=
1. Let h
→
0; then K
h
(
z
t
)
→
0 for any z
t
>
0. And let h
→
+
∞, the kernel K
h
(
z
t
)
→
1 for any z
t
>
0, which dues to the local linear model. For 0
<
h
<
+
∞, the kernel K
h
(
z
t
) changes from 0 to 1, which means that the weight at the same neighbor point changes from 0 to 1. The parameter a can adjust the weight and the parameter b can adjust the convergence speed of kernel K
h
(
z
t
) when h
→
0 or h
→
+
∞. In this study, we consider simulated chaotic systems and real life time series as examples to investigate a proper parameter h to achieve the best performance for modeling a LFAR model.It is also very crucial to choose a proper dimension of primary functions to achieve the best performance for modeling chaotic time series. Low dimension leads to bad simulation results, and high dimension increases the complexity of computation and leads to overfitting. Hence, our main purpose of this section is to determine the parameter set(
q
,
r
,
d
,
h
). Generally, the optimal dimension should be selected to minimize the mean squared error (MSE) or its improved versions. Here we use a simple and quick method which is proposed in [34]. It can be regarded as a modified generalized cross validation (
G
C
V
) criterion. Let k and Q be two given positive integers and satisfy n
>
k
Q. First we use Q subseries with sample size n
-
s
k
(
s
=
1,2
,
…
,
Q
) to estimate the unknown coefficient functions and then to compute the multistep forecasting errors of the next part with sample size m based on the estimated model. For example, the data with sample size n
-
4
k is used to get the estimated model, and the prediction errors for the next k data are computed. Then, the data with sample size n
-
3
k is used and so on. The average prediction error or the standard prediction errors [31] use the subseries which is given by(13)
SAPE
s
h
=
∑
i
=
n
-
s
k
+
1
n
-
s
k
+
k
x
i
-
x
^
i
2
∑
i
=
n
-
s
k
+
1
n
-
s
k
+
k
x
i
-
x
¯
2
,where x
¯
=
1
/
k
∑
i
=
n
-
s
k
+
1
n
-
s
k
+
k
x
i. The overall prediction error is given by(14)
SAPE
h
=
Q
-
1
∑
s
=
1
Q
SAPE
s
h
.Fan et al. [34] set k
=
[
0.1
n
] and Q
=
4, and Meng and Peng [28] set k
=
[
0.03
n
] and Q
=
3. The selected bandwidth does not critically depend on the choice of k and Q as long as k
Q is reasonably large. In practical implementations, we select Meng’s method.We select the proper bandwidth by minimizingS
A
P
E
(
h
). The function S
A
P
E
(
h
) is minimized by comparing its values in a finite set of scale parameters in a grid H
~
=
h
j
∣
h
j
=
h
0
l
j
,
j
=
1
,
…
,
J
,
l
>
1
,
h
0
>
0. We can obtain different S
A
P
E
(
h
) with different q, r, and d. It is also very important to choose the parameter set (
q
,
r
,
d
). Finally we can determinate the optimal parameter set (
q
o
p
t
,
r
o
p
t
,
d
o
p
t
,
h
o
p
t
).


### 2.4. Algorithm Description and Multistep Prediction

Select the optimal parameter set in executable LFAR algorithm that contains the process of the phase space reconstruction parameter set(
m
,
τ
), the number of nearest neighbor points q, the dimension of primary functions r, and the weight of neighbor observations parameter h. Here, we choose the number of nearest neighbor points from 2
m
(
r
+
1
)
+
1 to 2
m
(
r
+
1
)
+
50. To speed up the computation and avoid overfitting, we let r
<
8. We select J
=
19 and  l
=
1.2 and set h
0
=
0.05. The optimal parameters’ selection is described in detail in Algorithm 1.Algorithm 1:Searching algorithm for optimal parameter set.
Input: The relevant data, such as: chaotic time series, the initial parameters
Output: Results, such as: the optimal parameters
Scale chaotic time series between[
0,1
];
Form
∈
m
~ %  m
~
=
(
1
,
…
,
20
)
Forτ
∈
τ
~ %  τ
~
=
(
1
,
…
,
20
)
The autocorrelation function method, The Cao method.
End
End %Obtain the embedding dimension m
opt and the time delay τ
opt
Choose training set and select primary functions, Phase space reconstruction based on(
m
opt
,
τ
opt
),
Compute the Euclidean distancesd
(
⋅
) between phase points in the reconstructed phase space.
For(
h
,
q
,
r
,
d
)
∈
(
h
~
,
q
~
,
r
~
,
d
~
)
%h
0
=
0.05
,
J
=
19
,
q
~
=
2
m
r
+
1
+
1,2
m
r
+
1
+
6
,
…
,
2
m
r
+
1
+
51,  r
~
=
0,1
,
…
,
7,
d
~
=
τ
opt
,
2
τ
opt
,
…
,
m
opt
-
1
τ
opt
ComputeSAPE
(
h
,
q
,
r
,
d
).
End
Select the optimal parameter set:h
opt
,
q
opt
,
r
opt
,
d
opt
=
arg
⁡
min
h
,
q
,
r
,
d
⁡
SAPE
h
,
q
,
r
,
dIn Section2.2, the single-step prediction is given. For multistep forecasting, as we say in Section 1, there are two possible approaches. The iterative l-step prediction is to add x
^
t
+
1
,
x
^
t
+
2
,
…
,
x
^
t
+
l
-
1 to the training set and utilizes the LFAR model iteratively. The direct l-step prediction is to fit(15)
x
^
M
+
m
-
1
τ
+
1
+
l
=
∑
j
=
1
m
g
^
j
x
t
-
d
x
M
+
m
-
1
τ
-
j
-
1
τ
.The optimal parameter set (
r
o
p
t
,
d
o
p
t
,
h
o
p
t
) needs to be selected to minimize the S
A
P
E again. For multistep prediction in this paper, we use the iterative approach. The executable multistep prediction of LFAR model based on Algorithm 1 is described in detail in Algorithm 2.Algorithm 2:Multistep prediction algorithm.
Input: The relevant data, such as: chaotic time series, the optimal parameters
Output: Results, such as: the prediction values
Scale chaotic time series between[
0,1
];
Phase space reconstruction based on embedding dimensionm
opt and the time delay τ
opt;
Compute the Euclidean distancesd
(
⋅
) between phase points in the reconstructed phase space;
Fort
∈
t
~  %  t
~ is unknown data.
Compute prediction values by using the optimal parameters (
h
opt
,
q
opt
,
r
opt
,
d
opt
);
Add prediction values to the training set;
Compute Euclidean distances between the new phase point and others;
End


## 2.1. Phase Space Reconstruction for Chaotic Time Series

For a scalar chaotic time seriesx
1
,
x
2
,
…
,
x
n, the phase space can be reconstructed by Takes’ embedding theory and the reconstructed phase points are X
t
=
(
x
t
,
x
t
-
τ
,
…
,
x
t
-
(
m
-
1
)
τ
)
T, where t
=
1
+
(
m
-
1
)
τ
,
2
+
(
m
-
1
)
τ
,
…
,
n. The embedding dimension m and the time delay τ can be obtained by using Cao’s method [38]. Then, a continued vector mapping F
:
R
m
→
R
m or f
:
R
m
→
R can be described by the unknown evolution from X
t to X
t
+
1 or x
t
+
1. That is, X
t
+
1
=
F
(
X
t
) or x
t
+
1
=
f
(
X
t
).


## 2.2. The LFAR Model and Estimation Method

A chaotic time series prediction model for describing evolution fromX
t to x
t
+
1 can be written as(1)
x
t
+
1
=
f
X
t
=
f
x
t
,
x
t
-
τ
,
…
,
x
t
-
m
-
1
τ
.The continued vector mapping f is the best prediction function in the sense that f minimizes the expected prediction error:(2)
min
f
⁡
E
x
t
+
1
-
f
x
t
,
x
t
-
τ
,
…
,
x
t
-
m
-
1
τ
2
.The saturated nonparametric functionf
(
x
t
,
x
t
-
τ
,
…
,
x
t
-
(
m
-
1
)
τ
) cannot be estimated with reasonable accuracy due to the curse of dimensionality [30]. A LFAR model for chaotic reconstruction data is presented in m-dimensional reconstructed phase space. That is,(3)
x
t
+
1
=
∑
j
=
1
m
g
j
x
t
-
d
x
t
-
j
-
1
τ
,where d
∈
{
τ
,
2
τ
,
…
,
(
m
-
1
)
τ
} is the lag of the model-dependent variable x
t
-
d, m is the embedding dimension, τ is the time delay, and functional coefficients {
g
j
}
j
=
1
m are continued functions.The functional coefficients{
g
j
}
j
=
1
m are difficult to estimate because they are considered as nonparametric and do not have a conformed form. There are many nonlinear forms that can be used. Finding a good nonlinear form is hard by trying one model after another. Here the local nonparametric method is applied to obtain the estimations of unknown functional coefficients {
g
j
}
j
=
1
m. Using Taylor’s series expansion, g
j
(
u
) with r-order derivative near the point u
0 can be described as follows:(4)
g
j
u
=
∑
i
=
0
r
g
j
i
u
0
i
!
u
-
u
0
i
+
o
u
-
u
0
r
.Letf
i
(
u
)
=
(
u
-
u
0
)
i,  i
=
0,1
,
…
,
r,  u
∈
(
u
0
-
h
,
u
0
+
h
),  h
>
0 is the bandwidth for controlling the number of the nearest neighbor points, and a
j
i
=
g
j
(
i
)
(
u
0
)
/
i
!, i
=
0,1
,
…
,
r. Ignoring the higher order infinitesimal, we can approximate g
j
(
x
t
-
d
) by the r-dimensional primary functions {
f
0
,
f
1
,
…
,
f
r
} as follows:(5)
g
j
x
t
-
d
≈
∑
i
=
0
r
g
j
i
u
0
i
!
u
-
u
0
i
≡
∑
i
=
0
r
a
j
i
f
i
u
,where u is used to replace variable x
t
-
d and a
j
i represent the weight coefficients of primary functions.The LFAR model at the current state pointX
M in m-dimensional phase space can be described as follows:(6)
x
M
+
m
-
1
τ
+
1
=
∑
j
=
1
m
g
j
x
t
-
d
x
M
+
m
-
1
τ
-
j
-
1
τ
=
∑
j
=
1
m
∑
i
=
0
r
a
j
i
f
i
u
x
M
+
m
-
1
τ
-
j
-
1
τ
,where u
∈
(
x
M
+
(
m
-
1
)
τ
-
d
-
h
,
x
M
+
(
m
-
1
)
τ
-
d
+
h
), h
>
0. In order to obtain a LFAR model at the current state point X
M in the reconstructed phase space, we select q nearest neighbor points X
M
j
(
j
=
1,2
,
…
,
q
) by using the Euclidean distances d
(
i
)
=
X
i
-
X
M
(
i
=
1,2
,
…
,
M
-
1
). The estimations of the parameters θ
^
≡
(
a
^
00
,
…
,
a
^
0
r
,
…
,
a
^
m
0
,
…
,
a
^
m
r
)
(
m
+
1
)
(
r
+
1
)
×
1
T of the LFAR model at state point X
M can be obtained by solving the following weighted least squares (WLS) regression problem:(7)
θ
^
X
M
;
h
=
arg
⁡
min
θ
⁡
J
X
M
;
h
,where(8)
J
X
M
;
h
=
∑
t
=
1
q
x
M
t
+
m
-
1
τ
+
1
-
∑
j
=
1
m
∑
i
=
0
r
a
j
i
f
i
x
M
t
+
m
-
1
τ
-
d
x
M
t
+
j
-
1
τ
2
K
h
z
t
,and K
h
(
z
t
) is a kernel function and z
t
≡
d
(
t
)
-
min
⁡
d
1
,
d
2
,
…
,
d
q
(
o
r
z
t
≡
d
(
t
)
) is a nonnegative function which emphasizes neighbor observations around X
M. The parameter h is used to determine the weights of neighboring observations around X
M in estimating θ
^
(
X
M
).Letu
=
x
M
+
(
m
-
1
)
τ
-
d; we have g
^
j
(
x
t
-
d
)
=
a
^
j
0. Then we can obtain a prediction of next point:(9)
x
^
M
+
m
-
1
τ
+
1
=
∑
j
=
1
m
g
^
j
x
t
-
d
x
M
+
m
-
1
τ
-
j
-
1
τ
=
∑
j
=
1
m
a
^
j
0
x
M
+
m
-
1
τ
-
j
-
1
τ
.Set(10)
χ
=
U
1
x
M
1
U
1
x
M
1
+
τ
⋯
U
1
x
M
1
+
m
-
1
τ
U
2
x
M
2
U
2
x
M
2
+
τ
⋯
U
2
x
M
2
+
m
-
1
τ
⋮
⋮
⋱
⋮
U
q
x
M
q
U
q
x
M
q
+
τ
⋯
U
q
x
M
q
+
m
-
1
τ
q
×
m
r
+
1
,
Y
=
x
M
1
+
m
-
1
τ
+
1
,
…
,
x
M
q
+
m
-
1
τ
+
1
,
U
t
=
f
0
x
M
t
+
m
-
1
τ
-
d
,
f
1
x
M
t
+
m
-
1
τ
-
d
,
…
,
f
r
x
M
t
+
m
-
1
τ
-
d
,
t
=
1,2
,
…
,
q
.Then the WLS solution is(11)
θ
^
X
M
=
χ
T
W
χ
-
1
χ
T
W
Y
,where W is a q
×
q diagonal matrix with K
h
(
z
i
) as its ith diagonal element, which entails g
^
j
(
x
t
-
d
)
=
a
^
j
0.


## 2.3. Determination of Optimal Parameters

In the process of modeling LFAR, the parameter set(
m
,
τ
,
q
,
r
,
d
,
h
) needs to be estimated. The embedding dimension m and the time delay τ can be calculated by Cao’s method and the autocorrelation function method. For the remaining part of the parameter set, it is clear to see that the accuracy of prediction based on LFAR model is sensitive to the kernel function and the parameter h. Here we have the form of kernel function as follows:(12)
K
h
z
t
=
exp
⁡
-
z
t
a
h
b
.Conventionally, we have a
=
0.5 and b
=
1. Let h
→
0; then K
h
(
z
t
)
→
0 for any z
t
>
0. And let h
→
+
∞, the kernel K
h
(
z
t
)
→
1 for any z
t
>
0, which dues to the local linear model. For 0
<
h
<
+
∞, the kernel K
h
(
z
t
) changes from 0 to 1, which means that the weight at the same neighbor point changes from 0 to 1. The parameter a can adjust the weight and the parameter b can adjust the convergence speed of kernel K
h
(
z
t
) when h
→
0 or h
→
+
∞. In this study, we consider simulated chaotic systems and real life time series as examples to investigate a proper parameter h to achieve the best performance for modeling a LFAR model.It is also very crucial to choose a proper dimension of primary functions to achieve the best performance for modeling chaotic time series. Low dimension leads to bad simulation results, and high dimension increases the complexity of computation and leads to overfitting. Hence, our main purpose of this section is to determine the parameter set(
q
,
r
,
d
,
h
). Generally, the optimal dimension should be selected to minimize the mean squared error (MSE) or its improved versions. Here we use a simple and quick method which is proposed in [34]. It can be regarded as a modified generalized cross validation (
G
C
V
) criterion. Let k and Q be two given positive integers and satisfy n
>
k
Q. First we use Q subseries with sample size n
-
s
k
(
s
=
1,2
,
…
,
Q
) to estimate the unknown coefficient functions and then to compute the multistep forecasting errors of the next part with sample size m based on the estimated model. For example, the data with sample size n
-
4
k is used to get the estimated model, and the prediction errors for the next k data are computed. Then, the data with sample size n
-
3
k is used and so on. The average prediction error or the standard prediction errors [31] use the subseries which is given by(13)
SAPE
s
h
=
∑
i
=
n
-
s
k
+
1
n
-
s
k
+
k
x
i
-
x
^
i
2
∑
i
=
n
-
s
k
+
1
n
-
s
k
+
k
x
i
-
x
¯
2
,where x
¯
=
1
/
k
∑
i
=
n
-
s
k
+
1
n
-
s
k
+
k
x
i. The overall prediction error is given by(14)
SAPE
h
=
Q
-
1
∑
s
=
1
Q
SAPE
s
h
.Fan et al. [34] set k
=
[
0.1
n
] and Q
=
4, and Meng and Peng [28] set k
=
[
0.03
n
] and Q
=
3. The selected bandwidth does not critically depend on the choice of k and Q as long as k
Q is reasonably large. In practical implementations, we select Meng’s method.We select the proper bandwidth by minimizingS
A
P
E
(
h
). The function S
A
P
E
(
h
) is minimized by comparing its values in a finite set of scale parameters in a grid H
~
=
h
j
∣
h
j
=
h
0
l
j
,
j
=
1
,
…
,
J
,
l
>
1
,
h
0
>
0. We can obtain different S
A
P
E
(
h
) with different q, r, and d. It is also very important to choose the parameter set (
q
,
r
,
d
). Finally we can determinate the optimal parameter set (
q
o
p
t
,
r
o
p
t
,
d
o
p
t
,
h
o
p
t
).


## 2.4. Algorithm Description and Multistep Prediction

Select the optimal parameter set in executable LFAR algorithm that contains the process of the phase space reconstruction parameter set(
m
,
τ
), the number of nearest neighbor points q, the dimension of primary functions r, and the weight of neighbor observations parameter h. Here, we choose the number of nearest neighbor points from 2
m
(
r
+
1
)
+
1 to 2
m
(
r
+
1
)
+
50. To speed up the computation and avoid overfitting, we let r
<
8. We select J
=
19 and  l
=
1.2 and set h
0
=
0.05. The optimal parameters’ selection is described in detail in Algorithm 1.Algorithm 1:Searching algorithm for optimal parameter set.
Input: The relevant data, such as: chaotic time series, the initial parameters
Output: Results, such as: the optimal parameters
Scale chaotic time series between[
0,1
];
Form
∈
m
~ %  m
~
=
(
1
,
…
,
20
)
Forτ
∈
τ
~ %  τ
~
=
(
1
,
…
,
20
)
The autocorrelation function method, The Cao method.
End
End %Obtain the embedding dimension m
opt and the time delay τ
opt
Choose training set and select primary functions, Phase space reconstruction based on(
m
opt
,
τ
opt
),
Compute the Euclidean distancesd
(
⋅
) between phase points in the reconstructed phase space.
For(
h
,
q
,
r
,
d
)
∈
(
h
~
,
q
~
,
r
~
,
d
~
)
%h
0
=
0.05
,
J
=
19
,
q
~
=
2
m
r
+
1
+
1,2
m
r
+
1
+
6
,
…
,
2
m
r
+
1
+
51,  r
~
=
0,1
,
…
,
7,
d
~
=
τ
opt
,
2
τ
opt
,
…
,
m
opt
-
1
τ
opt
ComputeSAPE
(
h
,
q
,
r
,
d
).
End
Select the optimal parameter set:h
opt
,
q
opt
,
r
opt
,
d
opt
=
arg
⁡
min
h
,
q
,
r
,
d
⁡
SAPE
h
,
q
,
r
,
dIn Section2.2, the single-step prediction is given. For multistep forecasting, as we say in Section 1, there are two possible approaches. The iterative l-step prediction is to add x
^
t
+
1
,
x
^
t
+
2
,
…
,
x
^
t
+
l
-
1 to the training set and utilizes the LFAR model iteratively. The direct l-step prediction is to fit(15)
x
^
M
+
m
-
1
τ
+
1
+
l
=
∑
j
=
1
m
g
^
j
x
t
-
d
x
M
+
m
-
1
τ
-
j
-
1
τ
.The optimal parameter set (
r
o
p
t
,
d
o
p
t
,
h
o
p
t
) needs to be selected to minimize the S
A
P
E again. For multistep prediction in this paper, we use the iterative approach. The executable multistep prediction of LFAR model based on Algorithm 1 is described in detail in Algorithm 2.Algorithm 2:Multistep prediction algorithm.
Input: The relevant data, such as: chaotic time series, the optimal parameters
Output: Results, such as: the prediction values
Scale chaotic time series between[
0,1
];
Phase space reconstruction based on embedding dimensionm
opt and the time delay τ
opt;
Compute the Euclidean distancesd
(
⋅
) between phase points in the reconstructed phase space;
Fort
∈
t
~  %  t
~ is unknown data.
Compute prediction values by using the optimal parameters (
h
opt
,
q
opt
,
r
opt
,
d
opt
);
Add prediction values to the training set;
Compute Euclidean distances between the new phase point and others;
End


## 3. Numerical Experiments and Performance Evaluation

To compare our results with others, we looked for published research doing long-term prediction, testing it with Lorenz system, Mackey-Glass system, and the Sunspot time series. All these chaotic time series are scaled within[
0,1
] as follows:(16)
X
new
t
=
x
t
-
min
⁡
x
max
⁡
x
-
min
⁡
x
.


### 3.1. Prediction of Lorenz System

Lorenz time series can be produced as follows (Lorenz, 1963, see [1])(17)
x
˙
=
σ
y
-
x
,
y
˙
=
r
-
z
x
-
y
,
z
˙
=
x
y
-
b
z
,where σ, r, and b are commonly selected as σ
=
10, r
=
28, and b
=
8
/
3. The standard fourth-order Runge-Kutta method is used to get Lorenz time series, and the x-coordinate is used as observations. A time series with sample size 2300 is randomly generated. The first 1800 pieces of data are used for training. The data with sample size 1800
-
3
∗
50 is used to estimate model, and the prediction errors of the next 50 pieces of data are computed. Then, the data with sample size 1800
-
2
∗
50 is used, and the prediction errors are computed and so on. And the rest of data is treated as testing. The results of multistep prediction for Lorenz time series are shown in Figure 1.Figure 1
Results of Lorenz time series: (a) multistep prediction errors and (b) the original values and multistep prediction values, respectively.
(a)
(b)From Figure1, it can be seen that the LFAR model has small error values, and the values of multistep prediction are in good agreements with the real data. And the multistep prediction values of LFAR start diverging significantly from the 450th time step; this implies that the proposed model has a good performance of the multistep prediction of chaotic time series.The results of prediction are shown in Figures2, 3, and 4. Here we have four parameters, but we cannot show them in a 5-dimensional picture. So we sort the parameter set and show the S
A
P
E less than e
-
6 in Figure 2. From Figure 2, we can see that the LFAR model’s S
A
P
E with different parameters are similar. More details about the parameters are shown in Figures 3 and 4. Figure 3 shows that the functional coefficients’ order mainly takes four, the lag variable is selected from x
t to x
t
-
(
m
-
1
)
τ at different phase points, the weight parameter h changes with other parameters, and the number of nearest neighbor points is chosen in the vicinity of 2
m
(
r
+
1
)
+
41. We can obtain the optimal parameters from Figure 2. In Figure 4, we change one parameter and fix the other optimal parameters to investigate the influence of the changed parameter. From Figure 4, we can see that parameters can effectively influence the prediction error.Figure 2
Results of Lorenz time series’ multistep prediction: theS
A
P
E of parameter set.Figure 3
Results of Lorenz time series’ multistep prediction: (a) the functional coefficients’ order; (b) the delay of the lag variable (y-coordinate value y is equivalent to (
y
-
1
)
τ); (c) the weight parameter; (d) the number of nearest neighbor points (equivalent to number-2
m
(
r
+
1
)). Note: for Figure 4, the delay of the lag variable and the number of nearest neighbor points are the same as this figure, respectively.
(a)
(b)
(c)
(d)Figure 4
Results of Lorenz time series’ multistep prediction: (a) changed functional coefficients’ order; (b) changed delay of the lag variable; (c) changed number of nearest neighbor points; (d) changed radius of the neighbors.
(a)
(b)
(c)
(d)In Figure5, we show the functional coefficients’ estimation. From this figure, we can see that the estimated values are seasonal within 400 prediction steps, and the accuracy of prediction declines and does not follow precious law beyond 400 steps.Figure 5
The functional coefficients’ estimations: (a)–(f) the estimated values fromx
t
-
(
m
-
1
)
τ to x
t, respectively.
(a)
(b)
(c)
(d)
(e)
(f)


### 3.2. Prediction of Mackey-Glass System

Mackey-Glass system is used in literature as a benchmark model due to its chaotic characteristics [39]. Mackey-Glass time series is generated by the following discrete form:(18)
x
t
+
1
=
1
-
b
x
t
+
a
x
t
-
τ
1
+
x
t
-
τ
10
,where a
=
0.2, b
=
0.1, and τ
=
17 and initial conditions x
(
t
)
t
≤
0
=
1.2. Thus, we can obtain a scalar chaotic time series sample set with sample size of 3300. Then we choose 1800 pieces of data for training, and the remaining part is treated as testing data.Figure6 compares the real values with prediction values for the remaining part. From Figure 6, it can be seen that the LFAR model has small error values, and the values of prediction are in good agreements with the real values. The multistep prediction values of LFAR model start diverging significantly from the 900th time step, and this implies that the proposed model has a good performance for the multistep prediction of chaotic time series.Figure 6
Results of Mackey-Glass time series: (a) and (b) multistep prediction errors and values, respectively.
(a)
(b)The results of prediction are shown in Figures7, 8, and 9. From Figure 7, we can see that the LFAR model’s S
A
P
E with different parameters are similar. Figure 8 shows that the functional coefficients’ order mainly takes four, the lag variable is selected from x
t to x
t
-
(
m
-
1
)
τ at different phase points, the weight parameter h changes with other parameters, and the number of nearest neighbor points is chosen in the vicinity of 2
m
(
r
+
1
)
+
36. We can obtain the optimal parameters from Figure 7. From Figure 9, we can see that parameters can effectively influence the prediction error.Figure 7
Multistep prediction results of Mackey-Glass time series: theS
A
P
E of parameter set.Figure 8
Results of Mackey-Glass time series’ multistep prediction: (a) the order of functional coefficients; (b) the delay of the lag variable; (c) the weight parameter; (d) the number of nearest neighbor points.
(a)
(b)
(c)
(d)Figure 9
Multistep prediction results of Mackey-Glass time series: (a) the order of changed functional coefficients; (b) changed delay of the lag variable; (c) changed number of nearest neighbor points; (d) changed radius of the neighbors. Note: to speed up computation, we break off the computation when the prediction error values are large and setS
A
P
E
=
10.
(a)
(b)
(c)
(d)In Figure10, we show the functional coefficients’ estimations. From this figure, we can see that the estimated values are seasonal within 900 prediction steps, but the prediction accuracy declines and does not follow the precious law beyond 900 steps.Figure 10
The functional coefficients’ estimations: (a)–(f) the estimated values fromx
t
-
(
m
-
1
)
τ to x
t, respectively.
(a)
(b)
(c)
(d)
(e)
(f)


### 3.3. Prediction for Sunspot Time Series

Sunspot time series is a good indication of solar activity for solar cycles. The monthly smoothed Sunspot time series is obtained from the SIDC (World Data Center for the Sunspot Index). To compare the results with different models in the literature, data are selected in the same conditions reported by [8, 9]. Sunspot series from November 1834 to June 2001 (2000 points) are selected and scaled within [
0,1
]. The first 1000 samples of time series are selected for training and the remaining 1000 samples are kept to test the prediction models.Results are shown in Figures11 and 12. From Figure 11, it can be seen that the prediction values of LFAR model have accuracy of prediction.Figure 11
Results of Sunspot time series: (a) single-step prediction error and (b) the original and single-step prediction values, respectively.
(a)
(b)Figure 12
Results of Sunspot time series: (a) multistep prediction error and (b) the original and multistep prediction values, respectively.
(a)
(b)From Figure12, for the multistep prediction, we find that the prediction values of LFAR model start diverging significantly from the 60
t
h time step. This is because the real data has noise, and the noise affects the performance of prediction. Besides, the 60 time steps can still help us predict 5-year Sunspot data in the future.


### 3.4. Results and Discussion

We compare the proposed models with some of the models reported in the literature and measure performance with mean squared error (MSE), root mean squared error (RMSE) and normalized mean squared error (NMSE), relative error (p
e
r
r), and symmetric mean absolute percentage error (SMAPE), namely,(19)
MSE
=
1
n
∑
t
=
1
n
y
t
-
y
^
t
2
,
RMSE
=
1
n
∑
t
=
1
n
y
t
-
y
^
t
2
,
NMSE
=
∑
t
=
1
n
y
t
-
y
^
t
2
∑
t
=
1
n
y
t
-
y
¯
2
,
perr
=
∑
t
=
1
n
y
t
-
y
^
t
2
∑
t
=
1
n
y
t
2
,
SMAPE
=
1
n
∑
t
=
1
n
y
t
-
y
^
t
∑
t
=
1
n
1
/
2
y
t
+
y
^
t
.The comparative results are shown in Tables1–3. Data cannot be selected in the same conditions reported in the literature; thus the conclusions from Tables 1–3 have a few mistakes. In Table 1, we can see that the proposed models in this paper are better than some of the existing methods for predicting Lorenz time series. But the best results are of RBF and PSORBF [30]. This implies that the RBF can predict multistep results well, but the optimal parameters are difficult to obtain. The proposed model can obtain all optimal parameters. Table 2 presents that the proposed models are the best. For the real-world time series in Table 3, we can see that the proposed models are better than most of the existing methods only except the CERNN [29].Table 1
The comparative results of different models of Lorenz time series(
(
r
opt
,
d
opt
,
h
opt
,
q
opt
,
m
opt
,
τ
opt
)
=
(
4,5
,
1.1093,101,6
,
7
)
).
Multistep prediction
Prediction model
RMSE
perr
6-step prediction
CERNN (2008) [29]
1.25
E
-
04
The proposed LFAR (test data = 500)
1.66
E
-
04
6.12
E
-
08
100-step prediction
OOB-DMS (2014) [35]
3.16
The proposed LFAR (test data = 450)
4.9
E
-
03
5.63
E
-
05
300-step prediction
OOB-DMS (2014) [35]
5.48
RBF (2014) [30]
9.98
E
-
07
PSORBF (2014) [30]
6.92
E
-
08
The proposed LFAR (test data = 450)
2.51
E
-
02
1.5
E
-
03
450-step prediction
The proposed LFAR (test data = 450)
8.8
E
-
03
1.80
E
-
04Table 2
The comparative results of different models of Mackey-Glass time series(
(
r
opt
,
d
opt
,
h
opt
,
q
opt
,
m
opt
,
τ
opt
)
=
(
3,3
,
0.1792,54,6
,
6
)
).
Multistep prediction
Prediction model
RMSE
MSE
SMAPE
5-step prediction
EKF (2013) [33]
7.82
E
-
02
GEKF (2013) [33]
4.85
E
-
02
UKF (2013) [33]
1.89
E
-
02
GUKF (2013) [33]
2.6
E
-
03
GPF (2013) [33]
3.97
E
-
05
GGPF (2013) [33]
5.36
E
-
06
The proposed LFAR (test data = 500)
1.20
E
-
04
1.45
E
-
08
2.16
E
-
04
12-step prediction
CERNN (2008) [29]
3.79
E
-
03
The proposed LFAR (test data = 500)
1.71
E
-
04
2.91
E
-
08
2.91
E
-
04
30-step prediction
IGFC (2015) [36]
0.54
The proposed LFAR (test data = 500)
2.19
E
-
04
4.79
E
-
08
7.18
E
-
04
250-step prediction
N-best (2014) [37]
2.47
E
-
03
The proposed LFAR (test data = 500)
9.08
E
-
04
8.24
E
-
07
2.0
E
-
03
614-step prediction
HCNN (2011) [14]
0.5085
500-step prediction
The proposed LFAR (test data = 500)
1.3
E
-
03
1.72
E
-
06
2.9
E
-
03Table 3
The comparative results of different models of Sunspot time series(
(
r
opt
,
d
opt
,
h
opt
,
q
opt
,
m
opt
,
τ
opt
)
=
(
6,2
,
1.3312,139,7
,
13
)
).
Multistep prediction
Prediction model
RMSE
NMSE
1-step prediction
ERNN (2007) [8]
1.29
E
-
02
2.80
E
-
03
RBF-OLS (2009) [9]
4.60
E
-
02
HE-NARXNN (2010) [10]
1.19
E
-
02
5.90
E
-
04
ARMA-RESN (2013) [7]
15.7937
The proposed LFAR (test data = 500)
9.50
E
-
03 (1.38)
7.91
E
-
04
18-step prediction
CERNN (2008) [29]
8.23
E
-
03
MLPNN (2010) [31]
0.6915
FTLRNN (2010) [31]
0.3210
SOFM (2010) [31]
0.7085
The proposed LFAR (test data = 300)
9.59
E
-
02
0.1933
60-step prediction
The proposed LFAR (test data = 60)
0.1825
1.0680


## 3.1. Prediction of Lorenz System

Lorenz time series can be produced as follows (Lorenz, 1963, see [1])(17)
x
˙
=
σ
y
-
x
,
y
˙
=
r
-
z
x
-
y
,
z
˙
=
x
y
-
b
z
,where σ, r, and b are commonly selected as σ
=
10, r
=
28, and b
=
8
/
3. The standard fourth-order Runge-Kutta method is used to get Lorenz time series, and the x-coordinate is used as observations. A time series with sample size 2300 is randomly generated. The first 1800 pieces of data are used for training. The data with sample size 1800
-
3
∗
50 is used to estimate model, and the prediction errors of the next 50 pieces of data are computed. Then, the data with sample size 1800
-
2
∗
50 is used, and the prediction errors are computed and so on. And the rest of data is treated as testing. The results of multistep prediction for Lorenz time series are shown in Figure 1.Figure 1
Results of Lorenz time series: (a) multistep prediction errors and (b) the original values and multistep prediction values, respectively.
(a)
(b)From Figure1, it can be seen that the LFAR model has small error values, and the values of multistep prediction are in good agreements with the real data. And the multistep prediction values of LFAR start diverging significantly from the 450th time step; this implies that the proposed model has a good performance of the multistep prediction of chaotic time series.The results of prediction are shown in Figures2, 3, and 4. Here we have four parameters, but we cannot show them in a 5-dimensional picture. So we sort the parameter set and show the S
A
P
E less than e
-
6 in Figure 2. From Figure 2, we can see that the LFAR model’s S
A
P
E with different parameters are similar. More details about the parameters are shown in Figures 3 and 4. Figure 3 shows that the functional coefficients’ order mainly takes four, the lag variable is selected from x
t to x
t
-
(
m
-
1
)
τ at different phase points, the weight parameter h changes with other parameters, and the number of nearest neighbor points is chosen in the vicinity of 2
m
(
r
+
1
)
+
41. We can obtain the optimal parameters from Figure 2. In Figure 4, we change one parameter and fix the other optimal parameters to investigate the influence of the changed parameter. From Figure 4, we can see that parameters can effectively influence the prediction error.Figure 2
Results of Lorenz time series’ multistep prediction: theS
A
P
E of parameter set.Figure 3
Results of Lorenz time series’ multistep prediction: (a) the functional coefficients’ order; (b) the delay of the lag variable (y-coordinate value y is equivalent to (
y
-
1
)
τ); (c) the weight parameter; (d) the number of nearest neighbor points (equivalent to number-2
m
(
r
+
1
)). Note: for Figure 4, the delay of the lag variable and the number of nearest neighbor points are the same as this figure, respectively.
(a)
(b)
(c)
(d)Figure 4
Results of Lorenz time series’ multistep prediction: (a) changed functional coefficients’ order; (b) changed delay of the lag variable; (c) changed number of nearest neighbor points; (d) changed radius of the neighbors.
(a)
(b)
(c)
(d)In Figure5, we show the functional coefficients’ estimation. From this figure, we can see that the estimated values are seasonal within 400 prediction steps, and the accuracy of prediction declines and does not follow precious law beyond 400 steps.Figure 5
The functional coefficients’ estimations: (a)–(f) the estimated values fromx
t
-
(
m
-
1
)
τ to x
t, respectively.
(a)
(b)
(c)
(d)
(e)
(f)


## 3.2. Prediction of Mackey-Glass System

Mackey-Glass system is used in literature as a benchmark model due to its chaotic characteristics [39]. Mackey-Glass time series is generated by the following discrete form:(18)
x
t
+
1
=
1
-
b
x
t
+
a
x
t
-
τ
1
+
x
t
-
τ
10
,where a
=
0.2, b
=
0.1, and τ
=
17 and initial conditions x
(
t
)
t
≤
0
=
1.2. Thus, we can obtain a scalar chaotic time series sample set with sample size of 3300. Then we choose 1800 pieces of data for training, and the remaining part is treated as testing data.Figure6 compares the real values with prediction values for the remaining part. From Figure 6, it can be seen that the LFAR model has small error values, and the values of prediction are in good agreements with the real values. The multistep prediction values of LFAR model start diverging significantly from the 900th time step, and this implies that the proposed model has a good performance for the multistep prediction of chaotic time series.Figure 6
Results of Mackey-Glass time series: (a) and (b) multistep prediction errors and values, respectively.
(a)
(b)The results of prediction are shown in Figures7, 8, and 9. From Figure 7, we can see that the LFAR model’s S
A
P
E with different parameters are similar. Figure 8 shows that the functional coefficients’ order mainly takes four, the lag variable is selected from x
t to x
t
-
(
m
-
1
)
τ at different phase points, the weight parameter h changes with other parameters, and the number of nearest neighbor points is chosen in the vicinity of 2
m
(
r
+
1
)
+
36. We can obtain the optimal parameters from Figure 7. From Figure 9, we can see that parameters can effectively influence the prediction error.Figure 7
Multistep prediction results of Mackey-Glass time series: theS
A
P
E of parameter set.Figure 8
Results of Mackey-Glass time series’ multistep prediction: (a) the order of functional coefficients; (b) the delay of the lag variable; (c) the weight parameter; (d) the number of nearest neighbor points.
(a)
(b)
(c)
(d)Figure 9
Multistep prediction results of Mackey-Glass time series: (a) the order of changed functional coefficients; (b) changed delay of the lag variable; (c) changed number of nearest neighbor points; (d) changed radius of the neighbors. Note: to speed up computation, we break off the computation when the prediction error values are large and setS
A
P
E
=
10.
(a)
(b)
(c)
(d)In Figure10, we show the functional coefficients’ estimations. From this figure, we can see that the estimated values are seasonal within 900 prediction steps, but the prediction accuracy declines and does not follow the precious law beyond 900 steps.Figure 10
The functional coefficients’ estimations: (a)–(f) the estimated values fromx
t
-
(
m
-
1
)
τ to x
t, respectively.
(a)
(b)
(c)
(d)
(e)
(f)


## 3.3. Prediction for Sunspot Time Series

Sunspot time series is a good indication of solar activity for solar cycles. The monthly smoothed Sunspot time series is obtained from the SIDC (World Data Center for the Sunspot Index). To compare the results with different models in the literature, data are selected in the same conditions reported by [8, 9]. Sunspot series from November 1834 to June 2001 (2000 points) are selected and scaled within [
0,1
]. The first 1000 samples of time series are selected for training and the remaining 1000 samples are kept to test the prediction models.Results are shown in Figures11 and 12. From Figure 11, it can be seen that the prediction values of LFAR model have accuracy of prediction.Figure 11
Results of Sunspot time series: (a) single-step prediction error and (b) the original and single-step prediction values, respectively.
(a)
(b)Figure 12
Results of Sunspot time series: (a) multistep prediction error and (b) the original and multistep prediction values, respectively.
(a)
(b)From Figure12, for the multistep prediction, we find that the prediction values of LFAR model start diverging significantly from the 60
t
h time step. This is because the real data has noise, and the noise affects the performance of prediction. Besides, the 60 time steps can still help us predict 5-year Sunspot data in the future.


## 3.4. Results and Discussion

We compare the proposed models with some of the models reported in the literature and measure performance with mean squared error (MSE), root mean squared error (RMSE) and normalized mean squared error (NMSE), relative error (p
e
r
r), and symmetric mean absolute percentage error (SMAPE), namely,(19)
MSE
=
1
n
∑
t
=
1
n
y
t
-
y
^
t
2
,
RMSE
=
1
n
∑
t
=
1
n
y
t
-
y
^
t
2
,
NMSE
=
∑
t
=
1
n
y
t
-
y
^
t
2
∑
t
=
1
n
y
t
-
y
¯
2
,
perr
=
∑
t
=
1
n
y
t
-
y
^
t
2
∑
t
=
1
n
y
t
2
,
SMAPE
=
1
n
∑
t
=
1
n
y
t
-
y
^
t
∑
t
=
1
n
1
/
2
y
t
+
y
^
t
.The comparative results are shown in Tables1–3. Data cannot be selected in the same conditions reported in the literature; thus the conclusions from Tables 1–3 have a few mistakes. In Table 1, we can see that the proposed models in this paper are better than some of the existing methods for predicting Lorenz time series. But the best results are of RBF and PSORBF [30]. This implies that the RBF can predict multistep results well, but the optimal parameters are difficult to obtain. The proposed model can obtain all optimal parameters. Table 2 presents that the proposed models are the best. For the real-world time series in Table 3, we can see that the proposed models are better than most of the existing methods only except the CERNN [29].Table 1
The comparative results of different models of Lorenz time series(
(
r
opt
,
d
opt
,
h
opt
,
q
opt
,
m
opt
,
τ
opt
)
=
(
4,5
,
1.1093,101,6
,
7
)
).
Multistep prediction
Prediction model
RMSE
perr
6-step prediction
CERNN (2008) [29]
1.25
E
-
04
The proposed LFAR (test data = 500)
1.66
E
-
04
6.12
E
-
08
100-step prediction
OOB-DMS (2014) [35]
3.16
The proposed LFAR (test data = 450)
4.9
E
-
03
5.63
E
-
05
300-step prediction
OOB-DMS (2014) [35]
5.48
RBF (2014) [30]
9.98
E
-
07
PSORBF (2014) [30]
6.92
E
-
08
The proposed LFAR (test data = 450)
2.51
E
-
02
1.5
E
-
03
450-step prediction
The proposed LFAR (test data = 450)
8.8
E
-
03
1.80
E
-
04Table 2
The comparative results of different models of Mackey-Glass time series(
(
r
opt
,
d
opt
,
h
opt
,
q
opt
,
m
opt
,
τ
opt
)
=
(
3,3
,
0.1792,54,6
,
6
)
).
Multistep prediction
Prediction model
RMSE
MSE
SMAPE
5-step prediction
EKF (2013) [33]
7.82
E
-
02
GEKF (2013) [33]
4.85
E
-
02
UKF (2013) [33]
1.89
E
-
02
GUKF (2013) [33]
2.6
E
-
03
GPF (2013) [33]
3.97
E
-
05
GGPF (2013) [33]
5.36
E
-
06
The proposed LFAR (test data = 500)
1.20
E
-
04
1.45
E
-
08
2.16
E
-
04
12-step prediction
CERNN (2008) [29]
3.79
E
-
03
The proposed LFAR (test data = 500)
1.71
E
-
04
2.91
E
-
08
2.91
E
-
04
30-step prediction
IGFC (2015) [36]
0.54
The proposed LFAR (test data = 500)
2.19
E
-
04
4.79
E
-
08
7.18
E
-
04
250-step prediction
N-best (2014) [37]
2.47
E
-
03
The proposed LFAR (test data = 500)
9.08
E
-
04
8.24
E
-
07
2.0
E
-
03
614-step prediction
HCNN (2011) [14]
0.5085
500-step prediction
The proposed LFAR (test data = 500)
1.3
E
-
03
1.72
E
-
06
2.9
E
-
03Table 3
The comparative results of different models of Sunspot time series(
(
r
opt
,
d
opt
,
h
opt
,
q
opt
,
m
opt
,
τ
opt
)
=
(
6,2
,
1.3312,139,7
,
13
)
).
Multistep prediction
Prediction model
RMSE
NMSE
1-step prediction
ERNN (2007) [8]
1.29
E
-
02
2.80
E
-
03
RBF-OLS (2009) [9]
4.60
E
-
02
HE-NARXNN (2010) [10]
1.19
E
-
02
5.90
E
-
04
ARMA-RESN (2013) [7]
15.7937
The proposed LFAR (test data = 500)
9.50
E
-
03 (1.38)
7.91
E
-
04
18-step prediction
CERNN (2008) [29]
8.23
E
-
03
MLPNN (2010) [31]
0.6915
FTLRNN (2010) [31]
0.3210
SOFM (2010) [31]
0.7085
The proposed LFAR (test data = 300)
9.59
E
-
02
0.1933
60-step prediction
The proposed LFAR (test data = 60)
0.1825
1.0680


## 4. Conclusions

We propose a new methodology for forecasting chaotic time series based on FAR model, chaos theory, and local nonparametric technique. Firstly, the chaotic time series are reconstructed inm-dimensional phase space with a time delay τ by using chaos theory. Secondly, the neighbor points are selected by using local method in the m-dimensional phase space. Thirdly, we use the nearest neighbor points to identify a novel FAR model by using local least squares method. Finally, all parameters are calculated by the GCV criterion and judged by the sense of SAPE. The proposed functional coefficient autoregressive (FAR) model is used instead of local linear structure to approximate the local attractor in reconstructed phase space, which is a local nonparametric estimation of nonlinear dynamics.An algorithm based on the dynamic least squares criterion for estimation of local functional coefficients is proposed. For Mackey-Glass and Lorenz attractors, the parameters have been investigated carefully. For Sunspot time series, which is noise in data, have better one-step forecasting performance. And multistep prediction model can forecast accurately before 60th, and this can help us predict 5-year Sunspot data in the future. In these cases, we analyze and estimate the functional coefficients by using the proposed algorithm and examine the properties of iterative multistep prediction. By detailed investigation and comparing our results with published researches, we find that the LFAR model can effectively fit nonlinear characteristics of chaotic time series by using simple structure and has excellent performance for multistep forecasting.



---

*Source: 329487-2015-11-10.xml*