# Stability and Multistability of a Bounded Rational Mixed Duopoly Model with Homogeneous Product

Authors: Wei Zhou; Na Zhao; Tong Chu; Ying-Xiang Chang
Journal: Complexity
 (2020)

---

## Abstract

In this paper, a mixed duopoly dynamic model with bounded rationality is built, where a public-private joint venture and a private enterprise produce homogeneous products and compete in the same market. The purpose of this research is to study the stability and the multistability of the established model. The local stability of all the equilibrium points is discussed by using Jury condition, and the stability region of the Nash equilibrium point has been given. A special fractal structure called “hub of periodicity” has been found in the two-parameter space by numerical simulation. In addition, the phenomena of multistability (also called coexistence of multiple attractors) are also studied using basins of attraction and 1-D bifurcation diagrams with adiabatic initial conditions. We find that there are two different coexistences of multiple attractors. And, the fractal structure of the attracting basin is also analyzed, and the formation mechanisms of “holes” and “contact” bifurcation have been revealed. At last, the long-term profits of the enterprises are studied. We find that some enterprises can even make more profits under a chaotic circumstance.


---

## Body

## 1. Introduction

On account of the complexity of oligopoly market setting, many authors are committed to studying the structure of oligopoly market. Cournot is the first person who proposed a duopoly game model describing production competition in the oligopoly market in 1838; the built model by him is also called Cournot model by followers. But, in the early stage of the research of the classic Cournot model, the competitors are always assumed to have naïve expectations, where the players just guess the decision of the other party and will retain the current decision in the following period. Subsequently, great quantity literatures have discussed the stability and dynamical behavior of the Cournot model based on naïve expectation, such as literatures [1–6]. Under the assumption of naïve expectation, the players have complete market information and can accurately predict the response of their rivals, so such models are static in nature. However, the competitive market in which the enterprises are located is a complex system; the enterprises are unlikely to fully grasp market information due to the limitation of cognitive ability. So, the enterprises in reality always show bounded rationality when making decision. After that, a lot of researchers began to study the duopoly model with bounded rationality, in which the participants guess their opponent’s production decision in the next period based solely on their own marginal profit in the current period, as shown in literatures [7–11]. Thereafter, the complexity of competition among limited rational firms has attracted more and more attention of researchers. Some scholars have set out to study the stability of the equilibrium point and the stability condition, which is the set of parameters that they should be satisfied to ensure the stability of the equilibrium point. And then, some other scholars began to study the complex dynamical behaviors of the built model by using different methods. Agliari et al. [12] studied the changes in attractors and their basins of attraction after global bifurcation. Cavalli and Naimzada [13] analyzed the complex dynamical behaviors of the built nonlinear game model when the equilibrium point loses its stability, and the coexistence of multiple attractors is also found in their research. Bischi et al. [14] analyzed the global bifurcation of a nonlinear duopoly model with adaptive expectations using the critical curve, and the equilibrium selection problem arising in this model was also discussed by them. Zhou et al. [15] and Zhang et al. [16] studied the complex dynamical behaviors of the two-stage Cournot duopoly model with R&D spillover effect.However, enterprises are always assumed to be private enterprises in the above discussion. The private enterprises only focus on how to maximize their own profits without considering the social responsibilities. In fact, in the actual market environment, the enterprises need to consider not only their own profits but also their social responsibilities. Therefore, the analysis of mixed duopoly model has attracted the attention of more and more researchers in recent years. The authors studied that in the mixed duopoly market composed of public-private joint venture and private enterprise, partial privatization, or partial nationalization may be the best result in Matsumura [17]. The equilibrium of simultaneous decision-making in a mixed model composed of a semipublic enterprise and several private enterprises is analyzed in the model proposed by Bárcena-Ruiz and Garzón [18]. The welfare of partial privatization of state-owned companies in the differentiated mixed market model was discussed in Fujiwara [19]. It was shown by Fraja and Delbono [20] that the state-owned enterprises may obtain more benefits by maximizing their profits than by maximizing their benefits. Farrell and Shapiro [21] pointed out that, when costs are asymmetric, the merger of private enterprises will increase social welfare. However, most of the literatures on mixed oligopoly are analyzed based on the static game, and the dynamical phenomena of their models are not explained in detail.In this paper, the dynamical behavior of the mixed duopoly model composed of public-private joint venture and private enterprise is analyzed, in which the public-private joint venture takes the payoff of government into its objective function. We assume that these two enterprises produce homogeneous products. Next, we analyze the stability of the discrete time dynamic model and explain the global dynamic characteristics of the dynamic system through numerical simulation. In addition, the difference between the equilibrium profit and long-term average profit of these two enterprises is compared through the profit curve. We aim to discuss the influence of various parameters on the payoff of government and to analyze the phenomenon of coexistence of attractors and other dynamical behaviors of the mixed duopoly model in this paper.The framework of this paper is mainly composed of the following parts: in the second part, we establish a mixed Cournot duopoly model composed of a public-private joint venture and a purely private enterprise based on the limited rationality and the gradient adjustment mechanism. In the third part, the stability conditions of the equilibrium points are mainly discussed. And, in the fourth part, the complex dynamics of the built system is studied by means of numerical simulation. Furthermore, the long-term average profits of the duopoly firms are analyzed in the fifth part of this research. At last, the conclusion is given in Section6.


## 2. The Model

In this section, a mixed duopoly model with two enterprises, which are labeled byi (i = 1, 2) for convenience, is considered. In the built model, we assume that enterprise 1 is a public-private joint venture and enterprise 2 is a purely private enterprise. Assuming that the two enterprises produce homogeneous product, the prices and outputs of the produced goods of enterprise i are given by pi and qi, i=1,2 respectively. According to Dixit [22] and Singh and Vives [23], the inverse demand function of these two enterprises can be assumed as the form, which is given as(1)p=a−bQ,where Q=q1+q2 indicates the total output, a>0 denotes the maximum market price of the produced goods, and b>0 represents the reduction of the price when the quantity of the product increases by one unit. Assuming that the cost function of enterprise i has the form as Ciqi=ciqi,i=1,2, where ci>0 represents the marginal cost of enterprise i and a>ci always holds. Due to the difference of production technology between enterprise 1 and enterprise 2, we assume that c1≠c2. Thus, the profit function of enterprise i can be expressed as(2)πi=a−bq1+q2−ciqi,i=1,2.Because of the different nature of the two enterprises, they pursue different goals. For enterprise 1, it not only pursues the maximal profit but also considers corporate social responsibility, while enterprise 2 only pursues the maximization of its own profit. In this paper, we assume that enterprise 1 considers the payoff of government into its objective function. The payoff of governmentUG is a weighted average of the social welfare and the consumer surplus, where the social welfare W is the sum of consumer surplus and profits of these two enterprises, that is,(3)W=∫0QpQdQ−pQ+π1+π2=a−c1q1+a−c2q2−b2q1+q22.And, the payoff of governmentUG is given as(4)UG=W+β∫0QpQdQ−pQ=a−c1q1+a−c2q2−b21−βq1+q22,where β≥0 represents the weight of consumer surplus in the payoff of government UG. The government wants to make social welfare maximization when β=0, while if we choose β>0, that means the government is more concerned with consumer surplus than the corporate profits.According to the above hypotheses, the objective function of enterprise 1 and enterprise 2 can be written as(5)U1=αUG+1−απ1=αa−c1q1+a−c2q2−b21−βq1+q22+1−αa−bq1+q2−c1q1,U2=π2=a−bq1+q2−c2q2,where α∈0,1 is the proportion of state-owned shares in enterprise 1.The first-order conditions ofU1 and U2 can be written as follows:(6)∂U1∂q1=a−c1−2−α1+βbq1−1−αβbq2,∂U2∂q2=a−c2−bq1−2bq2.In reality, the enterprises are always facing a complex market environment, and they are constrained by factors such as computing ability and incomplete information. And thus, they cannot make decisions with full rationality. On the contrary, they usually make decisions under the premise of “bounded rationality.” Under the assumption of bounded rationality, each enterprise can only adjust its output in the next periodt+1 according to its own output qit and the marginal objective ∂Ui/∂qi,i=1,2 in the current period. That is, enterprise i will increase its output in period t+1, if the marginal objective of itself is positive in period t. Conversely, enterprise i will reduce its output in period t+1, if its marginal objective is negative in period t. According to Bischi et al. [24] and Elsadany and Awad [25], the adjustment mechanism is given as(7)q1t+1=q1t+V1q1t∂U1∂q1,q2t+1=q2t+V2q2t∂U2∂q2,where Viqi is a positive function which gives the extent of production variation of the ith enterprise following a given signal ∂Ui/∂qi,i=1,2. An adjustment mechanism of (7) has been proposed by some researchers [26] with a constant Viqi=vi. Instead, we assume that the function Viqi is a linear function here, i.e., Viqi=viqi, which is equivalent to the assumption that the relative production change is proportional to the estimated marginal objective [27], that is,(8)qit+1−qitqit=vi∂Ui∂qi,i=1,2,where vi>0 represents the speed of adjustment of enterprise i. By substituting the first-order conditions of U1 and U2 into (7), we can establish the following dynamic model:(9)T:q1t+1=q1t+v1q1ta−c1−2−α1+βbq1t−1−αβbq2t,q2t+1=q2t+v2q2ta−c2−bq1t−2bq2t.


