The System Identification and Prediction of the Social Earthquakes Burst in Human Society

The internal mechanism of the social “earthquakes”– the wars burst in human society, is investigated by nonlinear theory. We explore the scaling behavior in the human activity and the historical evolution of ancient China. Based on the war chronology and history chronology from B.C.2000 to A.D.1840, using statistical analysis together with the detrended fluctuation analysis, we find that the social evolution system shows self-organized critical behavior and there exists self-similar scaling behavior during the evolution system. Based on time-delay reconstructed phase space, a modified radial basis function neural network (modified RBF-NN) is used to predict the number of wars and the transformation of the historical stage. It shows that the prediction result matches quite well with the real value for the following two historical stages. Furthermore, by singular value decomposition and the sparse identification algorithm in time-delay coordinates, the coefficients of the candidate functions in identified system can be solved from a sparse regression problem. We extracted the governing equations of big events burst in social system. The identified system fits well with the real data, and it can capture the same topology with the original system of wars. In this paper, we investigated the wars burst in social system by statistic analysis, neural network as well as sparse identification. It sheds light on the study of wars or conflicts from perspective of data science and dynamical theory.


I. INTRODUCTION
Avariety of physical systems exhibit regular behavior in different time or space scales, such as the burst of the earthquakes on faults and crackling noise of paper crumpling [1], [2]. The behaviors are likely independent events, but there is universality for real systems which can be explained by simple models. If there is a slow driving force, the system crackles accompanied with various sizes of discrete events, such as the earthquakes burst at a critical state as the tectonic plates interacting with each other under a driving plate [3], [4]; We have explored the plastic deformation mechanism during compressive deformation for disordered The associate editor coordinating the review of this manuscript and approving it for publication was Junhua Li . materials such as bulk metallic glasses (BMGs) and high-entropy alloy and give stress-strain signals prediction [5]- [7]. The serrated flow in plastic dynamics manifests self-organization critical (SOC) state [8]- [12]. The intermittent serrated flow burst in the disordered materials during plastic deformation is induced by the interaction of the shear bands, which manifests as different size of sliding events [13]. The distribution for the sizes of the sliding events follows a power law. In the past decades, scientists have developed models and theories for understanding this kind of scale-invariant behavior in driven, nonlinear dynamical systems. The founder of the field cliodynamics, Turchin presented that we must transform history into an analytical, predictive science so that we can understand how human society evolves and how to develop a healthy society [14], [15].
Radial basis function neural network (RBF-NN) is a feed-forward neural network with excellent performance [16]. For example, it has good performance in learning and approximating nonlinear functions. So we use it for prediction based on the reconstruction of the time series. In this paper, we use a modified radial basis function neural network (modified RBF-NN), for which the input data is firstly reconstructed into a high dimension phase space. The prediction result could be related with the data in the dynamical system.
In this paper, we explore hidden order in the wars applying a scaling behavior analysis and the modified RBF-NN to study the dynamics of historical empires. Meanwhile, we extract the governing equation by system identification based on the collected data. To explore the hidden order in human activity of ancient China, we investigate the time series of the war data to characterize the dynamical behavior in the system of human activity which has not been revealed clearly by nonlinear theoretical analysis. The study on the complex dynamics of human activity is significant, which contributes to understand the history and predict the big events in the long run. The government management departments can make reasonable strategy to achieve stable development of the society. The war is not a single separated event. In fact, there is internal dynamical correlation between the units in the system. Richardson found that there exists power law distribution in the frequency versus size of wars [17]. Roberts reported that the power law is generated by self-organized criticality (SOC) [18]. Researchers dedicated to find out the dynamics rules in human activities by studying the timing of human activities, human dynamics [19], [20]. D. H. Tang et al. explored wars occurred in last two dynasties Ming and Qing (A.D.1367-A.D.1911), which shows that the distribution of recurrent time of wars obeys a stretched exponential form [21]. The hidden order in wars and the social system is still unclear and lack of quantitative and qualitative dynamic analysis. There are complex dynamical behaviors in the system of wars in human society, which has not been clearly investigated by nonlinear theory (for example, the dynamical structure is unclear).
The wars burst in human society are similar with the avalanches occurred in the sand pile model [3]. The society can be considered as a complex system composed of many interacting units (such as persons or small groups), and wars can be considered as the collapse of the dynamics, and we call the collapse as ''social earthquake''. There are some important questions such as the scaling behavior of the wars in long time-scale and the prediction of the burst of the big events during the evolution of the society, including the formation time of a new historical stage and the number of wars during a time interval. In order to study wars quantitatively, in this paper, we investigate the number of wars in each dynasty and the time interval for each new dynasty generated in ancient China during about four thousand of years (from B.C.2000 to A.D.1840), during which the country is almost a closed system. Statistical analysis show that the number of wars in different dynasties and the time intervals of dynasties display power law distribution, which means that the dynamics of wars is related with self-organized criticality. Meanwhile, the detrended fluctuation analysis suggests that there is scale-free behavior in the dynamics. It confirms the self-organization critical state. Furthermore, the prediction for the burst of wars in different dynasty and the formation time of a new dynasty are given by applying a modified radial basis function neural network and time-delay phase space reconstruction.
Most of the equations of the system are often established according to the physical laws or rules. There are often a finite number of state variables in the system equations considering some main factors. For the experimental data, most of the study is based on data analysis. There is still a gap between dynamical system and the data. The sparse identification in time-delay coordinates can extract the governing equation of the system, which is can be used to give further analysis on human social system.
In the current work, the novel aspect relates to three factors: (1) We investigate the power law and the scaling behavior in the human evolution system based on the number of wars burst in each dynasty and the time intervals of each dynasty.
(2) Using a modified RBF-NN based on time-delay phase space reconstruction, we predict the time interval for the each dynasty and the number of wars burst in each dynasty interval. The results are quite accordant with the real value.
(3) In time-delay coordinates, we apply the sparse identification of nonlinear dynamics method to extract the dynamics of the social system, and the identified system fits well with the system generated from real data.

