# A New Adaptive Local Linear Prediction Method and Its Application in Hydrological Time Series

Authors: Dunxian She; Xiaohua Yang
Journal: Mathematical Problems in Engineering
 (2010)

---

## Abstract

The embedding dimension and the number of nearest neighbors are very important parameters in the prediction of a chaotic time series. In order to reduce the uncertainties in the determination of the forgoing two parameters, a new adaptive local linear prediction method is proposed in this study. In the new method, the embedding dimension and the number of nearest neighbors are combined as a parameter set and change adaptively in the process of prediction. The generalized degree of freedom is used to help select the optimal parameters. Real hydrological time series are taken to examine the performance of the new method. The prediction results indicate that the new method can choose the optimal parameters of embedding dimension and the nearest neighbor number adaptively in the prediction process. And the nonlinear hydrological time series perhaps could be modeled better by the new method.


---

## Body

## 1. Introduction

The global and regional climates have already begun changing [1], and meteorological-driven processes have been studied by some researchers just as the study in the signal analysis and other fields [2–8]. Many hydrological processes, such as runoff, are usually nonlinear, complex, dynamic processes because of the involved physical process and the considerable spatial and temporal variability [9]. The simulation of the nonlinear hydrological time series was turned out to be very difficult with the traditional deterministic mathematic models. However, the emergence of chaos theory provides a new way to study this kind of highly complex system and makes it possible to extract deterministic regulation from the seemingly disordered hydrological phenomenon. Chaos theory is an important part of nonlinear science and some scientists have done research on its theories and applications [10, 11]. From the last decades, a series of theories and methods identifying chaotic essences in dynamic systems have gradually been established [12–17]. The first application of chaos theory in hydrological processes could be traced back to the analysis of a series of 1008 monthly rainfall values recorded in Nauru Island by Hense in 1987 [18], and then the chaos theory became a more and more reliable tool for the study of hydrological processes.With the development of chaos theory and research on its application technique, many methods have been proposed to predict chaotic time series, which can broadly be divided into two main categories: global method [19] and local method [20, 21]. The global method is made an attempt to approximate the whole time series on all attractors and seek a function valid at every point. It is obvious that this kind of method has the disadvantage that if new information is added into the model, the parameters may be changed and a lot of time will be wasted in the parameter estimation. However, the local method overcomes this drawback by building model only on the local attractor and utilizing only part of the past information. Farmer and Sidorowich [20] have already proved that the local prediction method is better than the global prediction method. On the other hand, by using some combining techniques, the forecast accuracy can be improved both in the global method and in the local method [21].In the process of the local prediction, the phase space reconstruction is the first step and three parameters (the embedding dimensionm, the time delay τ, and the number of the nearest neighbors q) should be determined. The studies of the nonlinear predictions indicate that the selection of delay time does not affect the reconstructed attractor reflecting the dynamics of the system unambiguously; so the key problem is the determination of the optimal m and q for prediction. (Usually, the number of the nearest neighbors q is taken greater than m in general condition.)Although there are many discussions on how to determine the optimal embedding dimension, three basic methods presented below are most usually used. The first kind of methods is based on calculating some geometrical invariants of the attractor, such as Grassberger-Procaccia method (G-P method) [16]. By increasing the embedding dimension, the minimum embedding dimension can be selected when the value of the invariant becomes saturated. The typical problems of this method are time-consuming for computation, unsuited for a short time series, and certainly subjective and sensitive to noise. The second kind of methods, such as the false nearest neighbors method (FNN) [22, 23] and Aleksić method [24], is based on the theory of false neighbors that the far apart points in the original phase space will moving closer together in the reconstruction phase space when a too low embedding dimension is selected. But the criterion is subjective in the judgment of the false neighbors, and Cao has improved the method [25]. The third kind of methods, such as principal component analysis (PCA) [26], singular-spectrum analysis (SSA) [27] and singular-value decomposition (SVD) [28], is based on a singular value decomposition proposed by Broomhead and King [26]. But the singular value decomposition is essentially a linear method based on the covariance matrix which reflects the linear dependence [29]. This approach is also subjective to some extent that the number of large singular values may depend on the details of the embedding process and the accuracy of data. Numerical experience has led several researchers [28, 30] to suspect the authenticity of this method in the analysis of nonlinear time series [31]. All of the methods mentioned above have some shortcomings like inapplicability to a short time series, more or less subjective, sensitivity to noise, and so forth [31].The uncertainties in the determination of the embedding dimension and the number of the nearest neighbors can affect the forecast precision in the process of the traditional local prediction. So as to improve the prediction accuracy, the problem that arises then is how to reduce the parameter uncertainties in the prediction method. As we know, the time series data are continually updated as the prediction process proceeded; the parameter estimated by the phase space reconstruction using the original time sequence may not be able to reconstruct the chaotic attractor of the new sequence. In order to reduce the uncertainties caused by the time sequence data updating, the embedding dimension and the number of nearest neighbors should be changed adaptively in the forecast procedure. For the purpose of getting a better prediction results, Jayawardena et al. [32] proposed a method using the generalized degrees of freedom (GDF) to determine the optimal number of the neighbors for the prediction. But in his method, the embedding dimensions take the same values in the whole forecast process like the traditional local prediction; so this method still needs to be improved.In this study, a parameter set(m,q) is established on the basis of the two uncertain parameters, that is, embedding dimension m and the number of nearest neighbors q. Then a new adaptive local prediction method is proposed, in which m and q are not fixed as certain values but are changed as the prediction steps developed. In the selection of the optimal parameter set (m,q) for each prediction step, the generalized degrees of freedom (GDF) are used and different error variances are calculated under different combinations of (m,q). The optimal parameter set (m,q) is chosen when the variance obtains the smallest value. In order to examine the validity of the new method, some real hydrological time series are used.


