Local Volterra multivariable chaotic time series multi-step prediction based on phase points clustering

1. Introduction

Chaotic time series prediction is an important application field and research hotspot of chaos theory [1]. In recent years, with the development of nonlinear science, multivariable chaotic time series prediction shows a greater advantage than univariate chaotic time series prediction, Multivariable chaotic time series prediction has wide application prospect in many aspects such as weather forecasting, runoff forecasting, economic forecasting, wind and electricity load forecast [2], stock forecasting [3] and target trajectory prediction. The research object of this paper is air pollutant PM 2.5. In general, air pollutant emissions are important indicators for evaluating environmental management and implementing new strategies. Fine particles (PM2.5 or PM10) contain particles no larger than 2.5 μm or 10 μm [4]. These particles are an important part of air pollution and are mainly caused by the combustion of fossil fuels and air pollutants such as NO₂, CO, and SO₂. High concentrations of PM2.5 can cause serious air pollution, causing a series of respiratory and pulmonary diseases, and have a serious impact on public health [5]. PM2.5 is influenced by multivariable factors such as NO₂, CO, humidity, wind speed and wind direction, and has non-linear relationship with air quality and meteorological data. In order to obtain higher prediction accuracy, researchers in recent years have made a lot of researches on the methods of predicting the concentration of PM2.5, PM10 and other air pollutants [6]. For the linear prediction method, Jian et al. used the ARIMA model to study the influence of meteorological factors on ultrafine particulate matter (UFP) and PM10 concentrations under heavy traffic conditions in Hangzhou [7]. Vlachogianni et al. using multivariate linear regression model to predict the NOx and PM10 concentrations in Athens, the NO, NO₂, CO and PM 2.5 concentrations were used as predictive variables [8]. The prediction models such as Arima and MLR have been widely used to predict PM concentration, but the linear method can provide reliable prediction result only when the time series data is linear or near linear, however, this kind of method still has wide application scope in the Prediction field [9]. For nonlinear predictive methods, for example, Cobourn, a nonlinear regression model is proposed to predict PM 2.5 air quality in Kentucky State Louisville metropolitan Area [10]. Artificial neural Network (ANN) model has the ability of detecting response and complex non-linear relationship between predictive factors, and can be trained by many effective training algorithms, therefore, artificial neural network model has aroused more and more attention in the field of air quality prediction [11]. For example, Pérez et al. uses the artificial neural network model to predict PM2.5 concentration and humidity by using wind speed and wind direction and PM2.5 concentration as predictive variables [12]. Fernando et al. show that the artificial neural network shows the advantages of fast calculation and high precision, can be used for the rapid development of air quality early-warning system, ANN prediction error changes according to the predicted range and different input variables. Although artificial neural network can produce better prediction results, the disadvantage of this method is that it requires a large number of training samples and results in unstable training results [13].

In recent years, with the progress and development of functional theory, Volterra filter has gained wide attention with its advantages of fast training speed, strong nonlinear approximation ability and high prediction accuracy [14]. At present, domestic and foreign scholars have conducted a series of researches on the multi-step prediction of high dimensional chaotic time series by Volterra functional technology, and its application is continuously deepening and refining. Literature [15] proposed the wind power prediction genetic algorithm-Volterra neural network (GA-VNN) model. This method combines the structural characteristics of the Volterra functional model and the BP neural network model, and uses the genetic algorithm to improve the combination model. In terms of global optimization capability, the predictive performance of wind power super short-term multi-step prediction is significantly higher than that of Volterra model and BP neural network model. However, this combination model has short prediction steps, large application limitations, and complex structures. Larger, less predictive efficiency. The above method only improves multi-step prediction performance from a single point of view, and does not fully consider the comprehensive factors such as prediction accuracy, prediction efficiency, multivariate information, and applicability. The literature [16] proposed a Volterra-based multivariate chaotic time series multi-step prediction model. This method expands the univariate Volterra filter into a multivariate, taking into account the effect of multivariate information on the performance of multi-step predictions. When the number is small, its multi-step prediction performance is better than that of the univariate Volterra model and the multivariable local polynomial model. However, because the model is based on global prediction and multi-step prediction, when the data volume is large, its multi-step prediction efficiency needs further improvement. It is one of the most efficient and common methods to predict high-dimensional chaotic complex systems through local modeling. The key to improving local multi-step prediction accuracy is to select a reasonable neighboring phase point. The current major criteria are: Distance, direction and distance comprehensive criteria, vector similarity and so on [17]. According to the theory of complex systems, the correlation information between multi-dimensional variables and observational data are all involved in the evolution of the system. The above criteria do not fully consider the association information between variables and multi-steps when selecting the optimal neighboring phase points. Comprehensive factors such as the evolutionary rule after the retrogression and the acting forces at different locations. At the same time, for the selection of the optimal neighboring phase points, the traditional local area prediction method compares the predicted phase points with all observed phase points in the phase space, and the amount of calculation is large. In particular, selecting the optimal neighboring phase point in observation points with a large amount of data has a long operation time, which greatly reduces the multi-step prediction efficiency. Currently, there is no effective solution to the related literature.