II. RELATED WORK
To reveal the hidden order in the time series, we have explored the scaling behavior in the serrated flow in bulk metallic glasses using the fractal, muti-fractal, the detrended fluctuation analysis [9], [11]. For the time series of the human social activity, it still needs clear nonlinear dynamical analysis and prediction. In recent years, the machine learning has been used to extract statistical or dynamical characteristics for prediction. Several artificial neural networks (ANNs) [22], [23], multi-input-multi-output network [24], [25], and deep learning method [26] have been used to extract the hidden information in time series data. Wang proposed a double-layer recurrent neural network to predict the PM2.5 value [27]. Wu established a novel framework with dedicated combination of data prediction, which is achieved by applying the Least Mean Square (LMS) dual prediction algorithm [28]. In data mining application, machine learning is a powerful method for the analysis of high dimensional data. He et al. proposed a decision graph based outlier detection (DGOD) method for social big data [29]. Wu reported an unsupervised linear dimension reduction algorithm named critical points preserving projection (CPPP) [30].
In [7], we presented a delay parameterized method (DPM) for the forecasting of low dimensional mid-term chaotic time VOLUME 8, 2020 series. The correlation function was introduced to bridge time series and hidden order of the system. We use the traversal algorithm, intelligent algorithm including particle swarm optimization and genetic algorithm to calculate the optimal parameters for prediction. The method was applied in Lorenz chaotic time series, stress-strain signals in material science and stock K-line maps, showing pretty good predictions. Based on the Takens embedding theorem [31], the time series can be reconstructed into a 2m + 1 dimensional phase space in the form of (t) = y(t), t(t + τ ), · · · , y(t + (2m)τ )), where τ is the time delay, y(t) is the time series. The prediction of the time series can be implied by tracing the trajectory in the reconstructed space. Then the correlation function is presented to predict for multi-dimensional coupling chaotic time series.
Take the linear correlation function for prediction as an example, let the time series < y 1 , y 2 , · · · , y m >, the phase space constructed from the multiple observations can be defined as where N is the length of each observation. The phase space reconstructed by the kth observation function y k , is where τ is seen as delay parameter. Y, Y k are m-dimensional phase spaces. According to the Theorem 1 in [7], the linear correlation where P is a m × m parameter matrix. The unknown data in Eq. (2) are P i,j and y k (N + k), k = 1, 2, · · · , N + (m − 1)τ , i.e., the state space Y. Define the first N − (m − 1)τ lines of Y k and Y asȲ k andȲ, the matrix P can be solved from, Then the prediction of {y k (N +k), k = 1, 2, , (m−1)τ } can be achieved by Y k = YP. The determination of the parameters is calculated by the least square method, where Y is the transpose of Y. Note that, Eq. (3) could be an over-determined equation, and (Ȳ Ȳ ) −1 is the inverse ofȲ Ȳ , andȲ Ȳ could take the pseudo inverse. The prediction by the nonlinear correlation function and multi-step prediction have been discussed in our previous work. The above prediction method is effective for the deterministic systems based on low dimensional mid-term time series. There is a question that can it predict the low dimensional short-term time series? Unfortunately this method is tough to make prediction for the short-term time series owing to the severely insufficient information.
In this paper, based on the data of wars, we reconstruct the time series into a phase space and predict applying RBF neural network. The model structure is shown in Fig. 1. The prediction is related with the points in the phase space, which relies on a dynamical system. Meanwhile, we use the sparse identification methods in the time delay coordinates to extract the dynamical system. There are several approaches to discover the dynamical system from data, such as symbolic regression [32], [33], sparse regression [34] and compressed sensing [35], [36]. Gaussian process regression is used for uncertain measurement and prediction of time series data [37]. Brunton et al. demonstrated the sparse identification algorithm on linear, nonlinear oscillators and the chaotic Lorenz system, as well as the high-dimensional PDE models. So far, there is no literature involving the dynamical equation for the burst of wars in human society. In the present work, we develop the sparse identification algorithm to discover the governing equations underlying a dynamical system simply from the data of wars. There are many critical data-driven problems, for example, the predicting and suppressing the spread of disease. Data-driven discovery of dynamics will continue to play an important role in human society.

