Generation mechanisms of strange nonchaotic attractors and multistable dynamics in a class of nonlinear economic systems

In this paper, we study strange nonchaotic attractors (SNAs) and multistable dynamics in a class of nonlinear economic systems. For quasiperiodically forced case, the generation and evolution mechanisms of SNAs are discussed. The fractal, Heagy–Hammel, torus doubling, and intermittency routes to SNAs are identified. The Lyapunov exponent, phase-sensitive function and power spectrum are used to characterize the dynamical and geometrical properties of SNAs. Moreover, when multistable phenomenon occur in the system, the boundaries of the basin of attraction may become intertwined, which leads to the economic unpredictability in the long run.


Introduction
In recent years, chaotic dynamics has made many applications in the field of economics. Benhabib and Day [1] first introduced the chaotic dynamics theory into economics. They showed that sequences of rational choices can be erratic when preferences depend on experience. Puu [2] used Cournot's duopoly theory to study the nonlinear dynamical behavior of two competing firms in the market. Under the assumption of isoelastic demand and constant unit production cost, the model shows repeated periodic and chaotic motions. Chiarella [3] considered the generalized nonlinear supply function in the traditional cobweb model and proved that there exists a route from period-doubling bifurcation to chaos in the locally unstable region. Matsumoto [4] investigated that it is better for the whole economy to fluctuate in chaos than to be stable in the long run.
The unemployment-inflation correlation is one of the most important correlations in economics. Phillips [5] reported his observations of rising wages and unemployment in 1958. As the first person to notice this, he pointed out that there is not only a correlation between inflation and unemployment, but also a negative correlation. In a sense, with the increase in wages, the unemployment rate will decline. In addition, inflation and unemployment can affect financial markets because it raises the level of uncertainty, which means an increase in financial market turbulence [6]. Hence, the dynam-ical behaviors of inflation and unemployment models are a problem worth studying.
In order to further explore the dynamical properties of a system, one usually applies external excitations, and these external excitations will make the system exhibit abundant dynamical phenomena [7][8][9], such as phase locking, attractor crises, and quasiperiodic motion. SNAs usually appear in some systems with quasiperiodically forced excitation [10][11][12][13]. On the one hand, this class of attractors exhibits the dynamical properties of regular attractors, that is, they are nonchaotic in the dynamics, since the maximum Lyapunov exponent is nonpositive. On the other hand, SNAs also show the geometric characteristics of chaotic attractors, that is, they have geometric fractal structure.
Over the years, SNAs have been widely studied experimentally and numerically in different dynamical systems. Pikovky and Feudal et al. [14] first proposed the method of rational number approximate frequency and phase sensitivity function to characterize the strangeness of SNAs. Ding et al. [15] studied a class of quasiperiodic excitation systems, in which the parameter region of SNAs has fractal structure, which lies between two critical curves, one of which marks the transition from quasiperiodic attractors to SNAs. The other marks the transition from SNAs to chaotic attractors. Romeiras et al. [16] identified SNAs with twofrequency quasiperiodic forcing in the damped pendulum equation and experimentally observed that SNAs have power spectral characteristics with different frequencies. Zhang et al. [17] proved that SNAs exist in the FHN neuron model under weak noise disturbance. Prasad et al. [18] found that SNAs can be generated by the collision of an invariant curve with itself. Venkatesan et al. [19] showed very many types of routes into chaos through SNAs in a straightforward quasiperiodic forcing cubic map and distinguished between the two classes of attractors by the phase diagram and finite time Lyapunov exponents. In particular, Linder et al. [20] used data collected by the Kepler space telescope to find strange nonchaotic stars, which further showed that strange nonchaotic phenomena are real in nature.
As a special type of attractor, many experts have devoted to the investigation of irregular dynamical transitions and mechanisms of SNAs. Several mechanisms and routes for forming SNAs are described in the literature, such as fractal route [21], torus collision route [22], intermittency route [23,24], bubbling route [25], and blowout bifurcation route [26]. The above litera-ture provides further reference for other routes [27][28][29][30] to form SNAs. In additional, a number of mathematically rigorous results [31][32][33][34] on the topic of SNAs have been reported in recent years.
We are interested in the dynamical behaviors of the unemployment-inflation system [35] under quasiperiodically forced excitations. The paper is organized as follows. In Sect. 2, we briefly describe a nonlinear economic system and introduce the phase sensitivity. Then, the generation mechanisms of SNAs are discussed in Sect. 3. The nonchaotic and strange properties of SNAs are analyzed via phase-sensitive function and power spectrum in Sect. 4. In Sec. 5, we uncover a dynamical phenomenon in which a 1T quasiperiodic attractor and 3T quasiperiodic attractor coexist in the system and obtain the basin of attraction of these coexisting attractors. The main results are summarized in Sect. 6.