## 2. The New Local Linear Prediction Method

### 2.1. The Phase Space Reconstruction

For a scalar time seriesx1,x2,…,xn, the multidimensional phase-space Yi can be reconstructed by the Takens embedding theory [15], according to(2.1)Yi=(xi,xi+τ,…,xi+(m-1)τ),
where i=1,2,…,M, M is the total points number of the phase-space, and M=n-(m-1)τ, τ is the delay time, m is the dimension of the vector Yi, called as embedding dimension, and n is the length of the time series.


#### 2.1.1. The Determination of Delay Time

Many methods have been developed in estimating the delay time. And the autocorrelation function [15, 16, 33, 34], the most widely used tool in determining the delay time, is employed in this study. The autocorrelation function can be described as(2.2)C(t)=1n-t∑i=1n-t(xi-μ)(xi+t-μ)σ2(x),
where C(t) is the autocorrelation coefficient, t is the lag time, and μ and σ are the mean and standard variation of the time series, respectively. The delay time τ is always selected when the autocorrelation coefficient has dropped to 1-1/e of its initial value (e is the base of natural logarithm).


#### 2.1.2. The Determination of Embedding Dimension

Among a number of embedding dimension calculation methods, the false nearest neighbors (FNN) method [22] is most widely used. Cao has proposed a new method on the basis of FNN method. The Cao method [25] can be described as follows. Assume that(2.3)a(i,m)=∥Yi(m+1)-Yj(m+1)∥∥Yi(m)-Yj(m)∥(i=1,2,…,n-mτ),
where ∥·∥ is some measurements of Euclidian distance,  Yi(m+1) is the ith phase points in the (m+1)-dimension reconstructed phase space, and Yj(m)(j=1,2,…,q) is the nearest neighbor of Yi(m). The mean value of all a(i,m) is (2.4)E(m)=1n-mτ∑i=1n-mτa(i,m).When the embedding dimension changes fromm to m+1, E1(m) can be defined as (2.5)E1(m)=E(m+1)E(m).IfE1(m) stops changing when m is greater than some value m0, then m0+1 can be taken as the minimum embedding dimension to reconstruct the phase space.


### 2.2. The Traditional Local Linear Prediction Method

Many local prediction methods have been developed during the last decades, and the local linear prediction method is considered in this study. The first step of the local linear prediction is to find the nearest neighbor points of the current phase pointYM in the reconstructed phase space. The Euclidean distances d(i)=∥Yi-YM∥(i=1,2,…,M-1) between the current vector YM and its M-1 preceding delay vectors Yi(i=1,2,…,M-1) will be calculated and then q nearest neighbor points YMk (k=1,2,…,q, q is taken value great than m in general condition) will be selected.The local linear prediction model in anm-dimensional reconstructed phase space is an autoregressive model and the prediction value, which is a linear superposition of the m elements in the delay vector YM, can be given as follows:(2.6)x̂(n+1)=AX1=a0+∑j=1majxM+(j-1)τ,
where A=[a0,a1,a2,…,am] is a coefficient vector that needs to be determined, and
(2.7)X1=[1,YM]T=[1,xM,xM+τ,…,xM+(m-1)τ]T.Deterministic predictions assume that if the phase space pointYt was similar to the current point YM, then the future point Yt+1 will also be close to the future point YM+1. The coefficient vector A can be estimated by the current phase vector YM and its q nearest neighbor points YMk(k=1,2,…,q) through the following equation:(2.8)AX=Y,
where Y=[xM1+(m-1)τ+1,xM2+(m-1)τ+1,…,xMq+(m-1)τ+1]T is the next series value of the points YMk(k=1,2,…,q), and(2.9)X=[11⋯1xM1xM2⋯xMq⋮⋮⋱⋮xM1+(m-1)τxM2+(m-1)τ⋯xMq+(m-1)τ].Because the value ofX and Y is known, so the estimation of coefficient vector A can be obtained by the least squares method as A=YXT(XXT)-1 from (2.8), and then by using (2.6), the prediction value x̂(n+1) can be calculated. The new prediction value x̂(n+1) is added to the original time series as the prediction steps developed, and the last phase point is YM+1 by the phase space reconstruction now. Following the same scheme in the forecast of x̂(n+1) and reestimating the coefficient vector A, x̂(n+2) can be computed.


### 2.3. The Determination of the Optimal Parameters(m,q)