III. SELF-ORGANIZED CRITICALITY AND THE SCALING BEHAVIOR ANALYSIS
Take China as an example considering the long history with abundant data, the data of wars in ancient China are recorded during the year B.C.2000 to A.D.1840 [38]. The number of wars in each year is shown in Fig. 2. The statistical distribution is investigated based on the number of wars δ N (i) in i-th dynasty period [ Fig. 3(a)]. Since the measurement of the social avalanche size is not possible to take, we use the number of wars occurred in i-th dynasty period, to reasonably represent the avalanche size. From the time series of δ N (i), [ Fig. 3(a)], every value of δ N (i) can be constructed to characterize the social energy relaxation process. The time series {y(i) = δ t (i), i = 1, 2, · · · , n − 1} is also investigated, which is the time interval between each historical turning point, i.e. the time of a new dynasty established. [ Fig. 3(b)].  We explore the distribution of avalanche size, which is an important fingerprint reflecting the avalanche of human dynamics. The statistical distribution of the normalized δ N (i) shows a power-law distribution: Fig. 4(a), insets], suggesting the dynamics manifests self-organized critical state (SOC). Furthermore, the statistical distribution of the normalized also shows a power-law  Fig. 4(b), insets], which implies the evolution of the SOC state in time scale. Self-organized criticality was first introduced by Bak et al. to explain the behavior of the cellular-automata model [39], [40]. In the model, particles are randomly dropped onto a square grid of boxes. Once the number of particles in a box exceeds the threshold value, the box will get released and all the particles are redistributed to the four nearest boxes or lost off at the edge of the grid. The avalanche comes to an end and the system evolves into a new stable configuration. Iterate above process, there will be different scales of avalanches. The ''avalanches'' in the model satisfied a power-law frequency-area distribution. Other cellular-automata models, such as the slider-block, forest-fire models, and earthquakes are also said to exhibit self-organized critical behavior [41]- [43].
For the war, it can be regarded as the avalanches in the system of human activities. The system is composed of many interacting units. As the system evolving, the social contradictions gradually increase, and this is the energy accumulation process. The establishment of a new emperor is the VOLUME 8, 2020 beginning of the energy accumulation. An avalanche will occur once the ''social energy'' exceeds the threshold, i.e. the war burst. Then the system will collapse in a short time, during which the energy will released and the system will go into a new stable state. The period of energy accumulation (usually for several years) is much longer than the collapse time (usually for several days) of the system. The overlap of the accumulation and relaxation of the energy can result in a hierarchy of spatiotemporal scales, and it manifests as the SOC behavior.
We further investigate the scaling behavior in the human activities based on the time series, {x(i) = δ N (i), i = 1, 2, · · · , n}. The detrended fluctuation analysis is applied to quantify the scaling behavior in the system [44]- [46]. The process of the detrended fluctuation analysis is described as following. Divide the signal {x(i), i = 1, 2, · · · , n} into N q (where N q = N /q) zones with each zone containing q elements. In the k-th zone, the local trend is defined as a linear function of {x k (j), j = 1, 2, · · · , q}, which is linearly fitted based on the original series, {x k (j), j = 1, 2, · · · , q}.
The detrended time series is, {x k (j) − x k (j), j = 1, 2, · · · , q}, with a mean-square error, and the root-mean-square is, zones. F q is a power function of the scale, q, F q ∼ q H , where H is the Hurst exponent. Different with the conventional detrended fluctuation analysis, in the present work, we focus on the scaling behavior in different scales of q, where q is defined as the temporal scale. Based on the number of wars in different dynasty periods, {x(i) = δ N (i), i = 1, 2, · · · , n}, and the time interval of each new dynasty formed, {y(i) = δ t (i), i = 1, 2, · · · , n − 1}, the scaling relation is plotted in Fig. 5. There exists linear relation on different temporal scales of q, which means that there is long range correlation and the self-similarity behavior in the dynamical system.