Nonlinear economic system
The unemployment and inflation rate are important factors that affect social stability and economic development. The dynamical analysis of unemployment and inflation models can deeply understand the current economic state and propose reasonable policies [36][37][38]. The unemployment rate refers to the proportion of surplus labor force in the whole labor force. Under the condition of high unemployment rate, the economy will have a downward trend. Inflation rate refers to the increase in average price in a continuous period, and consumer price index is generally used to describe inflation. There are many economic theories about inflation and unemployment. Among them, the more authoritative and effective theory is Phillips curve, which establishes a functional relationship between inflation rate and unemployment rate. Taking into account the influence of other relevant factors, we apply a small perturbation to the system of equations for the Phillips curve, and the form of the system is obtained as follows.
where a, b, c, m, ω, β 1 , and β 2 are the system parameters. In order to better understand the dynamical behav-iors of the nonlinear economic system, we take ε, b, and c as the control variables, and fix a = 1, m = 2, ω = ( √ 5 − 1)/2, β 1 = −2.5, and β 2 = 20 in what follows. Chaos in dynamical systems is characterized by its sensitive dependence on initial conditions, which is often called "butterfly effect." The sensitive dependence on small perturbations or initial conditions can be described by the Lyapunov exponent. In addition, the Lyapunov exponent can also show the average stretching or compression rate of each point in the moving orbit for a long time, which can describe the overall properties of the dynamics. Take a two-dimensional map as an example: According to the definition of Lyapunov exponent, we have If the Lyapunov exponent is nonpositive, then the system is nonchaotic, which means that there is no sensitive dependence on the initial conditions. The characterization of strange properties (geometric fractal structure) of SNAs is a difficult problem. Pikovsky and Feudel [14] proposed that strange and nonstrange attractors can be distinguished by tht phase sensitivity. From the map (2), we get the recurrence relation: Equation (4) can be further written as According to Eq. (2), the equation ∂ x m /∂θ m = ∂ x m /∂θ k = ∂ x m /∂θ 0 holds for any positive integer m and k. Therefore, the subscript of θ can be omitted and ∂ x n /∂θ can be called "derivative of variable with respect to phase." Hence, Eq. (5) can be rewritten as: where and R 0 = 1. The SNAs are nonchaotic, that is, the Lyapunov exponent λ is nonpositive. For sufficiently large n, the value of R n is very small, so as long as n = N is large enough, the derivative value can be approximated by Therefore, the strange properties of the attractor can be verified in terms of whether S N is bounded or not. To better characterize the strangeness of SNAs using this method, consider If the value of γ N increases with the number of iterations N , this means that |S N | tends to infinity, namely, the derivative of the state variable with respect to the phase is not a bounded value. Therefore, with the help of Eq. (9), the strange properties of the attractor can be verified. With the increase in the number of iterations, the growth rate of S N can represent the strange degree of the attractor, which gives a quantitative index of strangeness. Next, take m initial points (x 0 , θ 0 ) at random and compute the minimum value of γ N (x, θ) can more accurately determine whether the attractor is nonsmooth. In the process of calculation, if the attractor is nonsmooth, τ N will increase as N increases.
Assuming that this function increases by a power law as N → ∞ where the exponent μ is a quantitative characteristic of the strangeness of the attractor, which is called "phase sensitivity exponent." If μ ≈ 0, then the attractor is regular. while μ = 0, the value of the derivative of the state variable with respect to the phase increases with the increase in the number of iterations, and S N tends to infinity, so that the attractor is strange (nonsmooth).

Generation mechanisms of SNAs
In this section, we describe four routes for generating SNAs from quasiperiodic attractors and further discuss the generation mechanism of these SNAs.