In the regression analysis, the degrees of freedom are often used as a model complexity measure in various model selection criteria, such as MallowsCp, Akaike information criteria (AIC), and Bayesian information criterion (BIC). Yet these model selection criterion are asymptotic in nature and do not take into account the modeling procedure which can often be very complex [32]. So Ye [35] developed a concept of generalized degrees of freedom (GDF) that is applicable for evaluation of the model selection.The GDF is defined as the sum of the sensitivities of each fitted value of the model to perturbations in the corresponding observed value [32]. It is nonasymptotic in nature and thus is free of the sample-size constraint [35]. In the process of chaotic local linear prediction, the GDF can be viewed as the cost of the modeling process. Considering the uncertainty in the determination of the optimal embedding dimension and the number of nearest neighbors for the prediction, different unbiased estimates of the error variance can be obtained under different parameter sets (m,q). A better (m,q) for the model prediction is the one that has a smaller error variance. So the optimal parameters (m,q) can be selected by comparing the different error variances that are calculated. In every step of the prediction, the GDF can be calculated and the optimal parameters (m,q) can be selected; so the GDF can provide a novel way to guide the (m,q) changing adaptively.Using the local linear prediction method introduced earlier, the future values of the data series can be obtained. Here we will show how to choose the optimal(m,q) for the prediction by using the GDF method.From the foregoing description, the estimate of the coefficient vectorA in (2.8) can also be given as follows:(2.10)Y=AX+ε,
where Y, A, X are the same in (2.8) and ε is a row vector of error values.Using the least squares method,A and μY=AX (μY is the mean vector of Y) can be estimated as(2.11)A=YXT(XXT)-1,(2.12)μ̂Y=AX=YXT(XXT)-1X.The sum of squared residuals for a fixed embedding dimension RSS is then expressed as(2.13)RSS=(Y-μ̂Y)(Y-μ̂Y)T.Then an unbiased estimation of the error varianceσGDF2 is given as (2.14)σGDF2=RSSq-D=(Y-μ̂Y)(Y-μ̂Y)Tq-D,
where q is the number of nearest neighbors. D is the GDF and can be estimated by D=tr(H)=∑ihii, H=(hii)(m+1)×(m+1)=X(XTX)-1XT.This provides a tool to evaluate the goodness of the model with the chosen(m,q). With different values of the parameter set (m,q), the matrix X and fitted vector functions μ̂Y will be different, then different error variances can be obtained. An optimal (m,q) for the prediction is the one that has the smallest variance.


### 2.4. The New Local Linear Prediction Method

By comparing the estimations of error variancesσGDF2 for different (m,q), the optimal parameter set can be selected. It is then used for prediction. The new local linear method proposed in this study is described as follows.Step 1.
Phase space reconstruction. Determine the embedding dimensionmc and the time delay τ for the original time series.Step 2.
Calculate the error variances. Let the embedding dimensionm change from mmin to mmax, and the number of the nearest neighbor change from qmin to qmax (In this study, mmin and mmax are selected as 2 and mc, resp. This is not the only choice. qmin and qmax are taken as 2m+1 and 2m+10, the same as Jayawardena’s method [32].) (mmax-mmin+1)×(qmax-qmin+1) local linear models can be obtained. For each model, the error variance can be estimated using (2.13) and (2.14).Step 3.
Choose the optimal(mopt,qopt) which has the minimum variance.Step 4.
Use(mopt,qopt) to reconstruct the original time series and construct a local model to predict the next value x̂(n+1).Step 5.
Thenx̂(n+2) can be computed by the same scheme following the last point; now it is YM+1.Some real hydrological time series are chosen to examine the performance of the new local linear method (NLLP-2), and the new method is compared with the traditional local linear prediction method (TLLP) and Jayawardena’s method (NLLP-1) [32].


## 2.1. The Phase Space Reconstruction

For a scalar time seriesx1,x2,…,xn, the multidimensional phase-space Yi can be reconstructed by the Takens embedding theory [15], according to(2.1)Yi=(xi,xi+τ,…,xi+(m-1)τ),
where i=1,2,…,M, M is the total points number of the phase-space, and M=n-(m-1)τ, τ is the delay time, m is the dimension of the vector Yi, called as embedding dimension, and n is the length of the time series.


### 2.1.1. The Determination of Delay Time

Many methods have been developed in estimating the delay time. And the autocorrelation function [15, 16, 33, 34], the most widely used tool in determining the delay time, is employed in this study. The autocorrelation function can be described as(2.2)C(t)=1n-t∑i=1n-t(xi-μ)(xi+t-μ)σ2(x),
where C(t) is the autocorrelation coefficient, t is the lag time, and μ and σ are the mean and standard variation of the time series, respectively. The delay time τ is always selected when the autocorrelation coefficient has dropped to 1-1/e of its initial value (e is the base of natural logarithm).


### 2.1.2. The Determination of Embedding Dimension

Among a number of embedding dimension calculation methods, the false nearest neighbors (FNN) method [22] is most widely used. Cao has proposed a new method on the basis of FNN method. The Cao method [25] can be described as follows. Assume that(2.3)a(i,m)=∥Yi(m+1)-Yj(m+1)∥∥Yi(m)-Yj(m)∥(i=1,2,…,n-mτ),
where ∥·∥ is some measurements of Euclidian distance,  Yi(m+1) is the ith phase points in the (m+1)-dimension reconstructed phase space, and Yj(m)(j=1,2,…,q) is the nearest neighbor of Yi(m). The mean value of all a(i,m) is (2.4)E(m)=1n-mτ∑i=1n-mτa(i,m).When the embedding dimension changes fromm to m+1, E1(m) can be defined as (2.5)E1(m)=E(m+1)E(m).IfE1(m) stops changing when m is greater than some value m0, then m0+1 can be taken as the minimum embedding dimension to reconstruct the phase space.


