Reconstruction of Dynamical Forecasting Model between Western Pacific Subtropical High Area Index and Its Summer Monsoon Impact Factors Based on the Improved Self-Memorization Principle

1Research Center of Ocean Environment Numerical Simulation, Institute of Meteorology and Oceanography, PLA University of Science and Technology, P.O. Box 003, No. 60, Shuanglong Road, Nanjing 211101, China 2Department of Biological and Agricultural Engineering, Texas A & M University, College Station, TX 77843, USA 3Zachry Department of Civil Engineering, Texas A & M University, College Station, TX 77843, USA


Introduction
The western pacific subtropical high (WPSH) is one of the most important components of the East Asian Summer Monsoon (EASM) system [1].The intensity and position of WPSH show complex seasonal evolutions and the changes in them greatly affect the occurrence of rainy season in China, including floods, droughts, and heavy rains [1].For example, when WPSH reaches the northernmost position, especially in summer, it significantly influences rainfall over China and Japan [2].
Owing to its dominance on the East Asian climate, WPSH has become one of the leading topics of interest in atmospheric sciences [3][4][5].Over the past decades, much effort has gone into uncovering the forecast of the WPSH [6], especially the forecast of abnormal WPSH [7,8].Current forecasts for the WPSH can be divided into two categories: numerical forecasts and statistical forecasts.Numerical forecasts are widely used throughout the world; examples include the numerical forecast products of the European Centre for Medium-Range Weather Forecasts Model [9] and the Japanese FUFE502 numerical forecast products.However, numerical forecasts require field boundaries and the complex computations and low efficiency mean that the results are very unstable.Numerical forecast products have better results for large-scale weather systems, but for mesoscale weather systems, such as the WPSH, the results are less good and the forecast time is short.
Statistical methods, in contrast, have achieved some success in forecasting the WPSH.These methods can use historical data and the computation is simple.However, there are some inherent flaws in statistical methods.Using neural networks as an example, it is difficult to objectively determine the number of hidden layer neurons and the training process tends to predict a local optimum, which will limit the forecast accuracy [10].Moreover, the reliability of all these methods is gradually reduced with increasing forecast time, so the forecast results and credibility become very low after two weeks [11].Statistical forecasting products and numerical forecasting products both have some degrees of bias.In particular, error is more obvious in WPSH anomalies and longterm forecasting [3].So the prediction of unusual activities of WPSH within season and the long-term trend forecast of WPSH has become difficult problems in recent years.
A statistical-dynamical model of a weather system is reconstructed from actual data and can be used to describe the physical mechanisms of a complex weather system.Concerning the questions of local convergence of errors, Zhang et al. [12] introduced genetic algorithms (GA), which is widely used in a lot of fields [13,14], to improve the determination of root efficiency of model parameters.On that basis, Bai et al. [15,16] carried out research on the reconstruction of a nonlinear statistical-dynamical forecast model of the WPSH and achieved good results.
However, the dynamical prediction equations derived by Zhang et al. [12] and Bai et al. [15,16] greatly depend on the initial value, so the long-term forecast over 15 days diverged significantly and the results were not satisfactory.For the long-term forecast, the model should be improved.Cao [17] proposed the self-memorization principle, transforming the dynamical equation into memorization equation in a broader sense, named a differential-integral equation, wherein the memory coefficients can also be determined by actual data.This method has been widely used in prediction problems in meteorology, hydrology, and environmental field [18][19][20].Because this method avoids the problem of initial condition in differential equations, it can be introduced to improve the proposed dynamical forecast model.
The set of self-memory function is relatively simple [17] and is suitable for cyclical and linear systems.For nonlinear systems, especially chaotic systems, forecast results are unsatisfactory [20].As the atmosphere and ocean are nonlinear systems, the self-memory function is needed to be modified for nonlinear system modeling.The largest Lyapunov exponent is introduced to improve the traditional self-memorization function.Finally, the improved dynamical forecasting model of WPSH with a new self-memorization function is developed.The improved function not only takes into account the chaotic characteristics of the nonlinear system, but also absorbs the information of past observations.
In our study, we firstly define the WPSH area index (SI) as a measurement of the scope and form of the WPSH.The rest of this paper is organized as follows.Section 2 introduces data and factors.Three summer Monsoon factors which have higher correlation with SI are chosen.The dynamical model of SI and three factors is reconstructed in Section 3. In Section 4, self-memorization dynamics is introduced to improve the reconstructed model.The improved selfmemorization functions for chaotic systems are investigated and discussed in Section 5. Forecast experiment of 2010 and further forecast tests of an additional nine WPSH abnormal years are described in Section 6, and results are summarized in Section 7. Observed long-wave radiation (OLR) from May to October of the 30 years from 1982 to 2011 of National Oceanic and Atmospheric Administration (NOAA) satellites, with a resolution of 0.5