## 3. Stability Analysis

In this section, we will analyze the stability conditions of equilibrium points in the given model with varying parameters. To simplify the calculation, we definea−c1=A1>0, a−c2=A2>0, 2−α1+β=B1, and 1−αβ=B2; thus the system (9) can be simplified to(10)T:q1t+1=q1t+v1q1tA1−B1bq1t−B2bq2t,q2t+1=q2t+v2q2tA2−bq1t−2bq2t.In order to gain the equilibrium point, we letqit+1=qit in (10). By calculation, we can find four equilibrium points of system (10):(11)E0=0,0,E1=0,A22b,E2=A1B1b,0,E3=−C1C3b,C2C3b,where C1=2A1−A2B2, C2=A1−A2B1, and C3=B2−2B1.Because the equilibrium pointsE0, E1, and E2 are all on the coordinate axis, we call them boundary equilibrium points, and E3 is called Nash equilibrium point or interior equilibrium point. Obviously, the equilibrium point E3 exists if and only if B2−2B1≠0. The non-negativity of equilibrium point E2 is satisfied when B1>0. And, for the other two boundary equilibrium points, the non-negativity is always met when the selected parameter satisfies the economic meaning. So, based on the above analysis, we can draw a conclusion that all the boundary equilibrium points can be guaranteed to be non-negative if and only if B1>0. In order to make the equilibrium point E3 non-negative, the above defined parameters must meet one of the following conditions:(12)S1=A1>0,A2>0,B2≠2B1,B1>0,C3>0,C2>0,C1<0or(13)S2=A1>0,A2>0,B2≠2B1,B1>0,C3<0,C2<0,C1>0.


### 3.1. Stability of the Boundary Equilibria

The stability condition of equilibrium points of model (10) is generally discussed by solving the eigenvalues of Jacobian matrix. The Jacobian matrix of model (10) at any point q1,q2 is given as(14)Jq1,q2=1+v1A1−2B1bq1−B2bq2−B2bv1q1−bv2q21+v2A2−bq1−4bq2.Through analyzing the Jacobian matrix evaluated at the boundary equilibria, we can get the following propositions.Proposition 1.
E0 is an unstable node.Proof.
By calculation, we can get the Jacobian matrix atE0 as(15)J0,0=1+A1v1001+A2v2.
And, the eigenvalues ofJ0,0 can be obtained as λ1=1+A1v1 and λ2=1+A2v2 by solving the equation λJ−E=0. Because the conditions A1,A2>0 and v1,v2>0 are always satisfied, then both of the eigenvalues meet λ1,2>1. So, E0 is an unstable node.Proposition 2.
When the parameters meet the conditionsA1,A2,B1,B2,C1,C2,C3∈S2 and v2<2/A2, the boundary equilibrium point E1 is a saddle point. When v2<2/A2, the boundary equilibrium point E1 is an unstable node.Proof.
We can get the Jacobian matrix at the equilibrium pointE1 by calculation, which is given as(16)J0,A22b=1+A1v1−12A2B2v10−12A2v21−A2v2.
At the eigenvalues of Jacobian matrix evaluated at equilibrium point,E1 are given by λ1=1+1/22A1−A2B2v1 and λ2=1−A2v2, respectively. Since A1,A2,B1,B2,C1,C2,C3∈S2, it can be found that the eigenvalues meet λ1>1 and λ2<1 when v2<2/A2. And, hence the boundary equilibrium point E1 is a saddle point, when v2<2/A2, the eigenvalues satisfy λ1>1 and λ2>1, so the equilibrium point E1 is an unstable node.
Because the boundary equilibrium pointE2 is located on the coordinate axis as well, the same method can be used to analyze its stability and similar conclusions can be obtained. So, regarding the stability analysis of E2, we will not go into details here.


### 3.2. Stability of the Nash Equilibrium