The fractal route
The attractors get gradually wrinkled without any interaction with nearby periodic orbits in the fractal routes. The generation mechanism of such SNAs was described by Kaneko as early as 1984, but he did not refer to the emerging new invariant sets as SNAs. In contrast to other mechanisms for generating SNAs, this mechanism is difficult to relate to a precise bifurcation phenomenon. The fractal routes only change the structure of attractors gradually, and the number of attractors does not change obviously. Let us fix the parameter c = 0.75, ε = 0.05 and increase the value of b. For b = 1.35, there are two smooth branches in the phase diagram, which indicates that the system (1) is in a doubled torus state denoted as 2T (see Fig. 1a). As b is decreased further to 1.4, the 2T quasiperiodic attractor gets increasingly wrinkled and irregular, but the attractor remain continuous as shown in Fig. 1b. When b = 1.409, the attractor becomes extremely wrinkled and appears some discontinuous regions (see Fig. 1c). At such values, the attractor is strange, and the maximum Lyapunov exponent of the system is negative (λ max = − 0.0441 (see Fig. 2a)). Thus, the 2T quasiperiodic attractor transforms into an SNA. For b = 1.412, the SNA eventually evolved into a chaotic attractor (Fig. 1d) with a fractal structure and positive maximum Lyapunov exponen ( λ max = 0.0149 (see Fig. 2b)). Generally speaking, in the process of the attractor evolution, an SNA appears only after a smooth quasiperiodic attractor produces the wrinklelike geometry. Therefore, the wrinkling phenomenon can be regarded as a "prelude" to the generation of SNAs from fractal routes.

The Heagy-Hammel route
As the change of control parameters, two stable curve and unstable invariant curves collide, the stable invariant curves will lose stability, from a quasiperiodic attractor to an SNA, this route is known as Heagy-Hammel route. This mechanism may be closely related to homoclinic bifurcation.
To better illustrate the transition in the system (1), let c = 0.4, ε = 0.02 and b is to be taken as control parameter. For b = 1.15, there are four smooth invariant curves in the phase diagram, namely, a 4T quasiperiodic attractor (see Fig. 3a). As b is increased to 1.163, the 4T quasiperiodic attractor appears wrinkled and the distance between adjacent invariant curves starts to shrink (see Fig. 3b). As b increases further to 1.16575, the two adjacent stable invariant curves (blue lines) collide with the unstable invariant curve (red dashed lines), which causes the stable invariant curves become locally nonsmooth (see Fig. 3c). As b continues to increase to 1.166, the 4T quasiperiodic attractor evolves into an SNA (see Fig. 3d). Under this parameters, the maximum Lyapunov exponen of the system (1) is − 0.0177, so it can be verified that the attractor is nonchaotic (see Fig. 4). We find that the collision between the stable 4T quasiperiodic attractor and the unstable 2T quasiperiodic attractor leads to the creation of the SNA, and the attractor is discontinuous along the z-axis, which also further illustrates strangeness of the attractor.

The torus-doubling route
The torus-doubling route is that the torus-doubling bifurcation is interrupted by a subharmonic bifurcation, and the torus attractor loses stability and evolves into an SNA. The torus-doubling bifurcation and perioddoubling bifurcation have similarities and differences. The differences are as follows: (1) the dimension of attractor is different. (2) When period-doubling bifurcation occurs in a system, the Lyapunov exponent is equal to zero in one direction. When a torus-doubling bifurcation occurs, the Lyapunov exponent is equal to zero in two directions. (3) Period-doubling bifurcation is a typical route to chaos, and torus-doubling bifurcation may lead first to an SNA and then to a chaotic attractor, which is a typical route to SNAs.
In order to illustrate the transformation of torusdoubling route to SNAs, we can draw the phase diagrams in the (x, y) plane for c = 0.75, ε = 0.001. For b = 1.07, a 1T quasiperiodic attractor occurs in the (x, y) plane (see Fig. 5a). As the value of b is decreased b = 1.35, the 1T quasiperiodic attractor evolves into a 2T quasiperiodic attractor, which is gen-  Fig. 5b). As the value of b is decreased 1.4, the 2T quasiperiodic attractor again undergoes the torus doubling bifurcation and the corresponding period-4 quasiperiodic orbit is denoted as 4T (see Fig. 5c). In the generic case, the torus-doubling phenomena can continue to occur until a critical point is reached, beyond which the system will exhibit chaotic motion. However, in the present case, the torus-doubling phenomenon no longer occurs, instead the 8T quasiperiodic attractor becomes wrinkled and loses smoothness, which is shown in Fig. 5e. This is because the 8T quasiperiodic attractor collides with its unstable parent torus. For the attractor exhibited in Fig. 5e, the corresponding to the maximum Lyapunov exponent remains negative (λ max = − 0.0041) (see Fig. 6a). Therefore, the system is in the strange nonchaotic state for b = 1.43305. It is worth noting that the interval (1.43, 1.43305) is selected as the candidate interval for the existence of SNAs. In the next section, we will further verify that this interval is corresponding to SNAs through phase sensitivity function and power spectrum. Finally, the SNA evolves into a chaotic attractor with b = 1.44 (see Fig. 5f).