IV. THE PREDICTION BASED ON NEURAL NETWORK
In the following, we try to predict the number of wars in different dynasties and the time of a new historical stage applying a modified RBF-NN, for which the input data is firstly reconstructed into a high dimension phase space. The neural network contains the input layer, the hidden layer, and the output layer. It can be considered as a mapping between the input and output layer. For example, choose the Gaussian function as the radial basis function [16], where x is the input vector; c j is the center of the j-th RBF; σ j is the j-th RBF width; m is the number of nodes in the hidden layer. The symbol · indicates the Euclidean norm in this case. The input layer and the output layer represent the mapping from x to R j (x), and the linear mapping from R j (x) to y k respectively. The k-th output, y k , is a weighted summation of all RBFs, where K is the number of outputs, m is the number of nodes in the hidden layer, w jk is the weight of the j-th node of hidden layer to the k-th output layer. We predict the next three burst point based on the time interval in each dynasty, {x(i) = δ t (i), i = 1, 2, · · · , n − 1}, where n = 193. First of all, reconstruct the data into a m-dimension phase space, here m = 13, time delay τ = 1, we use the first 174 data to train and obtain the prediction value of y 1 (175). Use the data {x(2), x(3), · · · , x(174)} together with the predicted y 1 (175) in the red rectangle to predict y 2 (176). Repeat the above steps (See Fig. 6), we get the predicted value of y 3 (177). The initial training result and the following three predicted value are plotted in Fig. 7(a). It shows the following two steps prediction is quite close to the actual value. In addition, the prediction based on the number of wars burst in  i-th dynasty, {δ N (i), i = 1, 2, · · · , n} is also investigated, the result is shown in Fig. 7(b). Based on the above results, we obtain the prediction of the recurrent time of wars, i.e. the time of a new historical stage established, as well as the number of conflicts burst in the corresponding historical stage. We use 193 of the data (from the year B.C.922 to A.D.1862). The previous 174 data is used to train the modified RBF neural network, and we give the following 3 steps prediction. To compare model predictions, we add the prediction based on RBF-NN without reconstruction (See Fig. 7). To verify the improvement of the modified RBF-NN with phase space reconstruction, we compare the modified RBF-NN with existing RBF-NN. It is obvious that the prediction based on the modified RBF-NN is better for {δ N (i), i = 1, 2, · · · , n}. In addition, the Mean Squared Error (MSE) and Mean Absolute Error (MAE) are calculated based on the prediction of the time series, {δ t (i), i = 1, 2, · · · , n − 1}, which is shown in Table.1. It shows that the predicted value matches better with the real value in the following two steps by modified RBF-NN. We get the prediction on the time interval of each historical stage, that is to say, the beginning time point of a new historical stage can be forecasted.