For the stability of the Nash equilibrium point, we can derive the Proposition3 by solving the eigenvalues of the Jacobian matrix.Proposition 3.
If the system parameters belong to the setA1,A2,B1,B2,C1,C2,C3∈S2, the following conditions are always satisfied:(17)C24−C1v1v2>4C3+2B1C1v1,C22−C1v1v2<B1C1v1,where the Nash equilibrium point E3 will be locally stable.Proof.
First of all, if we calculate Jacobian matrix at Nash equilibrium pointE3, then we can get(18)JE3=1+C1B1C3v1C1B2C3v1−C2C3v21−2C2C3v2.
The characteristic polynomial of the Jacobian matrix evaluated at the Nash equilibrium pointE3 can be written as(19)Fλ=λ2−TrE3λ+DetE3,where TrE3 and DetE3 are the trace and determinant of the Jacobian matrix evaluated at E3, respectively. And, the specific mathematical expression of them can be obtained through a calculation, which are given as(20)Tr=2+C1B1v1−2C2v2C3,Det=C3+C1B1v1−2C2v2+C1C2v1v2C3.
The equilibrium point will be stable when the Jury condition [28] is met, i.e.,(i)
1+Tr+Det=4C3+2B1C1v1−4C2v2+C1C2v1v2/C3>0(ii)
1−Tr+Det=C1C2v1v2/C3>0(iii)
1−Det=2C2v2−B1C1v1−C1C2v1v2/C3>0
Obviously, when the parameters belong to the setA1,A2,B1,B2,C1,C2,C3∈S2, the condition ii always hold. And thus, the Nash equilibrium point will be stable if the other two conditions are satisfied. By a complicated calculation, the other two conditions can be rewritten as the following inequality system, which is given by(21)C24−C1v1v2>4C3+2B1C1v1,C22−C1v1v2<B1C1v1.
The above inequalities are so complicated that we have to discuss them under different situations in order to get a clearer stability condition.(i)
Case I: if the parameters meet the conditionsC1v1>4 and 4C3+2B1C1v1<0, then we can get v2>0>4C3+2B1C1v1/4−C1v1C2. While if the parameters meet C1v1>4 and 4C3+2B1C1v1>0, then we can get v2>4C3+2B1C1v1/4−C1v1C2>0.(ii)
Case II: if the parameters meet the following two conditions as2<C1v1<4 and 4C3+2B1C1v1/4−C1v1C2>B1C1v1/2−C1v1C2, then we can get 0<v2<B1C1v1/2−C1v1C2. While if the parameters meet the conditions as 2<C1v1<4 and 4C3+2B1C1v1/4−C1v1C2<B1C1v1/2−C1v1C2, then we can get 0<v2<4C3+2B1C1v1/4−C1v1C2.(iii)
Case III: if the conditionC1v1<2 holds, then we can get v2>0>B1C1v1/2−C1v1C2.
In particular, when we let one of the eigenvaluesλi=−1,i=1ori=2, that is, the following equation holds: 1+Tr+Det=4C3+2B1C1v1−4C2v2+C1C2v1v2/C3=0, then we can obtain a critical value of v2, i.e., v2∗=22C3+B1C1v1/4−C1v1C2. The Nash equilibrium E3 will lose its stability due to the occurrence of flip bifurcation when v2>v2∗. In other words, if the value of speed of adjustment v2 is greater than v2∗, the system will enter the chaotic state at a double speed. In addition, if we suppose 1−Det=0, then another critical value of v2 can be obtained, i.e., v2∗∗=B1C1v1/2−C1v1C2. And, the Nash equilibrium E3 would lose its stability due to the occurrence of Neimark-Sacker bifurcation when v2>v2∗∗.
However, since the discriminant of the characteristic polynomial is always positive, i.e.,Tr2E3−4DetE3=C1B1v1−2C2v22/C32−4C1C2v1v2/C3>0, all the eigenvalues of the Jacobian matrix evaluated at E3 cannot be a complex number. And, thus the stability of E3 cannot be destroyed by a Neimark-Sacker bifurcation. After the discussion of the local stability analysis of the equilibrium points, we will study the influence of varying parameters on the stability region and the path to chaos for the nonlinear system (9).
For a given fixed parameters’ setb=0.5,α=0.9,β=0.75,a−c2=7.99, Figure 1 shows the change of stable region in the parameter space v1,v2 for different values of a−c1. It can be seen that the Nash equilibrium point of system (9) will lose stability through two different routes when a−c1=4.57 (Figure 1(a)). The first route is the value of v1 or v2 in the stability domain pass through the arc M1N1 or arc N2M2, and a flip bifurcation will occur at such a situation. That is, the system will bifurcate from period-1 to period-2, then to period-4, period-8, and so on, and the system finally enters into the chaotic state. Another route is that the value(s) of v1 or/and v2 in the stability domain pass through the arc N1N2. According to the stability analysis given above, it is known that the Nash equilibrium of system (9) cannot lose stability through a Neimark-Sacker bifurcation. So, when the point in the stability domain crosses this arc N1N2, the system will directly enter into an escaping state. If we choose another value of a−c1 as a−c1=5.3, keep the values of other parameters unchanged; it can be seen from Figure 1(b) that points N1 and N2 have coincided. And, the stable Nash equilibrium point can only lose its stability through a flip bifurcation in such a situation. And, we can found that the stability region in Figure 1(b) becomes smaller than that of Figure 1(a). It can be seen from the exact expression of E3 that the value of E3 is independent of the speed of adjustment v1 and v2. However, the stability of E3 is easily affected by the values of v1 and v2 (Figure 1).Figure 1
The stable region in the parameter spacev1,v2 with different values of a−c1, where the parameters are fixed as b=0.5, α=0.9,β=0.75,anda−c2=7.99. (a) a−c1=4.57. (b) a−c1=5.3.
(a)(b)


## 3.1. Stability of the Boundary Equilibria

The stability condition of equilibrium points of model (10) is generally discussed by solving the eigenvalues of Jacobian matrix. The Jacobian matrix of model (10) at any point q1,q2 is given as(14)Jq1,q2=1+v1A1−2B1bq1−B2bq2−B2bv1q1−bv2q21+v2A2−bq1−4bq2.Through analyzing the Jacobian matrix evaluated at the boundary equilibria, we can get the following propositions.Proposition 1.
E0 is an unstable node.Proof.
By calculation, we can get the Jacobian matrix atE0 as(15)J0,0=1+A1v1001+A2v2.
And, the eigenvalues ofJ0,0 can be obtained as λ1=1+A1v1 and λ2=1+A2v2 by solving the equation λJ−E=0. Because the conditions A1,A2>0 and v1,v2>0 are always satisfied, then both of the eigenvalues meet λ1,2>1. So, E0 is an unstable node.Proposition 2.
When the parameters meet the conditionsA1,A2,B1,B2,C1,C2,C3∈S2 and v2<2/A2, the boundary equilibrium point E1 is a saddle point. When v2<2/A2, the boundary equilibrium point E1 is an unstable node.Proof.
We can get the Jacobian matrix at the equilibrium pointE1 by calculation, which is given as(16)J0,A22b=1+A1v1−12A2B2v10−12A2v21−A2v2.
At the eigenvalues of Jacobian matrix evaluated at equilibrium point,E1 are given by λ1=1+1/22A1−A2B2v1 and λ2=1−A2v2, respectively. Since A1,A2,B1,B2,C1,C2,C3∈S2, it can be found that the eigenvalues meet λ1>1 and λ2<1 when v2<2/A2. And, hence the boundary equilibrium point E1 is a saddle point, when v2<2/A2, the eigenvalues satisfy λ1>1 and λ2>1, so the equilibrium point E1 is an unstable node.
Because the boundary equilibrium pointE2 is located on the coordinate axis as well, the same method can be used to analyze its stability and similar conclusions can be obtained. So, regarding the stability analysis of E2, we will not go into details here.


## 3.2. Stability of the Nash Equilibrium

