# Local Prediction of Chaotic Time Series Based on Polynomial Coefficient Autoregressive Model

Authors: Liyun Su; Chenlong Li
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

We apply the polynomial function to approximate the functional coefficients of the state-dependent autoregressive model for chaotic time series prediction. We present a novel local nonlinear model called local polynomial coefficient autoregressive prediction (LPP) model based on the phase space reconstruction. The LPP model can effectively fit nonlinear characteristics of chaotic time series with simple structure and have excellent one-step forecasting performance. We have also proposed a kernel LPP (KLPP) model which applies the kernel technique for the LPP model to obtain better multistep forecasting performance. The proposed models are  flexible to analyze complex and multivariate nonlinear structures. Both simulated and real data examples are used for illustration.


---

## Body

## 1. Introduction

Chaos is widely encountered in nature and many scientific areas since Lorenz found chaotic motion in his research on meteorology [1]. Over the latest decades, researchers in many fields have paid much attention to chaos [2–13], and modeling and prediction of chaotic time series have become popular in the areas of meteorology [1, 2], medicine [3], economics [4], signal processing [5], traffic flow [6], power load [7], Sunspot prediction [8–11], and many others [12, 13].Although prediction of chaotic time series is very difficult, the chaos theory [14] provides a useful tool to predict chaotic time series. With the development of chaos theory and research on its application technique, the global [15] and local [16] methods are proposed as two main categories. The global method makes an attempt to approximate the whole time series on all attractors and seeks a function which is valid at every point. Then, we can get the future values from knowledge of all the previous elements of the time series. However, the parameters may be changed when new information is added into the model. The local method only use part of the past information to approximate the local attractor. The future values can be inferred by using the neighborhoods of the current points. Many prediction methods have gained popularity in practice during the last decades based on two main categories. Such as adaptive prediction [17–19], the support vector machine (SVM) [20–27], polynomial estimation [28–30], and neural network [8–11, 31–36].Recently, some researchers studies have shown that local methods can obtain generally better results than those obtained with global methods [18]. And some researchers found that the forecast accuracy can be improved by using some combining techniques both in the global and local method. By using those combining techniques, the parameters can be obtained faster and the analysis of residual can be underestimated if residuals are not randomness. Such as SVM combined with neurofuzzy model [21], SVM combined with PCA [26], ARMA combined with RESN [8], and neural network combined with neurofuzzy model [34]. Those methods can obtain generally better results than those obtained with single model, but they are complex, affected by personal experience and easy to be overfitted.To overcome these disadvantages, we propose to approximate the local attractor by using local polynomial coefficient autoregressive prediction (LPP) model. The functional coefficients of the state-dependent autoregressive model can effectively fit nonlinear characteristics of nonlinear time series. By using the polynomial function to approximate the functional coefficients, we can use LPP model to approximate the local attractor. The local attractor can be obtained by utilizing Takens embedding theorem [10]. The kernel prediction can improve the forecast accuracy [37, 38]. We combine the LPP with the kernel technique to develop a kernel LPP (KLPP) model to improve the multistep forecast accuracy. The KLPP model can be applied to choose proper radius of the nearest neighbors by using the kernel function. In this work, we compare our models with those reported in the literature to illustrate the effectiveness of our models in predicting chaotic time series.The paper is organized as follows. In Section2, the concept of the state-dependent autoregressive model, the LPP model, and the KLPP model are proposed, and the optimal parameter set is established by using GDF which is a tool to evaluate the goodness of the model with the chosen number of neighbors [12]. In Section 3, the prediction errors of the generated time series are computed to analyze the prediction effectiveness. Section 3 uses three simulated chaotic systems and one real life time series as examples to evaluate the proposed models. The conclusion of this paper is given in Section 4.


## 2. The LPP Model and KLPP Model

### 2.1. Data Structure in Phase Space Reconstruction

Suppose that we have a scalar chaotic time seriesx1,x2,…,xn. The first step of a local prediction is the phase space reconstruction. According to Takens embedding theory, the phase space can be reconstructed and the construct phase points Xt=(xt,xt+τ,…,xt+(m-1)τ)T, where t=1,2,…,M, and M=n-(m-1)τ. The embedding dimension m and the time delay τ can be obtained by the Cao method [39]. Then, a continued vector mapping F:Rm→Rm or f:Rm→R can be used to describe the unknown evolution from Xt to Xt+1. That is, Xt+1=F(Xt) and xt+1=f(xt).


### 2.2. The LPP Model

The state-dependent autoregressive model can be described as follows:(1)xt=g0Xt-1+∑j=1m‍gjXt-1xt-j,where Xt-1=(xt-1,xt-2,…,xt-m) and functional coefficient gj·, j=0,1,…,m is a continued vector mapping gj:Rm→R.On main purpose of this work is to estimate the mappingf by using this model. The state-dependent autoregressive model in an m-dimensional reconstructed phase space can be given as follows:(2)xt+1=g0Xt+∑j=1m‍gjXtxt-j-1τ,where Xt=(xt,xt-τ,…,xt-(m-1)τ), m is the embedding dimension, and τ is the time delay.We approximategj(Xt) by a polynomial function gj(Xt)=^aj(u)=∑i=0r‍ajiui, where the u as the lag d(∈(τ,2τ,…,(m-1)τ)) variable xt-d and j=0,1,…,m. Then, the LPP model in an m-dimensional phase space can be described as follows:(3)x^t+1=a0u+∑j=1m‍ajuxt-j-1τ,where aj(u)=∑i=0r‍ajiui, the u is the lag d(∈(τ,2τ,…,(m-1)τ)) variable xt-d, and r is the length of the polynomial function. Then, we have(4)x^t+1=∑i=0r‍a0ixt-di+∑j=1m‍∑i=0r‍ajixt-dixt-j-1r.The classical state-dependent autoregressive model has shortcomings such as low accuracy and low processing speed and it is affected by personal experience. The proposed method use the polynomial function to approximate the functional coefficients of the state-dependent autoregressive model, which makes it characterized by parameter model to reduce personal experience and improves processing speed. Also the proposed method is combined with the chaos theory, which makes it obtain high prediction accuracy.