## 2.1.1. The Determination of Delay Time

Many methods have been developed in estimating the delay time. And the autocorrelation function [15, 16, 33, 34], the most widely used tool in determining the delay time, is employed in this study. The autocorrelation function can be described as(2.2)C(t)=1n-t∑i=1n-t(xi-μ)(xi+t-μ)σ2(x),
where C(t) is the autocorrelation coefficient, t is the lag time, and μ and σ are the mean and standard variation of the time series, respectively. The delay time τ is always selected when the autocorrelation coefficient has dropped to 1-1/e of its initial value (e is the base of natural logarithm).


## 2.1.2. The Determination of Embedding Dimension

Among a number of embedding dimension calculation methods, the false nearest neighbors (FNN) method [22] is most widely used. Cao has proposed a new method on the basis of FNN method. The Cao method [25] can be described as follows. Assume that(2.3)a(i,m)=∥Yi(m+1)-Yj(m+1)∥∥Yi(m)-Yj(m)∥(i=1,2,…,n-mτ),
where ∥·∥ is some measurements of Euclidian distance,  Yi(m+1) is the ith phase points in the (m+1)-dimension reconstructed phase space, and Yj(m)(j=1,2,…,q) is the nearest neighbor of Yi(m). The mean value of all a(i,m) is (2.4)E(m)=1n-mτ∑i=1n-mτa(i,m).When the embedding dimension changes fromm to m+1, E1(m) can be defined as (2.5)E1(m)=E(m+1)E(m).IfE1(m) stops changing when m is greater than some value m0, then m0+1 can be taken as the minimum embedding dimension to reconstruct the phase space.


## 2.2. The Traditional Local Linear Prediction Method

Many local prediction methods have been developed during the last decades, and the local linear prediction method is considered in this study. The first step of the local linear prediction is to find the nearest neighbor points of the current phase pointYM in the reconstructed phase space. The Euclidean distances d(i)=∥Yi-YM∥(i=1,2,…,M-1) between the current vector YM and its M-1 preceding delay vectors Yi(i=1,2,…,M-1) will be calculated and then q nearest neighbor points YMk (k=1,2,…,q, q is taken value great than m in general condition) will be selected.The local linear prediction model in anm-dimensional reconstructed phase space is an autoregressive model and the prediction value, which is a linear superposition of the m elements in the delay vector YM, can be given as follows:(2.6)x̂(n+1)=AX1=a0+∑j=1majxM+(j-1)τ,
where A=[a0,a1,a2,…,am] is a coefficient vector that needs to be determined, and
(2.7)X1=[1,YM]T=[1,xM,xM+τ,…,xM+(m-1)τ]T.Deterministic predictions assume that if the phase space pointYt was similar to the current point YM, then the future point Yt+1 will also be close to the future point YM+1. The coefficient vector A can be estimated by the current phase vector YM and its q nearest neighbor points YMk(k=1,2,…,q) through the following equation:(2.8)AX=Y,
where Y=[xM1+(m-1)τ+1,xM2+(m-1)τ+1,…,xMq+(m-1)τ+1]T is the next series value of the points YMk(k=1,2,…,q), and(2.9)X=[11⋯1xM1xM2⋯xMq⋮⋮⋱⋮xM1+(m-1)τxM2+(m-1)τ⋯xMq+(m-1)τ].Because the value ofX and Y is known, so the estimation of coefficient vector A can be obtained by the least squares method as A=YXT(XXT)-1 from (2.8), and then by using (2.6), the prediction value x̂(n+1) can be calculated. The new prediction value x̂(n+1) is added to the original time series as the prediction steps developed, and the last phase point is YM+1 by the phase space reconstruction now. Following the same scheme in the forecast of x̂(n+1) and reestimating the coefficient vector A, x̂(n+2) can be computed.


## 2.3. The Determination of the Optimal Parameters(m,q)

In the regression analysis, the degrees of freedom are often used as a model complexity measure in various model selection criteria, such as MallowsCp, Akaike information criteria (AIC), and Bayesian information criterion (BIC). Yet these model selection criterion are asymptotic in nature and do not take into account the modeling procedure which can often be very complex [32]. So Ye [35] developed a concept of generalized degrees of freedom (GDF) that is applicable for evaluation of the model selection.The GDF is defined as the sum of the sensitivities of each fitted value of the model to perturbations in the corresponding observed value [32]. It is nonasymptotic in nature and thus is free of the sample-size constraint [35]. In the process of chaotic local linear prediction, the GDF can be viewed as the cost of the modeling process. Considering the uncertainty in the determination of the optimal embedding dimension and the number of nearest neighbors for the prediction, different unbiased estimates of the error variance can be obtained under different parameter sets (m,q). A better (m,q) for the model prediction is the one that has a smaller error variance. So the optimal parameters (m,q) can be selected by comparing the different error variances that are calculated. In every step of the prediction, the GDF can be calculated and the optimal parameters (m,q) can be selected; so the GDF can provide a novel way to guide the (m,q) changing adaptively.Using the local linear prediction method introduced earlier, the future values of the data series can be obtained. Here we will show how to choose the optimal(m,q) for the prediction by using the GDF method.From the foregoing description, the estimate of the coefficient vectorA in (2.8) can also be given as follows:(2.10)Y=AX+ε,
where Y, A, X are the same in (2.8) and ε is a row vector of error values.Using the least squares method,A and μY=AX (μY is the mean vector of Y) can be estimated as(2.11)A=YXT(XXT)-1,(2.12)μ̂Y=AX=YXT(XXT)-1X.The sum of squared residuals for a fixed embedding dimension RSS is then expressed as(2.13)RSS=(Y-μ̂Y)(Y-μ̂Y)T.Then an unbiased estimation of the error varianceσGDF2 is given as (2.14)σGDF2=RSSq-D=(Y-μ̂Y)(Y-μ̂Y)Tq-D,
where q is the number of nearest neighbors. D is the GDF and can be estimated by D=tr(H)=∑ihii, H=(hii)(m+1)×(m+1)=X(XTX)-1XT.This provides a tool to evaluate the goodness of the model with the chosen(m,q). With different values of the parameter set (m,q), the matrix X and fitted vector functions μ̂Y will be different, then different error variances can be obtained. An optimal (m,q) for the prediction is the one that has the smallest variance.