For the stability of the Nash equilibrium point, we can derive the Proposition3 by solving the eigenvalues of the Jacobian matrix.Proposition 3.
If the system parameters belong to the setA1,A2,B1,B2,C1,C2,C3∈S2, the following conditions are always satisfied:(17)C24−C1v1v2>4C3+2B1C1v1,C22−C1v1v2<B1C1v1,where the Nash equilibrium point E3 will be locally stable.Proof.
First of all, if we calculate Jacobian matrix at Nash equilibrium pointE3, then we can get(18)JE3=1+C1B1C3v1C1B2C3v1−C2C3v21−2C2C3v2.
The characteristic polynomial of the Jacobian matrix evaluated at the Nash equilibrium pointE3 can be written as(19)Fλ=λ2−TrE3λ+DetE3,where TrE3 and DetE3 are the trace and determinant of the Jacobian matrix evaluated at E3, respectively. And, the specific mathematical expression of them can be obtained through a calculation, which are given as(20)Tr=2+C1B1v1−2C2v2C3,Det=C3+C1B1v1−2C2v2+C1C2v1v2C3.
The equilibrium point will be stable when the Jury condition [28] is met, i.e.,(i)
1+Tr+Det=4C3+2B1C1v1−4C2v2+C1C2v1v2/C3>0(ii)
1−Tr+Det=C1C2v1v2/C3>0(iii)
1−Det=2C2v2−B1C1v1−C1C2v1v2/C3>0
Obviously, when the parameters belong to the setA1,A2,B1,B2,C1,C2,C3∈S2, the condition ii always hold. And thus, the Nash equilibrium point will be stable if the other two conditions are satisfied. By a complicated calculation, the other two conditions can be rewritten as the following inequality system, which is given by(21)C24−C1v1v2>4C3+2B1C1v1,C22−C1v1v2<B1C1v1.
The above inequalities are so complicated that we have to discuss them under different situations in order to get a clearer stability condition.(i)
Case I: if the parameters meet the conditionsC1v1>4 and 4C3+2B1C1v1<0, then we can get v2>0>4C3+2B1C1v1/4−C1v1C2. While if the parameters meet C1v1>4 and 4C3+2B1C1v1>0, then we can get v2>4C3+2B1C1v1/4−C1v1C2>0.(ii)
Case II: if the parameters meet the following two conditions as2<C1v1<4 and 4C3+2B1C1v1/4−C1v1C2>B1C1v1/2−C1v1C2, then we can get 0<v2<B1C1v1/2−C1v1C2. While if the parameters meet the conditions as 2<C1v1<4 and 4C3+2B1C1v1/4−C1v1C2<B1C1v1/2−C1v1C2, then we can get 0<v2<4C3+2B1C1v1/4−C1v1C2.(iii)
Case III: if the conditionC1v1<2 holds, then we can get v2>0>B1C1v1/2−C1v1C2.
In particular, when we let one of the eigenvaluesλi=−1,i=1ori=2, that is, the following equation holds: 1+Tr+Det=4C3+2B1C1v1−4C2v2+C1C2v1v2/C3=0, then we can obtain a critical value of v2, i.e., v2∗=22C3+B1C1v1/4−C1v1C2. The Nash equilibrium E3 will lose its stability due to the occurrence of flip bifurcation when v2>v2∗. In other words, if the value of speed of adjustment v2 is greater than v2∗, the system will enter the chaotic state at a double speed. In addition, if we suppose 1−Det=0, then another critical value of v2 can be obtained, i.e., v2∗∗=B1C1v1/2−C1v1C2. And, the Nash equilibrium E3 would lose its stability due to the occurrence of Neimark-Sacker bifurcation when v2>v2∗∗.
However, since the discriminant of the characteristic polynomial is always positive, i.e.,Tr2E3−4DetE3=C1B1v1−2C2v22/C32−4C1C2v1v2/C3>0, all the eigenvalues of the Jacobian matrix evaluated at E3 cannot be a complex number. And, thus the stability of E3 cannot be destroyed by a Neimark-Sacker bifurcation. After the discussion of the local stability analysis of the equilibrium points, we will study the influence of varying parameters on the stability region and the path to chaos for the nonlinear system (9).
For a given fixed parameters’ setb=0.5,α=0.9,β=0.75,a−c2=7.99, Figure 1 shows the change of stable region in the parameter space v1,v2 for different values of a−c1. It can be seen that the Nash equilibrium point of system (9) will lose stability through two different routes when a−c1=4.57 (Figure 1(a)). The first route is the value of v1 or v2 in the stability domain pass through the arc M1N1 or arc N2M2, and a flip bifurcation will occur at such a situation. That is, the system will bifurcate from period-1 to period-2, then to period-4, period-8, and so on, and the system finally enters into the chaotic state. Another route is that the value(s) of v1 or/and v2 in the stability domain pass through the arc N1N2. According to the stability analysis given above, it is known that the Nash equilibrium of system (9) cannot lose stability through a Neimark-Sacker bifurcation. So, when the point in the stability domain crosses this arc N1N2, the system will directly enter into an escaping state. If we choose another value of a−c1 as a−c1=5.3, keep the values of other parameters unchanged; it can be seen from Figure 1(b) that points N1 and N2 have coincided. And, the stable Nash equilibrium point can only lose its stability through a flip bifurcation in such a situation. And, we can found that the stability region in Figure 1(b) becomes smaller than that of Figure 1(a). It can be seen from the exact expression of E3 that the value of E3 is independent of the speed of adjustment v1 and v2. However, the stability of E3 is easily affected by the values of v1 and v2 (Figure 1).Figure 1
The stable region in the parameter spacev1,v2 with different values of a−c1, where the parameters are fixed as b=0.5, α=0.9,β=0.75,anda−c2=7.99. (a) a−c1=4.57. (b) a−c1=5.3.
(a)(b)


## 4. Numerical Simulations

We have studied the stability of the equilibrium points in Section3.2. In this section, the numerical simulation methods will be employed to study the complex dynamical behaviors of system (9). The main numerical methods used in this section are the chart of dynamic modes, largest Lyapunov exponent, bifurcation diagram, basin of attraction, critical curve, and so on.It is well known that the bifurcation diagram can provide an overview of the system dynamical behaviors by drawing one or two variables. With varying parameters, the dynamical behaviors of nonlinear system will suddenly change due to the bifurcation. Firstly, we fix the parameter asa−c1=3,a−c2=3.5,b=0.4,andβ=0.5 and choose the parameters v1, v2, and α as the bifurcation parameters, where α∈0,1 is the weight of state-owned shares in enterprise 1. In particular, we choose the parameters v1 and v2 as the main research objectives, as both firms can adjust their values at any time in the process of dynamic adjustment. Figure 2 shows the chart of dynamic modes and the corresponding chart of the largest Lyapunov exponent under this parameter set. Figure 2(a) is the chart of dynamic modes in the parameter plane v1,α with v2=0.8. It is well known that the chart of dynamic modes can reflect the strategic space of these two enterprises in the market. In Figure 2, different colors indicate different periodic modes. The readers can refer to the color bar, which is given in the right side of the chart of dynamic modes, for the detailed meaning of these colors. In the chart of dynamic modes, the black region represents the period number greater than 200, which may be a quasiperiodic state or a chaotic state. However, we cannot distinguish these states very well from the chart of dynamic modes, so they are all represented as black. It can be clearly seen from Figure 2(a) that the system bifurcates from period-1 to period-2, from period-2 to period-4, and so on and finally enters the chaotic state, as the value of parameter v1 increases. This phenomenon reminds us that if companies adjust their strategies by this route, the market will eventually enter chaos at a double speed. In Figure 2(a), it is important to note that when α∈0,0.12, the system goes straight into period-2, then into period-4, period-8, and so on and then into chaotic state. Figure 2(b) shows the chart of dynamic modes of parameter plane v1,v2, where α=0.09644163. In this figure, we can see two different paths into chaos for system (9). The first way is that the system enters chaos through a flip bifurcation. Another way is that the system first experiences flip bifurcation and then goes through a Neimark-Sacker bifurcation. That is, the stable Nash equilibrium point will lose its stability through a flip bifurcation and evolve into a stable period-2 cycle. And then, the stable period-2 cycle will lose stability through a Neimark-Sacker bifurcation, and quasiperiodic state arises. This phenomenon is not contradictory with the previous conclusion, as what we have stated is that the Nash equilibrium point will not lose stability through Neimark-Sacker bifurcation. And, the phenomenon observed here is that the period-2 cycle has lost its stability through a Neimark-Sacker bifurcation.Figure 2
The chart of dynamic modes and the corresponding chart of the largest Lyapunov exponent under parametersa−c1=3, a−c2=3.5, b=0.4, and β=0.5. (a) The chart of dynamic modes of parameters space v1,α, where v2=0.8. (b) The chart of dynamic modes of parameters space v1,v2, where α=0.09644163. (c) Partial enlargement of (a). (d) The corresponding chart of the largest Lyapunov exponent of (c).
(a)(b)(c)(d)Figure2(c) is a partial enlargement of Figure 2(a) in the region 0.87,0.915×0.962,0.97. From Figure 2(c), we can find very complex dynamical behavior of system (9). In order to further enclose the complex dynamical behaviors of system (9), two different methods are employed in the following discussion. The first method is the chart of the largest Lyapunov exponent, as shown in Figure 2(d). In Figure 2(d), the part of the largest Lyapunov exponent less than 0 is represented by the gradient from blue to white, and the system is in a stable periodic state, as shown in the corresponding chart of dynamic modes [29]. The boundary between white and black indicates that the largest Lyapunov exponent is equal to 0. That is, the system (9) is in a quasiperiodic state or a bifurcation is occurring. While if the largest Lyapunov exponent is greater than 0 and less than 2, the system will be in a chaotic state. However, the escaping trajectory will arise when the largest Lyapunov exponent is equal to 2. And, the largest Lyapunov exponent is represented by a gradient from black to gray, when the system is chaotic.Another method for analyzing the subtle dynamical behaviors in Figure2(c) is the 1-D bifurcation diagram and the largest Lyapunov exponent. Figure 3(a) shows the 1-D bifurcation diagram and the largest Lyapunov exponent corresponding to the red line in Figure 2(c) with v1=0.897, where α is taken as the variable parameter. It can be seen from Figure 3(a) that the system undergoes period-4, period-8, period-16, and so on. And, the system finally enters into chaotic state after the period-16 state. Meanwhile, we can also find two periodic windows in the middle of these chaotic states around α=0.96453 and α=0.96496. And, this phenomenon can also be detected from the largest Lyapunov exponent (Figure 3(c)) for details. In addition, we can also find a “jump” phenomenon near α=0.963. Then, a fractal structure called “hub of periodicity” ([29–32]) from α=0.9655 to α=0.967, which is a ring structure formed by the forward and reverse period-doubling bifurcation sequence, can be found in Figure 2(c). Through observing the bifurcation diagram Figure 3(a), we can confirm that the bifurcation sequence in the interval α∈0.9655,0.967 is period-12⟶period-24⟶chaos⟶period-24⟶period-12. Figure 3(b) is the 1-D bifurcation diagram corresponding to the white line in Figure 2(c) with α=0.9653, where v1 is taken as the bifurcation parameter. As it can be seen from Figure 3(b), the system undergoes the chaotic state, period-8 state, and then period-16, period-24, etc through a flip bifurcation. At last, the system enters into the chaotic state again. With a further increasement of v1, a reverse period-doubling bifurcation (also called period-halving bifurcation) occurs. Additionally, a “jump” phenomenon can also be found near v1=0.8758. And then, four multiperiodic windows appear, after that the system enters the chaotic state again. Furthermore, we can also find a quasiperiodic window from v1=0.90753 to v1=0.90985, see Figure 3(d) for more details.Figure 3
The 1-D bifurcation diagram and the corresponding largest Lyapunov exponent witha−c1=3, a−c2=3.5, b=0.4, β=0.5 and v2=0.8. (a) The 1-D bifurcation diagram with varying values of α when v1=0.897. (b) The 1-D bifurcation diagram with varying value of v1 when α=0.9653. (c) The corresponding largest Lyapunov exponent of Figure 3(a). (d) The corresponding largest Lyapunov exponent of Figure 3(b).
(a)(b)(c)(d)


