Stability andMultistability of a Bounded RationalMixedDuopoly Model with Homogeneous Product

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. *e purpose of this research is to study the stability and the multistability of the establishedmodel.*e 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.


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][2][3][4][5][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][8][9][10][11]. ereafter, 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 twostage Cournot duopoly model with R&D spillover effect.
However, enterprises are always assumed to be private enterprises in the above discussion. e 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. erefore, the analysis of mixed duopoly model has attracted the attention of more and more researchers in recent years. e 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]. e 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]. e welfare of partial privatization of stateowned 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 publicprivate 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. e 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 Section 6.

The Model
In this section, a mixed duopoly model with two enterprises, which are labeled by i (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 p i and q i , (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 where Q � q 1 + q 2 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 , where c i > 0 represents the marginal cost of enterprise i and a > c i always holds. Due to the difference of production technology between enterprise 1 and enterprise 2, we assume that c 1 ≠ c 2 . us, the profit function of enterprise i can be expressed as 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. e payoff of government U G 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, And, the payoff of government U G is given as where β ≥ 0 represents the weight of consumer surplus in the payoff of government U G . e 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 2 Complexity where α ∈ [0, 1] is the proportion of state-owned shares in enterprise 1.
e first-order conditions of U 1 and U 2 can be written as follows: 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 period t + 1 according to its own output q i (t) and the marginal objective zU i /zq i , (i � 1, 2) in the current period.
at 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 where V i (q i ) is a positive function which gives the extent of production variation of the ith enterprise following a given signal zU i /zq i , (i � 1, 2). An adjustment mechanism of (7) has been proposed by some researchers [26] with a constant V i (q i ) � v i . Instead, we assume that the function V i (q i ) is a linear function here, i.e., V i (q i ) � v i q i , which is equivalent to the assumption that the relative production change is proportional to the estimated marginal objective [27], that is, where v i > 0 represents the speed of adjustment of enterprise i. By substituting the first-order conditions of U 1 and U 2 into (7), we can establish the following dynamic model:

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 define a − c 1 � A 1 > 0, a − c 2 � A 2 > 0, 2 − α(1 + β) � B 1 , and 1 − αβ � B 2 ; thus the system (9) can be simplified to T: In order to gain the equilibrium point, we let q i (t + 1) � q i (t) in (10). By calculation, we can find four equilibrium points of system (10): where Because the equilibrium points E 0 , E 1 , and E 2 are all on the coordinate axis, we call them boundary equilibrium points, and E 3 is called Nash equilibrium point or interior equilibrium point. Obviously, the equilibrium point E 3 exists if and only if B 2 − 2B 1 ≠ 0. e non-negativity of equilibrium point E 2 is satisfied when B 1 > 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 B 1 > 0. In order to make the equilibrium point E 3 non-negative, the above defined parameters must meet one of the following conditions: e stability condition of equilibrium points of model (10) is generally discussed by solving the eigenvalues of Jacobian matrix. e Jacobian matrix of model (10) at any point (q 1 , q 2 ) is given as rough analyzing the Jacobian matrix evaluated at the boundary equilibria, we can get the following propositions.
Proof. By calculation, we can get the Jacobian matrix at E 0 as And, the eigenvalues of J(0, 0) can be obtained as are always satisfied, then both of the eigenvalues meet λ 1,2 > 1. So, E 0 is an unstable node.

Proposition 2. When the parameters meet the conditions
Proof. We can get the Jacobian matrix at the equilibrium point E 1 by calculation, which is given as At the eigenvalues of Jacobian matrix evaluated at equilibrium point, E 1 are given by it can be found that the eigenvalues meet λ 1 > 1 and λ 2 < 1 when v 2 < (2/A 2 ). And, hence the boundary equilibrium point E 1 is a saddle point, when v 2 < (2/A 2 ), the eigenvalues satisfy λ 1 > 1 and |λ 2 | > 1, so the equilibrium point E 1 is an unstable node.
Because the boundary equilibrium point E 2 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 E 2 , we will not go into details here.

Stability of the Nash Equilibrium.
For the stability of the Nash equilibrium point, we can derive the Proposition 3 by solving the eigenvalues of the Jacobian matrix.

Proposition 3. If the system parameters belong to the set
the following conditions are always satisfied: where the Nash equilibrium point E 3 will be locally stable.
Proof. First of all, if we calculate Jacobian matrix at Nash equilibrium point E 3 , then we can get e characteristic polynomial of the Jacobian matrix evaluated at the Nash equilibrium point E 3 can be written as where Tr(E 3 ) and Det(E 3 ) are the trace and determinant of the Jacobian matrix evaluated at E 3 , respectively. And, the specific mathematical expression of them can be obtained through a calculation, which are given as e equilibrium point will be stable when the Jury condition [28] is met, i.e., Obviously, when the parameters belong to the set A 1 , A 2 , B 1 , B 2 , C 1 , C 2 , C 3 ∈ S 2 , 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 4 Complexity e 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 conditions (ii) Case II: if the parameters meet the following two conditions as 2 In particular, when we let one of the eigenvalues , that is, the following equation holds: e Nash equilibrium E 3 will lose its stability due to the occurrence of flip bifurcation when v 2 > v * 2 . In other words, if the value of speed of adjustment v 2 is greater than v * 2 , the system will enter the chaotic state at a double speed. In addition, if we suppose 1 − Det � 0, then another critical value of v 2 can be obtained, i.e., v * * . And, the Nash equilibrium E 3 would lose its stability due to the occurrence of Neimark-Sacker bifurcation when v 2 > v * * 2 . However, since the discriminant of the characteristic polynomial is always positive, i.e., all the eigenvalues of the Jacobian matrix evaluated at E 3 cannot be a complex number. And, thus the stability of E 3 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' Figure 1 shows the change of stable region in the parameter space (v 1 , v 2 ) for different values of a − c 1 . It can be seen that the Nash equilibrium point of system (9) will lose stability through two different routes when a − c 1 � 4.57 (Figure 1(a)). e first route is the value of v 1 or v 2 in the stability domain pass through the arc M 1 N 1 or arc N 2 M 2 , and a flip bifurcation will occur at such a situation. at 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 v 1 or/and v 2 in the stability domain pass through the arc N 1 N 2 . 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 N 1 N 2 , the system will directly enter into an escaping state. If we choose another value of a − c 1 as a − c 1 � 5.3, keep the values of other parameters unchanged; it can be seen from Figure 1(b) that points N 1 and N 2 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 E 3 that the value of E 3 is independent of the speed of adjustment v 1 and v 2 . However, the stability of E 3 is easily affected by the values of v 1 and v 2 ( Figure 1).

Numerical Simulations
We have studied the stability of the equilibrium points in Section 3.2. In this section, the numerical simulation methods will be employed to study the complex dynamical behaviors of system (9). e 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 as a − c 1 � 3, a − c 2 � 3.5, b � 0.4, and β � 0.5 and choose the parameters v 1 , v 2 , and α as the bifurcation parameters, where α ∈ [0, 1] is the weight of state-owned shares in enterprise 1. In particular, we choose the parameters v 1 and v 2 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 (v 1 , α) with v 2 � 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. e 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 v 1 increases.
is 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 6 Complexity parameter plane (v 1 , v 2 ), where α � 0.09644163. In this figure, we can see two different paths into chaos for system (9). e 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. at 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.
is 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(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. e 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]. e boundary between white and black indicates that the largest Lyapunov exponent is equal to 0. at 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 Figure 2(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 v 1 � 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. en, a fractal structure called "hub of periodicity" ( [29][30][31][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). rough 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 v 1 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 v 1 , a reverse perioddoubling bifurcation (also called period-halving bifurcation) occurs. Additionally, a "jump" phenomenon can also be found near v 1 � 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 v 1 � 0.90753 to v 1 � 0.90985, see Figure 3(d) for more details. (9) is the multistability. In fact, the coexistence of multiple equilibrium points may arise in many economic systems. e 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 ∈ R n is invariant for a map T if it is mapped onto itself, i.e., T(A) � 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., lim n⟶+∞ T n (x) ∈ 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).

MultiStability. Another striking and interesting phenomenon in system
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 − c 1 � 6, a − c 2 � 7.6, b � 0.3, β � 0.25, v 1 � 0.44, and v 2 � 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. e 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). is can also be confirmed from the Lyapunov exponent spectrums of system (9) multiple attractors. In order to further enclose this interesting phenomenon, the basins of attraction related to Figure 4 are plotted. From Figure 5, 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. erefore, 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). e 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(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.

Global Dynamics, Critical Curve, and Noninvertible Map.
e 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  8 Complexity speaking, a noninvertible map T refers to those maps, which have multiple preimages. For example, the one-variable quadratic map y � x 2 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 Z k , (k � 0, 1, 2, . . .) according to the number of rank-1 preimages, and the points in region Z k always have k different rank-1 preimages. For example, the x-axis has divided the phase space of y � x 2 into two different regions: Z 0 and Z 2 . e points in the region Z 0 have no preimage, while the points in Z 2 have two different preimages. And, the borderline between different regions Z k is called critical curve, which is denoted by LC. e corresponding rank-1 preimage of LC under map T is denoted as LC − 1 , i.e., LC � T(LC − 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 � (q 1 , q 2 ) ∈ R 2 | De t � 0 . So, if we let the determinant of (14) to zero, then the following equation can be obtained: where M � 1 + A 1 v 1 and N � 1 + A 2 v 2 . For the built model (9), the curve LC − 1 has two branches, which are LC (a) − 1 and LC (b) − 1 . And, the critical line LC consists of LC (a) and LC (b) . e phase plane of system (9) is separated into three regions Z 0 , Z 2 , and Z 4 by these two branches. e regions Z 2 and Z 4 are divided by the branch LC (a) , while the regions Z 0 and Z 2 are divided by another branch LC (b) . e 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 Z 4 . at 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 − c 1 � 3, a − c 2 � 3.5, b � 0.4, β � 0.5, v 1 � 0.95, and v 2 � 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. is is due to the fact that the boundary of the attracting basin intersects with the critical curve LC (b) . And, the part enclosed by these two curves is labeled by H 0 , see Figure 6   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.
As the value of parameter α further increases to 0.2, we can find that the critical curve LC (b) moves to the left, and thus the intersecting part of boundary of attracting basin and LC (b) 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 LC (b) enter the yellow region completely. At such a situation, there is no intersection between the boundary of attracting basin and the curve LC (b) . And thus, there is no hole in the attracting domain. When α � 0.3, the curve LC (b) 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 as a − c 1 � 4.563, a − c 2 � 4.576, b � 0.011, α � 0.649, β � 0.057, and v 1 � 0.537 and choose the parameter v 2 as the bifurcation parameter, we will find different attractors and the corresponding basins of attraction with varying values of parameter v 2 . In Figure 7(a), the value of v 2 is chosen as v 2 � 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 Z 0 . And, the critical curve begins to intersect with the upper boundary of the attracting basin, as the value of v 2 increases to 0.606 (Figure 7(b)). e enclosed area surrounded by the basin's boundary and the critical curve LC (b) is labeled by H 0 , and it is very clear that the region H 0 , which is also part of the escaping basin, belongs to the region Z 2 . It means that the region H 0 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 v 2 further increases to v 2 � 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 v 2 , 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.
In addition, one of the holes in the attracting basin, which is labeled by I 0 , has been cut by the critical curve LC (a) , so a new hole I 1 will be created. However, since I 1 belongs to the region Z 0 , it will not produce new holes any more, which is different from H 0 . 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.

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 as a − c 1 � 1.8, a − c 2 � 1.76, b � 0.23, v 2 � 0.62, α � 0.07, and β � 0.4349294, and the parameter v 1 is chosen as the bifurcation parameter. Figure 8(a) and Figure 8(b) have shown the 1-D bifurcation diagram and the corresponding longrun 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 v 1 � 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 v 1 ≈ 1.78, these two enterprises have the same average profits. And then, with a further increase of v 1 , 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. erefore, we can draw a conclusion that chaos is not harmful to every enterprise; some enterprises can even earn more from chaos. ose 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).

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 parameter a − c 1 . 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.
irdly, 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. e 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. e 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.

Data Availability
No data were used to support the findings of this study.

Conflicts of Interest
e authors declare that they have no conflicts of interest. 14 Complexity