Research Data and Factors
These data are primarily used to calculate the SI and summer Monsoon impact factors in Section 2.3.

The Basic Facts of WPSH Activity in the Summer of 2010.
Change of WPSH within seasons in various years is very different compared with the average, especially the "abnormal" activities of the WPSH, which often lead to East Asia subtropical circulation anomalies and extreme weather events in China in some years, such as 1998, 2003, 2006, and 2010 [21, 22].
In order to better describe the changes in the WPSH, SI is used, which has been defined by the Central Meteorological Observatory [23] to characterize the WPSH range and intensity, that is, the average grid points of 500 hPa geopotential height greater than 588 gpm in the range [10 ∘ -90 ∘ N, 110 ∘ E-180 ∘ E].The higher the value is, the wider the range represented is, or the greater the intensity is.
SI of summer half-year (from May to October) from 1982 to 2011 can be calculated, and changes of SI in different years are shown in Figure 1.In Figure 1, the straight line represents the average SI in the observed years.SI varied widely in different years, representing irregular activities of the WPSH; particularly, in some years, SI showed dramatic changes.For example, in 2010 SI is 32, while in 1984 SI is just 6, which indicates that WPSH activities behaved abnormally in these particular years.If we use the average data of several years, these exceptions will not be apparent.So we should choose an abnormal SI year to study.
From Figure 1, 2010 is the most prominent year of WPSH anomalies as shown by the arrow.From May to October in 2010, SI is above the mean and it is the peak of nearly 30 years.The anomaly of WPSH strength resulted in a very unusual climate in China in 2010.Extreme heat and heavy precipitation events occurred frequently.The strength and wide range of these events are rarely seen in history.In particular, the most powerful WPSH happened from June to August and is to date the strongest one in the meteorological records (Figure 2).This most powerful WPSH directly resulted in rare flood disasters in the eastern part of south China, the Yangtze River, and northeast and northwest of China.Therefore, we selected the summer of 2010 as a typical case to analyze the relevance between the enhanced WPSH and the members of the Monsoon system.

Selection of Three Factors.
According to previous studies [21,22], there are many members of the summer Monsoon system, 21 of which are closely related to the WPSH.If all are used for modeling, the equation would be too complex.Therefore, the correlation analysis method is used to analyze the correlation between these factors and SI.(Specific definition of each factor can be seen in Xue et al. [24] and Yu et al. [22].Values of these factors can be calculated from CFSR data as in Section 2 and are not described in detail here.) Based on the above correlation analysis, three factors with the best correlation with the SI are filtered out for further study, which are The correlation coefficients between the above three factors and SI are 0.77, −0.80, and 0.86, all of which can reach above 0.75.The Mascarene high in the southern hemisphere enhances the WPSH early in the year, which is a positive correlation and a very close relationship, consistent with previous research [24].The close relationships among TH, J1V, and SI are basically consistent with previous research [22,25].