### 4.1. MultiStability

Another striking and interesting phenomenon in system (9) is the multistability. In fact, the coexistence of multiple equilibrium points may arise in many economic systems. The phenomenon of multistability is also called path-dependence in economics, as the final state of the system after long-run evolution would mainly depend on the initial/history conditions. Under the circumstance of same parameters and different initial conditions, the system maybe presents totally different bifurcation structures, which makes the dynamical behaviors of economical system even more complex. And, the selection of initial conditions plays a key role in analyzing the coexistence of multiple attractors. However, we must first introduce the concepts of attractor and attracting basin before the discussion of coexistence of multiple attractors. We recall that a set A∈Rn is invariant for a map T if it is mapped onto itself, i.e., TA=A. And, a closed invariant set A is an attractor if it is Lyapunov stable. And, the attracting basin of an attractor A is the set of all points x, which will eventually evolve to the attractor A as time goes on, i.e., limn⟶+∞Tnx∈A. In the following discussion, the bifurcation diagram with the so-called adiabatic initial conditions and basins of attraction are mainly employed to study the phenomena of multistability of system (9).Firstly, the influence of parameterα, which is the proportion of state-owned shares in enterprise 1, on the global dynamics of system (9) is discussed. With fixed parameters a−c1=6, a−c2=7.6, b=0.3, β=0.25, v1=0.44, and v2=0.38, it can be found that system (9) would exhibit different dynamical behaviors with adiabatic initial conditions, see Figures 4(a) and 4(b) for details. The 1-D bifurcation diagram is employed to explain this phenomenon, where the parameter α is selected as the bifurcation parameter. In Figures 4(a) and 4(b), the bifurcation diagrams are plotted with the monotonous increase and the monotonous decrease of parameter α, respectively. It can be seen from Figures 4(a) and 4(b) that the 1-D bifurcation diagrams have completely different shapes with adiabatic initial conditions, when the range of parameter α is chosen as α∈0.08,0.17 in Figures 4(a) and 4(b). This can also be confirmed from the Lyapunov exponent spectrums of system (9), see Figures 4(c) and 4(d) for details. Figures 4(c) and 4(d) are the corresponding Lyapunov exponent spectrums of Figures 4(a) and 4(b), respectively. The difference in shape between Figures 4(a) and 4(b) reflects the coexistence of multiple attractors. In order to further enclose this interesting phenomenon, the basins of attraction related to Figure 4 are plotted.Figure 4
The 1-D bifurcation diagrams with adiabatic initial conditions and the corresponding Lyapunov exponent spectrums, where the parameters are fixed asa−c1=6,a−c2=7.6,b=0.3,β=0.25,v1=0.44,andv2=0.38. The bifurcation diagram with (a) increasing parameter of α and (b) decreasing parameter of α.
(a)(b)(c)(d)From Figure5, we can see that the chaotic attractor can coexist not only with the multiperiodic attractor, but also with other chaotic attractors. In Figure 5(a), we can find that a multiperiodic attractor coexisted with a chaotic attractor, where the color of attracting domain of the multiperiod attractor is yellow and the attracting domain of the chaotic attractor is plotted using light blue, whereas the color of dark blue indicates escaping domain in Figure 5(a). In order to observe the shape of attractor more clearly, the blue attractor in Figure 5(a) is enlarged (Figure 5(b)). In Figure 5(b), we also changed the relative position of each pieces of attractor to pursue a high resolution of the attractor. Therefore, the coordinates of the figure are removed. Above methods to address figures can be suitably applied in Figures 5(d), 5(f), and 5(h). The multiperiodic attractor turns into a chaotic attractor through a Neimark-Sacker bifurcation, as the value of further increases, (Figures 5(c) and 5(e)). Compared with Figure 5(c), we can find that the structure of the chaotic attractor, which is shown in the yellow region, has changed in Figure 5(e). In Figure 5(c), the chaotic attractor in the yellow region is made up of 5 pieces, and each piece is completely enclosed. In Figure 5(c), the values of parameter α is chosen as α=0.128. However, each piece of the chaotic attractor has been broken down into 3 smaller pieces (Figure 5(e)), as the value of parameter α further increases to 0.135. The blue attractors in Figures 5(c) and 5(e) are also enlarged in Figures 5(d) and 5(f), respectively. In Figure 5(g), the coexistence of multiple attractors has once again become a chaotic attractor coexisted with a multiperiodic attractor (Figure 5(h) for enlargement), when the value of parameter α further increases to 0.15. From Figure 5, we can see that the yellow region is gradually shrinking and eventually disappears with the increase of α. It means that the phenomenon of coexistence of multiple attractors will disappear ultimately. It can also be found from Figures 4 and 5 that the parameter α not only affects the type of coexistence but also the size of attracting domain.Figure 5
Basins of attraction and partial enlargement of the attractor with different values of parameterα, where the rest of parameters are fixed as a−c1=6, a−c2=7.6, b=0.3, β=0.25, v1=0.44, and v2=0.38, and the coexistence of multiple attractors can be found. Basin of attraction with (a) α=0.087, (c) α=0.128, (e) α=0.135, and (g) α=0.15. (b) Sketch of enlargement of the blue attractor in (b) Figure 5(a), (d) Figure 5(c), (f) Figure 5(e), and (h) Figure 5(g).
(a)(b)(c)(d)(e)(f)(g)(h)


### 4.2. Global Dynamics, Critical Curve, and Noninvertible Map