## 2.4. The New Local Linear Prediction Method

By comparing the estimations of error variancesσGDF2 for different (m,q), the optimal parameter set can be selected. It is then used for prediction. The new local linear method proposed in this study is described as follows.Step 1.
Phase space reconstruction. Determine the embedding dimensionmc and the time delay τ for the original time series.Step 2.
Calculate the error variances. Let the embedding dimensionm change from mmin to mmax, and the number of the nearest neighbor change from qmin to qmax (In this study, mmin and mmax are selected as 2 and mc, resp. This is not the only choice. qmin and qmax are taken as 2m+1 and 2m+10, the same as Jayawardena’s method [32].) (mmax-mmin+1)×(qmax-qmin+1) local linear models can be obtained. For each model, the error variance can be estimated using (2.13) and (2.14).Step 3.
Choose the optimal(mopt,qopt) which has the minimum variance.Step 4.
Use(mopt,qopt) to reconstruct the original time series and construct a local model to predict the next value x̂(n+1).Step 5.
Thenx̂(n+2) can be computed by the same scheme following the last point; now it is YM+1.Some real hydrological time series are chosen to examine the performance of the new local linear method (NLLP-2), and the new method is compared with the traditional local linear prediction method (TLLP) and Jayawardena’s method (NLLP-1) [32].


## 3. Application in Hydrological Time Series

### 3.1. Study Area and Data Description

The real hydrological time series in this study are chosen as the daily discharges of Peschanaya at Tochil'noye (52.17°N and 85.17°E, basin area, 4720 km2) in Russian for the period January 1936–December 1980, and the daily discharges of Tim at Napas (59.85°N and 81.95°E, basin area, 24500 km2) in Russian for the period January 1953–December 1994. Some important statistical values are shown in Table 1.Table 1
Statistics indices of daily discharge time series at Peschanaya and Tim.
Statistical indicesPeschanayaTimTotal number of data1643715340Mean value (m3/s)32.3357192.5562Standard of deviation (m3/s)44.6106215.4076Maximum value (m3/s)661.00001360.0000Minimum value (m3/s)1.740028.4000Coefficient of variation1.37961.1187For the daily discharge from Peschanaya, the time series is divided into two data sets, the first training data set is the data during the period January 1936–December 1979, and the second prediction data set is the data during the period January 1, 1980–February 9, 1980. In the same way, the first training data set for Tim is the data during the period January 1953–December 1993, the second prediction data set is the data during the period January 1, 1994–February 9, 1994. The prediction lead time for both of the time series is 40 days.


### 3.2. Phase Space Reconstruction

In this study, the delay time and the embedding dimension of the above two daily discharge time series are determined by the autocorrelation function method and the Cao method, respectively. The changes of the autocorrelation functions for both of the time series are presented in Figure1 and the values of E1(m) are shown in Figure 2.Figure 1
Autocorrelation functions for daily discharge time series observed at Peschanaya and Tim. (The red line presents the situation that the autocorrelation function value is1-1/e of its initial value.)Figure 2
The values ofE1(m) for daily discharge time series observed at Peschanaya and Tim.The delay time is selected when the autocorrelation coefficient has dropped to1-1/e of its initial value. From Figure 1, the delay time for the Peschanaya and Tim can be determined as 9 and 20, respectively. (The red line in Figure 1 presents the situation that the autocorrelation function value is 1-1/e of its initial value.)From Figure2, the value of E1(m) becomes saturated from m0=5; so the embedding dimension can be determined as m0+1=6 for Peschanaya.  In the same way, the embedding dimension for Tim can be determined as 6 too.The third parameter in the local linear prediction of the daily discharge time series is the number of the nearest neighbors. In this study, only the traditional prediction method (TLLP) needs to give this parameter as a prior; the rest two methods (NLLP-1 and NLLP-2) can obtain the number of nearest neighbors adaptively in the prediction procedure.


### 3.3. The Prediction Results and Analysis