Reconstruction of the Dynamical
Model Based on GA Takens [26] set forth and tested the basic idea of reconstructing the dynamical system from the time series of observed data in his phase space reconstruction theory.Hence, nonlinear dynamics research has entered a new stage, where the evolution information of the dynamical system can be extracted from the time series, such as the calculation of fractal dimension and Lyapunov exponent [27].However, these studies help understand what a chaotic system is and what the periodic system is.Dynamical models for practical problems have not yet been developed fully.Huang and Yi [28] proposed a method of nonlinear dynamic model reconstruction from the actual data and tested it with the Lorenz system, the result of which was satisfactory.Therefore, we introduce this idea and improve it for reconstructing a dynamical model of the SI and three factors.
The principle of dynamical model reconstruction has been introduced in the previous studies [12,15,16] and is not described in detail here but can be found in Appendix A.
The existing parameter estimation methods (such as neighborhood search method and least square estimation) are mostly one-way search which needs to travel the entire parameter space, so the searching efficiency is low.Because of the limitation of the error gradient convergence and its dependence on initial solution, parameter estimation is prone to fall into local optimum, rather than global optimum.GA is a method extensively used as a global optimization method which has developed in recent years.GA has been found to be excellent in global search and parallel computing, and error convergence rate is greatly improved, so it is helpful in parameter optimization.
We can use a simplified second-order nonlinear dynamic model to depict the basic characteristics of the atmosphere and the ocean interactions [29].In our paper,  1 ,  2 ,  3 , and  4 are used to represent time coefficient series of SI, MH, TH, and J1V; then GA is introduced to reconstruct the dynamical model and for parameter optimization.
With  = ( − )  ( − ) minimum as a restriction, an optimization solution searching is made in the model parameter space by GA.
Suppose that the following nonlinear second-order ordinary differential equations are taken as the dynamical model of reconstruction and reconstruction; we choose the time coefficient series of SI, MH, TH, and J1V from May 1 to July 31 in 2010 as the expected data to optimize and retrieve model parameters: Suppose that the parameter matrix  = [ 1 ,  2 , . . .,  9 ;  1 ,  2 , . . .,  9 ;  1 ,  2 , . . .,  9 ] in the above equation is the population, the minimal residual sum of squares  = ( − )  ( − ) is the objective function, the individual fitness value is   = 1/, and the total fitness value is  = ∑  =1   .The idiographic operating steps include coding and creating of initial population, calculation of fitness, male individual choice, crossover, and variation.For a complete theoretical explanation, one can refer to the literature [30].During calculation, the step length is one day.After about 35 times of genetic operations and optimization search, it is possible to converge to the target adaptive value rapidly and retrieve each optimized parameter of the dynamical equations.After eliminating weak items with little dimension coefficient, the nonlinear dynamic model of SI and three factors are inverted as follows: The dynamical reconstruction model should be tested.So we chose the data of August 1 in 2010 of SI, MH, TH, and J1V which do not participate in the modeling as the forecast initial data.Then the Runge-Kutta method is used to carry out numerical integration of the above equations and every step of integration can be regarded as a day forecasting result.The forecast results of four time series within 25 days can be obtained.The correlation coefficients between forecast results and actual results of SI, MH, TH, and J1V within 25 days are, respectively, 0.5659, 0.5746, 0.6091, and 0.5023 and root mean square errors (RMSE) are 36.22%,34.54%, 32.71%, and 36.78%.This indicates that the results of long-term forecast within 25 days are unsatisfactory for the dynamical reconstruction model.That is because the integration results may diverge significantly with time, as well as our initial data for integration being relatively simple.Therefore, we should improve the model to get better long-term forecast results.
Following these analyses, correlation coefficients () were evaluated using the -test ( = 0.05).All resulting coefficients were found to be statistically significant at the 95% confidence level and were thus deemed to be statistically acceptable.