The critical curve plays an extremely crucial role in analyzing the structure of attracting basin and the global bifurcation phenomenon which leads to qualitative change of basins with varying parameters. In addition, the noninvertible map is a useful tool to enclose the formation mechanism of the global dynamics in nonlinear discrete systems. Generally speaking, a noninvertible mapT refers to those maps, which have multiple preimages. For example, the one-variable quadratic map y=x2 is a typical noninvertible map, as a single value of y always corresponds to two different values of x (here, x is called the preimage of y), when y>0. Moreover, the phase space can be divided into different regions Zk, k=0,1,2,… according to the number of rank-1 preimages, and the points in region Zk always have k different rank-1 preimages. For example, the x-axis has divided the phase space of y=x2 into two different regions: Z0 and Z2. The points in the region Z0 have no preimage, while the points in Z2 have two different preimages. And, the borderline between different regions Zk is called critical curve, which is denoted by LC. The corresponding rank-1 preimage of LC under map T is denoted as LC−1, i.e., LC=TLC−1. As shown by Bischi and Lamantia [33] and Zhou and Wang [34], the curve LC−1 is the locus of points where the Jacobian determinant is equal to 0, i.e., LC−1=q1,q2∈R2Det=0. So, if we let the determinant of (14) to zero, then the following equation can be obtained:(22)MN−Mv2+2NB1v1bq1−4Mv2+NB2v1bq2+2B1b2v1v2q12+4B2b2v1v2q22+8B1b2v1v2q1q2=0,where M=1+A1v1 and N=1+A2v2. For the built model (9), the curve LC−1 has two branches, which are LC−1a and LC−1b. And, the critical line LC consists of LCa and LCb. The phase plane of system (9) is separated into three regions Z0, Z2, and Z4 by these two branches. The regions Z2 and Z4 are divided by the branch LCa, while the regions Z0 and Z2 are divided by another branch LCb. The points in these three regions have 0, 2, and 4 rank-1 preimages, respectively. We can find that the origin is always in the region Z4. That is, the origin should have four different preimages.In order to reveal the complex global dynamics of system (9), the parameters are firstly fixed as a−c1=3, a−c2=3.5, b=0.4, β=0.5, v1=0.95, and v2=0.8, and the varying parameter α is chosen as the bifurcation parameter. Figure 6 gives the basins of attraction at these reference parameters. In Figure 6(a), the value of parameter α is chosen as α=0.1, and we can see that there are a lot of “holes” in the interior of the attracting basin. This is due to the fact that the boundary of the attracting basin intersects with the critical curve LCb. And, the part enclosed by these two curves is labeled by H0, see Figure 6(a) for details. It is very clear that H0 belongs to Z2, so the enclosed area H0 should have two rank-1 preimages, which are labeled by H−1a and H−1b. And, the holes in the interior of the attracting basin are made up of H−1a, H−1b and their higher rank preimages. In particular, we can see from Figure 6(a) that the regions H−1a and H−1b, which are actually a whole hole and are called the “main hole” in the attracting basin, are separated by the curve LC−1b. And, the rest of the holes are the higher rank preimages of the “main hole”. However, we can find that the holes will become smaller and smaller as the rank of preimages increases.Figure 6
Various basins of attraction under parametersa−c1=3, a−c2=3.5, b=0.4, β=0.5, v1=0.95, and v2=0.8 with different α values. (a) α=0.1; (b) α=0.2; (c) α=0.25; (d) α=0.3; (e) α=0.36; (f) α=0.372.
(a)(b)(c)(d)(e)(f)As the value of parameterα further increases to 0.2, we can find that the critical curve LCb moves to the left, and thus the intersecting part of boundary of attracting basin and LCb is reduced. It makes the holes inside the attracting basin become smaller, as shown in Figure 6(b). In Figure 6(c), a further increase of α makes the curve LCb enter the yellow region completely. At such a situation, there is no intersection between the boundary of attracting basin and the curve LCb. And thus, there is no hole in the attracting domain. When α=0.3, the curve LCb begins to intersect with the upper boundary of the attracting basin again (Figure 6(d)). And, a lot of “holes” have reappeared. However, these holes are mainly concentrated on the left side of the attracting domain rather than the lower part, and this phenomenon is different from Figures 6(a) and 6(b). As the value of α further increases, we can find these “holes” have also become larger (Figure 6(e)). In Figure 6(e), the value of parameter α is chosen as α=0.36. Additionally, a very interesting phenomenon can be found from Figure 6(e) that the chaotic attractor has become so big that it is almost in contact with the boundary of the attracting basin. In this situation, a global bifurcation called “contact bifurcation” will happen. When α=0.372, the attractor has been in contact with the boundary of the attracting basin, and the chaotic attractor has broken into “ghosts” [35].If we choose another set of parameters, we can find totally different fractal structures in the attracting basin. For example, if we fix the parameters asa−c1=4.563, a−c2=4.576, b=0.011, α=0.649, β=0.057, and v1=0.537 and choose the parameter v2 as the bifurcation parameter, we will find different attractors and the corresponding basins of attraction with varying values of parameter v2. In Figure 7(a), the value of v2 is chosen as v2=0.586. We can find from Figure 7(a) that there is no hole in the attracting basin (the yellow area in Figure 7(a)), as the boundary of the attracting basin belongs to the region Z0. And, the critical curve begins to intersect with the upper boundary of the attracting basin, as the value of v2 increases to 0.606 (Figure 7(b)). The enclosed area surrounded by the basin’s boundary and the critical curve LCb is labeled by H0, and it is very clear that the region H0, which is also part of the escaping basin, belongs to the region Z2. It means that the region H0 should have two rank-1 preimages and a lot of higher rank preimages, which are just the holes in the attracting domain. Moreover, the holes in the attracting domain will become larger and larger, as the value of parameter v2 further increases to v2=0.636 (Figure 7(c)). However, these “holes” do not break, as there are still two intersecting points between the critical curve and the boundary of attracting basin, as shown in Figure 7(c). With further increase of v2, it can be seen in Figures 7(d)–7(f) that there is only one intersection point between the critical curve and the basin of attraction, so the holes in the basins have been broken.Figure 7
Various basins of attraction under parametersa−c1=3,a−c2=3.5,b=0.4,β=0.5,v1=0.95,andv2=0.8 with different v2. (a) v2=0.586; (b) v2=0.606; (c) v2=0.636; (d) v2=0.666; (e) v2=0.696; (f) v2=0.796.
(a)(b)(c)(d)(e)(f)In addition, one of the holes in the attracting basin, which is labeled byI0, has been cut by the critical curve LCa, so a new hole I1 will be created. However, since I1 belongs to the region Z0, it will not produce new holes any more, which is different from H0. In Figure 7, we can see that the change of attractor is as follows: 4-cyclic, 8-cyclic, 4 pieces, and 2 pieces and finally turn into flower-like chaotic attractor through a flip bifurcation. It is worth pointing out that the type of attractor is independent of the structure of the attracting basin. It means that the fractal structure of attracting basin can be very complex when the attractor is simple. However, the attracting domain for the very complicated attractor can be very simple, too.


## 4.1. MultiStability

