The frequent financial critical states that occur in our world, during many centuries have attracted scientists from different areas. The impact of similar fluctuations continues to have a huge impact on the world economy, causing instability in it concerning normal and natural disturbances [1]. The anticipation, prediction, and identification of such phenomena remain a huge challenge. To be able to prevent such critical events, we focus our research on the chaotic properties of the stock market indices. During the discussion of the recent papers that have been devoted to the chaotic behavior and complexity in the financial system, we find that the Largest Lyapunov exponent and the spectrum of Lyapunov exponents can be evaluated to determine whether the system is completely deterministic, or chaotic. Accordingly, we give a theoretical background on the method for Lyapunov exponents estimation, specifically, we followed the methods proposed by J. P. Eckmann and Sano-Sawada to compute the spectrum of Lyapunov exponents. With Rosenstein's algorithm, we compute only the Largest (Maximal) Lyapunov exponents from an experimental time series, and we consider one of the measures from recurrence quantification analysis that in a similar way as the Largest Lyapunov exponent detects highly non-monotonic behavior. Along with the theoretical material, we present the empirical results which evidence that chaos theory and theory of complexity have a powerful toolkit for construction of indicators-precursors of crisis events in financial markets.
The frequent financial critical states that occur in our world, during many centuries
have attracted scientists from different areas. Such events appear to be the
embodiCopyright © 2020 for this paper by its authors. Use permitted under Creative Commons License Attribution 4.0 International (CC BY 4.0).
ment of chaos and chaotic behavior that has been the subject of research in various
fields, especially in economics and finance. The impact of similar fluctuations
continues to have a huge impact on the world economy, causing instability in it concerning
normal and natural disturbances [
Further events such as the Russian crisis in 1998, the Argentinian crisis in 2001,
and the global financial crisis in 2008 [
Similar conclusions were made based on many articles related to the paradigm of
nonlinear data analysis. For example, Zeman compared the chaotic behaviors of
Thailand in 1998 and Greece in 2013 in terms of economic indicators like GDP,
unemployment, exports, government debt, etc. without any further analysis [
The stock markets are a kind of a complex system with all kinds of interactions
that represent nonlinear characteristic in its dynamics. With the first contribution [
And yet the provided results show that the financial time series represent some little
evidence of chaotic dynamics. Other researchers [
In our contribution, we investigate changes in the dynamical properties of the financial datasets before a crisis event occurs using Lyapunov exponents (LE), which recommended themselves as a great tool for chaos quantification and indication. The key idea in this contribution, which we will ahead to throughout our research, is that the trajectory of the system and its complexity must diverge before the crisis state.
This divergence of the system should have the corresponding degree of chaotic patterns that our indicator can detect and monitor. Such an advantage allows us to use these instruments in predictive settings.
This paper is organized as follows: The brief list of literature devoted to chaos in finance and chaos detection in it with LE is given in Section 2. Definition of the LE methods that we use for its computation, and empirical results we present in Section 3. Some concluding remarks and future perspectives are given in Section 4.
With the high growth in computer science, computer simulations of complex and
chaotic systems become increasingly appreciated. For at least two decades, with
development in numerical computations and quantitative analysis, no doubt left that
chaos theory suggests the same unstable fluctuation that may be as common as the
extreme events and critical transitions in financial markets. For instance, Scheinkman
and LeBaron [
Lyapunov exponents are a natural first choice in exploring and indicating such
chaotic behaviors that occur in it. They do not only classify the system but also tell us
the limits of predictability of the chaotic system [
Sarkar S and Chadha V. [
Srinivasan S. and others in their paper [
As it is seen, chaos theory and its tools remain a huge challenge for researchers of
different fields of science and, namely, in the financial industry, and, as it was
suggested in [
Throughout our research, following the line of a growing body of literature, it was
noticed that despite the huge number of the research papers related to the topic of
Lyapunov exponents, there is a small number of papers devoted to the topic of their
construction as predictive indicators in the stock market. Therefore, relying on
methods and tools described in our previous papers [
Lyapunov exponent is a measure of the exponential rate of nearby trajectories in the phase-space of a dynamical system. In other words, it quantifies how fast converge or diverge trajectories that start close to each other, quantifying the strength of chaos in the system. The divergence of such trajectories can be defined as | (t) || (0) | et (1) where denotes the Lyapunov exponent; (t) is the distance between the reference point and its nearest neighbor after t iterations, (0) is the distance between the reference point and its nearest neighbor perturbed with small error at t 0 .
In such cases when our system is multi-dimensional, we have as many Lyapunov exponents as the dimensions in it. The existence of at least one positive Lyapunov exponent is generally seen as a strong indicator of chaos. Positive LE means that initially similar, phase space trajectories that are sensitive to initial conditions and diverge exponentially fast, characterize chaotic behavior of the system. Negative LE responds to the cases when trajectories remain close to each other, but it is not necessarily implied stability, and we have to examine our system in more detail. Zero or very close to zero exponents indicate that perturbations made along the trajectory neither diverge nor converge. Exactly the largest Lyapunov exponent is used to quantify the predictability of the systems, since exponential divergence means that in the system where the initial difference was infinitesimally small, start to rapidly lose its predictability, behaving differently. However, it should be noted that other exponents also contain important information about the stability of the system, including the directions of convergence and divergence of the trajectories [32].
With the great interest in LE, more and more methods and proposals for their calculating have appeared. Unfortunately, there has not been obtained accepted and universal method for estimating the whole spectrum of Lyapunov exponents from a time series data. One of the most common and popular algorithms have been applied by Wolf et al. [33], Sano and Sawada [34] and later improved by Eckmann et al. [35], Rosenstein et al. [36], Parlitz [37] and Balcerzak et al. [38]. Here, we followed the methods proposed by J. P. Eckmann and Sano-Sawada to compute the spectrum of Lyapunov exponents. With Rosenstein’s algorithm, we compute only the Largest (Maximal) Lyapunov exponents from an experimental time series. As again suggested by Eckmann et al [39] one of the measures from recurrence quantification analysis can be considered for estimation of the Largest Lyapunov exponent since it detects in a similar way highly non-monotonic behavior. 3.1
Eckmann et al. Method Firstly, according to the approach [35], we need to reconstruct attractor dynamics from a single time series {xt }tN1 . For this purpose, according to the delay embedding theorem of Takens [40], we need to choose embedding dimension dE . and after this, we construct dE - dimensional orbit representing the time evolution
X (ti ) [x(ti ), x(ti 1), x(ti 2), . . . , x(ti (dE 1))], i 1, 2, . . . , N dE 1 . Then we have to determine the most neighboring trajectories with X (ti ) : (2) (3) X (t j ) X (ti ) max { x(t j ) x(ti ) } .
0 dE 1 We sort the x(ti ) so that x(t(1) ) x(t(2) ) . . . x(t(N) ) and store the permutation and its inverse 1 . Then, we try to find the neighbors of x(ti ) in dimension 1 by looking at k 1(i) and scan the x(s) for s k 1, k 2, . . . and k 1, k 2 . . . until x(t(s) ) x(ti ) r . For chosen embedding dimension dE 1 , we select the value of s for which further condition is true
x(t(s) ) x(ti ) r, 0,1, 2, . . . dE 1.
Our next goal is to determine the dM dM matrix M i with a matrix dimension dM dE , which describes time evolution of small vectors that surround trajectory X (ti ) and how they map onto X (ti m) trajectory after m iterations. The dimension dM is chosen to avoid undetermined values in M i . Due to this, we have larger step size m and then, associate with X (ti ) a dM - dimensional vector such as
X (ti ) [x(ti ), x(ti m), . . . , x(ti (dM 1)m)].
Accordingly to the algorithm, it is assumed that dE (dM 1)m 1 , therefore, m (dE 1) / (dM 1). In the case when m 1, we need to estimate matrix M i which best satisfies
Mi (X (t j ) X (ti )) X (t j m) X (ti m). (4) The M i is then defined by the linear least-square method [41]. The last step of the algorithm is the classical QR matrix decomposition to find orthogonal matrices Qi and upper-triangular matrices Ri with non-negative diagonal elements such that
M1imQi Qi1Ri1, for i 0, 1, 2, . . .
In order to calculate dM Lyapunov exponents, the equation for the k th Lyapunov exponent with K number of points on the attractor, for which the Jacobian has been estimated, the diagonal eigenvalues of the matrix Ri , and the sampling step is given by: 1 K 1
ln((R(i) )kk ).
k mK i0 Thus, with linearizations by using the diagonal elements from the QR decomposition, we can calculate Lyapunov exponents.
The calculation results for the MLE on the example of index DAX are presented in Fig. 2.
(a) (b) Let us pay attention to the behavior of max at the moments of the known failures noted in the list of stock market crashes and bear markets. Definitely, we can see that in the pre-crisis period, the value of MLE decreases markedly, then increases in the post-crisis period. Rosenstein’s algorithm [36] uses the delay embedding method that reconstructs the most important features of a multi-dimensional attractor into a single one-dimensional time series of some finite size N . The reconstructed trajectory can be presented as (t1), (t2 ),. . . , (tM ) . For the time series {xi}1iN , each delay embedded vector (ti ) will be presented similarly to the vector of the form (2) with embedding dimension dE and time delay m . Then in the reconstructed trajectory we initialize searching for in the state space for the nearest neighbor (t j ) of the trajectory (ti ) : i (0)
min ji mean period (t j ) (ti ) , (5) where
is the Euclidian norm, and j i mean period denotes additional constraint that nearest neighbors have temporal separation greater than the mean period which can be calculated as the reciprocal of the mean frequency of the power spectrum, although authors of this method make a remark that they expect any comparable estimate. Such a condition gives us the possibility of considering each pair of neighbors as nearby initial conditions for different trajectories.
From (1), we have already known that the distance between states (ti ) and (t j ) will grow in time accordingly to a power law i (k) c ek where is a good approximation of the highest Lyapunov exponent. For further estimations, we look at the logarithm of the distance trajectory lni (k) (k t) ln ci , where i (k) is the distance between ith pair of the nearest neighbors defined in (5) after k time steps, ci is the initial separation of them and t is the time interval between measurements (sampling period of the time series).
Further result of this algorithm is not a numerical value, but a function of time: y(k, t) 1 1 M
ln i (k) , t M i1 where M N (dE 1) is the size of the reconstructed time series, and i (k) represents a set of approximately parallel lines, each with a slope roughly proportional to the maximal exponent. Then, it is proposed to be calculated as the angle of inclination of its most linear section. Finding such a section turns out to be a non-trivial task, and sometimes it is impossible to specify such a section at all. Despite this problem, Rosenstein’s method is easy for implementing and computing.
The MLE behavior for a window procedure with windows of different lengths is shown in Fig. 3.
(a) (b) It can be seen that, as before, MLE is also sensitive to the crisis conditions of the stock index. The method of Sano and Sawada [34] that is known in the literature as the Jacobian method, deals as previous ones with reconstructed trajectory as presented in (2). We will assume again that we are dealing with embedded vectors (ti ) [x(t1), x(t2), . . . , x(tidE 1)] in dE - dimensional space. Accordingly to (4), we denote by some operator that transfers trajectories from states (ti ) to (ti1) .
The goal is to choose a sphere with a sufficiently small radius in the phase space trajectory . After m number of iterations, some operator m transforms this sphere into an ellipsoid with a1, . . . , ap semiaxes. If the system has s number of the positive Lyapunov exponents, then the sphere will stretch along the axes a1, a2, . . . , as . Having an acceptable radius , operator m is going to be close enough to the sum of the shift operator and the linear operator which eigenvalues we need to estimate. Then, by averaging these eigenvalues over the entire attractor, we get an estimation of the spectrum of Lyapunov exponents.
For these purposes, suppose that we have a vector (ti ) for which we need to obtain the set of vectors {(tkj }1 jN that falls into the neighborhood of ith trajectory within the ball of some radius :
{y j} {(tkj ) (ti ) | (tkj ) (ti ) }, where y j is the displacement between (tkj ) (ti ) and Euclidean norm is defined as w (w12 w22 . . . wd2) for some vector w (w1, w2, . . . , wd ). After the evolution of a time interval m t , the trajectory (ti ) will be mapped to (tim) such as their neighbors. Then the displacement vector y j will proceed to
{z j} {(tkj m ) (tim ) | (tkj ) (ti ) }.
If the corresponding requirements were accomplished, then the evolution of y j to z j can be described by the operator i . For its optimal estimation, the least-square error algorithm can be a plausible procedure which minimizes the average of the squared error norm between z j and i yj with respect to all components of the matrix i as follows:
1 N 2 miin S miin N j0 z j i y j .
Denoting the component of matrix i by akl (i) , where k is a row, l is a column of the matrix and applying condition (6), we obtain dM dM system of equations of the form: iV C, (V)kl 1 N
y kj ylj , (C)kl N j1 1 N
zkj ylj , N j1 (6) (7) where V, C are dM dM - dimensional matrices, y kj is the k th component of the vector y j , and zkj is the k th component of the vector z j . If A is considered to be a solution of the equations above, then the spectrum of Lyapunov exponents can be calculated by the following formula: 1 n
ln Aieij , j nlim n i1 where A is the solve of the equations (7), and {ej} is the set of basis vectors in the tangent space (t j ) .
The MLE dynamics for the DAX index is presented in Fig. 4. Note that the presented method is the worst of those considered in the sense of sensitivity to crashes. Given that it requires a rather long time series to obtain a positive sign for MLE. Therefore, it can hardly be recommended for use as an indicator of crashes. Recurrence plots (RPs) have been introduced to study dynamics and recurrence states of complex systems [42, 43]. Similar to previous examples, a phase space trajectory (Fig. 5a) can be transformed from a time series into time-delay structures.
RP is a plot representation of those states which are recurrent (Fig. 5b). The recurrence matrix and the states are considered to be recurrent if the distance between them within the - radius. In this case, the recurrence plot is defined as:
Rij ( xi x j ), i, j 1, . . . , N, and
is a norm (representing the spatial distance between the states at times i and j ), is a predefined recurrence threshold, and is the Heaviside function (ensuring a binary R).
For the quantitative description of the system, the small-scale clusters such as diagonal and vertical lines can be used. The histograms of the lengths of these lines are the base of the recurrence quantification analysis [43].
Different elements of RP are distinguished and used, introducing different quantitative measures of complexity of recurrence diagrams. For our purposes, linear sections (lines) of the diagrams are important, which are consecutive sets of individual points. The black dots represent the recurrence of the dynamical process determined with a given resolution , and their organization characterizes the recurrence properties of the dynamics. A vertical line of length l starting from a dot (i, j) means that the trajectory starting from x j remains close to xi during l 1 time steps. A diagonal black line of length l starting from a dot (i, j) means that trajectories starting from xi and x j remain close during l 1 time steps, thus these lines are related to the divergence of the trajectory segments. The average diagonal line length
L lNlNllmmiinnlPP((ll)) is the average time that two segments of the trajectory are close to each other, and can be interpreted as the mean prediction time. Here P(l) is a histogram of diagonal lines of length l .
Another measure considers the length Lmax of the longest diagonal line found in the RP, or its inverse, the divergence,
Lmax max({li}iNl1), and Div 1 / Lmax, where Nl llmin P(l) is the total number of diagonal lines. These measures are related to the exponential divergence of the phase space trajectory. The faster the trajectory segments diverge, the shorter are the diagonal lines and the higher is the measure Div . Therefore, the measure of Div , according to Eckmann [39], can be used to estimate the largest positive Lyapunov exponent.
The comparative dynamics of the Div max measure and the DAX index are presented in Fig. 6.
(a) (b) A comparative analysis of the measures under consideration revealed an obvious advantage of the recursive measure. In addition to the smoothness of the measure itself, it can be calculated for windows of small sizes, which leads to inaccurate or incorrect results for other methods.
This research has made it clear that all three indices represent the behavior of deterministic and chaotic behavior in which the majority of crashes can be identified using the Maximal Lyapunov Exponent.
Fig. 7a presents the comparative dynamics of the daily values of the selected indices DAX, HSI, and S&P 500 with the considered dates of the main crashes. Fig. 7b illustrates the dynamics of the absolute values of the three highest Lyapunov Exponents for method Eckmann and the Maximal LE calculated by method Rosenstein ( R max ) and method Sano-Sawada, implemented in the well-known Tisean package ( T max ).
As can be noticed from Fig. 7b, for the method Rosenstein ( R max ) we have the slightly positive MLE values which can assure us of the chaotic and deterministic behavior which is peculiar to this market. The Sano-Sawada method gives negative MLE values, which most likely indicates method errors for short time series. Calculations show that positive MLE values stably indicate a chaotic picture only starting from windows larger than 1500. 4
The stability problems of financial markets in general and stock markets in particular, deserve more attention in order to ensure their stability and minimize losses as a result of critical changes. The methods of nonlinear dynamics make it possible to identify special states of complex dynamical systems, classify them, and indicate possible trajectories of motion. One of such universal tools for nonlinear dynamics is the spectrum of Lyapunov exponents, the largest of which determines the rate of spread of the trajectories of a dynamical system in phase space. If it is positive, the system is unstable and can assume chaotic states. The disadvantage of many classical methods for determining LE was the need to have a sufficiently long time series, otherwise, the results were irreproducible or incorrect. Moreover, it is interesting to observe the change in LE over time, identifying its characteristic changes. Comparing them with those for the initial series, we can try to predict the possible states of the system under study.
In this work, we have demonstrated the possibility of using LE as an indicator of stock market crashes. In the pre-crisis period, LE is markedly reduced, signaling a more predictable state. In a crisis, the growth of LE indicates the growth of the chaotic component of the market. Particularly promising are the relatively new methods for calculating LE based on the recurrent properties of the system and providing acceptable accuracy for short time series. Of interest is the scale-dependent version of the LE [44], to which we plan to devote a separate article. 31. Belinskyi, A., Soloviev, V., Semerikov, S., Solovieva, V.: Detecting stock crashes using Levy distribution. In: Kiv, A., Semerikov, S., Soloviev, V., Kibalnik, L., Danylchuk, H., Matviychuk, A. (eds) Proceedings of the Selected Papers of the 8th International Conference on Monitoring, Modeling & Management of Emergent Economy (M3E2EEMLPEED 2019) Odessa, Ukraine, May 22-24, 2019. CEUR Workshop Proceedings, vol. 2422, pp. 420–433. http://ceur-ws.org/Vol-2422/paper_34.pdf. 32. Nychka, D., Ellner, S., Gallant, A.R., McCaffrey, D.: Finding chaos in noisy system. Journal of the Royal Statistical Society. 54(2), 399-426 (1992). doi: 10.2307/2346135 33. Wolf, A., Swift, J.B., Swinney, H.L., Vastano, J.A.: Determining Lyapunov exponents from a time series. Phys. D Nonlinear Phenomena. 16(3), 285–317 (1985). doi: 10.1016/0167-2789(85)90011-9 34. Sano, M., Sawada, Y.: Measurement of the Lyapunov spectrum from a chaotic time series.
Phys. Rev. Lett. 55(10), 1082–1085 (1985). doi: 10.1103/PhysRevLett.55.1082 35. Eckmann, J. P., Kamphorst, S. O., Ruelle, D., Ciliberto, S.: Liapunov exponents from time series. Physical Review A, 34(6) (1986). doi: 10.1103/PhysRevA.34.4971 36. Rosenstein, M.T., Collins, J.J., De Luca, C.J.: A practical method for calculating largest Lyapunov exponents from small data sets. Phys. 65(1–2), 117–134 (1993). doi: 10.1016/0167-2789(93)90009-P 37. Parlitz, U.: Identification of true and spurious Lyapunov exponents from time series. Int. J.
Bifurc. Chaos. 2(1), 155–165 (1992). 38. Balcerzak, M., Pikunov, D., Dabrowski, A.: The fastest, simplified method of Lyapunov exponents spectrum estimation for continuous-time dynamical systems. Nonlinear Dyn. 94, 3053-3065 (2018). doi: 10.1007/s11071-018-4544-z 39. Eckmann, J.-P., Kamphorst, S.O., Ruelle, D.: Recurrence plots of dynamical systems, Europhys. Lett. 5, 973–977 (1987). 40. Takens F: Detecting strange attractors in turbulence. Springer, Coventry, England (1981). 41. Wilkinson, J. H., Reinsch, C.: Linear Algebra. Springer-Verlag, Berlin (1971). 42. Recurrence Plots and Their Quantifications: Expanding Horizons. Webber C.L., Marwan N.(eds.). Proceedings of the 6th International Symposium on Recurrence Plots, Grenoble, France, 17–19 June 2015, 180, 1-387. Springer, Heidelberg (2016). 43. Marwan, N., Romano, M.C., Theil, M, Kurths, J.: Recurrence plot for the analysis of complex systems. Phys. Rep. 438, 237–329 (2007). 44. Gao, J., Hu, J., Tung, W.-W., Zheng, Y.: Multiscale analysis of economic time series scale-dependent Lyapunov exponent. Quantitative Finance 13(2), 1-10 (2011).