#### 2.2.1. The LPP Estimation and Determination of the Optimal Parameters(q,r)

In order to obtain LPP model as the approximation of the mappingf at the current state point XM in the reconstructed phase space, we select q nearest neighbor points XMj(j=1,2,…,q) by using the Euclidean distances d(i)=∥Xi-XM∥2(i=1,2,…,M-1) and fit a LPP model to predict the future point xMj+(m-1)τ+1 of XMj(j=1,2,…,q). We can estimate the parameters a^ji through the least square equation(5)minaji∈R∑t=1q‍∑j=0m‍xMt+m-1τ+1-∑j=0m‍∑i=0r‍ajixMt+m-1τ-dixMt+j-1τ2,where j=0,1,…,m and i=0,1,…,r.Defineθ^≡(a^00,…,a^0r,…,a^m0,…,a^mr)(m+1)(r+1)×1T and set (6)χ=U1xM1-τU1xM1⋯U1xM1+m-1τU2xM2-τU2xM2⋯U2xM2+m-1τ⋮⋮⋱⋮UqxMq-τUqxMq⋯UqxMq+m-1τq×m+1r+1,Y=xM1+m-1τ+1,…,xMq+m-1τ+1T,Ut=1,xMt+m-1τ-d,…,xMt+m-1τ-dr,xMt-τ≡1t=1,2,…,q.It follows from least squares theory that(7)θ^=χTχ-1χTY.By using (4), the prediction value x^M+(m-1)τ+1 can be calculated. We add x^M+(m-1)τ+1 to the training set and utilize the prediction point X^M+1=(xM+1,xM+1+τ,…,x^M+1+(m-1)τ)T as the last phase point XM+1. Then, we can compute the multistep prediction x^M+(m-1)τ+2 by the same scheme in the forecast of xM+(m-1)τ+1.By using a concept of generalized degrees of freedom (GDF) [12, 13], the error variance can be measured. And the optimal parameters of the prediction model can be obtained by comparing the calculated GDF with different parameters. The unbiased estimators σ^GDF2 can be obtained in every step of the prediction with different parameters as follows:(8)σ^GDF2=Y-Y^Y-Y^Tq-D^,where Y^=χθ^=χ{χTχ}-1χTY, D=tr{H}=tr{χ{χTχ}-1χT}, and q is the number of nearest neighbor points.


### 2.3. The KLPP Model

We combine the LPP with the kernel technique to develop a kernel LPP (KLPP) model. The KLPP model can be applied to choose proper radius of the nearest neighbors by using the kernel function. We suppose that the spatial correlation between phase points could be measured by the kernel function. For the choice of the kernel function, we adopt the Epanechnikov kernel. It is a truncated function and is defined as follows:(9)Kzt=341-zt2+,where zt≡d(t)-min{d(1),d(2),…,d(q)}(orzt≡d(t)). Epanechnikov kernel is the optimal kernel for the kernel density estimation [40] and the local polynomial estimation [41] which minimizes the mean square error. Also, it is a truncated function which will help choose the neighborhood radius. So, we adopt the Epanechnikov kernel for prediction, although there are other kernel functions such as normal function.The estimatorsa^ji is the minimizer of the sum of weighted squares:(10)∑t=1q‍∑j=0m‍xMt+m-1τ+1-∑j=0m‍∑i=0r‍ajixMt+m-1τ-dixMt+j-1τ2Khzt,where j=0,1,…,m, i=0,1,…,r, and Kh(zt)≡(1/h)K(zt/h). The parameter h is used to control the boundary of the nonnegative kernel function which is introduced to choose radius of the nearest neighbors and emphasize neighboring observations around XM when estimating a^ji. It is clear to see that the prediction accuracy of KLPP model is sensitive to the radius of the neighbors. In this study, we choose the parameter h, in terms of GDF, and the optimal value of h is the minimizer of the estimators σ^GDF2.Thus, the estimatorsa^ji can be obtained from least squares theory that(11)θ^=χTWχ-1χTWY,where Y=(xM1+m-1τ+1,…,xMq+(m-1)τ+1)T and W is a q×q diagonal matrix with Kh(zi)≡(1/h)K(zi/h) as its ith diagonal element, χ is a q×[(m+1)(r+1)] matrix with [UixMi-τ,…,UixMi+(m-1)τ] as its ith row, and Ui=[1,xMi+m-1τ-d,…,xMi+(m-1)τ-dr]. And the GDF method chooses h to minimize(12)σ^GDF2h=Y-Y^Y-Y^Tq-D^,where Y^=χθ^=χ{χTWχ}-1χTWY, D=tr{H}=tr{χ{χTWχ}-1χTW}, and q is the number of nearest neighbor points with the Euclidean distance d(i)=∥Xi-XM∥(i=1,2,…,M-1) less than h.Now, we outline the algorithm based on the basic idea of our models.Step 1.
For a scalar time seriesx1,x2,…,xn, the phase space reconstruct with the embedding dimension m and the time delay τ by the autocorrelation function method and the Cao method; that is, Xt=(xt,xt+τ,…,xt+(m-1)τ)T.Step 2.
Compute the Euclidian distancesd(i)=∥Xi-XM∥2(i=1,2,…,M-1) and select the Epanechnikov kernel as the spatial correlation between phase points. Then, the kernel coefficient can be calculated by (9).Step 3.
Identify the parametersa^ji from (11).Step 4.
Calculate theσ^GDF2 in every step of the prediction with different parameters h, r, and a^ji and select the optimal parameters that make σ^GDF2 get minimum.Step 5.
Fit the KLPP with the optimal parameters and calculate the prediction valuex^M+(m-1)τ+1.Some additional remarks are now in order.(i) With the Epanechnikov kernel, letk be from 1 to 15 and hk=hmin+((hmax-hmin)/14)×(k-1) in Step 4, where hmin is the Euclidian distance between the 2(m+1)(r+1)th nearest neighbor and the reference point and hmax is the Euclidian distance between the [3(m+1)(r+1)+10]th nearest neighbor and the reference point.(ii) We acquire one-step and multistep prediction values by using real value and prediction value to the training set as the prediction steps developed, respectively.(iii) To speed up the computation and avoid overfitting, we may letr≤5.(iv) For the algorithm of LPP model, we calculate the parametersa^ji from (7) in Step 3.