Introduction of Self-Memorization Dynamics to Improve the Reconstructed Model
From the above discussion, we know that the accuracy of forecast results of the dynamical reconstruction model within 25 days is unsatisfactory.Literature suggests that introducing the principle of self-memorization into the mature model can improve long-term forecasting results [18,20].So the selfmemorization principle is introduced to improve the model.Following the study of Cao [17] and Feng et al. [19], the mathematical principle of self-memorization dynamics of systems is shown in Appendix B.
Using (2) as the dynamical kernel, the improved model based on the self-memorization equation (B.12) in Appendix B can be expressed as (3) Because we predict the value of SI, the first formula in (3) is used to get results of our final prediction.That is, If we can get the value of , , (4) can be used for prediction.Values of ,  are associated with the self-memorization function , which is defined in the next section.
We make a prediction with the self-memorization equation (4); the model uses the  values before  0 , therefore () in the self-memorization equation ( 4) in memorizing the +1 values of .This explains the reason why we call () the selfmemorization function.This is the mathematical basis for introducing the self-memorization principle.
The physical basis for introducing the self-memorization principle is that the thermodynamic equation is one of atmospheric motion equations, which implies that the atmospheric motion is an irreversible process.An important contribution to the study of irreversible process is to introduce the memory concept to physics.So the atmospheric development in the future is not only related to the state at the initial time, but also related to states in the past, which means the atmosphere does not forget its past.

Improved Self-Memorization Functions for Chaotic Systems
Cao [17] referred to the problem of self-memorization function assignment.Based on experience and tests, he found the memory level gradually decreased with the distance of initial prediction time, so the memory function can be expressed as It shows that the self-memorization function (5) indicated by Cao only considers the characteristics of memory level gradually decreasing with the distance of initial prediction time and ignores nonlinear characteristics of the selfmemorization function.So we modify the self-memorization function as Term  −(  −  ) represents the characteristics of memory level gradually decreasing with the distance of initial prediction time, while   (1−) reflects the nonlinear characteristics of selfmemorization function. and  are the parameters to be determined, and they are important for the specification of the assigned self-memorization function ().In Cao's study, he found no good way to determine parameters  and , only believing that they had powerful connection with the past observational data.
The Lyapunov exponent (LE) has long been used to study atmospheric and oceanic predictability [31].In chaotic dynamical systems theory, the largest Lyapunov exponent can be portrayed as a whole (long-term) average forecast error growth rate, which generally describes the divergence of nonlinear chaotic systems.So it is often introduced in the forecast study of chaotic systems such as atmospheric and oceanic systems [32,33].The traditional self-memorization function is useful in the linear periodic system forecast, but not appropriate in the forecast of nonlinear chaotic systems.The largest Lyapunov exponent (LLE) is one of the indexes which can represent the characteristics of chaotic systems.So here  can be considered as LLE of the system, which can be calculated from the reconstructed dynamical equation in Section 3.And  can be obtained through GA from historical data.

Determination of Parameter 𝑟.
Based on the above analysis, the solution of  requires the solution of the LLE in dynamical equation (2).

Definition of LLE.
The definition of the Lyapunov spectrum and the largest Lyapunov exponent (LLE) can be referred to the literature [34] and are not described in detail here.The LLE determines the notion of predictability for a dynamical system.A positive LLE is usually taken as an indication that the system is chaotic (provided some other conditions are met, e.g., phase space compactness).

The Principle of Lyapunov Exponent Spectrum Calculations.
If the dynamic differential equations are known, through theoretical derivation or using a numerical iterative algorithm to discretize differential equations, the exact Lyapunov exponents of known dynamical systems can be obtained.First, we work out the solution of the ordinary differential equations and get the Jacobi matrix.Then, we decompose the Jacobi matrix by QR while necessary orthogonal reformings in a multiple of small time intervals are performed.Finally, the Lyapunov spectrum of the system can be obtained by iterative calculations.The specific calculation principle of the Lyapunov exponents spectrum is not shown here, which can be referred to von Bremen et al. [35].

