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

A newmethodology, which combines nonparametricmethod based on local functional coefficient autoregressive (LFAR) formwith 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 LFARmodel can effectively fit nonlinear characteristics of chaotic time series by using simple structure and has excellent performance for multistep forecasting.


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][2][3][4][5][6][7][8][9][10][11][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][17][18][19][20], polynomial estimation [21][22][23][24], and neural network (NN) [25][26][27][28][29].In most of the published literature, singlestep 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][30][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.Section 2 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.

The LFAR Model and Estimation Method.
A chaotic time series prediction model for describing evolution from   to  +1 can be written as The continued vector mapping  is the best prediction function in the sense that  minimizes the expected prediction error: The saturated nonparametric function (  ,  − , . . .,  −(−1) ) cannot be estimated with reasonable accuracy due to the curse of dimensionality [30].A LFAR model for chaotic reconstruction data is presented in -dimensional reconstructed phase space.That is, where  ∈ {, 2, . . ., ( − 1)} is the lag of the modeldependent variable  − ,  is the embedding dimension,  is the time delay, and functional coefficients {  }  =1 are continued functions.
The functional coefficients {  }  =1 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 {  }  =1 .Using Taylor's series expansion,   () with -order derivative near the point  0 can be described as follows: Let   () = (− 0 )  ,  = 0, 1, . . ., ,  ∈ ( 0 −ℎ,  0 +ℎ), ℎ > 0 is the bandwidth for controlling the number of the nearest neighbor points, and   =  ()   ( 0 )/!,  = 0, 1, . . ., .Ignoring the higher order infinitesimal, we can approximate   ( − ) by the -dimensional primary functions { 0 ,  1 , . . .,   } as follows: where  is used to replace variable  − and   represent the weight coefficients of primary functions.

Determination of Optimal Parameters.
In the process of modeling LFAR, the parameter set (, , , , , ℎ) needs to be estimated.The embedding dimension  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 ℎ.
Here we have the form of kernel function as follows: Conventionally, we have  = 0.5 and  = 1.Let ℎ → 0; then  ℎ (  ) → 0 for any   > 0. And let ℎ → +∞, the kernel  ℎ (  ) → 1 for any   > 0, which dues to the local linear model.For 0 < ℎ < +∞, the kernel  ℎ (  ) changes from 0 to 1, which means that the weight at the same neighbor point changes from 0 to 1.The parameter  can adjust the weight and the parameter  can adjust the convergence speed of kernel  ℎ (  ) when ℎ → 0 or ℎ → +∞.In this study, we consider simulated chaotic systems and real life time series as examples to investigate a proper parameter ℎ 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 (, , , ℎ).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 (GCV) criterion.Let  and  be two given positive integers and satisfy  > .First we use  subseries with sample size  −  ( = 1, 2, . . ., ) to estimate the unknown coefficient functions and then to compute the multistep forecasting errors of the next part with sample size  based on the estimated model.For example, the data with sample size  − 4 is used to get the estimated model, and the prediction errors for the next  data are computed.Then, the data with sample size  − 3 is used and so on.The average prediction error or the standard prediction errors [31] use the subseries which is given by where  = (1/) ∑ −+ =−+1   .The overall prediction error is given by Fan et al. [34] set  = [0.1]and  = 4, and Meng and Peng [28] set  = [0.03]and  = 3.The selected bandwidth does not critically depend on the choice of  and  as long as  is reasonably large.In practical implementations, we select Meng's method.We select the proper bandwidth by minimizing SAPE(ℎ).The function SAPE(ℎ) is minimized by comparing its values in a finite set of scale parameters in a grid H = {ℎ  | ℎ  = ℎ 0   ,  = 1, . . ., ,  > 1, ℎ 0 > 0}.We can obtain different SAPE(ℎ) with different , , and .It is also very important to choose the parameter set (, , ).Finally we can determinate the optimal parameter set ( opt ,  opt ,  opt , ℎ opt ).
In Section 2.2, the single-step prediction is given.For multistep forecasting, as we say in Section 1, there are two possible approaches.The iterative -step prediction is to add x+1 , x+2 , . . ., x+−1 to the training set and utilizes the LFAR model iteratively.The direct -step prediction is to fit The optimal parameter set ( opt ,  opt , ℎ opt ) needs to be selected to minimize the SAPE 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.

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: 3.1.Prediction of Lorenz System.Lorenz time series can be produced as follows (Lorenz, 1963, see [1]) where , , and  are commonly selected as  = 10,  = 28, and  = 8/3.The standard fourth-order Runge-Kutta method is used to get Lorenz time series, and the -coordinate is  From Figure 1, 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 Figures 2, 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 SAPE less than  − 6 in Figure 2. From Figure 2, we can see that the LFAR model's SAPE 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   to  −(−1) at different phase points, the weight parameter ℎ changes with other parameters, and the number of nearest neighbor points is chosen in the vicinity of 2( + 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.
In Figure 5, 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.

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: where  = 0.2,  = 0.1, and  = 17 and initial conditions ()| ≤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.Figure 6 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.
The results of prediction are shown in Figures 7, 8,  and 9. From Figure 7, we can see that the LFAR model's SAPE with different parameters are similar.Figure 8 shows that the functional coefficients' order mainly takes four, the lag variable is selected from   to  −(−1) at different phase points, the weight parameter ℎ changes with other parameters, and the number of nearest neighbor points is chosen in the vicinity of 2( + 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.
In Figure 10, 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.

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.[30] 9.98 − 07 PSORBF (2014) [30] 6.92 − 08 The proposed LFAR (test data = 450) Results shown in 11 and 12. From Figure 11, it can be seen that the values LFAR model have accuracy of prediction.
From Figure 12, for the multistep prediction, we find that the prediction values of LFAR model start diverging significantly from the 60th 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.

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 The comparative results are shown in Tables 1-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     [33] 4.85 − 02 UKF (2013) [33] 1.89 − 02 GUKF (2013) [33] 2.6 − 03 GPF (2013) [33] 3.97 − 05 GGPF (2013) [33] 5  [31] 0.6915 FTLRNN (2010) [31] 0.3210 SOFM (2010) [31] 0.7085 The proposed LFAR (test data = 300) that the proposed models are better than most of the existing methods only except the CERNN [29].

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 in -dimensional phase space with a time delay  by using chaos theory.Secondly, the neighbor points are selected by using local method in the -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.

Figure 1 :
Figure 1: Results of Lorenz time series: (a) multistep prediction errors and (b) the original values and multistep prediction values, respectively.

Figure 2 :
Figure 2: Results of Lorenz time series' multistep prediction: the SAPE of parameter set.

Figure 3 :
Figure 3: Results of Lorenz time series' multistep prediction: (a) the functional coefficients' order; (b) the delay of the lag variable (coordinate value  is equivalent to ( − 1)); (c) the weight parameter; (d) the number of nearest neighbor points (equivalent to number-2( + 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.

Figure 4 :Figure 5 :
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.

Figure 8 :
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.

Figure 9 :
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 set SAPE = 10.

Figure 11 :
Figure 11: Results Sunspot time series: (a) single-step error and (b) the original and single-step prediction values, respectively.

Figure 12 :
Figure 12: Results of Sunspot time series: (a) multistep prediction error and (b) the original and multistep prediction values, respectively.
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 dimension  opt and the time delay  opt ; Compute the Euclidean distances (⋅) between phase points in the reconstructed phase space; For  ∈ t % { t is unknown data.}Compute prediction values by using the optimal parameters (ℎ opt ,  opt ,  opt ,  opt ); Add prediction values to the training set; Compute Euclidean distances between the new phase point and others; End Algorithm 2: Multistep prediction algorithm.