Both of the time series were predicted by TLLP, NLLP-1, and NLLP-2, respectively. In order to obtain the best forecast results with the TLLP, the number of nearest neighborsq for TLLP is changed from m+1 to m+19. By comparing the prediction results under different q, the optimal forecast results by TLLP in different q can be selected. The prediction results by TLLP were shown in Table 2.Table 2
The prediction results for the daily discharge time series observed at Peschanaya and Tim by TLLP method.
Number of nearest neighborsPeschanayaTimMAE (m3/s)RMSE (m3/s)CCMAE (m3/s)RMSE (m3/s)CCm+19.6882 × 1022.1893 × 1030.69661.6088 × 10104.5795 × 10100.2904m+24.2862 × 1051.4036 × 1060.554847.922478.86640.0600m+35.31516.20560.42938.1494 × 1041.4928 × 1050.5736m+41.47791.83100.95038.4746 × 10113.0884 × 10120.3439m+52.44154.28530.76962.0832 × 10179.1309 × 10170.2899m+60.84281.19070.83072.1107 × 10117.5627 × 10110.3491m+70.47810.68780.909222.463624.4395-0.3345m+80.53620.69170.889119.936520.9531-0.5286m+90.55990.76380.812719.111120.0202-0.5056m+100.49780.68150.923611.122314.21230.7822m+110.74450.96710.839918.263824.18650.8162m+122.62595.6787-0.656813.690517.09750.8609m+132.77426.0485-0.657411.015013.61990.8799m+141.95363.7148-0.684611.753514.73290.8702m+152.68775.8700-0.650511.726014.49790.8819m+161.83583.4853-0.66127.65179.47250.9089m+171.87633.3413-0.721611.948614.70640.8939m+182.51085.5012-0.639413.278116.13630.8916m+192.02664.0896-0.651614.820817.90760.8850The accuracy of prediction is evaluated by three measurement indices in this study, which are theMean absolute error (MAE), Root mean square error (RMSE), and Correlation coefficient (CC). The definitions of the three indices can be given as follows:(3.1)MAE=1n∑i=1n|yi-xi|,RMSE=1n∑i=1n(yi-xi)2,CC=∑i=1n(xi-x̅)(yi-y̅)∑i=1n(xi-x̅)2∑i=1n(yi-y̅)2,
where n is the number of the time series under investigation. The yi and xi are the predicted values and observed values in the time series, respectively, and x̅ and y̅ are the mean of observed values and predicted values, respectively.Generally, in the above measurement indices, lower values with the MAE and RMSE indicate better agreement, and higher positive values with the CC indicate better agreement between the observed values and predicted values.From Table2, it can be seen that the best forecast results for Peschanaya by TLLP method are obtained when the number of nearest neighbors is m+7, that is, 13, (with MAE = 0.4781, RMSE = 0.6878, CC = 0.9092). And the best forecast results for Tim by TLLP method are obtained when the number of nearest neighbors is m+16, that is, 22, (with MAE = 7.6517, RMSE = 9.4725, CC = 0.9089).The results of NLLP-1, NLLP-2, and the optimal results of the TLLP are showed in Figure3. The comparisons of the three methods between 20 lead time steps and 40 lead time steps are showen in Tables 3 and 4, respectively.Table 3
The comparison of the prediction results by TLLP (m=6, q=13), NLLP-1, and NLLP-2 for Peschanaya.
Method20 steps40 stepsMAE (m3/s)RMSE (m3/s)CCMAE (m3/s)RMSE (m3/s)CCTLLP0.47230.60050.87470.47810.68780.9092NLLP-10.23670.27670.96360.86131.13240.7158NLLP-20.35110.44550.69420.96581.18940.9516Table 4
The comparison of the prediction results by TLLP (m=6, q=22), NLLP-1, and NLLP-2 for Tim.
Method20 steps40 stepsMAE (m3/s)RMSE (m3/s)CCMAE (m3/s)RMSE (m3/s)CCTLLP2.57993.53690.95397.65179.47250.9089NLLP-17.44309.2437-0.53388.88729.8291-0.1647NLLP-22.60502.99280.91131.76002.24130.9054The prediction results of the hydrological time series. (a) The prediction results by TLLP (m=6, q=13), NLLP-1 and NLLP-2 for Peschanaya. (b) The prediction results by TLLP (m=6, q=22), NLLP-1 and NLLP-2 for Tim.
(a)(b)From Figure3(a) and Table 3, it can be seen that NLLP-1 gets the best prediction results at the beginning of the prediction process for Peschanaya (with MAE = 0.2367, RMSE = 0.2767 and CC = 0.9636 for NLLP-1 when the lead time step is 20), and under this condition, both of the results by NLLP-1 and NLLP-2 are better than those by TLLP. But as the prediction step increased, the forecast precision of NLLP-1 and NLLP-2 becomes worse than that of TLLP. When the lead time step becomes 40, the TLLP (with MAE = 0.4781, RMSE = 0.6878 and CC = 0.9092) obtains the best results.From Figure3(b) and Table 4, it can be seen that, not only at the beginning of the prediction procedure but also along the whole procedure when the lead time step is 40, the NLLP-2 gets the best prediction results (with MAE = 2.6050, RMSE = 2.9928, and CC = 0.9113 when the lead time step is 20 and MAE = 1.7600, RMSE = 2.2413, and CC = 0.9054 for TLLP when the lead time step is 40).From the above analysis, the prediction performance of the new method proposed in this study (NLLP-2) is better than that of TLLP at the beginning of prediction process for Peschanaya. For the daily discharge from Tim, the NLLP-2 is working superior than both TLLP and NLLP-1.The changes of the embedding dimension and the number of nearest neighbors in the foregoing three methods are shown in Figures4 and 5 (the lead time step is 20 days). It can be found that the embedding dimension and the number of nearest neighbors do change adaptively in the prediction processes.The changes of embedding dimension (a) and number of nearest neighbors (b) in the prediction process for Peschanaya (the lead time step is 20 days).
(a)(b)The changes of embedding dimension (a) and number of nearest neighbors (b) in the prediction process for Tim (the lead time step is 20 days).
(a)(b)