Lyapunov Spectrum of Dynamical Systems and Determination of 𝑟.
According to the above calculation method, the Lyapunov spectrum of the reconstruction dynamical model (see Figure 3) of SI and three factors can be obtained as in Figure 2 which shows that convergence speed and volatility are not too large and they are relatively stable.Final Lyapunov exponents are [0.3816,0.0016, −0.5715], containing both negative Lyapunov exponent and two positive Lyapunov The number of operations exponents, which demonstrate our dynamical system is indeed a chaotic system.
Finally,  is taken as the LLE, and  = 0.3816.

Determination of Parameters 𝑎.
Not only should the nonlinear characteristics of chaotic system be considered, but also the impact of the past actual data on self-memorization function should be taken into account.So parameter  should be optimized.Equation ( 4) is written as The initial value is (  ) =   ,  = ,  − 1, . . .,  − .
Our goal is to make the smallest error between the fitting values and the observed values of  periods before; that is, secondary quality index reaches a minimum: where x is the estimated value of the prediction model.x is the actual value.When further constraint is added, the selfmemorization function should be taken as Parameter  can be obtained by solving for the optimal solution of (8) in the constraint condition of (9).Setting () = ∑  =+1 ( x(  ) − x(  )) 2 as the objective function of GA, putting |(  )| ≤ 1 as a constraint, the specific process of GA is discussed in Section 3.
We can see that the value of parameter  is related with the  data before forecasting time .The forecasting value is preserved as the early data for the next prediction.With the change of data, the value of parameter  is varied.The first calculation is  1 = 0.244 through 39 iterations of GA, and then in next prediction step, parameter  is varied, which we do not list one by one here.In fact, the continuous adjustment of parameter  is more accurate.
Self-memorization function  is determined, and then it is substituted into (4), which can be used to do prediction.

Effective Steps and Maximum Effective Computation Time.
In order to further study the improved model, the algorithm to solve for effective steps and maximum effective computation time should be given, based on the dynamic core (2) of the improved model.Currently the optimal search method is used, which is achieved based on the size of the differences among many numerical solutions.Its principles can be seen in Kirkpatrick [36].
In step interval [ℎ min , ℎ max ], (ℎ min is a very small number),  steps are selected, ℎ  = ℎ min +  *   < ℎ max ( = 1, 2, . . ., ), and   is known as the time increment.The  steps are integral to  at the same time, getting the  number of numerical solutions, denoted by ỹ ().If their differences   () are less than the given tolerance limit , it means that they are relatively close to the true solution at  time.At that time the effective time step interval is [ℎ min , ℎ max ] and the width is  ℎ () = lg ℎ max − lg ℎ min .Otherwise, if   () > , that means values between some steps show larger deviations, so removing these steps from ỹ () will make the difference of the remaining  1 ( 1 < ) solutions less than .The optimum is to exclude the least number.
The effective step intervals which meet such conditions can be gotten.Given integral step, we continue to do integral calculation to  +   ; here,  is an integer.Integral operation continues until there are only two adjacent steps ℎ  ( 1 ), ℎ +1 ( 1 ) left.Then step interval [ℎ  ( 1 ), ℎ +1 ( 1 )] is divided into  parts for a another round of integral operation, until the difference of the solution in +1 step is smaller than .In this case, the maximum effective computation time of the dynamical kernel of the improved model ( 2) can be obtained, as shown in Figure 4.
From Figure 4, when  is certain, for our system, the maximum effective computation time increases with the increment of   and finally remains at 28.7, which is the maximum effective computation time (Figure 4(a)).When   remains constant, the maximum effective computation time oscillates with  (Figure 4(b)).

Prediction Experiment of 2010.
The period from August 2 to September 5 in 2010 is selected to do predictions (35 days in total).August contains abnormal activities of WPSH.Here, the forecasting result is obtained by the sum of (4), which is called as step-by-step forecast.The specific procedure is shown as follows.
Step-by-step forecasts are made after retrospective order  is fixed.That means when we forecast the SI value of August 2, we must get   based on the previous  + 1 time SI data and ( 1 ,  2 ,  3 ,  4 ) based on the previous  time SI, MH, TH, and J1V data (because  2 ,  3 , and  4 are MH, TH, and J1V).Then using them in (4) can obtain the SI value of August 2.Then, taking the SI value of August 2 as the initial value for the next prediction step, the SI value of August 3 can be generated, and so on.