## 2.1. Data Structure in Phase Space Reconstruction

Suppose that we have a scalar chaotic time seriesx1,x2,…,xn. The first step of a local prediction is the phase space reconstruction. According to Takens embedding theory, the phase space can be reconstructed and the construct phase points Xt=(xt,xt+τ,…,xt+(m-1)τ)T, where t=1,2,…,M, and M=n-(m-1)τ. The embedding dimension m and the time delay τ can be obtained by the Cao method [39]. Then, a continued vector mapping F:Rm→Rm or f:Rm→R can be used to describe the unknown evolution from Xt to Xt+1. That is, Xt+1=F(Xt) and xt+1=f(xt).


## 2.2. The LPP Model

The state-dependent autoregressive model can be described as follows:(1)xt=g0Xt-1+∑j=1m‍gjXt-1xt-j,where Xt-1=(xt-1,xt-2,…,xt-m) and functional coefficient gj·, j=0,1,…,m is a continued vector mapping gj:Rm→R.On main purpose of this work is to estimate the mappingf by using this model. The state-dependent autoregressive model in an m-dimensional reconstructed phase space can be given as follows:(2)xt+1=g0Xt+∑j=1m‍gjXtxt-j-1τ,where Xt=(xt,xt-τ,…,xt-(m-1)τ), m is the embedding dimension, and τ is the time delay.We approximategj(Xt) by a polynomial function gj(Xt)=^aj(u)=∑i=0r‍ajiui, where the u as the lag d(∈(τ,2τ,…,(m-1)τ)) variable xt-d and j=0,1,…,m. Then, the LPP model in an m-dimensional phase space can be described as follows:(3)x^t+1=a0u+∑j=1m‍ajuxt-j-1τ,where aj(u)=∑i=0r‍ajiui, the u is the lag d(∈(τ,2τ,…,(m-1)τ)) variable xt-d, and r is the length of the polynomial function. Then, we have(4)x^t+1=∑i=0r‍a0ixt-di+∑j=1m‍∑i=0r‍ajixt-dixt-j-1r.The classical state-dependent autoregressive model has shortcomings such as low accuracy and low processing speed and it is affected by personal experience. The proposed method use the polynomial function to approximate the functional coefficients of the state-dependent autoregressive model, which makes it characterized by parameter model to reduce personal experience and improves processing speed. Also the proposed method is combined with the chaos theory, which makes it obtain high prediction accuracy.


### 2.2.1. The LPP Estimation and Determination of the Optimal Parameters(q,r)

In order to obtain LPP model as the approximation of the mappingf at the current state point XM in the reconstructed phase space, we select q nearest neighbor points XMj(j=1,2,…,q) by using the Euclidean distances d(i)=∥Xi-XM∥2(i=1,2,…,M-1) and fit a LPP model to predict the future point xMj+(m-1)τ+1 of XMj(j=1,2,…,q). We can estimate the parameters a^ji through the least square equation(5)minaji∈R∑t=1q‍∑j=0m‍xMt+m-1τ+1-∑j=0m‍∑i=0r‍ajixMt+m-1τ-dixMt+j-1τ2,where j=0,1,…,m and i=0,1,…,r.Defineθ^≡(a^00,…,a^0r,…,a^m0,…,a^mr)(m+1)(r+1)×1T and set (6)χ=U1xM1-τU1xM1⋯U1xM1+m-1τU2xM2-τU2xM2⋯U2xM2+m-1τ⋮⋮⋱⋮UqxMq-τUqxMq⋯UqxMq+m-1τq×m+1r+1,Y=xM1+m-1τ+1,…,xMq+m-1τ+1T,Ut=1,xMt+m-1τ-d,…,xMt+m-1τ-dr,xMt-τ≡1t=1,2,…,q.It follows from least squares theory that(7)θ^=χTχ-1χTY.By using (4), the prediction value x^M+(m-1)τ+1 can be calculated. We add x^M+(m-1)τ+1 to the training set and utilize the prediction point X^M+1=(xM+1,xM+1+τ,…,x^M+1+(m-1)τ)T as the last phase point XM+1. Then, we can compute the multistep prediction x^M+(m-1)τ+2 by the same scheme in the forecast of xM+(m-1)τ+1.By using a concept of generalized degrees of freedom (GDF) [12, 13], the error variance can be measured. And the optimal parameters of the prediction model can be obtained by comparing the calculated GDF with different parameters. The unbiased estimators σ^GDF2 can be obtained in every step of the prediction with different parameters as follows:(8)σ^GDF2=Y-Y^Y-Y^Tq-D^,where Y^=χθ^=χ{χTχ}-1χTY, D=tr{H}=tr{χ{χTχ}-1χT}, and q is the number of nearest neighbor points.


## 2.2.1. The LPP Estimation and Determination of the Optimal Parameters(q,r)