## 3.1. Study Area and Data Description

The real hydrological time series in this study are chosen as the daily discharges of Peschanaya at Tochil'noye (52.17°N and 85.17°E, basin area, 4720 km2) in Russian for the period January 1936–December 1980, and the daily discharges of Tim at Napas (59.85°N and 81.95°E, basin area, 24500 km2) in Russian for the period January 1953–December 1994. Some important statistical values are shown in Table 1.Table 1
Statistics indices of daily discharge time series at Peschanaya and Tim.
Statistical indicesPeschanayaTimTotal number of data1643715340Mean value (m3/s)32.3357192.5562Standard of deviation (m3/s)44.6106215.4076Maximum value (m3/s)661.00001360.0000Minimum value (m3/s)1.740028.4000Coefficient of variation1.37961.1187For the daily discharge from Peschanaya, the time series is divided into two data sets, the first training data set is the data during the period January 1936–December 1979, and the second prediction data set is the data during the period January 1, 1980–February 9, 1980. In the same way, the first training data set for Tim is the data during the period January 1953–December 1993, the second prediction data set is the data during the period January 1, 1994–February 9, 1994. The prediction lead time for both of the time series is 40 days.


## 3.2. Phase Space Reconstruction

In this study, the delay time and the embedding dimension of the above two daily discharge time series are determined by the autocorrelation function method and the Cao method, respectively. The changes of the autocorrelation functions for both of the time series are presented in Figure1 and the values of E1(m) are shown in Figure 2.Figure 1
Autocorrelation functions for daily discharge time series observed at Peschanaya and Tim. (The red line presents the situation that the autocorrelation function value is1-1/e of its initial value.)Figure 2
The values ofE1(m) for daily discharge time series observed at Peschanaya and Tim.The delay time is selected when the autocorrelation coefficient has dropped to1-1/e of its initial value. From Figure 1, the delay time for the Peschanaya and Tim can be determined as 9 and 20, respectively. (The red line in Figure 1 presents the situation that the autocorrelation function value is 1-1/e of its initial value.)From Figure2, the value of E1(m) becomes saturated from m0=5; so the embedding dimension can be determined as m0+1=6 for Peschanaya.  In the same way, the embedding dimension for Tim can be determined as 6 too.The third parameter in the local linear prediction of the daily discharge time series is the number of the nearest neighbors. In this study, only the traditional prediction method (TLLP) needs to give this parameter as a prior; the rest two methods (NLLP-1 and NLLP-2) can obtain the number of nearest neighbors adaptively in the prediction procedure.


## 3.3. The Prediction Results and Analysis