V. THE SPARSE IDENTIFICATION OF THE SOCIAL SYSTEM
The above prediction using the modified RBF-NN is based on time-delay reconstructed phase space. In the following, we use the time-delay coordinates, and find the governing system of social system. Applying the Sparse Identification of Nonlinear Dynamics (SINDy) [34], we extract Return predicted values the governing system from the data of the wars, and it will bridge the gap between the data and the nonlinear dynamical system. We will describe the SINDy method and the sparse identification of the social system with time-delay coordinates in detail in the following. For a physical or biological system, the information of the system is hidden in the collected data or experimental data. For most of the system, there exist the governing equations, which can be written as a nonlinear differential system generally,ẋ where the vector x(t) = [x 1 (t)x 2 (t) · · · x n (t)] T ∈ R n is the state of the system at time t, and f is the dynamical constraints. The function f consists of only a few terms for many systems, such as some polynomial functions, or the trigonometric functions. Using x(t) at time t 1 , t 2 , · · · t m , anḋ x(t) can be calculated numerically. The matrix state X and its derivativeẊ (t) are as follows, . . .
It requires an assumption that there are only some active nonlinear functions. We solve the following sparse regression problem,Ẋ where = [ξ 1 , ξ 2 · · · ξ n ], each ξ k represents a sparse vector, which is corresponding to the coefficients of the candidate functions. It is the identified system extracted from data. The algorithm for sparse identification [34] is shown in Table 2.
For one of the row equationẋ k = f k (x) in Eq. (5), the sparse vector ξ k determines the active terms in the right hand side.
can be solved by sparse regression from Eq. (5), and the identified model with each row of the governing equations is constructed as follows, where (x T ) is a vector of symbolic functions of elements of x.
There is an important parameter λ in solving the sparse regression problem. If the component of the vector ξ k is less than λ, this component is set to be 0; If the vector ξ k is greater than λ, this component is preserved. The parameter λ is a threshold value, which can influence the number of active function terms.

Algorithm 2 Algorithm for Sparse Identification
Input: The time series of the number of wars in each year {N (i), i = 1, 2, · · · , n} reconstruct the matrix H ; get V by singular value decomposition.
Using sparse regression to solve: If (| | < λ ); (λ is the sparsification parameter). set the terms of to be zero; else Preserve the value in ; end Output: The coefficients of the candidate functions ξ 1 , ξ 2 and ξ 3 The identified social system applying SINDy in the time-delay coordinates  The Akaike Information Criterion (AIC) is used in the selection of parameter λ [47]. Given a range of the parameter λ and run the sparse identification algorithm, the model with the lowest AIC value is chosen as the VOLUME 8, 2020 identified system. AIC can be represented as where RSS denotes the residual sum of squares of the fitted model, k is the number of parameters, n is the length of the data.
In the following, we explore the system based on the time series of the number of wars δ N (i) in time-delay coordinates, which can synthesize additional dynamic variables from a single variable x i = δ N (i) [48]. In the time-delay coordinates, it produces a new attractor, and the identified system has the same topology with the original system based on Takens' theorem [31]. Constructing a Hankel matrix by stacking delayed time-series of x as rows: Applying the singular value decomposition (SVD), it gives where the columns of V can be considered as a hierarchical set of eigen-time-series. In the calculation, we use the data of x i from t = 1 to t = 194 with t = 0.003, and q = 10 in H . The value of q is possible to be a larger value which means it stacks more rows in H . Here we choose the first three dominant eigen-time-series given by the first three columns of V , and denote the coordinates as u, v, and w. Using fourth-order central difference, the derivativesu,v, andẇ can be compute numerically. Applying the SINDy algorithm in time-delay coordinates, we use the cubic order polynomial as the candidate nonlinear function. The system is identified in the space of polynomials in (u, v, w) up to three order. Set X = [u, v, w] T ,Ẋ = [uvẇ ] T , applying sparse regression, we seek the sparse solution to the system Eq. (5). The coefficients of the candidate functions ξ 1 , ξ 2 and ξ 3 are shown in Fig. 8. The identified social system applying SINDy in the time-delay coordinates are shown in Fig. 9. These coefficients are identified after normalizing the columns of (V ). In the time-delay coordinates, the identified system can capture the same topology with the original system. In the algorithm, the increasing of the polynomial order can results in over-fitting of the identified mode. Since we use the first three column of V to determine time-delay coordinates, there will be a small amount of information missing from the rest lower energy columns. The λ k here we used is 2. Generally, for the system with more variables, considering the variables may vary greatly, we solve the sparse regression with different parameters λ k . It contributes to get better fitting model. There is an advantage of this method is that it does not require large set of data. The identified system can capture the skeleton of the attractor of the original system. From Fig.10, we can see the identified system (black curve) fits well with the real data (red dashed curve).