In order to obtain LPP model as the approximation of the mappingf at the current state point XM in the reconstructed phase space, we select q nearest neighbor points XMj(j=1,2,…,q) by using the Euclidean distances d(i)=∥Xi-XM∥2(i=1,2,…,M-1) and fit a LPP model to predict the future point xMj+(m-1)τ+1 of XMj(j=1,2,…,q). We can estimate the parameters a^ji through the least square equation(5)minaji∈R∑t=1q‍∑j=0m‍xMt+m-1τ+1-∑j=0m‍∑i=0r‍ajixMt+m-1τ-dixMt+j-1τ2,where j=0,1,…,m and i=0,1,…,r.Defineθ^≡(a^00,…,a^0r,…,a^m0,…,a^mr)(m+1)(r+1)×1T and set (6)χ=U1xM1-τU1xM1⋯U1xM1+m-1τU2xM2-τU2xM2⋯U2xM2+m-1τ⋮⋮⋱⋮UqxMq-τUqxMq⋯UqxMq+m-1τq×m+1r+1,Y=xM1+m-1τ+1,…,xMq+m-1τ+1T,Ut=1,xMt+m-1τ-d,…,xMt+m-1τ-dr,xMt-τ≡1t=1,2,…,q.It follows from least squares theory that(7)θ^=χTχ-1χTY.By using (4), the prediction value x^M+(m-1)τ+1 can be calculated. We add x^M+(m-1)τ+1 to the training set and utilize the prediction point X^M+1=(xM+1,xM+1+τ,…,x^M+1+(m-1)τ)T as the last phase point XM+1. Then, we can compute the multistep prediction x^M+(m-1)τ+2 by the same scheme in the forecast of xM+(m-1)τ+1.By using a concept of generalized degrees of freedom (GDF) [12, 13], the error variance can be measured. And the optimal parameters of the prediction model can be obtained by comparing the calculated GDF with different parameters. The unbiased estimators σ^GDF2 can be obtained in every step of the prediction with different parameters as follows:(8)σ^GDF2=Y-Y^Y-Y^Tq-D^,where Y^=χθ^=χ{χTχ}-1χTY, D=tr{H}=tr{χ{χTχ}-1χT}, and q is the number of nearest neighbor points.


## 2.3. The KLPP Model

We combine the LPP with the kernel technique to develop a kernel LPP (KLPP) model. The KLPP model can be applied to choose proper radius of the nearest neighbors by using the kernel function. We suppose that the spatial correlation between phase points could be measured by the kernel function. For the choice of the kernel function, we adopt the Epanechnikov kernel. It is a truncated function and is defined as follows:(9)Kzt=341-zt2+,where zt≡d(t)-min{d(1),d(2),…,d(q)}(orzt≡d(t)). Epanechnikov kernel is the optimal kernel for the kernel density estimation [40] and the local polynomial estimation [41] which minimizes the mean square error. Also, it is a truncated function which will help choose the neighborhood radius. So, we adopt the Epanechnikov kernel for prediction, although there are other kernel functions such as normal function.The estimatorsa^ji is the minimizer of the sum of weighted squares:(10)∑t=1q‍∑j=0m‍xMt+m-1τ+1-∑j=0m‍∑i=0r‍ajixMt+m-1τ-dixMt+j-1τ2Khzt,where j=0,1,…,m, i=0,1,…,r, and Kh(zt)≡(1/h)K(zt/h). The parameter h is used to control the boundary of the nonnegative kernel function which is introduced to choose radius of the nearest neighbors and emphasize neighboring observations around XM when estimating a^ji. It is clear to see that the prediction accuracy of KLPP model is sensitive to the radius of the neighbors. In this study, we choose the parameter h, in terms of GDF, and the optimal value of h is the minimizer of the estimators σ^GDF2.Thus, the estimatorsa^ji can be obtained from least squares theory that(11)θ^=χTWχ-1χTWY,where Y=(xM1+m-1τ+1,…,xMq+(m-1)τ+1)T and W is a q×q diagonal matrix with Kh(zi)≡(1/h)K(zi/h) as its ith diagonal element, χ is a q×[(m+1)(r+1)] matrix with [UixMi-τ,…,UixMi+(m-1)τ] as its ith row, and Ui=[1,xMi+m-1τ-d,…,xMi+(m-1)τ-dr]. And the GDF method chooses h to minimize(12)σ^GDF2h=Y-Y^Y-Y^Tq-D^,where Y^=χθ^=χ{χTWχ}-1χTWY, D=tr{H}=tr{χ{χTWχ}-1χTW}, and q is the number of nearest neighbor points with the Euclidean distance d(i)=∥Xi-XM∥(i=1,2,…,M-1) less than h.Now, we outline the algorithm based on the basic idea of our models.Step 1.
For a scalar time seriesx1,x2,…,xn, the phase space reconstruct with the embedding dimension m and the time delay τ by the autocorrelation function method and the Cao method; that is, Xt=(xt,xt+τ,…,xt+(m-1)τ)T.Step 2.
Compute the Euclidian distancesd(i)=∥Xi-XM∥2(i=1,2,…,M-1) and select the Epanechnikov kernel as the spatial correlation between phase points. Then, the kernel coefficient can be calculated by (9).Step 3.
Identify the parametersa^ji from (11).Step 4.
Calculate theσ^GDF2 in every step of the prediction with different parameters h, r, and a^ji and select the optimal parameters that make σ^GDF2 get minimum.Step 5.
Fit the KLPP with the optimal parameters and calculate the prediction valuex^M+(m-1)τ+1.Some additional remarks are now in order.(i) With the Epanechnikov kernel, letk be from 1 to 15 and hk=hmin+((hmax-hmin)/14)×(k-1) in Step 4, where hmin is the Euclidian distance between the 2(m+1)(r+1)th nearest neighbor and the reference point and hmax is the Euclidian distance between the [3(m+1)(r+1)+10]th nearest neighbor and the reference point.(ii) We acquire one-step and multistep prediction values by using real value and prediction value to the training set as the prediction steps developed, respectively.(iii) To speed up the computation and avoid overfitting, we may letr≤5.(iv) For the algorithm of LPP model, we calculate the parametersa^ji from (7) in Step 3.


## 3. Numerical Experiments and Performance Evaluation

Here, we consider three simulated chaotic systems, Lorenz, Mackey-Glass, and Henon and one real life time series, Sunspot time series, as examples to evaluate the proposed models. Later, the results of proposed models are compared with the results reported in the literature for the above examples. For the multistep prediction of chaotic time series, it will produce data overflow when order of polynomial function is larger and normalization can avoid this phenomenon. Thus, those chaotic time series are selected and scaled between[0,1] as follows: (13)x-newt=xt-min⁡xmax⁡x-min⁡x.The prediction errors of the generated time series are computed to analyze the prediction effectiveness and compare the presented models with the results in the literature which are maximum absolute error (MAE), root mean squared error (RMSE), and normalized mean squared error (NMSE). Namely,(14)MAE=maxyt-y^t,RMSE=1n∑t=1n‍yt-y^t2,NMSE=∑t=1nyt-y^t2∑t=1nyt-y¯2.