Another striking and interesting phenomenon in system (9) is the multistability. In fact, the coexistence of multiple equilibrium points may arise in many economic systems. The phenomenon of multistability is also called path-dependence in economics, as the final state of the system after long-run evolution would mainly depend on the initial/history conditions. Under the circumstance of same parameters and different initial conditions, the system maybe presents totally different bifurcation structures, which makes the dynamical behaviors of economical system even more complex. And, the selection of initial conditions plays a key role in analyzing the coexistence of multiple attractors. However, we must first introduce the concepts of attractor and attracting basin before the discussion of coexistence of multiple attractors. We recall that a set A∈Rn is invariant for a map T if it is mapped onto itself, i.e., TA=A. And, a closed invariant set A is an attractor if it is Lyapunov stable. And, the attracting basin of an attractor A is the set of all points x, which will eventually evolve to the attractor A as time goes on, i.e., limn⟶+∞Tnx∈A. In the following discussion, the bifurcation diagram with the so-called adiabatic initial conditions and basins of attraction are mainly employed to study the phenomena of multistability of system (9).Firstly, the influence of parameterα, which is the proportion of state-owned shares in enterprise 1, on the global dynamics of system (9) is discussed. With fixed parameters a−c1=6, a−c2=7.6, b=0.3, β=0.25, v1=0.44, and v2=0.38, it can be found that system (9) would exhibit different dynamical behaviors with adiabatic initial conditions, see Figures 4(a) and 4(b) for details. The 1-D bifurcation diagram is employed to explain this phenomenon, where the parameter α is selected as the bifurcation parameter. In Figures 4(a) and 4(b), the bifurcation diagrams are plotted with the monotonous increase and the monotonous decrease of parameter α, respectively. It can be seen from Figures 4(a) and 4(b) that the 1-D bifurcation diagrams have completely different shapes with adiabatic initial conditions, when the range of parameter α is chosen as α∈0.08,0.17 in Figures 4(a) and 4(b). This can also be confirmed from the Lyapunov exponent spectrums of system (9), see Figures 4(c) and 4(d) for details. Figures 4(c) and 4(d) are the corresponding Lyapunov exponent spectrums of Figures 4(a) and 4(b), respectively. The difference in shape between Figures 4(a) and 4(b) reflects the coexistence of multiple attractors. In order to further enclose this interesting phenomenon, the basins of attraction related to Figure 4 are plotted.Figure 4
The 1-D bifurcation diagrams with adiabatic initial conditions and the corresponding Lyapunov exponent spectrums, where the parameters are fixed asa−c1=6,a−c2=7.6,b=0.3,β=0.25,v1=0.44,andv2=0.38. The bifurcation diagram with (a) increasing parameter of α and (b) decreasing parameter of α.
(a)(b)(c)(d)From Figure5, we can see that the chaotic attractor can coexist not only with the multiperiodic attractor, but also with other chaotic attractors. In Figure 5(a), we can find that a multiperiodic attractor coexisted with a chaotic attractor, where the color of attracting domain of the multiperiod attractor is yellow and the attracting domain of the chaotic attractor is plotted using light blue, whereas the color of dark blue indicates escaping domain in Figure 5(a). In order to observe the shape of attractor more clearly, the blue attractor in Figure 5(a) is enlarged (Figure 5(b)). In Figure 5(b), we also changed the relative position of each pieces of attractor to pursue a high resolution of the attractor. Therefore, the coordinates of the figure are removed. Above methods to address figures can be suitably applied in Figures 5(d), 5(f), and 5(h). The multiperiodic attractor turns into a chaotic attractor through a Neimark-Sacker bifurcation, as the value of further increases, (Figures 5(c) and 5(e)). Compared with Figure 5(c), we can find that the structure of the chaotic attractor, which is shown in the yellow region, has changed in Figure 5(e). In Figure 5(c), the chaotic attractor in the yellow region is made up of 5 pieces, and each piece is completely enclosed. In Figure 5(c), the values of parameter α is chosen as α=0.128. However, each piece of the chaotic attractor has been broken down into 3 smaller pieces (Figure 5(e)), as the value of parameter α further increases to 0.135. The blue attractors in Figures 5(c) and 5(e) are also enlarged in Figures 5(d) and 5(f), respectively. In Figure 5(g), the coexistence of multiple attractors has once again become a chaotic attractor coexisted with a multiperiodic attractor (Figure 5(h) for enlargement), when the value of parameter α further increases to 0.15. From Figure 5, we can see that the yellow region is gradually shrinking and eventually disappears with the increase of α. It means that the phenomenon of coexistence of multiple attractors will disappear ultimately. It can also be found from Figures 4 and 5 that the parameter α not only affects the type of coexistence but also the size of attracting domain.Figure 5
Basins of attraction and partial enlargement of the attractor with different values of parameterα, where the rest of parameters are fixed as a−c1=6, a−c2=7.6, b=0.3, β=0.25, v1=0.44, and v2=0.38, and the coexistence of multiple attractors can be found. Basin of attraction with (a) α=0.087, (c) α=0.128, (e) α=0.135, and (g) α=0.15. (b) Sketch of enlargement of the blue attractor in (b) Figure 5(a), (d) Figure 5(c), (f) Figure 5(e), and (h) Figure 5(g).
(a)(b)(c)(d)(e)(f)(g)(h)


## 4.2. Global Dynamics, Critical Curve, and Noninvertible Map

The critical curve plays an extremely crucial role in analyzing the structure of attracting basin and the global bifurcation phenomenon which leads to qualitative change of basins with varying parameters. In addition, the noninvertible map is a useful tool to enclose the formation mechanism of the global dynamics in nonlinear discrete systems. Generally speaking, a noninvertible mapT refers to those maps, which have multiple preimages. For example, the one-variable quadratic map y=x2 is a typical noninvertible map, as a single value of y always corresponds to two different values of x (here, x is called the preimage of y), when y>0. Moreover, the phase space can be divided into different regions Zk, k=0,1,2,… according to the number of rank-1 preimages, and the points in region Zk always have k different rank-1 preimages. For example, the x-axis has divided the phase space of y=x2 into two different regions: Z0 and Z2. The points in the region Z0 have no preimage, while the points in Z2 have two different preimages. And, the borderline between different regions Zk is called critical curve, which is denoted by LC. The corresponding rank-1 preimage of LC under map T is denoted as LC−1, i.e., LC=TLC−1. As shown by Bischi and Lamantia [33] and Zhou and Wang [34], the curve LC−1 is the locus of points where the Jacobian determinant is equal to 0, i.e., LC−1=q1,q2∈R2Det=0. So, if we let the determinant of (14) to zero, then the following equation can be obtained:(22)MN−Mv2+2NB1v1bq1−4Mv2+NB2v1bq2+2B1b2v1v2q12+4B2b2v1v2q22+8B1b2v1v2q1q2=0,where M=1+A1v1 and N=1+A2v2. For the built model (9), the curve LC−1 has two branches, which are LC−1a and LC−1b. And, the critical line LC consists of LCa and LCb. The phase plane of system (9) is separated into three regions Z0, Z2, and Z4 by these two branches. The regions Z2 and Z4 are divided by the branch LCa, while the regions Z0 and Z2 are divided by another branch LCb. The points in these three regions have 0, 2, and 4 rank-1 preimages, respectively. We can find that the origin is always in the region Z4. That is, the origin should have four different preimages.In order to reveal the complex global dynamics of system (9), the parameters are firstly fixed as a−c1=3, a−c2=3.5, b=0.4, β=0.5, v1=0.95, and v2=0.8, and the varying parameter α is chosen as the bifurcation parameter. Figure 6 gives the basins of attraction at these reference parameters. In Figure 6(a), the value of parameter α is chosen as α=0.1, and we can see that there are a lot of “holes” in the interior of the attracting basin. This is due to the fact that the boundary of the attracting basin intersects with the critical curve LCb. And, the part enclosed by these two curves is labeled by H0, see Figure 6(a) for details. It is very clear that H0 belongs to Z2, so the enclosed area H0 should have two rank-1 preimages, which are labeled by H−1a and H−1b. And, the holes in the interior of the attracting basin are made up of H−1a, H−1b and their higher rank preimages. In particular, we can see from Figure 6(a) that the regions H−1a and H−1b, which are actually a whole hole and are called the “main hole” in the attracting basin, are separated by the curve LC−1b. And, the rest of the holes are the higher rank preimages of the “main hole”. However, we can find that the holes will become smaller and smaller as the rank of preimages increases.Figure 6
Various basins of attraction under parametersa−c1=3, a−c2=3.5, b=0.4, β=0.5, v1=0.95, and v2=0.8 with different α values. (a) α=0.1; (b) α=0.2; (c) α=0.25; (d) α=0.3; (e) α=0.36; (f) α=0.372.
(a)(b)(c)(d)(e)(f)As the value of parameterα further increases to 0.2, we can find that the critical curve LCb moves to the left, and thus the intersecting part of boundary of attracting basin and LCb is reduced. It makes the holes inside the attracting basin become smaller, as shown in Figure 6(b). In Figure 6(c), a further increase of α makes the curve LCb enter the yellow region completely. At such a situation, there is no intersection between the boundary of attracting basin and the curve LCb. And thus, there is no hole in the attracting domain. When α=0.3, the curve LCb begins to intersect with the upper boundary of the attracting basin again (Figure 6(d)). And, a lot of “holes” have reappeared. However, these holes are mainly concentrated on the left side of the attracting domain rather than the lower part, and this phenomenon is different from Figures 6(a) and 6(b). As the value of α further increases, we can find these “holes” have also become larger (Figure 6(e)). In Figure 6(e), the value of parameter α is chosen as α=0.36. Additionally, a very interesting phenomenon can be found from Figure 6(e) that the chaotic attractor has become so big that it is almost in contact with the boundary of the attracting basin. In this situation, a global bifurcation called “contact bifurcation” will happen. When α=0.372, the attractor has been in contact with the boundary of the attracting basin, and the chaotic attractor has broken into “ghosts” [35].If we choose another set of parameters, we can find totally different fractal structures in the attracting basin. For example, if we fix the parameters asa−c1=4.563, a−c2=4.576, b=0.011, α=0.649, β=0.057, and v1=0.537 and choose the parameter v2 as the bifurcation parameter, we will find different attractors and the corresponding basins of attraction with varying values of parameter v2. In Figure 7(a), the value of v2 is chosen as v2=0.586. We can find from Figure 7(a) that there is no hole in the attracting basin (the yellow area in Figure 7(a)), as the boundary of the attracting basin belongs to the region Z0. And, the critical curve begins to intersect with the upper boundary of the attracting basin, as the value of v2 increases to 0.606 (Figure 7(b)). The enclosed area surrounded by the basin’s boundary and the critical curve LCb is labeled by H0, and it is very clear that the region H0, which is also part of the escaping basin, belongs to the region Z2. It means that the region H0 should have two rank-1 preimages and a lot of higher rank preimages, which are just the holes in the attracting domain. Moreover, the holes in the attracting domain will become larger and larger, as the value of parameter v2 further increases to v2=0.636 (Figure 7(c)). However, these “holes” do not break, as there are still two intersecting points between the critical curve and the boundary of attracting basin, as shown in Figure 7(c). With further increase of v2, it can be seen in Figures 7(d)–7(f) that there is only one intersection point between the critical curve and the basin of attraction, so the holes in the basins have been broken.Figure 7
Various basins of attraction under parametersa−c1=3,a−c2=3.5,b=0.4,β=0.5,v1=0.95,andv2=0.8 with different v2. (a) v2=0.586; (b) v2=0.606; (c) v2=0.636; (d) v2=0.666; (e) v2=0.696; (f) v2=0.796.
(a)(b)(c)(d)(e)(f)In addition, one of the holes in the attracting basin, which is labeled byI0, has been cut by the critical curve LCa, so a new hole I1 will be created. However, since I1 belongs to the region Z0, it will not produce new holes any more, which is different from H0. In Figure 7, we can see that the change of attractor is as follows: 4-cyclic, 8-cyclic, 4 pieces, and 2 pieces and finally turn into flower-like chaotic attractor through a flip bifurcation. It is worth pointing out that the type of attractor is independent of the structure of the attracting basin. It means that the fractal structure of attracting basin can be very complex when the attractor is simple. However, the attracting domain for the very complicated attractor can be very simple, too.


