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

Cyclicity 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).


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 socio-economic model [1][2][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 boom-bust 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 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 offen, 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 generated a lot of discussion regarding the views on the sources of irregular fluctuations (see, e.g. [4][5][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 and Ueda [7,8]. Many famous economists (see, e.g. [4,[9][10][11][12][13][14][15][16][17][18][19][20][21][22][23][24][25][26][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][29][30][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 fine-tuning 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][33][34][35][36][37][38][39][40]. The reviewed literature shows 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][42][43][44][45][46][47][48][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][53][54][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, two-dimensional 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][59][60][61][62][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 low-dimensional 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 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.

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 Sharovalov model proposed in [65] which describes the behavior of a mid-size firm where coefficients α, β, σ, δ, µ, γ at variables (x, y, z) 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: where X, Y , and Z are 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 using the following coordinate transformation 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 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] using coordinate substitutions 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 < 2c (Fig. 1 , and η is a chosen parameter. 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][69][70][71][72][73][74][75][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)  ± 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 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 exhibits irregular behavior, at which a chaotic attractor can be reveal. As an example, Shapovalov et al. [65,77], and Gurina and Dorofeev [78] show that system (1) exhibits chaotic behavior for some values of parameters.
Localization of global attractor and furthest calculation of the limit values of the finite-time Lyapunov exponents and the finite-time Lyapunov dimension along various trjectories of this attractor are nontrivial tasks. While trivial attractors (stable equilibrium) can be easily found analytically or numerically, the search for periodic and 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 of 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 deduсed analytical formula for the Lyapunov dimension.

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 summarised in Appendix.
Rewrite system (2) in a formu where f is a continuously differentiable vector-function. Let u(t, u 0 ) be any solution of (8) such that u(0, u 0 ) = u 0 ∈ B exists for t ∈ [0, ∞), it is unique and stays in the absorbing set B. 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 (U, || · ||): {ϕ t } t≥0 , (U ⊆ 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, I is 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 3 where j(t, u) = max{m : The finite-time Lyapunov dimension is defined as: 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 4 sup u∈A 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 < lim inf Thus, the best estimation (11) takes the form [79] dim 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 exact Lyapunov dimension in [82]. A conjecture on the Lyapunov dimension of self-excited 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 Lyapunov dimension is hardly possible. However, they can be estimated by the eigenvalues of the symmetrized Jacobian matrix [80,83]. The Kaplan-Yorke formula with respect to the ordered set of eigenvalues ν i (J(u)) = ν i (u), ν 1 (u) ≥ ν 2 (u) ≥ ν 3 (u), 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: Generally speaking, one cannot get the same values of ) 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. Another approach is based on the Leonov method of analytical estimation of the Lyapunov dimension 5 whereV 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.

Main result
Using an effective analytical approach, proposed by Leonov [79,84], which is based on a combination of the Douady-Oesterlé approach with the direct Lyapunov method 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 then Proof. Consider system (2) with the Jacobian matrix under the conditions (15) and (16). We apply the transformation (3)  . Then the symmetrized Jacobian matrix of this system 1 2 (SJS −1 + (SJS −1 ) * ) 6 has the following eigenvalues The inequalities imply λ 1 ≥ λ 2 ≥ λ 3 . From (21) following [84] we get the ratio where s ∈ [0, 1] is a real number. Using the famous inequality √ k + l ≤ √ k + l 2 √ k , ∀k > 0, l ≥ 0, we obtain an estimate We introduce the function V (x, y, z) = θ(x,y,z) P and Q i (i = 0, 2) are some positive real parameters. Then 6 Symbol * denotes the transposition of matrix.
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, 2) of the function θ(x, y, z) such that Substituting W (x, y, z) andθ in (27), we get where , Since RHS of the second inequality in (31) is positive, we obtain 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 ) = 3P +2b+2 8P (b−2) > 0, . Inequalities (34) mean that a positive Q 2 exists such that where . Since the denominator of fraction (35) is positive, we obtain required condition (17) This completes the proof.
We obtain a formula of the exact Lyapunov dimension of the global attractor for certain region of the parameters (b, c) (Fig. 2) of system (2). 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 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 . We show the grid of points  τ k (k = 1, . . . , N ), N = 1000 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 Note that if conditions on dissipativiness of system (2) are not satisfied 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 dimension (e.g. if there are trajectories with initial conditions in the cuboid, which tend to infinity). 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 attractors. 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 equilibrium 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 finite-time 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 into a chaotic attractor. The "skeleton" of a chaotic attractor for this system comprises embedded UPOs. 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.
In general, the closeness of the real trajectory u(t, u 0 ) and the corresponding pseudo-trajectorỹ 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ũ(t, u upo 1 0 ) gives 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 r = 51, b = 5.7, c = 18.3. At the same time we get formula of the exact Lyapunov dimension of the global attractor for certain region of the parameters (b, c, r) (15) and (16) of system (2). Thus, we get the following relations

Conclusion
In this paper, we studied the irregular behavior (chaotic attractor, unstable limit cycles) of the mid-size firm model, assuming the deterministic endogenous mechanism for generating these fluctuations. Using an analytical approach, we calculated quantitative characteristics of irregular dynamics, such as the Lyapunov dimension and topological entropy, and demonstrated the complexity and ambiguity of using numerical procedures for calculating these indicators. We have obtained a number of new results. First, we proved a theorem about the exact formula for the Lyapunov dimension of the global attractor in the model. Similar way used for getting the formula of 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.