### 3.1. Performance Analysis in Lorenz Equations

The Lorenz time series can be produced as follows [1]:(15)x˙=σy-x,y˙=r-zx-y,z˙=xy-bz,where σ, r, and b are dimensionless parameters and most commonly selected to be σ=10, r=28, and b=8/3. The standard fourth-order Runge-Kutta method is used to get the Lorenz time series and the x-coordinate is used for prediction. A time series with a length of 2000 is randomly generated. 1800 samples are used for training, and the rest are treated as testing. The results of one-step and multistep predictions for Lorenz time series are shown in Figure 1, in which we use LPP model and KLPP model.Figure 1
Results of the Lorenz time series: (a) and (b) are one-step prediction value and error, respectively; (c) and (d) are multistep prediction value and error, respectively.
(a)
(b)
(c)
(d)From Figure1, it can be seen that, for the one-step prediction, both the prediction values of LPP and KLPP models have small error values. For the multistep prediction, we find out that KLPP model has smaller error values than LPP model and the values of prediction are in good agreement with the real time series.The results of prediction are shown in Figures2 and 3. From Figure 2, we can see that the LPP model parameters are in good agreement with the KLPP model parameters. The polynomial orders mainly take one or two, which avoids overfitting. The lag variables are selected from xt-τ to xt-(m-1)τ at different phase points. The number of nearest neighbor points of both LPP model and KLPP model is chosen in the vicinity of 50, respectively. The radius of the neighbors changes from 0.05 to 0.3. From Figures 1 and 2, for one-step prediction, we can see that both LPP model and KLPP model have similar prediction performance. For the multistep prediction, the prediction values of LPP start diverging significantly from the 120th time step, and the error values are larger than KLPP model. In Figure 3, it can be seen that the parameters almost do not change in account with the parameters in Figure 2. This implies that the kernel can improve the performance of the multistep prediction of chaotic time series.Figure 2
Results of the Lorenz time series one-step prediction: (a) is the polynomial order; (b) is the delay of the lag variable; (c) is the number of nearest neighbor points; (d) is the radius of the neighbors.
(a)
(b)
(c)
(d)Figure 3
Results of the Lorenz time series multistep prediction: (a) is the polynomial order; (b) is the delay of the lag variable; (c) is the number of nearest neighbor points; (d) is the radius of the neighbors.
(a)
(b)
(c)
(d)


### 3.2. Performance Analysis in Mackey-Glass Equations

The Mackey-Glass model has been used in literature as a benchmark model due to its chaotic characteristics [42]. Mackey-Glass time series is generated by the following discrete form:(16)xt+1=1-bxt+axt-τ1+xt-τ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 samples set with length of 2300 from (16). A segment of 1800 samples is used for training, while the remaining part is treated as testing data. Figure 3 compares the real value with prediction value for the rest samples for testing.From Figure4, we can see that both the one-step prediction and the multistep prediction values of LPP and KLPP models have small error values. And the multistep prediction values come up with the real time series with a length of 500. The multistep prediction values are of more accuracy than Lorenz time series, which means the proposed models have different performances in different chaotic time series.Figure 4
Results of Mackey-Glass time series: (a) and (b) are one-step prediction value and error, respectively; (c) and (d) are multistep prediction value and error, respectively.
(a)
(b)
(c)
(d)From Figures5 and 6, we can see that the parameters are almost the same. The polynomial orders mainly take one or two. The number of nearest neighbor points is chosen in the vicinity of 50 and the radius of the neighbors changes from 0.05 to 0.3.Figure 5
Results of Mackey-Glass time series one-step prediction: (a) is the polynomial order; (b) is the delay of the lag variable; (c) is the number of nearest neighbor points; (d) is the radius of the neighbors.
(a)
(b)
(c)
(d)Figure 6
Results of Mackey-Glass time series multistep prediction: (a) is the polynomial order; (b) is the delay of the lag variable; (c) is the number of nearest neighbor points; (d) is the radius of the neighbors.
(a)
(b)
(c)
(d)


### 3.3. Performance Analysis in Henon Equations

Henon mapping is an evident dynamic system [43]. Its equations are written by(17)xt+1=1-axt2+yt,yt+1=bxt,where a=1.4 and b=0.3. The x values with a length of 1800 are generated to reconstruct the state space as training data, and 500 samples are used for testing. Results are shown in Figures 7, 8, and 9.Figure 7
Results of Henon time series: (a) and (b) are one-step prediction value and error, respectively; (c) and (d) are multistep prediction value and error, respectively.
(a)
(b)
(c)
(d)Figure 8
Results of Henon time series one-step prediction: (a) is the polynomial order; (b) is the delay of the lag variable; (c) is the number of nearest neighbor points; (d) is the radius of the neighbors.
(a)
(b)
(c)
(d)Figure 9
Results of Henon time series multistep prediction: (a) is the polynomial order; (b) is the delay of the lag variable; (c) is the number of nearest neighbor points; (d) is the radius of the neighbors.
(a)
(b)
(c)
(d)From Figure7, we can see that the one-step prediction values of LPP and KLPP models have small error values. The multistep prediction values of LPP and KLPP models start diverging significantly from the 25th time step and 35th step, respectively. From Figures 8 and 9, we can see that the polynomial orders mainly take 3, 4, and 5; the number of nearest neighbor points is chosen from 50 to 100. This implies that the proposed models are overfitting.


### 3.4. Performance Analysis in Sunspot Time Series