VI. CONCLUSION AND DISCUSSION
There is hidden order in the evolution of human society although the number of wars and the interval of each historical period seem irregular. We investigate the internal scaling behavior based on the war data of the ancient China, the time interval between each historical turning point, {δ t (i), i = 1, 2, · · · , n − 1}, and the number of wars occurred in i-th historical period, {δ N (i), i = 1, 2, · · · , n}. We find that there is obvious power law in the two time series, which suggests the self-organized critical state. The distribution of war frequency results in a power-law spectrum. The rule is similar to the earthquake in nature, and we have investigated the complex dynamical behaviors in a spring-block model for earthquakes. The sliding speed of the blocks shows a power-law distribution suggesting the self-organized critical state [4], which is a feature in the real observed earthquakes. For the self-organized critical state, the system is in a state relatively controllable when the accumulated energy is less than the threshold. While the system will collapse once the energy accumulated exceeds a certain threshold. The wars or conflicts in human system can be considered as the ''social earthquakes''. Similarly, for the wars or conflicts, once the social contradiction (energy) accumulates to a certain threshold, it will burst wars or big events, which is accompanied by energy release. The system will go into a stable state for a period and evolves to the next new critical state.
Meanwhile, we investigate the scaling behavior applying detrended fluctuation analysis based on the time series, {x(i) = δ N (i), i = 1, 2, · · · , n}. It shows that there is long range correlation in the system. In addition, we predicted when the new historical stage begins and how many conflicts or wars burst in the time interval. The prediction provides early warning, according which the politicians could make effective policy to reduce the contradictions in the social system. It will postpone the self-organization of the system to a critical state, or weaken the social instability. This work provides a theoretical basis for the evolution of human social system, which is conducive to the long-term stability of the society. The multi-step prediction is based on RBF neural network method applying phase space reconstruction, which can be used in other time series analysis, such as the disease warning, natural disaster warning.
Using time-delay coordinates, we introduce the sparse identification of nonlinear dynamics to the social system and extract the governing equation of the system which can capture the same topology with the original system. This is a critical important for the further study of the social system by dynamical system such as the bifurcation theory and stability analysis. For the algorithm in IV, the input data is reconstructed into a phase space and then apply RBF neural network to predict the data. The modified RBF-NN is based on phase space reconstruction and dynamical set. For the algorithm in V, we use time-delay coordinates, the input data is also reconstructed into a phase space. This is something in common in IV and V. The algorithm in IV is focus on the prediction based on the modified RBF-NN. The algorithm in V is mainly to extract the dynamical system which can capture the same topological structure with the original war system. In sparse identification, the coefficient of the candidate function could be time dependent, and it will be helpful to extract more accurate system. In addition, the sparse identification of nonlinear system with control is another issue which needs further consideration.
In this paper we have explored the scaling behavior and given prediction on when the war burst, which is helpful to illuminate the intrinsic mechanism of the human society. We transform the history into an analytical and predictive science applying statistical methods and nonlinear analysis, which will contribute to the study of the social system and benefit to develop a healthy society.
CUN CHEN received the Ph.D. degree in applied mathematics from Zhengzhou University, in 2015. She is currently an Associate Professor with the Henan Academy of Big Data, School of Mathematics and Statistics, Zhengzhou University. Her research interests include differential equation, bifurcation and chaos, machine learning, statistical analysis, and data science. XIAOXIANG GUO received the master's degree from Zhengzhou University, in 2017. He is currently pursuing the Ph.D. degree with the School of Mathematics and Statistics, Zhengzhou University, China. His current research interests include the machine learning and data science.
JINGLI REN received the Ph.D. degree in applied mathematics from the Beijing Institute of Technology, in 2004. She became a Professor at the School of Mathematics and Statistics, Zhengzhou University, since 2006. She is currently the Deputy Dean of the Henan Academy of Big Data. Her research interests include differential equation, bifurcation and chaos, percolation theory, machine learning, and data science. VOLUME 8, 2020