Both of the time series were predicted by TLLP, NLLP-1, and NLLP-2, respectively. In order to obtain the best forecast results with the TLLP, the number of nearest neighborsq for TLLP is changed from m+1 to m+19. By comparing the prediction results under different q, the optimal forecast results by TLLP in different q can be selected. The prediction results by TLLP were shown in Table 2.Table 2
The prediction results for the daily discharge time series observed at Peschanaya and Tim by TLLP method.
Number of nearest neighborsPeschanayaTimMAE (m3/s)RMSE (m3/s)CCMAE (m3/s)RMSE (m3/s)CCm+19.6882 × 1022.1893 × 1030.69661.6088 × 10104.5795 × 10100.2904m+24.2862 × 1051.4036 × 1060.554847.922478.86640.0600m+35.31516.20560.42938.1494 × 1041.4928 × 1050.5736m+41.47791.83100.95038.4746 × 10113.0884 × 10120.3439m+52.44154.28530.76962.0832 × 10179.1309 × 10170.2899m+60.84281.19070.83072.1107 × 10117.5627 × 10110.3491m+70.47810.68780.909222.463624.4395-0.3345m+80.53620.69170.889119.936520.9531-0.5286m+90.55990.76380.812719.111120.0202-0.5056m+100.49780.68150.923611.122314.21230.7822m+110.74450.96710.839918.263824.18650.8162m+122.62595.6787-0.656813.690517.09750.8609m+132.77426.0485-0.657411.015013.61990.8799m+141.95363.7148-0.684611.753514.73290.8702m+152.68775.8700-0.650511.726014.49790.8819m+161.83583.4853-0.66127.65179.47250.9089m+171.87633.3413-0.721611.948614.70640.8939m+182.51085.5012-0.639413.278116.13630.8916m+192.02664.0896-0.651614.820817.90760.8850The accuracy of prediction is evaluated by three measurement indices in this study, which are theMean absolute error (MAE), Root mean square error (RMSE), and Correlation coefficient (CC). The definitions of the three indices can be given as follows:(3.1)MAE=1n∑i=1n|yi-xi|,RMSE=1n∑i=1n(yi-xi)2,CC=∑i=1n(xi-x̅)(yi-y̅)∑i=1n(xi-x̅)2∑i=1n(yi-y̅)2,
where n is the number of the time series under investigation. The yi and xi are the predicted values and observed values in the time series, respectively, and x̅ and y̅ are the mean of observed values and predicted values, respectively.Generally, in the above measurement indices, lower values with the MAE and RMSE indicate better agreement, and higher positive values with the CC indicate better agreement between the observed values and predicted values.From Table2, it can be seen that the best forecast results for Peschanaya by TLLP method are obtained when the number of nearest neighbors is m+7, that is, 13, (with MAE = 0.4781, RMSE = 0.6878, CC = 0.9092). And the best forecast results for Tim by TLLP method are obtained when the number of nearest neighbors is m+16, that is, 22, (with MAE = 7.6517, RMSE = 9.4725, CC = 0.9089).The results of NLLP-1, NLLP-2, and the optimal results of the TLLP are showed in Figure3. The comparisons of the three methods between 20 lead time steps and 40 lead time steps are showen in Tables 3 and 4, respectively.Table 3
The comparison of the prediction results by TLLP (m=6, q=13), NLLP-1, and NLLP-2 for Peschanaya.
Method20 steps40 stepsMAE (m3/s)RMSE (m3/s)CCMAE (m3/s)RMSE (m3/s)CCTLLP0.47230.60050.87470.47810.68780.9092NLLP-10.23670.27670.96360.86131.13240.7158NLLP-20.35110.44550.69420.96581.18940.9516Table 4
The comparison of the prediction results by TLLP (m=6, q=22), NLLP-1, and NLLP-2 for Tim.
Method20 steps40 stepsMAE (m3/s)RMSE (m3/s)CCMAE (m3/s)RMSE (m3/s)CCTLLP2.57993.53690.95397.65179.47250.9089NLLP-17.44309.2437-0.53388.88729.8291-0.1647NLLP-22.60502.99280.91131.76002.24130.9054The prediction results of the hydrological time series. (a) The prediction results by TLLP (m=6, q=13), NLLP-1 and NLLP-2 for Peschanaya. (b) The prediction results by TLLP (m=6, q=22), NLLP-1 and NLLP-2 for Tim.
(a)(b)From Figure3(a) and Table 3, it can be seen that NLLP-1 gets the best prediction results at the beginning of the prediction process for Peschanaya (with MAE = 0.2367, RMSE = 0.2767 and CC = 0.9636 for NLLP-1 when the lead time step is 20), and under this condition, both of the results by NLLP-1 and NLLP-2 are better than those by TLLP. But as the prediction step increased, the forecast precision of NLLP-1 and NLLP-2 becomes worse than that of TLLP. When the lead time step becomes 40, the TLLP (with MAE = 0.4781, RMSE = 0.6878 and CC = 0.9092) obtains the best results.From Figure3(b) and Table 4, it can be seen that, not only at the beginning of the prediction procedure but also along the whole procedure when the lead time step is 40, the NLLP-2 gets the best prediction results (with MAE = 2.6050, RMSE = 2.9928, and CC = 0.9113 when the lead time step is 20 and MAE = 1.7600, RMSE = 2.2413, and CC = 0.9054 for TLLP when the lead time step is 40).From the above analysis, the prediction performance of the new method proposed in this study (NLLP-2) is better than that of TLLP at the beginning of prediction process for Peschanaya. For the daily discharge from Tim, the NLLP-2 is working superior than both TLLP and NLLP-1.The changes of the embedding dimension and the number of nearest neighbors in the foregoing three methods are shown in Figures4 and 5 (the lead time step is 20 days). It can be found that the embedding dimension and the number of nearest neighbors do change adaptively in the prediction processes.The changes of embedding dimension (a) and number of nearest neighbors (b) in the prediction process for Peschanaya (the lead time step is 20 days).
(a)(b)The changes of embedding dimension (a) and number of nearest neighbors (b) in the prediction process for Tim (the lead time step is 20 days).
(a)(b)


## 4. Conclusion and Prospect

To obtain the optimum prediction results, a new adaptive local linear prediction method is proposed in which the embedding dimension and the number of nearest neighbors are combined as a parameter set(m,q) and change adaptively as the forecast steps increased in this study. The main results are given as follows.(i)
The optimal parameters(m,q) can be gotten by using the GDF method. In the process of the selection of the optimal parameters (m,q) for each prediction step, the generalized degree of freedom (GDF) is used and different error variances are calculated under different combinations of (m,q). The optimal parameters (m,q) are chosen when the variance is the smallest. Then the optimal parameters are used to reconstruct the phase space and predict the next value of the time series.(ii)
For the sake of comparing the performance of the different prediction method, real nonlinear hydrological time series are chosen to be examined. The results from the used example indicate that the new adaptive local linear prediction method proposed in this study can choose the embedding dimension and the number of the nearest neighbors adaptively in the prediction process. Compared with the methods of TLLP and NLLP-1, the new method is better than TLLP at the beginning of the prediction steps for both of the time series, and NLLP-2 is better than TLLP and NLLP-1 for the time series from Tim during the whole forecast period.(iii)
The new adaptive local linear prediction method can be used in predicting other nonlinear time series in the future and its theory will be further studied.



---

*Source: 205438-2010-06-28.xml*