The Sunspot time series is a good indication of solar activity for solar cycles. The monthly smoothed Sunspot time series has been 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 [9, 10]. Sunspot series from November 1834 to June 2001 (2000 points) are selected and scaled between [0,1]. The first 1000 samples of time series are selected to train the prediction models and the remainder 1000 samples are kept to test the prediction models.From Figure10, it can be seen that, for the one-step prediction, both prediction values of LPP and KLPP models have accuracy prediction. For the multistep prediction, we find out that KLPP model has smaller error values than LPP model, and the values of prediction are in good agreement with the real time series. The prediction values of LPP and KLPP models start diverging significantly from the 10th and 320th time step, respectively. This implies that the kernel can improve the performance of the multistep prediction of chaotic time series.Figure 10
Results of Sunspot time series: (a) and (b) are one-step prediction value and error, respectively; (c) and (d) are multistep prediction value and error, respectively.
(a)
(b)
(c)
(d)From Figures11 and 12, we can see that the polynomial orders mainly take one or two. The number of nearest neighbor points of both LPP and KLPP models is chosen in the vicinity of 20, respectively. The radius of the neighbors changes from 0.05 to 0.3 except for a few phase points. For one-step prediction, the LPP and KLPP models have better prediction performance for Sunspot time series. And for the multistep prediction of Sunspot time series, the KLPP model has more accuracy prediction than LPP model.Figure 11
Results of Sunspot time series one-step prediction: (a) is the polynomial order; (b) is the delay of the lag variable; (c) is the number of nearest neighbor points; (d) is the radius of the neighbors.
(a)
(b)
(c)
(d)Figure 12
Results of Sunspot time series multistep prediction: (a) is the polynomial order; (b) is the delay of the lag variable; (c) is the number of nearest neighbor points; (d) is the radius of the neighbors.
(a)
(b)
(c)
(d)


### 3.5. Results and Discussion

We compare the proposed models with some of the models reported in the literature and the results are shown in Tables1–4. In Table 1, we can see that the proposed models in this paper are better than some of the existing methods. But the best results are of ERNN [9], ARMA-RESN [8], and IEC-LSSVM [20]. This is because these methods have optimized the values for the embedding dimensions for the phase space reconstruction. They also have the advantage of the architectural properties of different models in residual analysis or some models combined with different models which further improve the results. These are also seen for the rest of time series in Tables 2–4. And the data cannot be selected in the same conditions reported in the literature; thus, the conclusions from Tables 1–4 have a few mistakes.Table 1
The comparative results of different models of Lorenz time series.
Prediction model
RMSE
MAE
NMSE
RBF-OLS (2006) [34]
1.41
E
-
09
ERNN (2007) [9]
8.79
E
-
06
9.90
E
-
10
PSO-LSSVM (2009) [25]
1.80
E
-
04
FBMNN (2009) [10]
2.05
E
-
05
HE-NARXNN (2010) [11]
1.08
E
-
04
1.98
E
-
10
IEOP (2010) [33]
1.05
E
-
05
CSAA-SVR (2012) [23]
8.76
E
-
04
COA-SVR (2012) [24]
3.03
E
-
03
ESN (2012) [35]
2.50
E
-
03
HR (2012) [35]
1.50
E
-
03
WESN (2012) [31]
9.58
E
-
04
ARMA-RESN (2013) [8]
3.03
E
-
06
RBFNN (2013) [10]
5.00
E
-
03
LNF (2013) [21]
6.40
E
-
06
IEC-LSSVM (2014) [20]
1.18
E
-
07
4.41
E
-
07
The proposed LPP
2.03
E
-
05
1.04
E
-
04
1.05
E
-
08
The proposed KLPP
2.03
E
-
05
1.17
E
-
04
1.04
E
-
08Table 2
The comparative results of different models of Mackey-Glass time series.
Prediction model
RMSE
MAE
NMSE
RBF-OLS (2006) [34]
1.02
E
-
03
ERNN (2007) [9]
4.20
E
-
05
3.15
E
-
08
PSO-LSSVM (2009) [25]
2.80
E
-
03
FBMNN (2009) [10]
3.02
E
-
04
HE-NARXNN (2010) [11]
3.72
E
-
05
2.70
E
-
08
IEOP (2010) [33]
1.94
E
-
05
ODCF-SVM (2012) [22]
1.60
E
-
02
LNF (2013) [21]
7.90
E
-
04
LSSVR (2013) [27]
5.70
E
-
03
IEC-LSSVM (2014) [20]
2.82
E
-
06
8.32
E
-
06
The proposed LPP
2.71
E
-
05
2.45
E
-
04
1.27
E
-
08
The proposed KLPP
2.46
E
-
05
2.06
E
-
04
1.04
E
-
08Table 3
The comparative results of different models of Henon time series.
Prediction model
RMSE
MAE
NMSE
LSSVM (2005) [10]
4.40
E
-
03
RBF-OLS (2009) [10]
4.13
E
-
02
FBMNN (2009) [10]
2.70
E
-
05
Incremental Algorithm(2012) [36]
7.85
E
-
02
K2 Algorithm (2012) [36]
1.14
E
-
02
The proposed LPP
1.20
E
-
06
1.46
E
-
05
2.70
E
-
12
The proposed KLPP
8.44
E
-
07
6.14
E
-
06
1.33
E
-
12Table 4
The comparative results of different models of Sunspot time series.
Prediction model
RMSE
MAE
NMSE
ERNN (2007) [9]
1.29
E
-
02
2.80
E
-
03
RBF-OLS (2009) [10]
4.60
E
-
02
HE-NARXNN (2010) [11]
1.19
E
-
02
5.90
E
-
04
ARMA-RESN (2013) [8]
15.7937
The proposed LPP
2.70
E
-
02
(
3.93
)
0.48
6.40
E
-
03
The proposed KLPP
3.30
E
-
02
(
4.80
)
0.52
9.50
E
-
03Table2 presents that the proposed models are better than most of the existing methods except for the IEC-LSSVM [20]. In Table 3, the proposed models have the best results. For the real-world time series in Table 4, we can see that the proposed models are better than some of the existing methods, which is similar to the best results reported in the literature.


## 3.1. Performance Analysis in Lorenz Equations