The type-I intermittency route
In type-I intermittency routes, the torus is eventually replaced by a strange nonchaotic attractor through the saddle-node bifurcation. When the Floquet multiplier of unperturbed system passes through the unit circle, the attractor loses stability and the phenomenon is called "type-I intermittency." Saddle-node bifurcation is a necessary but not sufficient condition for type-I intermittency occurrence. Another necessary condition for type-I intermittency occurrence is that the orbit must repeatedly enter the neighborhood of the original periodic orbit.
In order to describe SNAs in terms of torus intermittency, let us fix the parameters c = 0.75, ε = 0.00314 and b varies from 1.42 to 1.43. For b = 1.42, a 4T quasiperiodic attractor with four smooth branches appears in the (θ, x) plane (see Fig. 7a). As b is increased to 1.43, the 4T quasiperiodic attractor suddenly appears many disordered points near the orbit, which are the characteristics of type-I intermittency (see Fig. 7c and d). From the enlarged figure in Fig. 7d, it can be observed that the attractor has the geo-  metric fractal structure, and the attractor is no longer smooth. Therefore, the attractor is strange. The nonchaotic property is represented by the maximum Lyapunov exponent λ max = − 0.0358 (see Fig. 8).

Determining the strange properties of SNAs
In Sect. 2, we introduce the phase sensitivity function to verify the strange properties of SNAs. Here we take the SNA generated by the torus-doubling route as an example (Fig. 5e) to describe the strange properties of this attractor. The phase sensitivity exponent μ and the maximum derivative value τ N of the state variable with respect to the phase z are calculated. For b = 1.4, the phase sensitivity exponent μ = 0 and the value of τ N tends to a small bounded value as the number of iterations increases, which means the attractor is regular (see Fig. 9a). In sharp contrast, for b = 1.43305, the phase sensitivity exponent μ = 2.3561 and the value of τ N tends to infinity as the number of iterations increases; namely, the attractor is strange; see Fig. 9b.
The power spectrum (Fourier amplitude spectrum) corresponding to periodic attractors or quasiperiodic attractors is discrete, and the discrete power spectrum has some δ-peaks. The power spectrum corresponding to the chaotic attractor is continuous and has no δ-peaks. However, SNAs usually appear in the transition region from quasiperiodic attractors to chaotic attractors, the power spectrum corresponding to SNAs is a special spectrum between discrete and continuous. This particular spectrum is called singular continuous spectrum, which exhibits a combination of continuous and discrete components, and has many δ-peaks [39]. We take the Fourier transform of the time series {x n } at frequency ω as {x n } at frequency ω as x n e i2π nω .
Hence, the power spectrum of the attractor is defined as [40] Therefore, we can use singular continuous spectrum to examine strange properties of the attractor. For the SNA shown in Fig. 5e, The power spectrum is continuous with many δ-peaks, indicating that the attractor is strange (see Fig. 10).

Multistable analysis
In complex dynamical systems coexistence of attractors is called multistability or multistable dynamics. The final state of multistable systems is closely related to the choice of initial conditions. Small changes in initial conditions may lead to changes in attractor types. In this section we extend our analysis to determining basins of attraction in parameter space. For the system parameters c = 0.75, ε = 0.01 and b = 0.75, a 1T quasiperiodic attractor and 3T quasiperiodic attractor coexist in the system under different initial values (2.5, 3.5, 0) and (2, 2, 0), as shown in Fig. 11. The attracting domain corresponding to the 1T quasiperiodic attractor is green region and the attracting domain corresponding to white 3T quasiperiodic attractor is red region. The blue region is infeasible region (see Fig. 12). The attracting domain of the 3T quasiperiodic attractor is nested within the 1T quasiperiodic attractor, and the coexisting basins of attraction are highly intertwined implying that small uncertainties in the specification of the parameters result in qualitatively different types of behavior.

Conclusions
In this work, the dynamics of a class of nonlinear economic systems under quasiperiodic excitation is considered. We identify four types of transitions to strange nonchaotic attractors as the system parameters are varied: fractal route, Heagy-Hammel route, torus-doubling route, and intermittent route. The methods used to numerical investigations are maximal Lyapunov exponent, Fourier spectrum, and phase sensitivity function. In additional, a novel dynamical behavior that the 1T quasiperiodic attractor coexists with 3T quasiperiodic attractors is uncovered. In such cases, basin of attraction is also obtained to better understand the long-term unpredictability of the economic system. The results of this work offer ideas for the study of strange nonchaotic dynamics in other systems as well as provide support for the economic theory and policy research.