Study of irregular dynamics in an economic model : attractor localization and Lyapunov exponents

This is a self-archived version of an original article. This version may differ from the original in pagination and typographic details. Author(s): Title: Year: Version: Copyright: Rights: Rights url: Please cite the original version: CC BY 4.0 https://creativecommons.org/licenses/by/4.0/ Study of irregular dynamics in an economic model : attractor localization and Lyapunov exponents © 2021 The Author(s). Published by Elsevier Ltd. Published version Alexeeva, Tatyana A.; Kuznetsov, Nikolay V.; Mokaev, Timur N. Alexeeva, T. A., Kuznetsov, N. V., & Mokaev, T. N. (2021). Study of irregular dynamics in an economic model : attractor localization and Lyapunov exponents. Chaos, Solitons and Fractals, 152, Article 111365. https://doi.org/10.1016/j.chaos.2021.111365 2021 Chaos, Solitons and Fractals 152 (2021) 111365 Contents lists available at ScienceDirect Chaos, Solitons and Fractals Nonlinear Science, and Nonequilibrium and Complex Phenomena journal homepage: www.elsevier.com/locate/chaos Study of irregular dynamics in an economic model: attractor localization and Lyapunov exponents Tatyana A. Alexeeva a , Nikolay V. Kuznetsov b , c , d , ∗, Timur N. Mokaev b a St. Petersburg School of Physics, Mathematics, and Computer Science, HSE University, St. Petersburg 194100, Russia b Faculty of Mathematics and Mechanics, St. Petersburg State University, St. Petersburg 198504, Russia c Faculty of Information Technology, University of Jyväskylä, Jyväskylä 40014, Finland d Institute for Problems in Mechanical Engineering RAS, St. Petersburg 199178, Russia a r t i c l e i n f o Article history: Received 28 July 2021 Accepted 18 August 2021 Keywords: Lyapunov exponents Lyapunov dimension Chaos Unstable periodic orbit Absorbing set Mid-size firm model a b s t r a c t Cyclicality and instability inherent in the economy can manifest themselves in irregular fluctuations, including chaotic ones, which significantly reduces the accuracy of forecasting the dynamics of the economic system in the long run. We focus on an approach, associated with the identification of a deterministic endogenous mechanism of irregular fluctuations in the economy. Using of a mid-size firm model as an example, we demonstrate the use of effective analytical and numerical procedures for calculating the quantitative characteristics of its irregular limiting dynamics based on Lyapunov exponents, such as dimension and entropy. We use an analytical approach for localization of a global attractor and study limiting dynamics of the model. We estimate the Lyapunov exponents and get the exact formula for the Lyapunov dimension of the global attractor of this model analytically. With the help of delayed feedback control (DFC), the possibility of transition from irregular limiting dynamics to regular periodic dynamics is shown to solve the problem of reliable forecasting. At the same time, we demonstrate the complexity and ambiguity of applying numerical procedures to calculate the Lyapunov dimension along different trajectories of the global attractor, including unstable periodic orbits (UPOs). ©2021 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ ) 1. Introduction Increasing uncertainty, unpredictability, and instability in the world, nature cataclysms, a series of economic crises, self-fulfilling expectations which give rise to bubbles and crashes, as well as rapid development and implementation of digital technologies in everyday life have posed a number of new challenges for scientists, governments, and policy makers: to study, understand and interpret the behavior of complex dynamical systems, including socioeconomic models [1–3] . An inherent component of observed economic processes is cyclicality, which is manifested through the occurrence of various types of fluctuations in the economic system under consideration. In particular, regular fluctuations could be either periodic boombust phenomena associated with predictable changes in some elements of the economic system that reappear at fairly constant time intervals, or seasonal fluctuations that are permanent in nature. Regular and stable periodic oscillations lead to the predictable ∗Corresponding author at: Faculty of Information Technology, University of Jyväskylä, 40014 Jyväskylä, Finland. dynamics of the process’ model and are quite simple to describe mathematically. A number of straightforward quantitative measures, such as phase-frequency characteristics and amplitude, can be calculated for them. However more often, economic systems exhibit irregular (including chaotic) behavior. The role of irregular oscillatory dynamics for forecasting and stabilization of economic processes significantly depends on the source and nature of these fluctuations. On the one hand, irregular economic fluctuations could be the result of unusual events such as large bankruptcies, oil and currency shocks, floods, strikes, civil unrest, epidemics, etc. These events could be thought of as initiated by exogenous shocks. On the other hand, irregular fluctuations could be generated by endogenous mechanisms inherent in the very nature of economic systems. Thus, there are two ways to examine of irregularity in the economy. First approach takes into account random processes that are considered in the model as exogenous shocks. Second one is based on identification of a deterministic endogenous mechanism of occurrence of irregular fluctuations, which may also be chaotic. These two approaches were developed in economics literature in parallel and https://doi.org/10.1016/j.chaos.2021.111365 0960-0779/© 2021 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license ( http://creativecommons.org/licenses/by/4.0/ ) T.A. Alexeeva, N.V. Kuznetsov and T.N. Mokaev Chaos, Solitons and Fractals 152 (2021) 111365 generated a lot of discussion regarding the views on the sources of irregular fluctuations (see, e.g. [4–6] ). Since the 1970s, there has been keen interest in the study of deterministic chaotic dynamics in economic models within the framework of the second approach. This research was stimulated by the discovery of chaos in dynamical systems by Lorenz [7] , Ueda et al. [8] . Many famous economists (see, e.g. [4,9–27] ) have suggested numerous examples of economic models in which qualitatively and quantitatively reasonable irregular fluctuations might occur in purely deterministic settings. For instance, the larger literature [9,13,14,19,20,28–31] examines the endogenous cycles and irregular chaotic dynamics which could be generated by deterministic, equilibrium models of the economy. The models often exhibit complex dynamics characterized by both chaotic behavior and instability. Such combination suggests a nonlinear dynamical system, somewhat unstable at the core, but effectively contained further out. The contribution of these models has been to demonstrate the compatibility of endogenous irregular fluctuations with equilibrium dynamics in economics. At the same time, theoretical tools were developed for effective chaos control, which, by small finetuning the parameters of system, made it possible to stabilize selected orbits embedded in a chaotic attractor and nudge the dynamics toward a desired trajectory. Examples applications of these tools can be found in [32–40] . The reviewed literature showed the relevance of chaos for economic models and contributed to development of advanced mathematical tools for study of complex nonlinear dynamical systems in economics, which continues up to now. During the last few years, highly influential authors published a number of significant papers (see, e.g. [41–49] ). The studies of models with irregular dynamics have received a new impetus and spread into many subfields of economic theory. Especially, such models offer important contributions in macroeconomics, dynamical game theory, theory of rational inattention, finance, environmental economics, and industrial organization (for survey of the literature, see [50] ). To understand, describe and make measurable the properties of irregular dynamics it is important to calculate its quantitative characteristics. Indicators based on Lyapunov exponents, including such as entropy and dimension, naturally arise in economics [51] . In economic models these characteristics could be considered as indicators of irregular (primarily, chaotic) behavior, as the growth rate of the value of some economic variable (for instance, technology level), or as a measure of costs of making decisions by a rationally inattentive agent who acquires information about the values of alternatives through a limited-capacity channel (see, e.g. [52–55] ). In this paradigm important results and arguments were presented which provide novel support for the idea that business cycles may be largely driven by endogenous deterministic cyclical forces (see, e.g. [6,56,57] ). There are two main approaches in studying this topic. The first approach is based on the possibility of obtaining analytical results for low-dimensional nonlinear models (in the literature, twodimensional dynamical systems are most often studied). The second one is based on the ability to study complex irregular dynamics using numerical procedures. However, the possibility of obtaining reliable results using them is significantly limited due to the necessity of performing calculations only over finite time intervals, rounding-off errors in numerical methods, and the unbounded space of initial data sets [58–63] . It should be noted that the sensitivity to small changes in the initial data, inherent in irregular (chaotic) dynamics, can cause significant forecasting errors. This, on the one hand, can explain some of the difficulties associated with forecasting behavior of the models, and on the other hand could be interpreted as unpredictability in real world problems (see, e.g. [6] ). Trajectories in models of such processes may be attracted not to a stationary point or a periodic cycle, but to an irregular invariant set, including chaotic attractor. Additional complexity of the dynamics can be also associated with various unstable orbits embedded into the chaotic attractor of the dynamical system. Stabilization of unstable orbits makes it possible to improve the forecasting of the model dynamics [63] . Analytical methods allow overcoming these limitations at least for some lowdimensional models (see, e.g. [62,64] ) and are able to mitigate the influence of computer errors. Thus, this is capable of making reliable forecasts of model dynamics and of getting its exact qualitative and quantitative characteristics. We continue the line of research on the limiting dynamics for a mid-size firm model, which began in [62,63] , where we have obtained conditions for the global stability. In this paper we focus on a different approach, associated with the identification of deterministic endogenous mechanisms of irregular fluctuations in economic systems. We use an analytical approach for localization of a global attractor and study limiting dynamics of the model. We estimate the Lyapunov exponents and get the exact formula for the Lyapunov dimension of the global attractor of this model analytically. With the help of DFC, the possibility of transition from irregular limiting dynamics to regular periodic dynamics is shown to solve the problem of reliable forecasting. At the same time, we demonstrate the complexity and ambiguity of applying numerical procedures to calculate the Lyapunov dimension along different trajectories of the global attractor, including UPOs. 2. Problem statement For understanding and reliable predicting the behavior of economic models in continuous time the study of its limit oscillations is an important task. This task could be solved by an analytical localization of the global attractor (whenever applicable) for the corresponding system of ODE, i.e., constructing a bounded closed positively invariant region (an absorbing set). On this attractor, along with the corresponding solution for the system we obtain some estimates of irregular (including chaotic) dynamics. This allows us to calculate various quantitative characteristics based on the Lyapunov exponents such as the Lyapunov dimension of the attractor and entropy. Consider the Shapovalov model proposed in [65] which describes the behavior of a mid-size firm  ˙ x = −σx + δy, ˙ y = μx + μy −βxz , ˙ z = −γz + αxy , (1) where coefficients α, β, σ, δ, μ, γat variables (x, y, z) ∈ R 3 are positive control parameters with the economic meaning. We define this model in terms of the differences between actual levels of the variables X, Y , and Z, denoted the growth of three main factors of production: the loan amount X, fixed capital Y and the number of employees Z(as an increase in human capital), and its potential (natural) levels x p , y p , and z p respectively 1 . Thus, we consider the gap between the actual and potential levels of factors of production: x = X −x p , y = Y −y p , and z = Z −z p , where X, Y , and Zare nonnegative. Note that system (1) describes the behavior of a mid-size firm correctly when the global attractor or its absorbing set lays in the domain x ≥−x p , y ≥−y p , and z ≥−z p . System (1) can be reduced to a Lorenz-like system  ˙ x = −cx + cy , ˙ y = rx + y −xz , where b = γ μ, c = σ μ, r = δ σ, ˙ z = −bz + xy , (2) 1 We assume that the potential (natural) levels of factors of production correspond to the production possibilities of a mid-size firm as a whole, reflecting its natural, technological, and institutional constraints. 2 T.A. Alexeeva, N.V. Kuznetsov and T.N. Mokaev Chaos, Solitons and Fractals 152 (2021) 111365 Fig. 1. Analytical localization of the chaotic attractor of system (2) with parameters set at b = 5 . 7 , c = 18 . 3 , r = 51 by the global absorbing set B = B R  1 , where B R is the ellipsoid (gray), 1 is the parabolic cylinder (brown). Here M = 1 2 1 c + b 2 c = 0 . 1052 , A = 0 . 1111 , and η= 29651 . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) using the following coordinate transformation (x, y, z) →  μ  αβ x, μσ δ αβ y, μσ δβ z  , t → t μ. (3) System (2) in crucial respect differs from the classical Lorenz system [7] in the sign of the coefficient at y in the second equation, which is 1 here and -1 in the Lorenz system. Accordingly, the inverse transformation (x, y, z) →   αβ μx, r  αβ μy, rβ μz  , t → μt (4) reduces system (2) to system (1) with coefficients σ= cμ, δ= rcμ, γ= bμ2 . In addition, system (1) with parameters satisfying the relations σ2 / (σ−δ) = μand δ< σ< μcan be reduced to the well-known Chen system [67]  ˙ x = −dx + dy , ˙ y = ( c −d ) x + cy −xz , with b = γ, c = σ2 σ−δ= μ, d = σ, d < c, ˙ z = −bz + xy , (5) using coordinate substitutions (x, y, z) →  1  αβ x, σ δ αβ y, σ δβ z  . (6) The possibility of reducing system (1) to the Chen system (5) under the above conditions shows the complexity of studying a mid-size firm model. The problem of analytical calculation of the dimension of the attractor for the Chen system remains an issue [66] . It was shown in [62] that for system (2) the global absorbing set B = 1  B R can be constructed under conditions 2 < b < 2 c( Fig. 1 ), where 1 =  (x, y, z) ∈ R 3 | z ≥x 2 2 c  is the parabolic cylinder, B R = {(x, y, z) ∈ R 3 | 1 2 [ Ax 2 −2Mxy + y 2 + 2 Transformations (3) and (4) do not change the direction of time, which is essential for analysis of the Lyapunov dimension and Lyapunov exponents [66] . (z −(r + (A + M ) c −M )) 2 ] ≤η} is the ellipsoid, M = 1 2 1 c + b 2 c , A > M 2 , and η= η(b, c, r, A ) > 0 . The presence of an absorbing set implies the existence of a global attractor A glob , which contains all local self-excited and hidden attractors [68–76] and a stationary set. In the interior of the global absorbing set model (1) can show both regular and irregular limit dynamics depending upon values of model’s parameters [62] . In case of the global stability we observe regular dynamics when all trajectories of system (2) tend to the stationary set { S 0 , S ±} , where S 0 =(0 , 0 , 0) , S ±=(± b(r + 1) , ± b(r + 1) , r + 1) are equilibria of system (2) . As it was shown in [62] , the system is globally stable in the following parameter domain ( b + 1 ) b c −1 < r < b c + 1 ( b −1 ) , 2 < b < 2 c. (7) Thus, in [62] the regular dynamics of system (2) was studied and the conditions of global stability were obtained. On the other hand, if condition (7) is violated, the system may exhibit irregular behavior, at which a chaotic attractor can be reveal. As an example, Shapovalov et al. [65] , Shapovalov and Kazakov [77] , and Gurina and Dorofeev [78] show that system (1) exhibits chaotic behavior for some values of parameters. Localization of a global attractor and furthest calculation of the limit values of the finite-time Lyapunov exponents and the finitetime Lyapunov dimension along various trajectories of this attractor are nontrivial tasks. While trivial attractors (stable equilibrium) can be easily found analytically or numerically, the search for periodic or chaotic attractors can be a challenging problem. For numerical localization of the attractor, one needs to choose an initial point in its basin of attraction. After a transient process, a trajectory, starting in a neighborhood of an unstable equilibrium, is attracted to the state of oscillation and then traces it. Next, the computations are being performed for a grid of points in vicinity of the state of oscillation to explore the basin of attraction and improve the visualization of the attractor. However, for an arbitrary system possessing a transient chaotic set, the time of transient process depends strongly on the choice of initial data in the phase space and also on the parameters of numerical solvers to integrate a trajectory (e.g., order of the method, step of integration, relative and absolute tolerances). This complicates the task of distinguishing a transient chaotic set from a sustained chaotic set (attractor) in numerical experiments. Since the “lifetime” of a transient chaotic process can be extremely long and in view of the limitations of reliable integration of chaotic ODEs, even long-time numerical computation of the finite-time Lyapunov exponents and the finite-time Lyapunov dimension does not guarantee a relevant approximation of the Lyapunov exponents and the Lyapunov dimension [59,61,63] . In this paper, we obtain analytical formula for the exact Lyapunov dimension for global attractor of system (2) . We demonstrate difficulties in numerical computation of the finite-time Lyapunov exponents and the finite-time Lyapunov dimension along one randomly chosen trajectory over a long time interval which are caused by finite precision numerical integration of ODE, UPOs embedded into the attractor, and choice of various initial data. This confirms the significance of the deduced analytical formula for the Lyapunov dimension. 3. Analytical estimation of finite-time Lyapunov dimension and exact Lyapunov dimension In this section, we give the main definitions and explanations. Some definitions, proofs and technical parts used from now onwards in this section are summarized in Appendix. 3 T.A. Alexeeva, N.V. Kuznetsov and T.N. Mokaev Chaos, Solitons and Fractals 152 (2021) 111365 Rewrite system (2) in a form ˙ u = f(u ) , f : R 3 → R 3 , (8) where fis a continuously differentiable vector-function. Let u (t, u 0 ) be any solution of (8) such that u (0 , u 0 ) = u 0 ∈ R 3 exists for t ∈ [0 , ∞ ) . For system (8) the evolutionary operator ϕ t (u 0 ) = u (t, u 0 ) defines a smooth dynamical system { ϕ t } t≥0 in the phase space (R 3 , || ·|| ) : { ϕ t } t≥0 , (R 3 , || ·|| ) , with Euclidean norm. We consider fundamental matrix Dϕ t (u ) = y 1 (t) , y 2 (t) , y 3 (t) , Dϕ 0 (u ) = I, with cocycle property, where { y i (t) } 3 i =1 are linearly independent solutions of the linearized system, Iis the unit 3 ×3 matrix. The finite-time local Lyapunov dimension [59,79] can be defined via an analog of the Kaplan-Yorke formula with respect to the set of ordered finite-time Lyapunov exponents . { LE i (Dϕ t (u )) = LE i (t, u ) } 3 i =1 at the point u : dim L ( t, u ) = d KY { LE i ( t, u ) } 3 i =1 = j ( t, u ) + LE 1 ( t, u ) + ···+ LE j ( ,u ) ( t, u ) LE j ( t,u ) +1 ( t, u ) | , (9) where j(t, u ) = max { m :  m i =1 LE i (t, u ) ≥0 } , dim L (t, u ) = 3 for j(t, u ) = 3 , or t = 0 . If j(t, u ) ∈ { 1 , 2 } , then  j(t,u ) i =1 LE i (t, u ) ≥0 , LE j(t,u )+1 (t, u ) < 0 and dim L ( t, u ) = j ( t, u ) + s ( t, u ) : j ( t,u )  i =1 LE i ( t, u ) + s ( t, u ) LE j ( t,u ) +1 ( t, u ) = 0 . (10) The finite-time Lyapunov dimension is defined as: dim L (t, A ) = sup u ∈A dim L (t, u ) , (11) where A is a compact invariant set. The Douady–Oesterlé theorem [80] implies that for any fixed t > 0 the finite-time Lyapunov dimension on set A , defined by (11) , is an upper estimate of the Hausdorff dimension: dim H A ≤dim L (t, A ) . By the Horn inequality [81, p.50] , cocycle property, and invariance of A we have 3 sup u ∈A (  j 1 LE i (kt, u ) + s LE j+1 (kt, u )) ≤sup u ∈A (  j 1 LE i (t, u ) + s LE j+1 (t, u )) for j ∈ { 1 , 2 } , s ∈ [0 , 1] and any integer k > 0 . The infimum is achieved at infinity, otherwise for d : 0 < dim (T , A ) < d < lim inf k → + ∞ dim (kT , A ) from (10) and the Horn inequality one gets a contradiction: 0 < liminf k → + ∞ sup u ∈A  d 1 LE i (Dϕ kT (u )) ≤liminf k → + ∞ sup u ∈A  d 1 LE i (Dϕ T (u )) < 0 . Thus, the best estimation (11) takes the form Kuznetsov [79] dim L A = inf t> 0 sup u ∈A dim L (t, u ) = lim inf t→ + ∞ sup u ∈A dim L (t, u ) (12) and is called the Lyapunov dimension . If the supremum of finite-time local Lyapunov dimensions on set A is achieved at such an equilibrium point u eq ≡ϕ t (u eq ) ∈ A : dim L A = dim L u eq , then the Lyapunov dimension can be represented in analytical form and it is called the exact Lyapunov dimension in [82] . A conjecture on the Lyapunov dimension of selfexcited attractor [59,61,79] is that for a typical system, the Lyapunov dimension of a self-excited attractor does not exceed the Lyapunov dimension of one of the unstable equilibria, the unstable manifold of which intersects with the basin of attraction and visualizes the attractor. In a general case, analytical computation of the Lyapunov exponents and the Lyapunov dimension is hardly possible. However, they can be estimated by the eigenvalues of the symmetrized Jacobian matrix [80,83] . The KaplanYorke formula with respect to the ordered set of eigenvalues νi (J(u )) = νi (u ) , ν1 (u ) ≥ν2 (u ) ≥ν3 (u ) , 3 see Appendix. i = 1 , 2 , 3 , of the symmetrized Jacobian matrix 1 2 (J(u ) + J(u ) ∗) , J(u ) = Df(u ) [79] gives an upper estimation of the Lyapunov dimension of an attractor A : dim L A = in f t> 0 sup u ∈A d KY { LE i ( t, u ) } 3 i =1 ≤sup u ∈A d KY { νi ( u ) } 3 i =1 . (13) Generally speaking, one cannot get the same values of { νi (u ) } 3 i =1 at different points u ; thus, the supremum of d KY ({ νi (u ) } 3 i =1 ) on A has to be computed. To obtain estimate (13) , it is not necessary to integrate the solutions of the system; however, the analytical estimation of { νi (u ) } 3 i =1 on the attractor may be a challenging task. At the same time, an effective analytical estimation of the Lyapunov dimension via (13) can be obtained by the Leonov method . 4 The inequality dim H A ≤dim L A < j + s holds, if sup u ∈A ν1 ( u, S ) + ···+ νj ( u, S ) + sνj+1 ( u, S ) + ˙ V ( u ) < 0 , (14) where ˙ V (u ) = ( grad (V )) ∗f(u ) , V : R 3 → R 1 is a differentiable scalar function, Sis a nonsingular 3 ×3 matrix, νi (u, S) = νi (SJ(u ) S −1 ) is the ordered set of eigenvalues ν1 (u, S) ≥ν2 (u, S) ≥ν3 (u, S) , i = 1 , 2 , 3 , of the symmetrized Jacobian matrix 1 2 (SJ(u ) S −1 + (SJ(u ) S −1 ) ∗) , j ∈ { 1 , 2 } is an integer number, and s ∈ [0 , 1] is a real number. 4. Main result Using the Leonov method [79,84] we estimate the Lyapunov exponents and obtain the Lyapunov dimension for the global attractor in system (2) . Theorem 1. If for parameters of system (2) the following relations hold 2 < b < 2 c, (15) r > b c + 1 ( b −1 ) , (16) ( b + 1 ) ( b −2 ) b 2 + 6 bc −3 c 2 + b + c ( 2 c −b )  −c b 2 + b −c ( 8 −b ) r ≤0 , (17) then dim L A glob = 3 −2(b + c −1) c −1 +  (c + 1) 2 + 4 cr . (18) Proof. Consider system (2) with the Jacobian matrix J =  −c c 0 r −z 1 −x y x −b  (19) under the conditions (15) and (16) . We apply the transformation (3) with a nonsingular matrix S =  −1 a 0 0 −b+1 c 1 0 0 0 1  (20) to this system, where a = c √ ( 1+ b ) ( c−b ) + rc . Then the symmetrized Jacobian matrix of this system 1 2 SJS −1 + (SJS −1 ) ∗5 has the following eigenvalues λ2 = −b, λ1 , 3 =−c −1 2 ±1 2  ( 2 b + 1 −c ) 2 + a 2 b + 1 c x + y 2 + az −2 b a 2  1 2 . (21) 4 see Appendix. 5 Symbol ∗denotes the transposition of matrix. 4 T.A. Alexeeva, N.V. Kuznetsov and T.N. Mokaev Chaos, Solitons and Fractals 152 (2021) 111365 The inequalities 2 λj −λj+1 ≥− ( −1 ) j ( 2 b + 1 −c ) + | 2 b + 1 −c | ≥0 , j = 1 , 2 , (22) imply λ1 ≥λ2 ≥λ3 . From (21) following [84] we get the ratio 2 ( λ1 + λ2 + sλ3 ) = −( s + 1 ) ( c −1 ) −2 b + ( 1 −s ) ( 2 b + 1 −c ) 2 + a 2 b+1 c x + y 2 + az −2 c a 2 (23) where s ∈ [0 , 1] is a real number. Using the inequality √ k + l ≤ √ k + l 2 √ k , ∀ k > 0 , l ≥0 , we obtain an estimate 2 ( λ1 + λ2 + sλ3 ) ≤−( c −1 + 2 b ) −s ( c −1 ) + ( 1 −s ) ( c + 1 ) 2 + 4 cr 1 2 + 2 ( 1 −s ) [ ( c+1 ) 2 +4 cr ] 1 2  −cz + a 2 z 2 4 + a 2 4 b+1 c x + y 2  . (24) We introduce the function V (x, y, z) = θ(x,y,z) [ (c+1) 2 +4 cr] 1 2 , where θ(x, y, z) = a 2 Q 0 x 2 + a 2 (−c Q 1 + Q 2 ) y 2 + a 2 Q 2 z 2 + a 2 4 c Q 1 x 4 −a 2 Q 1 x 2 z −a 2 P Q 1 xy −c b z, (25) P and Q i (i = 0 , 1 , 2) are some positive real parameters. Then 2 ( λ1 + λ2 + sλ3 ) + 2 ˙ V ≤− ( c −1 + 2 b ) −s ( c −1 ) + ( 1 −s ) ( c + 1 ) 2 + 4 cr 1 2 + 2 ( 1 −s ) [ ( c+1 ) 2 +4 cr ] 1 2 W ( x, y, z ) + ˙ θ, (26) where W (x, y, z) = −cz + a 2 z 2 4 + a 2 4 b+1 c x + y 2 . Choose the parameters P and Q i (i = 0 , 1 , 2) of the function θ(x, y, z) such that F = W (x, y, z) + ˙ θ≤0 , ∀ x, y, z ≥x 2 2 c . (27) Substituting W (x, y, z) and ˙ θin (27) , we get F = A 0 z 2 + A 1 x 2 + A 2 xy + A 3 y 2 , (28) where A 0 = a 2 2 c ( b + P ) Q 1 −2 bQ 2 + 1 4 , A 1 = a 2 ( b+1 ) 2 4 c 2 −rP Q 1 , A 2 = a 2 ( ( c −1 ) P −2 c ) Q 1 + 2 rQ 2 + b+1 2 c −c ba 2 , A 3 = a 2 1 4 + 2 Q 2 −c ( 2 + P ) Q 1 . (29) Then A 0 ≤0 A 3 ≤0 4 A 1 A 3 −A 2 2 ≥0  ⇒ F ≤0 , ∀ x , y , z ≥x 2 2c , (30) ⇔ ⎧ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎩ Q 1 ≤b c ( b + P ) Q 2 −1 8 c ( b + P ) , Q 1 ≥2 c ( 2 + P ) Q 2 + 1 4 c ( 2 + P ) , Q 1 ≥2 2 c + P Q 2 + ( b + c + 1 ) 2 ba 2 −4 c 3 4 a 2 bc 2 ( r + 1 ) ( 2 c + P ) . (31) Since RHS of the second inequality in (31) is positive, we obtain ⎧ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎩ b c ( b + P ) Q 2 −1 8 c ( b + P ) −2 c ( 2 + P ) Q 2 + 1 4 c ( 2 + P ) ≥0 , b c ( b + P ) Q 2 −1 8 c ( b + P ) −2 2 c + P Q 2 + ( b + c + 1 ) 2 ba 2 −4 c 3 4 a 2 bc 2 ( r + 1 ) ( 2 c + P ) ≥0 , (32) Fig. 2. Parameters of system (2) complying with the conditions (15) and (17) . ⇔ ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ( b −2 ) P c ( b + P ) ( 2 + P ) Q 2 −3 P + 2 b + 2 8 c ( b + P ) ( 2 + P ) ≥0 , −( 2 c −b ) P c ( b + P ) ( 2 c + P ) Q 2 + c ( 8 c −b ) r −2 b 3 −4 ( 3 c + 1 ) b 2 −2 + 13 c −6 c 2 b + 8 c 2 P+ 8 bc 2 ( b + P ) ( 2 c + P ) ( r + 1 ) +6 bc 2 r −2 b ( b + 1 ) b 2 + b + 6 bc −3 c 2 ≥0 . (33) It follows from condition (15) that the coefficient at Q 2 in the first inequality of (33) is positive and the coefficient at Q 2 in the second inequality of (33) is negative. Hence, we can reduce (33) to the following inequalities L (b, c, r, P ) ≤Q 2 ≤R (b, c, r, P ) , (34) where L (b, c, r, P ) = 3 P+2 b+2 8 P(b−2) > 0 , R (b, c, r, P )= (c(8 c−b) r−2 b 3 −4(3 c+1) b 2 −(2+13 c−6 c 2 ) b+8 c 2 ) P+ 8 bc (2 c−b)(r+1) P +6 bc 2 r−2 b(b+1)(b 2 + b+6 bc −3 c 2 ) . Inequalities (34) mean that a positive Q 2 exists such that R ( b, c, r, P ) −L ( b, c, r, P ) = −( b + P ) ( k 1 r + k 0 ) 4 bc ( b −2 ) ( 2 c −b ) ( r + 1 ) P ≥0 , (35) where k 1 = −c b 2 + b −c(8 −b) , k 0 = (b + 1) (b −2)(b 2 + 6 bc − 3 c 2 + b) + c(2 c −b) . Since the denominator of fraction (35) is positive, we obtain required condition (17) k 1 r + k 0 ≤0 . (36) This completes the proof.  We obtain a formula for the exact Lyapunov dimension of the global attractor for certain region D of the parameters (b, c) of system (2) ( Fig. 2 ). Here D is the region such that band cin D satisfy (15) and (17) , and ris such that conditions (16) and (17) are held. The same approach allows one to estimate of the topological entropy of the global attractor [60,81,85,86] . To demonstrate significance of this analytical result we compare it with numerical simulations. We discuss the difficulties of numerical procedures for reliable estimation of the Lyapunov dimension and Lyapunov exponents along one randomly chosen trajectory over a long time interval. A natural way to get reliable estimation of the Lyapunov dimension of an attractor A is to localize the attractor A ⊂C, to consider a grid of points C grid on C, and to find the maximum of the corresponding finite-time local Lyapunov dimensions for a certain time t = T . In Fig. 3 is shown the grid of points C grid filling the basin of attraction: the grid of points fills 5 T.A. Alexeeva, N.V. Kuznetsov and T.N. Mokaev Chaos, Solitons and Fractals 152 (2021) 111365 Fig. 3. Numerical localization of the chaotic attractor of system (2) with parameters set at b = 5 . 7 , c = 18 . 3 , r = 51 by the cuboid Cand the corresponding grid of points C grid . cuboid C = [ −27 , 27] ×[ −65 , 65] ×[3 , 95] (containing the attractor) rotated by 45 degrees around the z-axis, with the distance between points equal to 0.5. The time interval considered is [0 , T = 500] at the time points t = t k = τk (k = 1 , . . . , N) , N = 10 0 0 according to the time step τ= t k −t k −1 = 0 . 5 , and the integration method is MATLAB ode45 with predefined parameters. The infimum on the time interval is computed at the points { t k } N 1 . For system (2) with parameters under consideration, we use a MATLAB realization of the adaptive algorithm of the finite-time Lyapunov dimension and Lyapunov exponents computation [59] and obtain the maximum of the finite-time local Lyapunov dimensions at the grid of points ( max u ∈C grid dim L (t, u ) is computed for trajectories of system (2) using MATLAB ode45 integration method with predefined parameters and with threshold parameter ξ= 0 . 01 for adaptively adjusting the number of SVD approximations). For parameters b = 5 . 7 , c = 18 . 3 , r = 51 we get max u ∈C grid dim L ( 100 , u ) =2 . 0808 , max u ∈C grid dim L ( 500 , u ) =2 . 0792 . (37) Note that if for a certain time, t = t k , the computed trajectory is out of the cuboid, the corresponding value of the finite-time local Lyapunov dimension is not taken into account in the computation of the maximum of the finite-time local Lyapunov dimensions. If the maximum of local Lyapunov dimensions on the global attractor, which involves all equilibria, is achieved at an equilibrium point: dim L (u cr eq ) = max u ∈A dim L (u ) , then this allows one to get analytical formula for the exact Lyapunov dimension [82] . The exact Lyapunov dimension dim L A glob = dim L S 0 = 2 . 4347 > dim L A ≈max u ∈C grid dim L (t k , u ) ≈2 . 0808 (see (37) ) obtained by formula (18) and the estimation (37) are consistent with the hypothesis on the Lyapunov dimension of self-excited attractor. Using Theorem 1 we can get the value of the exact Lyapunov dimension on the global attractor, which coincides with the Lyapunov dimension at a stationary (zero) point. This result is nontrivial since to compute reliably numerically the dimensions on the trajectories of the global attractor is extremely difficult. We demonstrate challenging nature of this task by the following examples. Choosing the initial data somewhere in the phase space, we can obtain the values of the dimensions along the various trajectories by a numerical procedure. Generally speaking, these values of the dimensions will also be different. For instance, system (2) has the analytical solution u (t) = (0 , 0 , z 0 e −bt ) which tends to the equilibFig. 4. Period-1 UPO u upo 1 (t) (red, period τ1 = 0 . 69804 ) stabilized using the UDFC method, and pseudo-trajectory ˜ u (t, u upo 1 0 ) (blue, t ∈ [0 , 100] ) in system (2) with parameters set at b = 5 . 7 , c = 18 . 3 , r = 51 . rium S 0 = (0 , 0 , 0) from any initial point (0 , 0 , z 0 ) ∈ R 3 . The existence of such solutions in the phase space complicates the procedure of visualization of a chaotic attractor (pseudo-attractor) by one pseudo-trajectory with arbitrary initial data computed for a sufficiently large time interval. In particular, the numerical computation of finite-time local Lyapunov exponents along this trajectory during any time interval does not lead to averaging of these values across the attractor, but to tending of these values to the finitetime local Lyapunov exponents of S 0 . The challenges of the finite-time Lyapunov dimension computation along the trajectories over large time intervals is connected with the existence of UPOs embedded in a chaotic attractor. Along with the existence of the analytical solution u (t) = (0 , 0 , z 0 e −bt ) the global attractor of system (2) contains a period-1 UPO. Consider system (8) . Let u upo (t, u upo 1 0 ) be its UPO with period τ> 0 , u upo (t −τ, u upo 1 0 ) = u upo (t, u upo 1 0 ) , and initial condition u upo 1 0 = u upo (0 , u upo 1 0 ) . To compute the UPO, we add the unstable delayed feedback control (UDFC) [87] in the following form: ˙ u ( t ) = f ( u ( t ) ) −KB [ F N ( t ) + w ( t ) ] , ˙ w ( t ) = λ0 c w ( t ) + λ0 c −λ∞ c F N ( t ) , F N ( t ) = C ∗u ( t ) −( 1 −R ) N  k =1 R k −1 C ∗u ( t −kT ) , (38) where 0 ≤R < 1 is an extended DFC parameter, N = 1 , 2 , . . . , ∞ defines the number of previous states involved in delayed feedback function F N (t) , λ0 c > 0 , and λ∞ c < 0 are are additional UDFC parameters, B, Care vectors and K > 0 is a feedback gain. For the initial condition u upo 1 0 and T = τwe have F N (t) ≡0 , w (t) ≡0 , and, thus, the solution of system (38) coincides with the periodic solution of initial system (8) . For system (2) with parameters b = 5 . 7 , c = 18 . 3 , r = 51 , using (38) with B ∗= ( 0 , 1 , 0 ) , C ∗= ( 0 , 1 , 0 ) , R = 0 . 7 , N = 100 , K = 10 , λ0 c = 0 . 1 , λ∞ c = −5 , one can stabilize a period-1 UPO u upo 1 (t, u 0 ) with period τ1 = 0 . 69804 from the initial point u 0 = (0 . 1 , 0 . 1 , 0 . 1) , w 0 = 0 on the time interval [0 , 100] (see Fig. 4 ). We use the Pyragas procedure [87,88] for numerical stabilization and visualization of UPOs. For the initial point u upo 1 0 ≈(29 . 6688 , 26 . 1650 , 73 . 8221) on the UPO u upo 1 (t) = u (t, u upo 1 0 ) we numerically compute the trajectory of system (38) without the stabilization (i.e. with K = 0 ) on the time interval [0 , T = 100] (see Fig. 4 ). We denote it by ˜ u (t, u upo 1 0 ) to distinguish this pseudo-trajectory from the periodic orbit u (t, u upo 1 0 ) . On the initial small time interval [0 , T 1 ≈2 τ1 ] , 6 T.A. Alexeeva, N.V. Kuznetsov and T.N. Mokaev Chaos, Solitons and Fractals 152 (2021) 111365 Fig. 5. Period-1 UPO u upo 1 (t) (red, period τ1 = 0 . 69804 ) stabilized using the UDFC method, pseudo-trajectory ˜ u (t, u upo 1 0 ) (blue), and the analytical value LE 1 (u upo 1 0 ) (green) for t ∈ [0 , 100] in system (2) with parameters set at b = 5 . 7 , c = 18 . 3 , r = 51 . (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) even without the control, the obtained trajectory ˜ u (t, u upo 1 0 ) approximately traces the ”true” trajectory (periodic orbit) u (t, u upo 1 0 ) . But for t > T 1 , without a control, the pseudo-trajectory ˜ u (t, u upo 1 0 ) diverges from u (t, u upo 1 0 ) and visualize a local chaotic attractor A . In general, the closeness of the real trajectory u (t, u 0 ) and the corresponding pseudo-trajectory ˜ u (t, u 0 ) calculated numerically can be guaranteed on a limited short time interval only. The obtained values of the largest finite-time Lyapunov exponent LE 1 (t, u upo 1 0 ) computed along the stabilized UPO u (t, u upo 1 0 ) and the trajectory without stabilization ˜ u (t, u upo 1 0 ) give us the following results. On the initial part of the time interval [0 , T 1 ≈2 τ1 ] , one can indicate the coincidence of these values with a sufficiently high accuracy. After t > T 2 ≈10 the difference in values becomes significant and the corresponding graphs diverge in such a way that the graph corresponding to the unstabilized trajectory is higher than the parts of the graphs corresponding to the UPO and the analytical value largest Lyapunov exponent: LE 1 (u upo 1 0 ) = 1 . 80401 , computed via Floquet multipliers (see Fig. 5 ). Using numerical experiments, we analyze the chaotic dynamics of system (2) and visualize a self-excited attractor for values of parameters b = 5 . 7 , c = 18 . 3 , r = 51 . At the same time, we get formula for the exact Lyapunov dimension of the global attractor for certain region of the parameters (b, c, r) (15) and (16) of system (2) by the analytical way. Thus, we get the following relations dim L A glob = dim L S 0 = 2 . 4347 > dim L A ≈max u ∈C grid dim L ( t k , u ) ≈2 . 0808 ≥dim L u up o 1 ≈2 . 0738 . (39) 5. Conclusion In this paper, we studied the irregular behavior (including chaotic attractor) of the mid-size firm model, assuming the deterministic endogenous mechanism for generating these fluctuations in the economic system. Using an analytical approach, we calculated quantitative characteristics of irregular dynamics, such as the Lyapunov dimension, and demonstrated the complexity and ambiguity of using numerical procedures for calculating these indicators. First, we proved a theorem about the exact formula for the Lyapunov dimension of the global attractor in the model. Similar way could be used for getting the formula for the topological entropy. Second, we identified an UPO for the model and stabilized it using the Pyragas control procedure. Third, we numerically calculated the finite-time Lyapunov dimension along the trajectories of the global attractor, including UPO, thereby providing support for arguments about difficulties of application of the numerical procedures and importance of the obtained exact formula for the Lyapunov dimension of the global attractor. We believe that expanding our knowledge of the role, sources, as well as qualitative and quantitative characteristics of irregular oscillatory dynamics may diminish researchers’ reliance on unrealistically large shocks to explain economic data. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgments This paper was prepared with the support by the Leading Scientific Schools of Russia: project NSh-2624.2020.1 (sections 3, 4). Authors from the St.Petersburg State University acknowledge support from St.Petersburg State University grant Pure ID 75207094 (section 1,2). This work was motivated by research conducted at the Institute for Nonlinear Dynamical Inference at the International Center for Emerging Markets Research ( http://icemr.ru/ institutefornonlineardynamicalinference/ ) and Financial Research Institute of the Ministry of Finance of the Russian Federation, a number of whose employees the authors thank for helpful suggestions and comments. Especially, we thank William A. Barnett, with whom the authors started to collaborate in the direction considered above [62] , for his extremely valuable comments and support. Appendix Let { ϕ t } t≥0 denote a smooth dynamical system with continuous time, and let set A be its compact invariant set. Fundamental matrix Dϕ t (u ) = y 1 (t) , y 2 (t) , y 3 (t) , Dϕ 0 (u ) = Iconsists of linearly independent solutions { y i (t) } 3 i =1 of the linearized system, where I is the unit 3 ×3 matrix, with the following cocycle property: Dϕ t+ s ( u ) = Dϕ t ( ϕ s ( u ) ) Dϕ s ( u ) , ∀ t, s ≥0 , ∀ u ∈ R 3 . (40) Let LE i (·) = t −1 ln σi (·) for t > 0 , where σi (Dϕ t (u )) = σi (t, u ) , i = 1 , 2 , 3 , be the singular values of Dϕ t (u ) (i.e. σi (t, u ) > 0 and σi (t, u ) 2 are the eigenvalues of the symmetric matrix Dϕ t (u ) ∗Dϕ t (u ) with respect to their algebraic multiplicity), ordered so that σ1 (t, u ) ≥σ2 (t, u ) ≥σ3 (t, u ) > 0 for any u ∈ R 3 , t ≥ 0 . Consider a set of finite-time Lyapunov exponents { LE i (Dϕ t (u )) = LE i (t, u ) } 3 i =1 at the point u : LE i (t, u ) = 1 t ln σi (t, u ) , t > 0 , i = 1 , 2 , 3 , (41) ordered by decreasing for all t > 0 . We can introduce the following concepts –the finite-time local Lyapunov dimension (of map ϕ t at point u ): dim L (t, u ) = dim L (ϕ t , u ) , the finite-time Lyapunov dimension (of map ϕ t with respect to set A ): dim L (t, A ) = dim L (ϕ t , A ) , and the Lyapunov dimension (of dynamical system { ϕ t } t≥0 with respect to set A ): dim L A = dim L ({ ϕ t } t≥0 , A ) . Consider the dynamical system { ϕ t } t≥0 , (R 3 , || ·|| ) under the change of coordinates w = h (u ) , where h : R 3 → R 3 is a diffeomorphism. In this case the dynamical system { ϕ t } t≥0 , (R 3 , || ·|| ) is transformed to the dynamical system { ϕ t h } t≥0 , and the compact 7 T.A. Alexeeva, N.V. Kuznetsov and T.N. Mokaev Chaos, Solitons and Fractals 152 (2021) 111365 set A ⊂R 3 invariant with respect to { ϕ t } t≥0 is mapped to the compact set h (A ) ⊂R 3 . Here Dϕ t h (w ) = Dh (ϕ t (u )) Dϕ t (u ) Dh (u ) −1 . (42) Proposition 1. (see, e.g. [79 , 89] ) For any diffeomorphism h : R 3 → R 3 the Lyapunov dimension is invariant with respect to diffeomorphism, i.e. dim L ({ ϕ t } t≥0 , A ) = dim L ({ ϕ t h } t≥0 , h (A )) . (43) The proof of this proposition uses the Horn inequality for (42) and the fact that singular values of Dh (ϕ t (u )) and (Dh (ϕ t (u ))) −1 are uniformly bounded in ton A . Moreover, instead of Dh one can consider any 3 ×3 matrix H(u ) such that all its elements are scalar continuous functions of u and det H(u )  = 0 for all u ∈ A , and get 6 lim t→ + ∞ LE i H(ϕ t (u )) Dϕ t (u ) H(u ) −1 −LE i Dϕ t (u ) = 0 , i = 1 , 2 , 3 , dim L ({ ϕ t } t≥0 , A ) = lim inf t→ + ∞ sup u ∈A d KY { LE i H(ϕ t (u )) Dϕ t (u ) H(u ) −1 } 3 1 . (44) If an equilibrium u eq ≡ϕ(u eq ) ∈ A has simple real eigenvalues, then a nonsingular 3 ×3 matrix Sexists such that the linearization takes the form SD f (u eq ) S −1 = diag λ1 (u eq ) , ···, λ3 (u eq ) , where λj (u eq ) ≥λj+1 (u eq ) , i = 1 , 2 . Then, by the linear change of variables w = h (u ) = Su and the invariance we get lim t→ + ∞ LE i (t, u eq ) = λi (u eq ) and dim L u eq = d KY ({ λi (u eq )) } 3 i =1 . For analytical estimation of the Lyapunov dimension via the eigenvalues of the symmetrized Jacobian matrix we use the generalized Liouville’s relation (see, e.g., [83] , [81, p.68] ) and get, ∀ t > 0 , u ∈ A , the following: j  i =1 LE i ϕ t ( u ) + s LE j+1 ϕ t ( u )  ≤1 t t  0 j  i =1 νi ( ϕ τ( u ) ) +sνj+1 ( ϕ τ( u ) ) dτ ≤sup u ∈A j  i =1 νi ( u ) + sνj+1 ( u ) . (45) From (45) we obtain the upper estimation of the Lyapunov dimension (13) . The Leonov method of analytical estimation of the Lyapunov dimension is based on (44) and (13) . Following [84,90,91] , we consider H(u ) = p(u ) S, where p : R 3 → R 1 is a continuous scalar function, Sis a nonsingular 3 ×3 matrix. Then we compute the Lyapunov dimension by (44) : dim L A = lim inf t→ + ∞ sup u ∈A d KY { LE i p(ϕ t (u )) p(u ) −1 SDϕ t (u ) S −1 } 3 1 , and estimate it by (13) . For that by (41) and (45) we get the estimation:  j i =1 LE i p ϕ t ( u ) p ( u ) −1 SD ϕ t ( u ) S −1  ≤j 1 t ln p ϕ t ( u ) p ( u ) −1 + 1 t  t 0  j i =1 νi SJ ( u ) S −1 dτ. (46) In general, while under the diffeomorphism h (u ) = Su the Lyapunov dimension is invariant and J(u ) → SJ(u ) S −1 , the values νi (SJ(u ) S −1 ) = νi (u, S) are not invariant and, thus, Stogether with p(u ) may be used to simplify their computation (the idea with S was introduced in [90, Eq.(8)] and p(u ) was introduced in [84] ). 6 By the Horn inequality for the matrices D H (ϕ t (u )) = H(ϕ t (u )) Dϕ t (u ) H(u ) −1 and Dϕ t (u ) = H(ϕ t (u )) −1 D H (ϕ t (u )) H(u ) . The scalar multiplier of the type p(ϕ t (u ))(p(u )) −1 can be interpreted as the changes of Riemannian metrics [92] (see, also [81] ). The following theorem is a reformulation of the results from Leonov [91] , 93 ] (see also [79,81] ). Theorem 2. If there exist an integer j ∈ { 1 , 2 } , a real s ∈ [0 , 1] , a differentiable scalar function V : R 3 → R 1 , and a nonsingular 3 ×3 matrix Ssuch that condition (14) , i.e. sup u ∈A ν1 ( u, S ) + ···+ νj ( u, S ) + sνj+1 ( u, S ) + ˙ V ( u ) < 0 , is satisfied, where ˙ V (u ) = ( grad (V )) ∗f(u ) , then dim H A ≤dim L A < j + s. Proof. Let p(u ) = e V (u )(j+ s ) −1 . Then (j + s ) 1 t ln (p(ϕ t (u )) p (u ) −1 )= 1 t (  t 0 ˙ V (ϕ τ(u )) dτ) . Thus by invariance of Aand (45) from (46) we get j  i =1 LE i (SDϕ t (u ) S −1 ) + s LE j+1 (SDϕ t (u ) S −1 ) +(j + s ) 1 t ln p(ϕ t (u )) p(u ) −1 ≤ ≤sup u ∈A j  i =1 νi (u, S) + sνj+1 (u, S) + ˙ V (u ) < 0 . (47) Since lim t→ + ∞ (j + s ) 1 t ln p(ϕ t (u )) p(u ) −1 = 0 for any u ∈ A there exists T > 0 such that j  i =1 LE i SD ϕ t ( u ) S −1 + s LE j+1 SD ϕ t ( u ) S −1 < 0 , ∀ t > T , u ∈ A . (48) Thus, taking into account (10) , dim L A < j + s .  References [1] Scheffer M , et al. Anticipating critical transitions. Science 2012;338:344–8 . [2] Battiston S . Complexity theory and financial regulation economic policy needs interdisciplinary network analysis and behavioral modeling. Science 2016;351(6275):818–19 . [3] Aliber RZ , Kindleberge CP . Manias, panics, and crashes: a history of financial crises. Palgrave Macmillan UK; 2015 . [4] Benhabib J . Chaotic dynamics in economics. In: The new palgrave dictionary of economics. London: Palgrave Macmillan UK; 2016. p. 1–4 . [5] Akhmet M , Akhmetova Z , Fen M .Chaos in economic models with exogenous shocks. J Econ Behav Organ 2014;106:95–108 . [6] Beaudry P , Galizia D , Portier F . Putting the cycle back into business cycle analysis. Am Econ Rev 2020;110(1):1–47 . [7] Lorenz E . Deterministic nonperiodic flow. J Atmos Sci 1963;20(2):130–41 . [8] Ueda Y , Akamatsu N , Hayashi C . Computer simulations and non-periodic oscillations. Trans IEICE Jpn 1973;56A(4):218–55 . [9] Benhabib J , Nishimura K . The Hopf bifurcation and the existence and stability of closed orbits in multi sector models of optimal economic growth. J Econ Theory 1979;21:421–44 . [10] Benhabib J , Nishimura K . Competitive equilibrium cycles. J Econ Theory 1985;35(2):284–306 . [11] Day RH . Irregular growth cycles. Am Econ Rev 1982;72(3):406–14 . [12] Day RH . The emergence of chaos from classical economic growth. Q. J. Econ. 1983;98(2):201–13 . [13] Grandmont JM . On endogenous competitive business cycles. Econometrica 1985;53(5):995–1045 . [14] Boldrin M , Montrucchio L . On the indeterminacy of capital accumulation paths. J Econ Theory 1986;40:26–39 . [15] Boldrin M , Woodford M . Equilibrium models displaying endogenous fluctuations and chaos: a survey. J Monet Econ 1990;25(2):189–222 . [16] Day RH , Shafer W . Ergodic fluctuations in deterministic economic models. J Econ Behav Organ 1987;8(3):339–61 . [17] Baumol W , Benhabib J . Chaos: significance, mechanism, and economic applications. J Econ Perspect 1989;3(1):77–105 . [18] Medio A . Chaotic dynamics: theory and applications to economics. Cambridge University Press; 1992 . [19] Sorger G . On the minimum rate of impatience for complicated optimal growth paths. J Econ Theory 1992;56:160–79 . [20] Mitra T . An exact discount factor restriction for period three cycles in dynamic optimization models. J Econ Theory 1996;69:281–305 . [21] Brock WA , Hommes CH . A rational route to randomness. Econometrica 1997;65:1059–95 . [22] Brock W , Hommes C . Heterogeneous beliefs and routes to chaos in a simple asset pricing model. J Econ Dyn Control 1998;22:1235–74 . 8