Determination of 𝑝.
When the principle of selfmemorization is introduced, the retrospective order  is related with self-memorization of the system [17].If the system "forgets" slowly, which means parameters  and  are smaller, a high level of  should be used.WPSH abnormal activities are on ten-day scale [22,24], which is a slow process in contrast to the large-scale atmospheric motion.So parameters  and  are smaller, and generally  is in the range of 5 to 15.
The correlation coefficients and root mean square errors (RMSE) between forecast value and real value can be calculated with different retrospective order, as shown in Table 1.
From the table, we can see when  = 8, the correlation coefficient is the largest, and the root mean square error is the smallest.So we choose  = 8.
After  is determined, the improved self-memorization equation ( 4) can be used for prediction experiments.Short-, medium-, and long-term integration forecasts of 15 days, 25 days, and 35 days are carried out.Retrospective order  = 8 means that earlier observational data ( + 1 = 9) are used for integration to begin.The integral result per day is preserved as preliminary information and will be used for integration in the next period.

Prediction Results of 2010.
The forecast results within 35 days are shown in Figure 5, which show that the forecast result of SI is better.
As can be seen in Figure 5, forecast performance of the first 15 days is better.The correlation coefficient has reached 0.9542 and the root mean square errors (RMSE) are 2.45%.Two peaks and one valley of SI are also forecast very accurately.The forecast time series from 15 days to 25 days gradually diverges, but the trend is accurate.The correlation coefficient reaches 0.9254 and the RMSE is 6.37%.Forecast rate is still accurate.The peak of SI in August 18 is also forecast accurately.After nearly 25 days, the RMSE begins to increase: the forecast trend is still accurate with correlation coefficient of 0.8136.Forecast curve from 25 days to 35 days is accurate compared with the actual situation.But the peaks and valleys are not forecast accurately.In particular after 28 days, the error has started to increase significantly.For example, the forecast value of August 28 is nearly 100 more than the actual value, resulting in a false peak.The root mean square error increases to 19.18%, but it is still controlled within 20%. Figure 5 shows that the forecast results after 28 days are obviously unsatisfactory with the occurrence of greater oscillations.The effective forecast period (within 28 days) happens to be in accordance with the maximum integration time (28.7 units) calculated in Section 6.1.
In brief, from Figure 5, we can see that the short-term forecast within 28 days is accurate, and the average RMSEs are less than 8%, which indicates that the reconstruction model can perform accurate predictions of the changing trends of indices.

More Forecasting Experiments of WPSH Abnormal Years.
To further test the forecasting performance of the improved model, cross testing of more experiments is performed.From Figure 2 discussed in Section 2.2, we choose another four years (1998, 2006, 1987, and 1983) in which WPSH intensity is abnormally strong (bigger SI) and five years (1984, 2000, 1994, 1999, and 1985)