PM2.5 concentration chaotic time series data is affected by multiple variables such as temperature and wind speed in the prediction. To improve the prediction accuracy, a multi-step local model Volterra multivariable chaotic time series prediction model based on phase clustering is proposed. The experimental data was collected from UCI Machine Learning Repository’s US Embassy in Beijing within five years of weather reporting and pollution levels. Modeling the PM2.5 concentration, based on the multi-variable local Volterra multi-step prediction model, combined with the evolution law of multi-step backtracking of adjacent phase points and the influence of the correlation information between variables on the evolution trajectory of phase points, A new multivariate evolution trajectory similarity synthetical criterion is used to improve the precision of multi-step prediction of PM2.5 concentration by selecting a reasonable neighboring phase point. According to the characteristics of the similarity between observation points in phase space, the observation phase points are first clustered, and then the similarity degree is compared between the phase points and the nearest cluster phase point to select the optimal neighboring phase point. Reduce the number of calculation iterations and reduce multi-step prediction time. The idea of phase cluster analysis is used to improve the multi-step prediction efficiency of PM2.5 concentration, and ultimately achieve the purpose of comprehensively improving multi-step prediction performance of multi-variable local Volterra model. It has broad application prospect and value in air quality prediction.

2. Local Volterra multivariable chaotic time series prediction

3. Select the optimal adjacent phase point

The selection of traditional adjacent phase points is mostly based on the similarity of the one-step evolution of the chaotic orbit, which lacks the consideration of the overall evolution of the system, and ignores the influence of the correlation information between variables on the evolutionary trajectory. When the prediction step is longer, the prediction trajectory and the original trajectory deviate greatly. In order to improve the learning performance of the model, avoid increasing the complexity of the model, while avoiding the influence of the relative weak phase point on the model precision, the selection of the adjacent phase point should not only select the phase points similar to the prediction point evolution trajectory in many phase points, but also control the number of adjacent points. Therefore, by defining the mutual information relative contribution rate, the multi-step evolution distance similarity and multi-step evolution direction to measure the similarity of the multivariate evolutionary trajectory, this paper proposes a selecting optimal adjacent phase points synthesis criterion.

Mutual information can better reflect the statistical dependence of variables, and mutual information can be used for multivariate correlation analysis [20]. In order to avoid the direct estimate of high dimensional probability density function, based on the estimation of $k$ adjacent mutual information [21], the relative contribution rate of the input variable to the nearest neighbor information value is used to represent the correlation information between variables. The basic method of $k$ adjacent mutual information estimation is: if $\epsilon \left(i\right)/2$ represents point $e_i=(x_i,y_i)$ to its $k$th neighbor distance, the distance between the $X$-axis and the $Y$-axis is respectively, and the mutual information estimation of variables $X$ and $Y$ is:

where $\psi \left(x\right)$is the digamma function, which is generally written as $\psi \left(x\right)={\mathrm{\Gamma }}^{\mathrm{\text{'}}}\left(x\right)/\mathrm{\Gamma }\left(x\right)$; $n_x$ and $n_y$ are the number of data points of the conditions of $‖x_i-x_j‖\le {\epsilon }_x\left(i\right)/2$and $‖y_i-y_j‖\le {\epsilon }_y\left(i\right)/2\mathrm{}\mathrm{}\mathrm{}$; $⟨\cdot \cdot \cdot ⟩$ represents the average of all of these variables.

For multivariate, the mutual information can be extended by Eq. (10) to:

where $M$ is the dimension of the variable. Based on the relative contribution of different input variables to the overall mutual information value, the correlation between multivariate variables is analyzed. Let $S=\left\{X_1,X_2,\dots ,X_M\right\}$ be the subset of input variables, $S_i$ means to remove a subset of $X_i$ variables in the $S$ set ($i=\mathrm{1,2},\dots ,M$), and $Y$ is the output variable.

The relative contribution rate of the variable $X$ to the mutual information is defined as:

Criterion 1: The distance similarity between the predicted point $X\left(p\right)$ and the phase point $X\left(t\right)$ is as follows:

where $d_p=‖X\left(p\right)-X(p-q)‖$, $\mathrm{}{d}_t=‖X\left(t\right)-X(t-q)‖$, $\delta (p,t,q)\in \left[\mathrm{0,1}\right]$, and $\delta (p,t,q)$ is smaller, the closer the evolution distance of the $q$ step of the adjacent phase point $X\left(t\right)$ is to the $q$ step evolution of $X\left(p\right)$. Considering the smaller phases of the backtracking step, the greater the role of the prediction point, the weighted processing of the action of the multi-step backtracking point. Define the back step distance weight vector $\gamma \left(j\right)$, $(j=\mathrm{1,2},\dots ,q)$, and satisfy the condition: $\gamma \left(1\right)\ge \gamma \left(2\right)\ge \cdot \cdot \cdot \ge \gamma \left(q\right)$, $\gamma \left(j\right)=2(q-j+1)/q(q+1)$.

The weighted treatment of the similarity of multi-step evolution distance is:

2) Criterion 2: The direction distance similarity between the predicted point $X\left(p\right)$ and the phase point $X\left(t\right)$ is as follows:

where $F\left(p,q\right)=X\left(p\right)-X\left(p-q\right)$, $F\left(t,q\right)=X\left(t\right)-X\left(t-q\right)$, $\left[\right]$ represents the inner product of vectors, $‖‖$ represents the distance norm of the vector, $\left|\right|$ represents the absolute value, and the smaller the value of $\mathrm{c}\mathrm{o}\mathrm{s}\theta \left(p,t,q\right)$, the closer the $q$ step evolution direction of the adjacent phase point $X\left(t\right)$ is to the $q$ step of $X\left(p\right)$. Define backward step weight vector $\phi \left(j\right)\left(j=\mathrm{1,2},\dots ,q\right)$, meet the condition: $\phi \left(1\right)\ge \phi \left(2\right)\ge \cdot \cdot \cdot \ge \phi \left(q\right)$, take $\phi \left(j\right)=2\left(q-j+2\right)/q\left(q+3\right)$. The weighted processing of the similarity of multi-step evolution direction is:

3) Comprehensive criterion multivariate evolutionary trajectory similarity.

The evolutionary trajectory similarity between the prediction point $X_i\left(p\right)$ and the phase point $X_i\left(t\right)$ of the single variable $X_i$ is:

where ${\mu }_1$, ${\mu }_2$ are the weight values of the distance index and the direction index, respectively, and ${\mu }_1+{\mu }_2=$ 1, $i=$ 1, 2, …, $M$ represents the number of variables. The value of ${\lambda }_i\left(p,t\right)$ is smaller, indicating that the prediction point $X_i\left(p\right)$ of a single variable $X_i$ is similar to that of phase point $X_i\left(t\right)$. Based on the relative contribution of input variables to the mutual information value, the similarity of the single variable evolution trajectory are extended to the multivariable, and the similarity of the evolution trajectory of the predicted point $S_i\left(p\right)$ and the phase point $S_i\left(t\right)$ of the input variable set $S=\left\{X_1,X_2,\dots ,X_M\right\}$ is as follows:

Finally, according to the multivariate evolution similarity index selected the optimal point data, carrying on the training model.

4. Adjacent points clustering analysis

For the selection of the optimal adjacent point, the traditional local prediction method is to compare the predicted phase point with the global observation point, the calculation time is longer and the prediction efficiency is lower. For this reason, based on the $K$-means clustering algorithm, this paper proposes a new scheme based on adjacent phase point clustering analysis. Firstly, the observation phase points are analyzed by clustering, and then the predicted phase points are compared with the nearest cluster phase points, and the optimal adjacent phase points are selected to reduce the prediction time.

Due to the complexity of multi-variable observation point, the cluster with different sizes and densities may be formed. The $K$-means algorithm is not very good for the data clustering of heterogeneous distribution. In order to prevent the density of the cluster from being too high and influence the prediction accuracy, the weighted processing of the sum of the errors in the cluster is carried out, and the criterion function of the clustering is defined as:

where $K$ is the number of cluster centers, $N$ represents the total number of elements of the data set, $W_j$ stands for category $j$, $N_j$ represents the number of elements in the $j$th class, ${\stackrel{-}{X}}_{w_j}$ represents the clustering center vector of class $j$, and $X$ represents the element vector belonging to class $j$.

The clustering analysis of adjacent phase points takes the predicted phase points as the reference standard of the evolutionary trajectory similarity, and then clustering the evolutionary trajectory similarity of all the observation points. By the Eq. (17), the similarity vector of the multivariate evolution trajectory is defined as:

By Eq. (18), the similarity value of multi-variable evolution trajectory of $l$ observation phase points is calculated, which is ranked as $\left[{\xi }_1,{\xi }_2,\cdot \cdot \cdot ,{\xi }_l\right]\left({\xi }_1\le {\xi }_2\le \cdot \cdot \cdot \le {\xi }_l\right)$ in ascending order and $\left[l/K\right]$ (integer) as the interval. In turn, the evolution path similarity vector $\stackrel{-}{\lambda }\left(1\right),\stackrel{-}{\lambda }\left(2\right),\dots ,\stackrel{-}{\lambda }\left(K\right)$ of $K$ corresponding phase point is selected as the initial cluster center, and the initial cluster center can be dispersed as far as possible to close to the final clustering center and reduce the time of clustering.

The update formula of clustering center is:

Based on the relative contribution rate of input variables to the mutual information value, by Eq. (12), the vector of the multi-variable mutual information contribution rate is:

Define a new operator symbol:

According to the multidimensional evolutionary trajectory similarity clustering feature of the phase point, from Eqs. (20-23), (19) can be expressed as:

5. The implementation step of multivariable multi-step prediction

For the multi-step prediction problem of high-dimensional chaotic system, a lot of practice shows that most nonlinear dynamical systems can be characterized by Volterra series [22-24], so we can construct multi-step chaotic time series multi-step prediction model with Volterra series expansion. Therefore, in order to make the multivariate adjacent phase points selection more reasonable and more efficient, this paper proposes a local Volterra multivariable chaotic time series multi-step prediction based on phase points clustering. The concrete steps are as follows:

Selection of input variables and multivariate phase space reconstruction;

Taking the predicted phase points as the reference standard of evolution trajectory similarity, the multi variable trajectory similarity vector $\lambda$ of all observation points is calculated, which is used as the data set to be clustered, and determined the number of clustering centers $k$;

The prediction phase point is used as the reference standard for the similarity of evolutionary trajectory. In order to form $K$ clusters, we calculate the similarity value $\xi$ of multivariate evolutionary trajectory of all the observed phase points, and take the value of ascending arrangement and equal interval, then select $K$ corresponding phase point's evolution trajectory similarity vector $\stackrel{-}{\lambda }\left(1\right),\stackrel{-}{\lambda }\left(2\right),\dots ,\stackrel{-}{\lambda }\left(K\right)$ as the initial cluster center.

By redistributing the phase point, the value of the criterion function is minimized, and $k$ new clusters are formed and the cluster centers are updated.

Cycle step 4, until the cluster centers are not changed, clustering is completed;

The evolution trajectory similarity vector $\lambda$ of the predicted phase point is calculated, and calculate the distance to each cluster center $\stackrel{-}{\lambda }\left(j\right)$, choose the nearest cluster, and the prediction phase point is compared with the similarity of the evolution trajectory of all the observation points in the cluster, and the d optimal adjacent phase points are selected from this.

The optimal adjacent phase points are used to train the multivariable Volterra model. and the multivariable coefficient vector $H\left(n\right)$ is obtained, and the input vector $U\left(n\right)$ is generated by the prediction phase point, and the one-step prediction value $\widehat{x}\left(n+1\right)$ is given.

The predictive value is entered into the observational data, then iterated, and the function of multi-step prediction is finally realized.

The multi-step prediction performance of this model depends on the rationality and efficiency of the selection of the optimal adjacent phase points.

The realization step of multi-step prediction. According to the local Volterra prediction model proposed in this paper, we can get the first step prediction value. In the second step prediction, the first step prediction value is added as the known value to the historical data to calculate the second step Prediction value, when calculating the prediction value of the third step, the prediction value of the second step of the first step is added as the known value to the historical data to calculate the prediction value of the third step. According to this method, according to the model Continuous iteration for multi-step prediction.

6. Simulation and result analysis

6.2. The judgment of chaotic characteristics

In this paper, the chaotic characteristics of PM2.5 concentration, temperature and wind speed are identified by the saturation correlation dimension method. Crassberger and Procaccia [25, 26] proposed the correlation dimension method to investigate the chaotic characteristics of time series. For random sequences, the correlation dimension increases with the increase of embedding dimension, and does not reach saturation. For the chaotic sequence, correlation dimension tends to saturation gradually along with the increase of the embedding dimension. Therefore, we can distinguish the chaotic sequence from the random sequence according to whether the correlation dimension is saturated or not.

Consider the vector set near the attractor $\left\{Y_i=\mathrm{1,2},\dots ,N\right\}$, the associated integral $C\left(r\right)$ is defined as the number of points that are less than r for any given r, and the ratio of the number of points $\left(Y_i,Y_j\right)$ to all possible points, that is:

25

$C\left(r\right)=\frac{2}{N(N-1)}{\sum }_{i=1}^N{\sum }_{j=i+1}^N\mathrm{\Theta }\left(r-\left|Y_i-\right|Y_j\right),$

where:

$\mathrm{\Theta }\left(x\right)=\left\{\begin{array}{l}0,\mathrm{}\mathrm{}\mathrm{}\mathrm{}x\le 0,\\ 1,\mathrm{}\mathrm{}\mathrm{}\mathrm{}x>0,\end{array}\right.$

when $N\to \mathrm{\infty }$, for any small $r$, we can follow the law of the exponential power change of $C$, which is $C\left(r\right)\subset r^D$, and we can calculate $D=d\left(\mathrm{l}\mathrm{n}\right[C\left(r\right)\left]\right)/d\left(\mathrm{l}\mathrm{n}\right(r\left)\right)$ by the derivative of $\mathrm{l}\mathrm{n}\left(C\right(r\left)\right)$ with $\mathrm{l}\mathrm{n}\left(r\right)$. In this paper, the G-P algorithm [27] is used to calculate the saturation correlation dimension of the time series. Figs. 1-3 respectively show the $\mathrm{l}\mathrm{n}\left(C\right(r\left)\right)$ and $\mathrm{l}\mathrm{n}\left(r\right)$ curves of PM2.5, temperature and wind speed time series, where $C\left(r\right)$ is the correlation function and $r$ is the radius.

Fig. 1PM2.5 time series lnC(r)-lnr graph

Fig. 2temperature time series lnC(r)-lnr graph

Fig. 3wind time series lnC(r)-lnr graph

As shown in Figs. 1-3, the $\mathrm{l}\mathrm{n}\left(C\right(r\left)\right)$ and $\mathrm{l}\mathrm{n}\left(r\right)$ curves gradually become parallel with the embedding dimension increasing, that is, the correlation dimension gradually reaches saturation. According to Eq. (8), the attractor dimension of PM2.5 in temperature and wind speed time series is respectively $D\left(p\right)=$ 2.3576, $D\left(t\right)=$ 0.7144 and $D\left(w\right)=$ 1.2835. The attractor dimension is fractional Form, indicating that there are chaos characteristics in time series of PM2.5, temperature and wind speed.

6.5. Select the number of cluster centers

The 5000 groups of training data produced by the Lorenz model for multivariate phase space reconstruction, and a total of 4976 phase points are generated, the 4976th phase point as the prediction phase point, the first 4975 phase points as the observation phase point, and the prediction phase point as the reference standard of the multivariate evolution trajectory similarity, The multivariate evolutionary trajectory similarity vectors of all the observed phase points are calculated as data sets for clustering. The evolution trajectory similarity vector of the first 4975 observation points is shown in Table 3.

In order to study the influence of the number of clustering centers on the prediction time and precision, the root mean square error of the ten-step prediction is used as the evaluation criterion, and the time of consumption is calculated by the 50-step prediction range. The prediction model uses the model proposed in this paper. the experimental results are shown in Table 4.

Table 3The evolution trajectory similarity vector of the observation point

X component	Y component	Z component
0.0574	0.0613	0.2837
0.0385	0.0572	0.2635
0.0176	0.0245	0.1439
⋅⋅⋅	⋅⋅⋅	⋅⋅⋅
0.0016	0.0047	0.0075

Table 4The influence of the number of clustering centers on predicting time and accuracy

Number of cluster centers	Prediction error	Total time-consuming/s
2	4.123×10-1	7.153
4	1.734×10-1	5.306
6	5.167×10-2	2.527
8	8.109×10-3	1.128
10	3.285×10-2	2.637
12	6.497×10-2	4.417
14	8.536×10-2	6.973

As can be seen from Table 4, with the increasing number of clustering centers, the root mean square error of the ten-step prediction decreases first and then increases, and the total time-consuming of 50 steps is reduced and finally stabilizes at about 1.128 seconds. When the number of cluster centers is 8, the multi-variable evolution path similarity vector of 4976 observation phase points is clustered, and the projection of the X and Y direction is shown in Fig. 4. The clustering center and phase point distribution are shown in Table 5.

Due to the complexity and uncertainty of the multivariate evolutionary trajectory similarity vector of the observation phase points, the uneven distribution of the data sets may lead to the formation of clusters with different sizes and densities. It is shown from Fig. 4 and Table 5 that the $K$-means algorithm clustering with weighted processing avoids the high density of elements in the cluster, makes the clustering result relatively balanced, and greatly reduces the effect on the multi-step prediction precision of the model.

Fig. 4The XY direction projection of the adjacent phase clustering

Table 5Distribution of cluster centers and phase points

Clustering center	X component	Y component	Z component	Number of phase points
1	0.1909	0.1626	0.3952	601
2	0.8377	0.5395	0.0271	576
3	0.5249	0.3048	0.2138	1398
4	0.0976	0.4185	0.5216	421
5	0.0976	0.1096	0.1646	763
6	0.3931	0.1464	0.2385	453
7	0.0262	0.0310	0.0602	354
8	0.2909	0.2666	0.6347	409

6.6. Multi-step prediction performance comparison of each model

6.6.1. The local multivariable linear autoregressive model

The AR model is a commonly used time series prediction model, which has a relationship between variables of time series, and the current variable can be represented by the previous time [28, 29]. The local multivariable linear autoregressive model function can be written as:

27

$x_{1,n+1}=a+{\sum }_{i=1`}^M{\sum }_{j=1}^{m_i}{b}_{i,j}{x}_{i,n-\left(j-1\right){\mathrm{\tau }}_i}=a+b_{\mathrm{1,1}}{x}_{1,n}$
$+\cdots +b_{1,m_1}{x}_{1,n-\left(m_i-1\right){\tau }_1}+\cdots +b_{M,m_M}{x}_{M,n-\left(m_M-1\right){\tau }_M},$

where $a$, $b_1$, 1, …, $b_{M,m}$ is the fitting parameter. Combined with the evolutionary trajectory comprehensive criterion, the parameters are estimated by the least square method to calculate the predicted value.

6.6.2. The local multivariable least learning machine (ELM) prediction model

For a multivariable time series training set:

$S=\left\{\left.({\mathrm{x}}_{\mathrm{i}},y_{\mathrm{i}})|x_i\in X\in R^n,\mathrm{}\mathrm{}\mathrm{}{y}_i\in Y\in R,\mathrm{}\mathrm{}\mathrm{}i=1,\dots ,N\right\}\right.,$

with $N$ samples, where $x_i$ is the input variable and$y_i$is the output variable, the objective function of ELM [30] is: $min\frac{1}{2}{‖\beta ‖}^2+\frac{C}{2}{\sum }_{i=1}^N{e_i}^2$, $y_i=h\left(x_i\right)\beta +e_i$, $(i=1,\dots ,N)$. Where $x$ is the output vector of the hidden layer; $L$ is the number of hidden nodes, $W_j$ is the input weight of the $j$th hidden node; $b_j$ is the deviation of the $j$th hidden node; $\beta$ is the output weight of the hidden layer; And for the $j$th hidden node excitation function. According to the Lagrange multiplier method and KKT optimal conditions, we have: $\beta =H^T\alpha =H^T[H{H}^T+diag\left(\frac{1}{C},\dots ,\frac{1}{C}\right){]}^{-1}Y$:

28

$H=\left[\begin{array}{c}h\left(x_1\right)\\ ⋮\\ h\left(x_N\right)\end{array}\right]={\left[\begin{array}{ccc}g(w_1,b_1,x_1)& \cdots & g(w_L,b_L,x_1)\\ ⋮& & ⋮\\ g(w_1,b_1,x_N)& \cdots & g(w_L,b_L,x_N)\end{array}\right]}_{N\times L}.$

So as to obtain the prediction model $f\left(x\right)=h\left(x\right)\beta$.

6.6.3. The local multivariable RBF neural network multi-step prediction model

In the RBF neural network structure as shown in Fig. 5, the input vector of the network is the transpose of the input of m-dimensional vector $X_n$, then the number of network input layer nodes is $m={\sum }_l^M{m}_l$. Let the nodes of the hidden layer of the network have radial basis vectors $H={\left[h_1,h_2,\dots ,h_j,\dots ,h_k\right]}^T$ and $h_j$ as Gauss basis functions [31]:

29

$h_j=\mathrm{e}\mathrm{x}\mathrm{p}\left(-\frac{‖X-C_j‖}{2b_j}\right),\mathrm{}\mathrm{}\mathrm{}j=\mathrm{1,2},\dots ,k.$

In this hidden layer, the center vector of Gaussian basis function of the $j$th node is: $C_j={\left[c_{j1},c_{j2},\dots ,c_{ji},\dots ,c_{jm}\right]}^T$, $i=$ 1, 2, …,$m$. The same dimensionality of $C_j$ and $X_n$ can determine the quality of the network to a large extent. Let the base width vector of hidden layer node of network be: $B={\left[b_1,b_2,\dots ,b_j,\dots ,b_k\right]}^T$. Where $b_j$ is the Gaussian basis function radius of node $j$, the size of which determines the complexity of the network, and the appropriate value can be selected through experiments and error information; $‖•‖$ is the vector norm, which is generally the Euclidean norm.

The weight vector from the hidden layer to the output layer is:

30

$W={\left[w_1,w_2,\dots ,w_j,\dots ,w_k\right]}^T.$

The output function of the network is:

31

$x_{l,n+1}={\widehat{G}}_l\left(X_n\right)={\sum }_{j=1}^k{w}_j{h}_j={\sum }_{j=1}^k{w}_j\exp \left(-\frac{{‖X-C_j‖}^2}{2b_j}\right).$

In order to verify the superiority of the proposed model algorithm, the experimental performance is compared with the multi-variable local linear AR model, multivariable local RBF neural network multi-step prediction and multivariable ELM multi-step prediction model. The experimental results are shown in Figs. 6-9.

Fig. 5RBF network structure

Fig. 6The multi-step prediction model proposed in this paper

a)

b)

Fig. 7multi-step prediction of local linear AR model

a)

b)

Fig. 8Multi-step prediction of local RBF model model

a)

b)

Fig. 9Multi-step prediction of local ELM model

a)

b)

In order to verify the multi-step prediction performance of the proposed model, compared with other multi-step prediction models, the effective step size of multi-step prediction is 50 steps, as shown in Figs. 6- 9. It can be seen from Fig. 6 that the prediction accuracy of the model proposed in the previous 12 steps is very high. Although the predicted trajectory deviates from the original trajectory after 12 steps, the deviation is relatively small, indicating that this model predicts a longer prediction step Long and small prediction error. As can be seen from Fig. 7, the multivariate local linear AR model has only 4 steps for multi-step prediction. After 4 steps, the predicted trajectory deviates from the original trajectory, the prediction error fluctuates greatly, and the prediction effect is not good. As can be seen from Fig. 8 and Fig. 9, the multivariable local RBF neural network model and Multivariable Local Least Learning Machine (ELM) model have five-step and seven-step effective prediction steps. Although the prediction error is small in the effective step, the same deviation from the predicted trajectory also limits the effective step of the prediction and the prediction error is larger.

In order to verify the high efficiency of the model, the root mean square error is used as the evaluation criterion, as shown in Eq. (26), to estimate the time consumption of the 50 step prediction step. The prediction error and total time-consuming of each model are shown in Table 6. It can be seen from Table 6 that the total time-consuming forecast for multi-step prediction is only 0.835 s, which greatly improves the prediction efficiency. In summary, the local Volterra multivariable chaotic time series multi-step prediction based on phase point clustering proposed in this paper has better multi-step prediction performance.

Table 6Predictive errors and total time-consuming comparisons for each model

Multi-step prediction model	One step prediction error	Ten steps prediction error	Total time-consuming / s
Multivariate local linear AR model	6.317×10-2	8.165×10-1	3.725
Multivariable local RBF model	2.054×10-2	5.483×10-1	9.043
Multivariate local ELM model	8.523×10-2	2.692×10-1	5.416
The model presented in this paper	4.564×10-3	8.627×10-2	0.835

6.6.4. Cram’er-Rao lower bound performance evaluation

Cram’er-Rao Lower Bound (CRLB) determines a lower bound for the variance of any unbiased estimator, and the estimate for reaching this lower bound is called the minimum variance unbiased (MVU) estimator [32]. Let $\theta$ be the parameter to be estimated, estimate it from the observed value $x$, they satisfy $a$ probability density function $f(x;\theta )$, the lower bound of the variance of any unbiased estimate is the reciprocal of Fisher’s information:

32

$\mathrm{v}\mathrm{a}\mathrm{r}\left(\widehat{\theta }\right)\ge \frac{1}{I\left(\theta \right)}.$

The Fisher information is defined as:

33

$I\left(\theta \right)=E[\left(\frac{\partial \mathrm{l}\mathrm{n}f\left(x;\theta \right)}{\partial \theta }{)}^2\right]=-E\left[\left(\frac{{\partial }^2\mathrm{l}\mathrm{n}f\left(x;\theta \right)}{\partial {\theta }^2}\right)\right].$

If and only if:

$\frac{\partial \mathrm{l}\mathrm{n}f\left(\mathbf{x};\theta \right)}{\partial \theta }=I\left(\theta \right)\left(g\left(\mathbf{x}\right)-\theta \right).$

The unbiased estimate for all $\theta$ to reach the lower bound can be calculated as $\widehat{\theta }=g\left(x\right)$, which is the MVU estimator and the minimum variance is $1/I\left(\theta \right)$.

We use the parameters of the multivariate Volterra model as an estimator, and use CRLB to evaluate the performance of the prediction result. The time series is $x\left(1\right),x\left(2\right),...,x\left(N\right)$. For multivariate estimated parameter $\theta =[{\theta }_1,{\theta }_2,\dots ,{\theta }_N{]}^T$, its probability density function is $f(x;\theta )$, which satisfies two regularization conditions [33]. The Fisher information matrix is an $N\times N$ matrix.

$I(\theta {)}_{m,k}$ is defined as:

$I(\theta {)}_{m\times k}=E\left[\left(\frac{\partial \mathrm{l}\mathrm{n}f\left(x;\theta \right)}{\partial {\theta }_m}\right)\left(\frac{\partial \mathrm{l}\mathrm{n}f\left(x;\theta \right)}{\partial {\theta }_k}\right)\right]=-E\left[\left(\frac{{\partial }^2\mathrm{l}\mathrm{n}f\left(x;\theta \right)}{\partial {\mathrm{\theta }}_m\partial {\theta }_k}\right)\right].$

Let$T\left(x\right)$ be an estimate of a parameter vector, $T\left(x\right)=\left(T_1\right(x),T_2(x),\dots ,T_k(x){)}^T$, and remember that its expected vector $E\left(T\right(x\left)\right)$ is $\varphi \left(\theta \right)$. The Cramer-Rao lower bound holds that the covariance matrix of $T\left(x\right)$ satisfies:

$Co{v}_{\theta }\left(T\right(x\left)\right)\ge \frac{\partial \varphi \left(\theta \right)}{\partial \theta }[I\left(\theta {)}^{-1}\right]{\left(\frac{\partial \varphi \left(\theta \right)}{\partial \theta }\right)}^T,$

where $A\ge B$ means that $A-B$ is a positive semidefinite matrix; $\partial \varphi \left(\theta \right)/\partial \theta$ is a Jacobian matrix whose $ij$th element is $\partial \varphi (\theta {)}_i/\partial {\theta }_j$.

Fig. 10Plot of CRLB performance evaluation

In order to better illustrate the effectiveness of the proposed algorithm, CRLB is used to estimate the lower bound of the MSE of the predicted value in a one-step prediction. Fig. 10. shows the MSE of the lower limit assessment of CRLB. The method presented in this paper also clearly shows that MSE has the same trend as CRLB and can track CRLB’s MSE well, which verifies the effectiveness and superiority of the proposed method.

7. Conclusions

Studies have shown that PM2.5 is strongly linked to other health problems such as respiratory diseases and heart and lung diseases. Protecting the environment to protect humans from air pollution, developing effective and reliable PM2.5 prediction models, and assisting environmental management departments in formulating effective environmental protection measures are urgent and indispensable requirements. PM2.5 concentration chaotic time series data is influenced by multiple variables such as temperature and wind speed in prediction, and in order to improve the prediction accuracy, a multi-step local Volterra multivariable chaotic time series prediction model based on phase clustering is proposed. Through reasonable selection of forecasting model input variables, and combining with clustering phase points, effectively solve the traditional local prediction model selection in the adjacent point problems of unreasonable and time-consuming long, and the proposed model has greatly increased the multi-step prediction step length and precision effectively, reduce the error of cumulative effect, and greatly reduce the multi-step prediction time, has broad application prospects and value in air quality forecast.