The Lorenz time series can be produced as follows [1]:(15)x˙=σy-x,y˙=r-zx-y,z˙=xy-bz,where σ, r, and b are dimensionless parameters and most commonly selected to be σ=10, r=28, and b=8/3. The standard fourth-order Runge-Kutta method is used to get the Lorenz time series and the x-coordinate is used for prediction. A time series with a length of 2000 is randomly generated. 1800 samples are used for training, and the rest are treated as testing. The results of one-step and multistep predictions for Lorenz time series are shown in Figure 1, in which we use LPP model and KLPP model.Figure 1
Results of the Lorenz time series: (a) and (b) are one-step prediction value and error, respectively; (c) and (d) are multistep prediction value and error, respectively.
(a)
(b)
(c)
(d)From Figure1, it can be seen that, for the one-step prediction, both the prediction values of LPP and KLPP models have small error values. For the multistep prediction, we find out that KLPP model has smaller error values than LPP model and the values of prediction are in good agreement with the real time series.The results of prediction are shown in Figures2 and 3. From Figure 2, we can see that the LPP model parameters are in good agreement with the KLPP model parameters. The polynomial orders mainly take one or two, which avoids overfitting. The lag variables are selected from xt-τ to xt-(m-1)τ at different phase points. The number of nearest neighbor points of both LPP model and KLPP model is chosen in the vicinity of 50, respectively. The radius of the neighbors changes from 0.05 to 0.3. From Figures 1 and 2, for one-step prediction, we can see that both LPP model and KLPP model have similar prediction performance. For the multistep prediction, the prediction values of LPP start diverging significantly from the 120th time step, and the error values are larger than KLPP model. In Figure 3, it can be seen that the parameters almost do not change in account with the parameters in Figure 2. This implies that the kernel can improve the performance of the multistep prediction of chaotic time series.Figure 2
Results of the Lorenz time series one-step prediction: (a) is the polynomial order; (b) is the delay of the lag variable; (c) is the number of nearest neighbor points; (d) is the radius of the neighbors.
(a)
(b)
(c)
(d)Figure 3
Results of the Lorenz time series multistep prediction: (a) is the polynomial order; (b) is the delay of the lag variable; (c) is the number of nearest neighbor points; (d) is the radius of the neighbors.
(a)
(b)
(c)
(d)


## 3.2. Performance Analysis in Mackey-Glass Equations

The Mackey-Glass model has been used in literature as a benchmark model due to its chaotic characteristics [42]. Mackey-Glass time series is generated by the following discrete form:(16)xt+1=1-bxt+axt-τ1+xt-τ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 samples set with length of 2300 from (16). A segment of 1800 samples is used for training, while the remaining part is treated as testing data. Figure 3 compares the real value with prediction value for the rest samples for testing.From Figure4, we can see that both the one-step prediction and the multistep prediction values of LPP and KLPP models have small error values. And the multistep prediction values come up with the real time series with a length of 500. The multistep prediction values are of more accuracy than Lorenz time series, which means the proposed models have different performances in different chaotic time series.Figure 4
Results of Mackey-Glass time series: (a) and (b) are one-step prediction value and error, respectively; (c) and (d) are multistep prediction value and error, respectively.
(a)
(b)
(c)
(d)From Figures5 and 6, we can see that the parameters are almost the same. The polynomial orders mainly take one or two. The number of nearest neighbor points is chosen in the vicinity of 50 and the radius of the neighbors changes from 0.05 to 0.3.Figure 5
Results of Mackey-Glass time series one-step prediction: (a) is the polynomial order; (b) is the delay of the lag variable; (c) is the number of nearest neighbor points; (d) is the radius of the neighbors.
(a)
(b)
(c)
(d)Figure 6
Results of Mackey-Glass time series multistep prediction: (a) is the polynomial order; (b) is the delay of the lag variable; (c) is the number of nearest neighbor points; (d) is the radius of the neighbors.
(a)
(b)
(c)
(d)


## 3.3. Performance Analysis in Henon Equations

Henon mapping is an evident dynamic system [43]. Its equations are written by(17)xt+1=1-axt2+yt,yt+1=bxt,where a=1.4 and b=0.3. The x values with a length of 1800 are generated to reconstruct the state space as training data, and 500 samples are used for testing. Results are shown in Figures 7, 8, and 9.Figure 7
Results of Henon time series: (a) and (b) are one-step prediction value and error, respectively; (c) and (d) are multistep prediction value and error, respectively.
(a)
(b)
(c)
(d)Figure 8
Results of Henon time series one-step prediction: (a) is the polynomial order; (b) is the delay of the lag variable; (c) is the number of nearest neighbor points; (d) is the radius of the neighbors.
(a)
(b)
(c)
(d)Figure 9
Results of Henon time series multistep prediction: (a) is the polynomial order; (b) is the delay of the lag variable; (c) is the number of nearest neighbor points; (d) is the radius of the neighbors.
(a)
(b)
(c)
(d)From Figure7, we can see that the one-step prediction values of LPP and KLPP models have small error values. The multistep prediction values of LPP and KLPP models start diverging significantly from the 25th time step and 35th step, respectively. From Figures 8 and 9, we can see that the polynomial orders mainly take 3, 4, and 5; the number of nearest neighbor points is chosen from 50 to 100. This implies that the proposed models are overfitting.


## 3.4. Performance Analysis in Sunspot Time Series

The Sunspot time series is a good indication of solar activity for solar cycles. The monthly smoothed Sunspot time series has been 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 [9, 10]. Sunspot series from November 1834 to June 2001 (2000 points) are selected and scaled between [0,1]. The first 1000 samples of time series are selected to train the prediction models and the remainder 1000 samples are kept to test the prediction models.From Figure10, it can be seen that, for the one-step prediction, both prediction values of LPP and KLPP models have accuracy prediction. For the multistep prediction, we find out that KLPP model has smaller error values than LPP model, and the values of prediction are in good agreement with the real time series. The prediction values of LPP and KLPP models start diverging significantly from the 10th and 320th time step, respectively. This implies that the kernel can improve the performance of the multistep prediction of chaotic time series.Figure 10
Results of Sunspot time series: (a) and (b) are one-step prediction value and error, respectively; (c) and (d) are multistep prediction value and error, respectively.
(a)
(b)
(c)
(d)From Figures11 and 12, we can see that the polynomial orders mainly take one or two. The number of nearest neighbor points of both LPP and KLPP models is chosen in the vicinity of 20, respectively. The radius of the neighbors changes from 0.05 to 0.3 except for a few phase points. For one-step prediction, the LPP and KLPP models have better prediction performance for Sunspot time series. And for the multistep prediction of Sunspot time series, the KLPP model has more accuracy prediction than LPP model.Figure 11
Results of Sunspot time series one-step prediction: (a) is the polynomial order; (b) is the delay of the lag variable; (c) is the number of nearest neighbor points; (d) is the radius of the neighbors.
(a)
(b)
(c)
(d)Figure 12
Results of Sunspot time series multistep prediction: (a) is the polynomial order; (b) is the delay of the lag variable; (c) is the number of nearest neighbor points; (d) is the radius of the neighbors.
(a)
(b)
(c)
(d)