## 5. Long-Term Average Profit Analysis

In the previous discussion, the dynamical behaviors of system (9) have been analyzed using different methods. However, all the enterprises want to know more about their individual profits under a very complex competitive environment. In addition, they also want to know under what conditions they can make more profits. It is generally well known that chaos is harmful to the economics. However, it was pointed out by Wu et al. [36] and Matsumoto [37, 38] that chaos may bring higher profits in a long-term competition. In this section, the long-run profits of the participants in the built Cournot duopoly model will be discussed through the 1-D bifurcation diagram and the corresponding profit curve.In the following discussion, the parameters are fixed asa−c1=1.8, a−c2=1.76, b=0.23, v2=0.62, α=0.07, and β=0.4349294, and the parameter v1 is chosen as the bifurcation parameter. Figure 8(a) and Figure 8(b) have shown the 1-D bifurcation diagram and the corresponding long-run profits. By comparing Figures 8(a) and 8(b), we can find that these two figures have very good consistency. It can be seen from Figure 8 that the average profits equal the equilibrium profits when the system is in a stable period-1 state, and the profit of enterprise 1 is higher than that of enterprise 2. When v1=1.3782, the first flip bifurcation occurs. At such a situation, the average profit of enterprise 1 starts to decline, while the average profit of enterprise 2 begins to increase. When v1≈1.78, these two enterprises have the same average profits. And then, with a further increase of v1, the average profit of enterprise 2 will be higher than that of enterprise 1. When the system enters into the chaotic state, we can clearly see the average profit of enterprise 1 has dropped rapidly, but the average profit of enterprise 2 has increased sharply. Therefore, we can draw a conclusion that chaos is not harmful to every enterprise; some enterprises can even earn more from chaos. Those enterprises are just the players who fish in the trouble water. After that, the system enters a periodic state, and the average profit of enterprise 1 continues to decline. And, the sustainable growth in the average profit of enterprise 2 can also be found in Figure 8(b).Figure 8
The parameters are chosen asa−c1=1.8,a−c2=1.76,b=0.23,v2=0.62,α=0.07,andβ=4.349294e−01. (a) The 1-D bifurcation diagram. (b) The profit curves of the two enterprises under the given parameters, where blue represents the profit curve of enterprise 1 and red represents the profit curve of enterprise 2.
(a)(b)


## 6. Conclusions

In this paper, a Cournot model with homogeneous products is built. In the built model, enterprise 1 wants to make the weight of its own profit and the payoff of government maximization, while enterprise 2 just wants to make its individual profit maximization. And then, the local stability and the complex dynamical behaviors of the given system are discussed through the analytical method and numerical method. Firstly, the local stability of all the equilibrium points is studied, and we find that the Neimark-Sacker bifurcation cannot occur in our model. It is generally known that Neimark-Sacker bifurcation often means that the economic system will oscillate and become unstable. And, this kind of oscillation will also cause certain damage to the normal operation of the economic system. In addition, we also find that the stability region of Nash equilibrium point is not only related to the controlling interest of state-owned enterprise but also greatly affected by the parametera−c1.Secondly, the dynamical behaviors of the model are analyzed by means of numerical simulation. A lot of numerical methods including bifurcation diagram, largest Lyapunov exponent, basin of attraction, and critical curve are employed to reveal the complex dynamical behaviors of system (9). Furthermore, the chart of dynamic modes with different parameter’s combination is used to analyze the complex dynamical behaviors in this article. And, a very special fractal structure called “eye of chaos,” which is a ring structure composed of forward and reverse period-doubling bifurcation sequences, has been found and discussed in our research.Thirdly, the multistability and the coexistence of attractors are explained by using the 1-D bifurcation diagrams with different initial conditions and basins of attraction. Moreover, the evolution of attractors and basins of attraction with varying parameters are also analyzed. We find that there are two different coexistences of multiple attractors, which are a multiperiodic attractor coexisting with a chaotic attractors and a chaotic attractor coexisting with another chaotic attractor. And, the phenomenon of coexistence of multiple attractors will disappear when the proportion of state-owned shares in enterprise 1 is large enough. In addition, the fractal structure of attracting basin is also analyzed using the critical curve. And, we find that the “holes” might appear in the attracting domain. The formation mechanism of these “holes” has been revealed, and we find that the appearance of these “holes” is related to the intersection between the boundary of attracting basin and the critical curves. In the process of studying multistability of the given model, a very important global bifurcation called “contact bifurcation” has also been found in our research. The contact bifurcation means that the attractor has come into contact with the boundary of the attracting basin in which the attractor is located. And, both the attractor and the attracting basin will be destroyed after the contact bifurcation, and only their “ghosts” remain in the phase space. At last, the long-term profits of enterprises are analyzed in our research. And, we find that some enterprises can even earn more from chaos.



---

*Source: 4580415-2020-08-12.xml*