Summary and Discussion
7.1.Summary.Although in recent years the WPSH has become one of the leading topics of atmospheric sciences, it is still difficult to be forecast.Combining dynamical system reconstruction idea with the improved self-memorization principle, a new dynamical forecasting model of SI index is developed.The approach of this paper consists of the following steps.
(1) We use correlation analysis method to discuss contributions and impact of the key members of EASM to the abnormal WPSH in 2010 and identify the three most significant factors: the Mascarene high, the Tibetan high (eastern type), and Monsoon circulation index at Bay of Bengal.
(2) Take the SI, MH, XZ, and J1V time series and consider them as trajectories of a set of 4 coupled quadratic differential equations based on the dynamic system reconstruction idea.Parameters of this dynamical model are estimated by a genetic algorithm (GA).
(3) Apply the "principle of self-memorization" in order to improve the forecasting results of the above dynamical model.This involves a manipulation of the time series using an exponentially decaying (in time) function (), which also depends on the largest Lyapunov exponent of the above dynamical system.A second free parameter of () is determined by minimizing the distance of the modified reconstructed time series.
(4) The improved model is used to forecast the SI index.Through 2010 experiments and other 9 experiments of WPSH abnormal years, we find that forecast results within 25 days are good, which proves that our improved model has better long-term forecasting results.Given the complexity of the mechanism of the WPSH [3], the new dynamical forecasting model has some scientific significance and practical value.
Based on the discussion in Section 1, we know that the forecast results and credibility of the previous methods are very low after two weeks [11].So our improved statisticaldynamical model represents an exploration of and supplement to the traditional numerical forecasting and statistical forecasting methods and also extends the prediction time.
7.2.Discussion.Why are the forecasting results of our improved model especially good?(1) Previous studies have showed that the long-term forecasting results of dynamical system reconstructed model are unsatisfactory; thus, we introduce the self-memorization principle to improve the long-term forecasting results.Previous studies have proved that this approach is feasible.For example, Gu [18] used the self-memorization principle to improve the traditional T42 model.And Wang et al. [25] also carried out dynamical prediction of building subsidence deformation with selfmemorization model.Long-term forecasting results of their models were very good.Our study also shows that introducing the self-memorization principle to improve a mature model can bring better long-term forecasting results.(2) The largest Lyapunov exponent is firstly introduced to improve the traditional self-memorization function, making it more appropriate to describe the chaotic systems, such as WPSH.
(3) In developing the improved model, the parameters are obtained from the historical data, which contain sufficient information of WPSH abnormal process.Statistical data combining with improved dynamical model makes forecasting results more reliable.
Although the forecast results of improved model are better, there are still some issues which are needed for further research.
(1) The physical meaning of the factors in the model is not clear, so its dynamical characteristics should be further analyzed.
(2) The forecast accuracy has a great relationship with self-memorization functions.Maybe we can find a better self-memorization function to improve forecast accuracy for long term in the future.
These two items are the focus of our next work.
Coefficients of the above generalized unknown equation can be identified through inverting the observed data.Given a vector , vector  can be solved to satisfy the above equation.It is a nonlinear system with respect to ; however it is a linear system with respect to  (assume  is unknown).So the classical least square method can be introduced to estimate the equation and the regular equation    =    can be derived by making the residual sum of squares  = ( − )  ( − ) minimum.
Based on the above approach, coefficients of the nonlinear dynamical systems can be determined and the nonlinear dynamical equations of observed data can be established.

B. Mathematical Principle of Self-Memorization Dynamics of Systems
In general, dynamical equations of a system can be written as

Figure 1 :
Figure 1: SI of summer half-year (from May to October) from 1982 to 2011.

Figure 2 :
Figure 2: The average monthly geopotential height field at 500 hPa of July and August in 2010.

Figure 3 :
Figure 3: The Lyapunov spectrum of the reconstructed dynamical model.

Figure 4 :
Figure 4: The maximum effective computation time changing with time increment   (a) and tolerance limit  (b).

Figure 5 :
Figure 5: The 35-day forecast results of the subtropical high area index.

Table 1 :
The correlation coefficient (C.C.) and root mean square errors (RMSE) when the retrospective order  is different.

Table 2 :
The correlation coefficients (C.C.) and root mean square errors (RMSE) between forecast value and real value of different events of WPSH abnormal years.
in which WPSH intensity is abnormally weak (smaller SI) to carry out integral forecast experiments of SI.In accordance with the previous idea, the common 2010 model which is reconstructed in Section 4 is used to carry out the cross testing.Forecast results of different time periods (1-15 days as short term, 16-25 days as medium term, and 26-35 days as long term) are compared with the actual situation.Cross test results are displayed in Table2.From the table, we can see the forecast results of the short term and the medium term are accurate, and those of the long term (>26 days) can also be accepted.The results of MH, FLH, and TH are similar to those of SI.