## 3.5. Results and Discussion

We compare the proposed models with some of the models reported in the literature and the results are shown in Tables1–4. In Table 1, we can see that the proposed models in this paper are better than some of the existing methods. But the best results are of ERNN [9], ARMA-RESN [8], and IEC-LSSVM [20]. This is because these methods have optimized the values for the embedding dimensions for the phase space reconstruction. They also have the advantage of the architectural properties of different models in residual analysis or some models combined with different models which further improve the results. These are also seen for the rest of time series in Tables 2–4. And the data cannot be selected in the same conditions reported in the literature; thus, the conclusions from Tables 1–4 have a few mistakes.Table 1
The comparative results of different models of Lorenz time series.
Prediction model
RMSE
MAE
NMSE
RBF-OLS (2006) [34]
1.41
E
-
09
ERNN (2007) [9]
8.79
E
-
06
9.90
E
-
10
PSO-LSSVM (2009) [25]
1.80
E
-
04
FBMNN (2009) [10]
2.05
E
-
05
HE-NARXNN (2010) [11]
1.08
E
-
04
1.98
E
-
10
IEOP (2010) [33]
1.05
E
-
05
CSAA-SVR (2012) [23]
8.76
E
-
04
COA-SVR (2012) [24]
3.03
E
-
03
ESN (2012) [35]
2.50
E
-
03
HR (2012) [35]
1.50
E
-
03
WESN (2012) [31]
9.58
E
-
04
ARMA-RESN (2013) [8]
3.03
E
-
06
RBFNN (2013) [10]
5.00
E
-
03
LNF (2013) [21]
6.40
E
-
06
IEC-LSSVM (2014) [20]
1.18
E
-
07
4.41
E
-
07
The proposed LPP
2.03
E
-
05
1.04
E
-
04
1.05
E
-
08
The proposed KLPP
2.03
E
-
05
1.17
E
-
04
1.04
E
-
08Table 2
The comparative results of different models of Mackey-Glass time series.
Prediction model
RMSE
MAE
NMSE
RBF-OLS (2006) [34]
1.02
E
-
03
ERNN (2007) [9]
4.20
E
-
05
3.15
E
-
08
PSO-LSSVM (2009) [25]
2.80
E
-
03
FBMNN (2009) [10]
3.02
E
-
04
HE-NARXNN (2010) [11]
3.72
E
-
05
2.70
E
-
08
IEOP (2010) [33]
1.94
E
-
05
ODCF-SVM (2012) [22]
1.60
E
-
02
LNF (2013) [21]
7.90
E
-
04
LSSVR (2013) [27]
5.70
E
-
03
IEC-LSSVM (2014) [20]
2.82
E
-
06
8.32
E
-
06
The proposed LPP
2.71
E
-
05
2.45
E
-
04
1.27
E
-
08
The proposed KLPP
2.46
E
-
05
2.06
E
-
04
1.04
E
-
08Table 3
The comparative results of different models of Henon time series.
Prediction model
RMSE
MAE
NMSE
LSSVM (2005) [10]
4.40
E
-
03
RBF-OLS (2009) [10]
4.13
E
-
02
FBMNN (2009) [10]
2.70
E
-
05
Incremental Algorithm(2012) [36]
7.85
E
-
02
K2 Algorithm (2012) [36]
1.14
E
-
02
The proposed LPP
1.20
E
-
06
1.46
E
-
05
2.70
E
-
12
The proposed KLPP
8.44
E
-
07
6.14
E
-
06
1.33
E
-
12Table 4
The comparative results of different models of Sunspot time series.
Prediction model
RMSE
MAE
NMSE
ERNN (2007) [9]
1.29
E
-
02
2.80
E
-
03
RBF-OLS (2009) [10]
4.60
E
-
02
HE-NARXNN (2010) [11]
1.19
E
-
02
5.90
E
-
04
ARMA-RESN (2013) [8]
15.7937
The proposed LPP
2.70
E
-
02
(
3.93
)
0.48
6.40
E
-
03
The proposed KLPP
3.30
E
-
02
(
4.80
)
0.52
9.50
E
-
03Table2 presents that the proposed models are better than most of the existing methods except for the IEC-LSSVM [20]. In Table 3, the proposed models have the best results. For the real-world time series in Table 4, we can see that the proposed models are better than some of the existing methods, which is similar to the best results reported in the literature.


## 4. Conclusions

In this paper, we apply the polynomial function to approximate the functional coefficients of the state-dependent autoregressive model. Based on the phase space reconstruction, we present a novel local nonlinear model called LPP model for chaotic time series prediction. The LPP model can effectively fit nonlinear characteristics of chaotic time series with simple structure and have better one-step forecasting performance. But, we find out that the LPP model have bad multistep forecasting performance. Then, we propose a kernel LPP model, which applies the kernel technique for the LPP model, so that it may have better multistep forecasting performance than LPP model.The LPP and KLPP models are simple and are not affected by personal experience. Simulated and real data examples illustrate that the proposed models are flexible to analyze complex and multivariate nonlinear structures. The numerical experiments show the KLPP model has more accuracy prediction than LPP model for multistep prediction. We compare LPP and KLPP models with different models and the results show that our models are feasible.



---

*Source: 901807-2015-07-21.xml*