# Stability, Global Dynamics, and Social Welfare of a Two-Stage Game under R&D Spillovers

Authors: Wei Zhou; Tong Chu; Xiao-Xue Wang
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

In this paper, a repeated two-stage oligopoly game where two boundedly rational firms produce homogeneous product and apply gradient adjustment mechanism to decide their individual R&D investment is considered. Results concerning the equilibrium in the built model and the stability are discussed. The effects of system parameters on the complex dynamical behaviors of the built game are analyzed. We find that the system can lose stability through a flip bifurcation or a Neimark–Sacker bifurcation. In addition, the coexistence of multiattractors is also discussed using the so-called basin of attraction. At the end of this research, the social welfare of the given duopoly game is also studied.


---

## Body

## 1. Introduction

Technological innovation is the first driving force for the development of social progress and economy. Research and development (denoted as R&D for short) activities are the main driving force for enterprises, which are the microeconomic foundation to carry out technological innovation and have become an important way for enterprises to gain greater competitive advantage. Therefore, the research about R&D has attracted more and more attention of both theoretical and applied researchers.There are many tools to analyze the R&D behavior of enterprises, among which the industrial organization theory mainly uses two or three stages of noncooperative game model to investigate R&D competition. The problems studied by industrial organization theory are more in line with the actual situation of enterprises, and the results obtained from the research of industrial organization theory also have stronger explanatory power for R&D problems.Along this line, D’Aspremont and Jacquemin [1] first put forward a famous AJ model of two-stage duopoly game, which laid a foundation for later scholars to study R&D problems by using game theory. Since then, Kamien et al. [2] have extended the AJ model to many enterprises and analyzed four different R&D forms, which are R&D competition, R&D cartel, RJV competition, and RJV cartel. Amir et al. [3] analyzed the impact of endogenous spillovers on cooperative and noncooperative game and compared them. On the research of R&D competition among enterprises, using multistage game theory can be further referenced by Qiu [4]; Brod et al. [5]; Matsumura [6]; and so on.However, the above literature usually assumes that every enterprise has complete information of profit function and complete rationality, so the established models are often static game models. However, in the real market economy, enterprises cannot grasp enough complete decision-making information. Considering that the decision-makers of enterprises are limited by objective conditions such as perception ability, the decision-making of enterprises cannot be completely rational, but they can be only boundedly rational. They often adopt a simpler dynamic adjustment process in the decision-making process. In the process of dynamic adjustment, the competition among enterprises will converge to Nash equilibrium or never converge to Nash equilibrium through repeated games, and even chaos will occur.As one of the three main research topics of nonlinear science, chaos has become an important research object of the economic system since its inception. In recent years, more and more scholars began to apply chaos theory to the economic system so as to explore the law of development in economy and reveal the complex economic phenomena [7–9]. Oligopoly game, as a very important economic model, has also become a hot issue for scholars. Firstly, some scholars prove the existence of chaos in the oligopoly game model by mathematical methods. Li et al. [10] proved the Kato chaos in the duopoly game model. Pireddu [11] assumed that the three oligarchs in the same market are heterogeneous and the existence of chaos in the built game is proved by using a topological approach. In addition, some scholars have studied the stability of equilibrium in the competitions of multi-oligarchs. For example, Elsadany [12] considered the Nash equilibrium stability of the duopoly Cournot model with relative profit maximization, where he supposed the cost function has external effects. Tramontana et al. [13] discussed the effect of increasing numbers of competitors on the stability of Cournot–Nash equilibrium and found that the Nash equilibrium would become unstable when the number of competitors increased. Matsumoto et al. [14] constructed an oligopoly model with bounded rationality and discussed the existence of the unique equilibrium state of the model with simple price and cost functions. Based on isoelastic demand function, Snyder et al. [15] established a continuous Cournot model of oligopoly competition and discussed the stability of the given model. Askar [16] established a duopoly Cournot model and analyzed the stability of duopoly competition in the case of concave demand function and uncertain cost function. In addition, more scholars have focused on the complex dynamic behavior of multi-oligarch game according to the degree of product differentiation (see [17, 18]), the difference of dynamic adjustment strategy (see [19–21]), the existence of delayed decision-making (see [22, 23]), and so on. Finally, a model with time varying delays and a multiscale time approach is proposed by Cavalli and Naimzada [24] for a monopoly and then generalized to a Cournotian competition in Cavalli et al. [25].In recent years, more and more scholars began to use the multistage dynamic game model to analyze the complexity and formation mechanism of R&D competition among enterprises in order to reveal the decision-making rules of R&D competition. Bischi and Lamantia [26] used a two-stage game model to simulate the complexity of the R&D network. Matsumura et al. [6] established a two-stage model and discussed the impact of cooperation degree between enterprises on rival profits. Shibata [27] extended Matsumura’s model and analyzed the impact of R&D spillovers on the profits of enterprises. By assuming that firms cooperate in the production stage and compete in the R&D stage, Zhang et al. [28] built a duopoly game model with semicollusion and discussed the complex dynamical behaviors in the built model. Zhou et al. [29] first established a two-stage R&D game model based on product differentiation. And then they discussed the intermittent chaos, bifurcation, and coexisting attractors using basin of attraction.However, the complex dynamical behavior of the two-stage duopoly game considering R&D spillover, where the duopoly firms produce homogeneous products and compete in both stages, has not been studied yet. The main purpose of this research is to establish a two-stage duopoly Cournot model and investigate the effects of parameters on the dynamical behaviors of the built game. In this paper, we suppose that the boundedly rational duopoly firms produce homogeneous products in a market and conduct R&D activities to reduce the production cost. Therefore, the model we built here would be a two-stage game. At the first stage, the duopoly firms compete with R&D investments, and we also allow the existence of R&D spillovers in order to better simulate the real situation. While at the second stage, we presume that the firms choose the outputs as their decision variables.The rest of this research is organized as follows: the two-stage duopoly model will be built in Section2. In Section 3, the local stability of the equilibrium points and the stability region of Nash equilibrium will be discussed. While in Section 4, a series of numerical simulations will be conducted in order to disclose the complex dynamical behaviors of the built model. The effects of parameters on the social welfare will be analyzed in Section 5. At last, this research is summarized in the final section.


## 2. The Duopoly Model

In this section, we first introduce a Cournot duopoly model. We assume that there are two firms producing homogeneous products in an oligopolistic industry, and the firms are labeled byi, (i = 1, 2) for convenience. Both firms conduct R&D to reduce their respective product costs and improve the quality of their products. Furthermore, the competition between these two firms can be simulated by a two-stage duopoly game. At the first stage, these two firms compete in R&D level in order to reduce the product cost and further maximize their respective profits. Since the limited access to market information, these two firms are considered as boundedly rational. While at the second stage, the two firms conduct quantity competition and will share their information so as to maximize their own profits after the selection of R&D investment in the first stage. Following D’Aspremont and Jacquemin [1], the market, in which the quantity-setting firms operate, can be characterized by a linear inverse demand function (the inverse demand function is the inverse function of the demand function), which can be given as(1)pqi,qj=a−bqi+qj,i=1,2,i≠j,where qi≥0 is the production output of firm i (i = 1, 2) and a,b>0 are positive constants. Parameter a represents the reservation price of produced goods, i.e., the maximum possible price that consumers are willing to pay. The production cost of firm i can be given by(2)cixi,xjqi=c−xi−βxjqi,i=1,2,i≠j,where c∈0,a is some constant, which means the marginal costs of these two firms (it is supposed that firm 1 and firm 2 have the same marginal costs in this research). xi≥0 is the R&D investment of firm i, and β∈0,1 measures the R&D spillovers between firm i and its rival. Here, β=0 means that the protection of intellectual rights is perfect. That is, firm i’s knowledge acquired through R&D can only affect its own product cost, whereas β=1 means that the two firms share their R&D achievements completely. Or it can be understood as the R&D spillover is perfect. According to D’Aspremont and Jacquemin [1], the R&D cost of firm i can be represented as(3)gxi=12γxi2,i=1,2,where γ>0 is a positive parameter, which can be used to measure the R&D efficiency of these two firms. In view of the built game is two-stage game, we know that the competitive strategies of these two firms should be composed of a level of R&D investment and a subsequent quantity competition strategies based on the R&D choice. We shall now analyze the duopoly game in which both firms do not cooperate in R&D level and output level. It is obvious that the firm i’s profit can be given as(4)πi=a−bqi+qjqi−c−xi−βxjqi−12γxi2,i=1,2,i≠j.It is clear that equation (4) can be regarded as the residue that subtracts the production cost and the R&D cost from sales volume. The marginal profit for firm i at any given point can be given by(5)∂πi∂qi=a−2bqi−bqj−c+xi+βxj,i=1,2,i≠j.It is obvious that the maximum profit of firmi can be solved by setting the derivative of πi with respected to qi to zero. That is,(6)∂π1∂q1=a−2bq1−bq2−c+x1+βx2=0,∂π2∂q2=a−2bq2−bq1−c+x2+βx1=0.Equation (6) can be regarded as a linear system of equations about unknowns q1 and q2. Through solving the above two equations simultaneously, we can obtain the optimal output of firm i, which is given by(7)qi=a−c+2−βxi+2β−1xj3b,i=1,2,i≠j.Substituting equation (7) into equation (4), then the expected profit of firm i can be written as(8)πiex=19ba−c+2−βxi+2β−1xj2−12γxi2,i=1,2,i≠j.It is worth noting that equation (8) can be regarded as a function with variables xi and xj. All the firms will maximize their own profits by choosing the proper R&D investment xi. Then, an optimization problem can be obtained. Differentiating the function πiex with respect to the variable xi, then we can get(9)∂πiex∂xi=29b2−βa−c+2−βxi+2β−1xj−γxi,i=1,2,i≠j.It is necessary to presume that the firms are boundedly rational owing to the limited market information. That is to say, all the firms can only make decisions on the basis of the market information they have. Furthermore, the process of decision-making is also dynamic for the boundedly rational firms. Namely, ifxit is used to represent the R&D investment decision of firm i at the period t, then the R&D investment of firm i at the period t + 1 (i.e., xit+1) will depend on the marginal profit of firm i at period t. According to Bischi and Naimzada [30], the dynamical adjustment mechanism employed here is xit+1=xit+αixiΦixi, i=1,2.WhereΦixi=∂πiex/∂xi is marginal relative profit of firm i, the specific form of it has been given in equation (9). And αi(i=1,2) are the adjustment speeds (or speeds of adjustment by some researchers) of firm i’s investment in R&D. This parameter reflects the response speed of firm i to its marginal profit signal. In period t, if the estimated marginal profit firm iΦixi is positive/negative, then firm i will increase/decrease its respective R&D investment at period t + 1 with a rate of αi. According the above analysis, the final dynamic game can be expressed as(10)x1t+1=x1t+α1x1t29b2−βa−c+2−βx1+2β−1x2−γx1,x2t+1=x2t+α2x2t29b2−βa−c+2−βx2+2β−1x1−γx2.


## 3. Stability Analysis for the Duopoly Model

Obviously, previous duopoly model (10) consists of two nonlinear difference equations. The dynamical behaviors of it can be very complex. If we want to study the complex dynamical behaviors of system (10), the equilibrium points (also called fixed points) should be solved in the first place. Therefore, if we set xit+1=xit in equation (10), then a set of nonlinear algebraic equations can be obtained. Through solving the nonlinear algebraic equations, four equilibrium points can be found, which are E0=0,0, E1=2a−c2−β/9bγ−22−β2,0, E2=0,2a−c2−β/9bγ−22−β2, and E∗=x1∗,x2∗, where x1∗=x2∗=x∗=2a−c2−β/9bγ−22−β1+β.It is clear that the equilibrium pointsE0, E1, and E2are on the coordinate axes of phase plane x1,x2, so they are often called boundary equilibrium points. In addition, E∗ is usually called internal equilibrium point as it is located in the interior of phase plane x1,x2. We stress that it could be shown that E∗ corresponds to the Nash equilibrium of the two-stage static game. The stability of all these equilibrium points and the complex dynamical behaviors of model (10) after these equilibrium points lose their stability are our major concern in the future discussion. Both the analytical method and numerical method would be employed in order to analyze the stability of all these equilibrium points and the complex dynamical behavior of system (10). However, the nonnegativity of all these equilibrium points should be ensured before the discussion of these issues, as negative equilibrium has no economic significance. According to economic meaning of the parameters, we know that a, b, c, γ, and β are all positive, so it can be inferred that the boundary equilibrium points E1, E2 and the unique Nash equilibrium point E∗ are positive if and only if the inequalities given below hold:(11)a>c,2>β,9bγ>22−β2,9bγ>22−β1+β.Otherwise, there will be at least one firm out of the game. In this case, the duopoly market will become a monopoly market. The following discussion will mainly focus on the local stability analysis of all the boundary equilibrium points. It is generally known that the Jacobian matrix is a primary tool to analyze the local stability of these boundary equilibrium points. For simplicity purposes, we use the symbolJx1,x2 to represent the Jacobian matrix of game (10) at any given point x1,x2, which has a specific form as(12)Jx1,x2=1+α129b2−βa−c+22−βx1+2β−1x2−2γx1,29bα1x12−β2β−1,29bα2x22−β2β−1,1+α229b2−βa−c+22−βx2+2β−1x1−2γx2.The characteristic equation ofJx1,x2 can be given as(13)Pλ=λ2−Trλ+Det=0,where Tr and Det are the trace and the determinant of Jx1,x2, respectively.According to the theory of nonlinear dynamics, we know that any fixed pointEx¯1,x¯2 will be locally asymptotically stable if all the eigenvalues λii=1,2,… of the Jacobian matrix evaluated at Ex¯1,x¯2 are located inside the unit circle of the complex plane. That is, all the eigenvalues need to satisfy λi<1i=1,2,….Proposition 1.
The trivial equilibrium pointE0 is a repelling node.Proof.
Obviously, the Jacobian matrix evaluated at the trivial equilibrium pointE0 can be rewritten as(14)JE0=1+α129b2−βa−c,0,0,1+α229b2−βa−c.
Quite evidently, the eigenvalues ofJE0 are λ1=1+α12/9b2−βa−c and λ2=1+α22/9b2−βa−c. Since all the parameters are positive constants and the conditions a>c, b>0, αi>0i=1,2, and 0<β<1 always meet according to the economic meaning, then we have λ1>1 and λ2>1. Therefore, we can easily draw a conclusion that E0 is a repelling node.In the same way, the similar results about the boundary equilibrium pointsE1 and E2 can be obtained as follows.Proposition 2.
The boundary equilibrium pointsE1 and E2 are saddle points or attracting nodes.Proof.
The Jacobian matrix evaluated at the equilibrium pointE1 is(15)JE1=1−α129b2−βa−c,49bα1a−c2−β22β−1,0,1+α229b2−βa−c9bγ−62−β1−β9bγ−22−β2.
It is clear that the Jacobian matrixJE1 given above is an upper triangular matrix. The eigenvalues of this matrix can be easily gained, which are λ1=1−α12/9b2−βa−c and λ2=1+α22/9b2−βa−c9bγ−62−β1−β/9bγ−22−β2. According to inequalities (11) and the nonnegative conditions α1>0, α2>0, we know that λ1<1 and λ2>1 if and only if the inequalities 1/2≤β≤1 and 9bγ>62−β1−β meet. While if the condition 0<β<1/2 and 22−β2<9bγ<62−β1−β meet, then the modules of the eigenvalues will be less than one, i.e., λ1<1 and λ2<1. Therefore, the fixed point E1 would be an attracting node in such a situation.
Similarly, we can also prove that the boundary equilibrium pointE2 is a saddle point or an attracting node.


### 3.1. Local Stability Analysis of the Internal Equilibrium Point

Another topic of this research is the local stability of the internal equilibrium pointE∗=x1∗,x2∗, (x1∗=x2∗=x∗). Actually, the local stability analysis of the internal equilibrium point is more complicated than that of the boundary equilibrium points. The Jacobian matrix at E∗ is(16)JE∗=1+α129b2−βa−c+3x∗−2γx∗,29bα1x∗2−β2β−1,29bα2x∗2−β2β−1,1+α229b2−βa−c+3x∗−2γx∗.The specific forms of trace and determinant ofJE∗ are given as(17)TrJE∗=2+α1+α2Δ1,DetJE∗=1+α1+α2Δ1+α1α2Δ12−Δ22,where Δ1=2/9b2−βa−c+3x∗−2γx∗ and Δ2=2/9bx∗2−β2β−1.According to the expression of Jury discriminant in stability analysis [31], we know that the internal equilibrium point is asymptotically stable if the following nonlinear inequalities hold. That is,(18)P1=1−TrJE∗+DetJE∗>0,P−1=1+TrJE∗+DetJE∗>0,DetJE∗<1.The dual inequalities of these three conditions (18) correspond to three different ways that one of the eigenvalues may cross the unit circle. If P1<0, then one of the real eigenvalues of JE∗ is larger than 1. While if P−1<0, then one of the real eigenvalues ofJE∗ will cross the unit circle of complex plane from −1,0. Finally, if Det>1, then JE∗ has a complex conjugate pair of eigenvalues lying outside the unit circle of complex plane. It can be obtained by calculation that the first condition of (18) is(19)1−TrJE∗+DetJE∗=α1α222−βa−c9b29bγ−62−β1−β9bγ−22−β1+β>0.The equation (19) can be expressed equivalently as bγ>2/32−β1−β. In addition, condition (19) would be always fulfilled according to the nonnegative conditions. It means that all the eigenvalues cannot cross the unit circle in the complex plane through the point 1,0. In other words, the pitchfork bifurcation will never happen. Nevertheless, the other two conditions of (18) are equivalent to(20)54b2+x∗α1α22−βa−c3bγ−22−β1−β−3bx∗α1+α29br−22−β2>0,(21)x∗27b22α1α22−βa−c3bγ−22−β1−β−3bα1+α29br−22−β2<0.The above two inequalities together constitute the stability region of system (10) in the parameter space. That is, the internal equilibrium point would be stable if the parameters meet these two conditions. However, we can draw a conclusion that the internal equilibrium point will lose its stability, when one (or both) of conditions (20) and (21) is violated. When one of the system parameters crosses this stability region, the system will undergo a flip bifurcation or a Neimark–Sacker bifurcation. As we will mainly discuss the influence of speed of adjustment αi on the system’s dynamical behavior in this research, the bifurcation curve of flip bifurcation is defined by(22)α1f=3bα29br−22−β2−27b29bγ−22−β1+β2−βa−cα22−βa−c3bγ−22−β1−β−3b9br−22−β2,and the Neimark–Sacker bifurcation curve is defined by α1=α1ns, where α1ns is shown as(23)α1ns=3bα29br−22−β22α22−βa−c3bγ−22−β1−β−3b9br−22−β2,where α1f and α1ns are defined by the system parameters except for α1. They give the flip bifurcation curve and the Neimark–Sacker bifurcation curve, respectively. If the values of parameters cross the bifurcation curves, then a bifurcation will arise and the stable/unstable internal equilibrium point will become unstable/stable. Based on the above discussions, we can get the following proposition.Proposition 3.
The internal equilibrium pointE∗ is a(i)
sink ifα1f<α1<α1ns(ii)
source ifα1>α1f and α1>α1ns(iii)
saddle ifα1>α1f and α1>α1ns or α1<α1f and α1>α1ns(iv)
nonhyperbolic point if eitherα1=α1f or α1=α1ns
The specific mathematical expressions of the stability region in the parameter space have been studied analytically in the above discussion. However, the stability region of the internal equilibrium point can also be analyzed numerically.In the following discussion, we will analyze the influence of each parameter on the stability region of system (10). As an example, the numerical stability region of adjustment speeds α1,α2 for the internal equilibrium point is shown in Figure 1(a), where the parameters are fixed as a=50, c=42, r=0.8, b=1, and β=0.7, which are also the reference parameters set of Figure 1. From Figure 1(a), we can see that the stability region is formed by coordinate axes and a curve, and the stability region is symmetric about the diagonal of plane α1,α2. It means that the position of parameters α1 and α2 is equivalent, which can also be discovered from system (10). Furthermore, we can also find that an increasing of α1 and/or α2 can bring the internal equilibrium point out of the stability region, and then a flip (or Neimark–Sacker) bifurcation may occur.Moreover, the influence of other parameters on the stability region is also worth studying. For example, the value of the parameterb is positively related to the size of the stability region. If we choose the values of the parameters as b=1.2, a=50, c=42, γ=0.8, and β=0.7, where only the parameter b is different compared with the reference parameters set, the stability region can be seen from Figure 1(b). Comparing with Figure 1(a), we can find that the size of stability region in Figure 1(b) has increased. Similar argument holds true if the value of the parameter c increases; when all the other parameters are fixed, the region of stability will become larger (see Figure 1(d)). In Figure 1(d), the parameters values are chosen as c=45, a=50, γ=0.8, b=1, and β=0.7. However, the value of the parameter β is negatively correlated with the stability region. If the parameters set are taken as β=0.4, a=50, c=42, b=1, and γ=0.8, we will find that the stability region will be become larger as the value of β decreases (see Figure 1(c)). While when the parameter γ decreases, the region of stability will become smaller than Figure 1(a), as shown in Figure 1(e), where the parameters set are fixed as γ=0.6, a=50, c=42, b=1, and β=0.7. If we fix the values of α1, α2, b, c, γ, and β and increase the value of parameter a, we will find that the region of stability will become smaller than Figure 1(a) (see Figure 1(f)). The values of the parameters in Figure 1(f) are given as a=55, c=42, b=1, β=0.7, and γ=0.8.Figure 1
The stability region of system (10) in the α1,α2 parameter space with different sets of parameters. (a) The parameter set is chosen as a=50, c=42, r=0.8, b=1, and β=0.7. (b) The parameter set is chosen as b=1.2, a=50, c=42, γ=0.8, and β=0.7. (c) The parameter set is chosen as β=0.4, a=50, c=42, b=1, and γ=0.8. (d) The parameter set is chosen as c=45, a=50, γ=0.8, b=1, and β=0.7. (e) The parameter set is chosen as γ=0.6, a=50, c=42, b=1, and β=0.7. (f) The parameter set is chosen as a=55, c=42, b=1, β=0.7, and γ=0.8.
(a)(b)(c)(d)(e)(f)


## 3.1. Local Stability Analysis of the Internal Equilibrium Point

Another topic of this research is the local stability of the internal equilibrium pointE∗=x1∗,x2∗, (x1∗=x2∗=x∗). Actually, the local stability analysis of the internal equilibrium point is more complicated than that of the boundary equilibrium points. The Jacobian matrix at E∗ is(16)JE∗=1+α129b2−βa−c+3x∗−2γx∗,29bα1x∗2−β2β−1,29bα2x∗2−β2β−1,1+α229b2−βa−c+3x∗−2γx∗.The specific forms of trace and determinant ofJE∗ are given as(17)TrJE∗=2+α1+α2Δ1,DetJE∗=1+α1+α2Δ1+α1α2Δ12−Δ22,where Δ1=2/9b2−βa−c+3x∗−2γx∗ and Δ2=2/9bx∗2−β2β−1.According to the expression of Jury discriminant in stability analysis [31], we know that the internal equilibrium point is asymptotically stable if the following nonlinear inequalities hold. That is,(18)P1=1−TrJE∗+DetJE∗>0,P−1=1+TrJE∗+DetJE∗>0,DetJE∗<1.The dual inequalities of these three conditions (18) correspond to three different ways that one of the eigenvalues may cross the unit circle. If P1<0, then one of the real eigenvalues of JE∗ is larger than 1. While if P−1<0, then one of the real eigenvalues ofJE∗ will cross the unit circle of complex plane from −1,0. Finally, if Det>1, then JE∗ has a complex conjugate pair of eigenvalues lying outside the unit circle of complex plane. It can be obtained by calculation that the first condition of (18) is(19)1−TrJE∗+DetJE∗=α1α222−βa−c9b29bγ−62−β1−β9bγ−22−β1+β>0.The equation (19) can be expressed equivalently as bγ>2/32−β1−β. In addition, condition (19) would be always fulfilled according to the nonnegative conditions. It means that all the eigenvalues cannot cross the unit circle in the complex plane through the point 1,0. In other words, the pitchfork bifurcation will never happen. Nevertheless, the other two conditions of (18) are equivalent to(20)54b2+x∗α1α22−βa−c3bγ−22−β1−β−3bx∗α1+α29br−22−β2>0,(21)x∗27b22α1α22−βa−c3bγ−22−β1−β−3bα1+α29br−22−β2<0.The above two inequalities together constitute the stability region of system (10) in the parameter space. That is, the internal equilibrium point would be stable if the parameters meet these two conditions. However, we can draw a conclusion that the internal equilibrium point will lose its stability, when one (or both) of conditions (20) and (21) is violated. When one of the system parameters crosses this stability region, the system will undergo a flip bifurcation or a Neimark–Sacker bifurcation. As we will mainly discuss the influence of speed of adjustment αi on the system’s dynamical behavior in this research, the bifurcation curve of flip bifurcation is defined by(22)α1f=3bα29br−22−β2−27b29bγ−22−β1+β2−βa−cα22−βa−c3bγ−22−β1−β−3b9br−22−β2,and the Neimark–Sacker bifurcation curve is defined by α1=α1ns, where α1ns is shown as(23)α1ns=3bα29br−22−β22α22−βa−c3bγ−22−β1−β−3b9br−22−β2,where α1f and α1ns are defined by the system parameters except for α1. They give the flip bifurcation curve and the Neimark–Sacker bifurcation curve, respectively. If the values of parameters cross the bifurcation curves, then a bifurcation will arise and the stable/unstable internal equilibrium point will become unstable/stable. Based on the above discussions, we can get the following proposition.Proposition 3.
The internal equilibrium pointE∗ is a(i)
sink ifα1f<α1<α1ns(ii)
source ifα1>α1f and α1>α1ns(iii)
saddle ifα1>α1f and α1>α1ns or α1<α1f and α1>α1ns(iv)
nonhyperbolic point if eitherα1=α1f or α1=α1ns
The specific mathematical expressions of the stability region in the parameter space have been studied analytically in the above discussion. However, the stability region of the internal equilibrium point can also be analyzed numerically.In the following discussion, we will analyze the influence of each parameter on the stability region of system (10). As an example, the numerical stability region of adjustment speeds α1,α2 for the internal equilibrium point is shown in Figure 1(a), where the parameters are fixed as a=50, c=42, r=0.8, b=1, and β=0.7, which are also the reference parameters set of Figure 1. From Figure 1(a), we can see that the stability region is formed by coordinate axes and a curve, and the stability region is symmetric about the diagonal of plane α1,α2. It means that the position of parameters α1 and α2 is equivalent, which can also be discovered from system (10). Furthermore, we can also find that an increasing of α1 and/or α2 can bring the internal equilibrium point out of the stability region, and then a flip (or Neimark–Sacker) bifurcation may occur.Moreover, the influence of other parameters on the stability region is also worth studying. For example, the value of the parameterb is positively related to the size of the stability region. If we choose the values of the parameters as b=1.2, a=50, c=42, γ=0.8, and β=0.7, where only the parameter b is different compared with the reference parameters set, the stability region can be seen from Figure 1(b). Comparing with Figure 1(a), we can find that the size of stability region in Figure 1(b) has increased. Similar argument holds true if the value of the parameter c increases; when all the other parameters are fixed, the region of stability will become larger (see Figure 1(d)). In Figure 1(d), the parameters values are chosen as c=45, a=50, γ=0.8, b=1, and β=0.7. However, the value of the parameter β is negatively correlated with the stability region. If the parameters set are taken as β=0.4, a=50, c=42, b=1, and γ=0.8, we will find that the stability region will be become larger as the value of β decreases (see Figure 1(c)). While when the parameter γ decreases, the region of stability will become smaller than Figure 1(a), as shown in Figure 1(e), where the parameters set are fixed as γ=0.6, a=50, c=42, b=1, and β=0.7. If we fix the values of α1, α2, b, c, γ, and β and increase the value of parameter a, we will find that the region of stability will become smaller than Figure 1(a) (see Figure 1(f)). The values of the parameters in Figure 1(f) are given as a=55, c=42, b=1, β=0.7, and γ=0.8.Figure 1
The stability region of system (10) in the α1,α2 parameter space with different sets of parameters. (a) The parameter set is chosen as a=50, c=42, r=0.8, b=1, and β=0.7. (b) The parameter set is chosen as b=1.2, a=50, c=42, γ=0.8, and β=0.7. (c) The parameter set is chosen as β=0.4, a=50, c=42, b=1, and γ=0.8. (d) The parameter set is chosen as c=45, a=50, γ=0.8, b=1, and β=0.7. (e) The parameter set is chosen as γ=0.6, a=50, c=42, b=1, and β=0.7. (f) The parameter set is chosen as a=55, c=42, b=1, β=0.7, and γ=0.8.
(a)(b)(c)(d)(e)(f)


## 4. Numerical Illustrations of Dynamical Economy Model

It is well known that many complex dynamical behaviors in nonlinear dynamical systems can only be studied by numerical methods. Therefore, a numerical analysis is employed in this section to analyze the effects of varying parameters on the stability of the Nash equilibrium point of nonlinear system (10), as well as the irregular dynamical behaviors after the Nash equilibrium point loses its stability.In the following discussion, a lot of numerical tools, such as 1-D bifurcation diagram, 2-D bifurcation diagram, phase portrait, and the largest Lyapunov exponent plot, are employed to reveal the complex dynamical behaviors of system (10). It is worth emphasizing that the 2-D bifurcation diagram is a more effective tool [32] (see Figure 2) in the numerical analysis of nonlinear dynamics. Hence, the 2-D bifurcation diagram is chosen as the principal tool in our analysis. In the next discussion, let us start with a set of parameters, which are chosen as a=50, c=42, r=0.8, b=1, and β=0.7. In fact, this set of parameters is the same with that of Figure 1(a). Based on this set of parameters, the unique Nash equilibrium point is reached after some iterations, if the values of α1 and α2 are chosen in the green area of Figure 2(a). However, an increasement of α1 and/or α2 will bring the point α1,α2 out the stability region and result in the instability of the Nash equilibrium point. As a consequence, a stable period-2 cycle arises. This phenomenon is mainly caused by the so-called flip bifurcation (see the boundary between the green area and the light blue area in Figure 2(a)). From Figure 2(a), we can also obtain the direct effect of parameter α1 and α2on the local stability of the Nash equilibrium point. In particular, a further increasement of parameter αi(i=1,2), with the other parameters are held fixed, may turn the period-2 points unstable through a flip bifurcation or a Neimark–Sacker bifurcation. This research result has very important economic significance and wide applications. This destabilizing effect owing to larger adjustment speed has attracted a lot of attention and has already been proved by Flam [33]. Figure 2(b) presents the 2-D bifurcation diagram in the parameter plane of α1,α2, in which different colors represent different periods. The dark green indicates stable period-1 state, green for stable period-2 cycles, yellow for period-4, little green for period-8, and so on. However, the states of quasiperiodic, chaos, and divergent trajectories are all expressed in dark black, as the restriction of color numbers we can used here.Figure 2
(a) The stability region of the internal equilibrium point and the period-2 cycle with the parametersa=50, c=42, r=0.8, b=1, and β=0.7. (b) The 2-D bifurcation diagram in the α1,α2 plane with the same parameter selection in Figure 2(a).
(a)(b)


### 4.1. Global Dynamical Behavior Analysis

The local stability analysis given above mainly focuses on the dynamical behavior of system (10) near the equilibrium point. It means that the initial conditions should be chosen in the neighborhood of equilibrium point, if we want to analyze the local stability of an equilibrium point. However, the initial state of a real economic system often does not belong to the neighborhood of an equilibrium point, so the necessity of a global dynamical behavior analysis is highlighted. Moreover, the global dynamical analysis of a nonlinear system can also help us to know well the long run behavior of the system [18].For the sake of better understanding of some global dynamic characteristics that may arise in the preset quantity-setting game, the 2-D bifurcation diagram is mainly employed. To clarify the problem more clearly, the parameters in Figure3(a) are chosen as the same with that of in Figure 2(a). Firstly, we analyze the dynamical behavior around the line I3, which is shown in Figure 2(a). In this region, the period-2 cycle will lose stability through a Neimark–Sacker bifurcation and the system will enter quasiperiodic state after this kind of bifurcation. However, if the parameter combination α1,α2 passes through the line I1 or I2, the period-2 cycle will develop into a state of period-4 cycle, and the system will ultimately enter into chaotic state through a flip bifurcation. In Figure 3(a), we can find that different colors are associated to different periodic cycles and the detailed corresponding relationship between colors and period number can be found from the colorbar, which is given in the rightmost of Figure 3(a). Furthermore, we can also find a lot of areas of different color overlap; that is, some parts of light green area are overlapped by the pink area. This phenomenon is due to the coexistence of several attractors, which will be studied in the next section. It should be pointed out that the color of dark black may correspond to four different states, which may be chaotic state, quasiperiodic state, large periodic state, and unfeasible trajectories (or divergent trajectories by some researchers). To distinguish the states in dark black area of Figure 3(a), the 2-D largest Lyapunov exponent diagram is employed (see Figure 3(b)). The gradient colors are shown in Figure 3(b). We can obtain from Figure 3(b) that different colors represent different values of the largest Lyapunov exponent. The dark gray, which represents the relatively small values (negative), is for the stable state. As the color gradually fades, the value of the largest Lyapunov exponent also increases little by little. However, the maximum is set to 2 because of the limitation of the algorithm, and it corresponds to the white area in Figure 3(a). In addition, it is worth pointing out that the largest Lyapunov exponent is zero at the bifurcation point. By comparing the differences between these two figures, we can find the divergent trajectories very clearly, which is shown in the white region of Figure 3(b). Furthermore, we can find that Figures 3(a) and 3(b) have similar fractal structures.Figure 3
(a) The 2-D bifurcation diagram in theα1,α2 parameter space with parameters a=50, c=42, r=0.8, b=1, and β=0.7. (b) The 2-D largest Lyapunov exponent diagram corresponding to Figure 3(a).
(a)(b)Similarly, there is also a symmetric fractal structure in Figure4(a). All the parameters of Figure 4(a) are the same with the parameters in Figure 3(a) except for β. In Figure 4(a), the value of β is set as 0.75, and it is little bigger than that of in Figure 3(a). Figure 4(b) gives the largest Lyapunov exponent diagram against the parameter α1,α2, which corresponds to the 2-D bifurcation diagram of Figure 4(a). Comparing with Figure 3(b), we can find that in addition to the reduction of the stability region of E∗, the divergent area has also increased.Figure 4
(a) The 2-D bifurcation diagram against parametersα1 and α2 with a=50, c=42, r=0.8, b=1, and β=0.75. (b) The corresponding 2-D largest Lyapunov exponent diagram.
(a)(b)Through the above discussion, we can see that when other parameters are fixed, an increase of parametersαi(i=1,2) will lead to the instability of Nash equilibrium through a flip bifurcation or a Neimark–Sacker bifurcation (see Figure 2(b)). For example, if the parameter set is chosen as a=50, c=42, r=0.8, b=1, and β=0.75 and the parameters α1 and α2 are chosen as the bifurcation parameter, we can find that the Nash equilibrium point is stable until a flip bifurcation occurs at which a set of stable period-2 points arise as the parameter α1 grows. Then, with the speed of adjustment α1 further increases, a series of period-doubling bifurcations occurs, and the cycles with higher periods and chaotic state will arise eventually. However, if the value of α2 meets 0.476<α2<0.800, the period-2 cycle can also lose its stability through a Neimark–Sacker bifurcation as the value of α1 increases.Figure5(a) gives a series of period-doubling bifurcations that cause instability, where the set of parameters are chosen as α2=0.15, a=50, c=42, r=0.8, b=1, and β=0.7. We can notice that the Nash equilibrium point is stable for 0≤α1≤0.6105, and the Nash equilibrium point will lose stability through a flip bifurcation at α1≈0.6105. Successively, a series of different state of period-2, period-4, period-8, etc., will appear constantly. At last, the system will enter into chaotic state when α1>0.8012.Figure 5
(a) The 1-D bifurcation diagram ofα1, where the reset of parameters is chosen as α2=0.15, a=50, c=42, r=0.8, b=1, and β=0.7. (b) The largest Lyapunov exponent of α1corresponding to Figure 5(a). (c) The 1-D bifurcation diagram of α1, where the reset of parameters is chosen as α2=0.35, a=50, c=42, r=0.8, b=1, and β=0.7. (d) The largest Lyapunov exponent of α1corresponding to Figure 5(c). (e) The 1-D bifurcation diagram of α1, where the reset of parameters is chosen as α2=0.50, a=50, c=42, r=0.8, b=1, and β=0.7. (f) The largest Lyapunov exponent of α1 corresponding to Figure 5(e). (g) The 1-D bifurcation diagram of α1, where the reset of parameters is chosen as α2=0.76, a=50, c=42, r=0.8, b=1, and β=0.7. (h) The largest Lyapunov exponent of α1 corresponding to Figure 5(g).
(a)(b)(c)(d)(e)(f)(g)(h)Another useful tool to study the complex dynamical behaviors of nonlinear systems is the so-called Lyapunov exponent. The Lyapunov exponent describes the convergence or divergence of adjacent orbits in nonlinear dynamical systems, which is suitable for distinguishing the regular attractors and the strange attractors. If there is no positive Lyapunov exponent, we can make sure that the attractor is a regular attractor. If there is at least one positive Lyapunov exponent, then we can induce that a strange attractor (or chaotic attractor) exists. Zero Lyapunov exponent corresponds to bifurcation point or quasiperiodic orbit, while the especially large Lyapunov exponent always implies that there is an unfeasible trajectory that would evolve to infinity. Therefore, the largest Lyapunov exponent will be employed to distinguish the regular orbits, the quasiperiodic orbits, the chaotic orbit, and the unfeasible orbits (or divergent trajectories). The largest Lyapunov exponent with parameterα1 is calculated and plotted in Figure 5(b), which also corresponds to Figure 5(a). It can be seen from Figure 5(b) that the largest Lyapunov exponent is negative when the bifurcation parameter is chosen as α1<0.6105. When α1≈0.6105, the first period-doubling bifurcation takes place. As the value of parameter α1 further increases, the occurrence of period-doubling bifurcation leads to chaos, eventually.Similarly, the bifurcation diagram of system (10) with α1 under the parameter set α2=0.35, a=50, c=42, r=0.8, b=1, and β=0.7 is shown in Figure 5(c). The largest Lyapunov exponent associated with Figure 5(c) is disposed in Figure 5(d). The bifurcation diagram and the corresponding largest Lyapunov exponents of system (10) with the parameters set α2=0.5, a=50, c=42, r=0.8, b=1, and β=0.7 are shown in Figures 5(e) and 5(f), respectively. From these two figures, we can see that the Nash equilibrium is locally stable when α1<0.4851. While if α1≈0.4581, a period-doubling bifurcation occurs. Then, the period-2 cycle loses stability through a Neimark–Sacker bifurcation at α1=0.715. Subsequently, an attracting invariant circle appears when the parameter α1 exceeds 0.715 and finally the chaotic attractor arises. The Neimark–Sacker bifurcation occurs when the modulus of a pair of complex eigenvalues crosses the unit circle. When a Neimark–Sacker bifurcation appears, the dynamics of market suddenly becomes quasiperiodic, and this is more difficult to deal with for a boundedly rational firm. The time series plot of a quasiperiodic trajectory is sometimes hardly distinguishable from a chaotic state or even a random state.According to the previous method, more complex phenomena can also be revealed. If we chose the parameters asα2=0.76, a=50, c=42, r=0.8, b=1, and β=0.7, the bifurcation diagram of system (10) with respect to α1 and the corresponding largest Lyapunov exponent can be plotted, which are shown in Figures 5(g) and 5(h), respectively. When α1=0.2628, a period-4 cycle appears and then it will lose stability through Neimark–Sacker bifurcation at α1=0.4716. In Figure 5(h), the positive largest Lyapunov exponent corresponds to chaotic state. This means that the dynamics of market becomes unstable and easily accesses to the chaotic state for large value of α1.In the following discussion, we will discuss the impact of other parameters on the complex dynamical behaviors of system (10). It can be seen from Figure 6(a) that the internal equilibrium point E∗ is stable for small value of parameter a. However, as the value of parameter a increases, the internal equilibrium point becomes unstable, suddenly. Gradually, the irregular dynamical behaviors arise, including high-periodic cycles, chaotic states, and so on. In Figure 6(b), the bifurcation diagram about parameter b is plotted. We can see that there is a mixed bifurcation process with period-doubling bifurcation and inverse Neimark–Sacker bifurcation merged together. In this process of bifurcation, there is also a phenomenon called “chaotic bubbles.” The “chaotic bubbles” begin in a multiperiodic state and end in a multiperiodic state, but there exists chaotic state between these two multiperiodic states. As the parameter value of b goes further, there exists a period-halving bifurcation. So, we can conclude that it is an effective way to make system (10) more stable by improving the product differentiation. Or we can say that the greater the degree of product differentiation, the more stable the market is.Figure 6
The bifurcation diagrams and the largest Lyapunov exponent diagram of system (10) with fixed values of α1 and α2, where α1=0.63 and α2=0.76. (a) The bifurcation diagram of parametera, where the rest parameters are b=1, c=42, β=0.7, and γ=0.8. (b) The bifurcation diagram of parameterb, where the rest parameters area=50, c=42, β=0.7, andγ=0.8. (c) The bifurcation diagram of parameter c, where the rest parameters are a=50, b=1, β=0.7, and γ=0.8. (d) The bifurcation diagram of parameter γ, and the rest parameters area=50, b=1, c=42, and β=0.7. (e) The bifurcation diagram of parameter β, where the rest parameters are a=50, b=1, c=42, and γ=0.8. (f) The largest Lyapunov exponents of β corresponding to Figure 6(e).
(a)(b)(c)(d)(e)(f)Similarly, from Figures6(c) and 6(d), we can see that the system is chaotic if the marginal cost c or the efficient measure of R&D investment cost γ is small enough. As c or γ increases, there exist mixed bifurcation processes, too. The bifurcation diagram of investment spillovers parameter β about the R&D investment is even more complicated (see Figure 6(e)). We can see from Figure 6(e) that the stable period-2 cycle loses its stability via a series of period-doubling bifurcations and Neimark–Sacker bifurcation, successively. Irregular states appear ultimately. Moreover, we can also observe “chaotic bubbles” in this bifurcation process. These numerical results fit very well with Proposition 3. The largest Lyapunov exponent plot corresponding to Figure 6(e) is displayed in Figure 6(f). In addition, a lot of irregular points can be found from Figure 6(e). These irregular points correspond to the coexistence of multiple attractors, which is the key issue of our research and will be discussed in the coming discussion.


### 4.2. Basins of Attraction and Multistability

A dramatic change in the long run dynamics occurs if we take different parameters and chooseα1 as the bifurcation parameter. By assuming α2=0.76, a=50, c=42, r=0.8, b=1, and β=0.7, some interesting events take place when the bifurcation parameter α1 belongs to the interval 0.2628,0.6491 (see Figure 5(g)). We have shown in Figure 5(g) that a supercritical Neimark–Sacker bifurcation happens after the period-2 stable cycle loses its stability. Then, the annular attractor, which consists of two closed curves, arises in the phase space, as shown in Figure 7(a). If the parameter α1 further increases, different phase-locking situations will take place and the periodic-7 cycle undergoes a period-doubling sequence that gives rise to a periodic-14 attractor as shown in Figure 7(b). Then, a homoclinic bifurcation of the two period-14 saddle cycle happens, and it leads to the emergence of a “cyclical chaotic attractor,” see Figure 7(c). The enlargement of one piece of the attractor allows us to appreciate that each of the closed curve exhibits loops and self-intersections. It means that the cyclical chaotic attractor is made up of two weakly chaotic rings, which merge into a unique annular chaotic attractor, when α1 is further increased and the homoclinic bifurcation arises again, as shown in Figure 7(d). In addition, the basins of attraction of these attractors are also given in Figure 7. In Figure 7, the white region represents the attracting domain, while the red region indicates the divergent area. We can find that the shape of basin of attraction has scarcely changed in the process of variation of α1. However, the attractors are getting bigger and bigger. In Figure 7(d), the attractor has almost connected with the boundary of the basin of attraction. Then, a global bifurcation called “contact bifurcation” will happen if the value of α1 further increases. The chaotic attractor will be destroyed and the so-called “ghost” will occupy the entire basin of attraction after the “contact bifurcation” happens.Figure 7
The basins of attraction and the attractors with varyingα1, where the rest parameters are chosen as α2=0.76, a=50, c=42, r=0.8, b=1, and β=0.7. In these figures, the red region is the divergent area and the white region is the basin of attractors. (a) α1=0.48, the period-2 limit cycles (two attracting closed orbits) arise after the Neimark–Sacker bifurcation. (b) α1=0.49, the phenomenon of phase-locking appears, where we can find a period-14 cycle. (c) α1=0.50, an attractor made up by 2 weakly chaotic rings comes up. (d) A chaotic attractor and its basin of attraction when α1=0.63.
(a)(b)(c)(d)Another interesting phenomenon may arise in system (10) is the coexistence of multiple attractors. When the parameter α1 is chosen as the bifurcation parameter, and the values of other parameters are fixed as α2=0.8, a=50, c=42, r=0.8, b=1, and β=0.7, the basins of attraction at this circumstances can be plotted numerically, see Figure 8. In these two pictures, we can find the emergence of very complex basins of attraction with fractal boundaries. The value of parameter α1 in Figure 8(a) is chosen as α1=0.4952. From Figure 8(a), we can see that there is a period-2 cycle coexisting with a chaotic attractor, where the purple region is the attracting area of the period-2, while the yellow region is the attracting area of the chaotic attractor, and the red region is the divergent area. However, the purple region and the yellow region are intertwined, and it is difficult to separate them. It means that the initial conditions have a great influence on the final state of system (10). This fact also tells us to be very cautious in choosing the initial state of the system. As the value of α1 further increases to α1=0.4996, the yellow region is replaced by numerous yellow spots, and this is due to the contact between the chaotic attractor and the boundary of its basin of attraction. Hence, only the “ghost” of the chaotic attractor survives in the yellow region, see Figure 8(b). If the value of α1 increases further, we will find that all the yellow spots will eventually disappear, and these yellow spots will also be replaced by purple. This means that the coexistence of multiattractors has also disappeared.Figure 8
The basins of attraction with varyingα1 in the phase plane with the rest parameters are chosen as α2=0.8, a=50, c=42, r=0.8, b=1, and β=0.7. (a) α1=0.4952, a period-2 cycles coexist with a chaotic attractors. (b)α1=0.4996, the basin becomes very complex after the disappearance of the chaotic attractor.
(a)(b)


## 4.1. Global Dynamical Behavior Analysis

The local stability analysis given above mainly focuses on the dynamical behavior of system (10) near the equilibrium point. It means that the initial conditions should be chosen in the neighborhood of equilibrium point, if we want to analyze the local stability of an equilibrium point. However, the initial state of a real economic system often does not belong to the neighborhood of an equilibrium point, so the necessity of a global dynamical behavior analysis is highlighted. Moreover, the global dynamical analysis of a nonlinear system can also help us to know well the long run behavior of the system [18].For the sake of better understanding of some global dynamic characteristics that may arise in the preset quantity-setting game, the 2-D bifurcation diagram is mainly employed. To clarify the problem more clearly, the parameters in Figure3(a) are chosen as the same with that of in Figure 2(a). Firstly, we analyze the dynamical behavior around the line I3, which is shown in Figure 2(a). In this region, the period-2 cycle will lose stability through a Neimark–Sacker bifurcation and the system will enter quasiperiodic state after this kind of bifurcation. However, if the parameter combination α1,α2 passes through the line I1 or I2, the period-2 cycle will develop into a state of period-4 cycle, and the system will ultimately enter into chaotic state through a flip bifurcation. In Figure 3(a), we can find that different colors are associated to different periodic cycles and the detailed corresponding relationship between colors and period number can be found from the colorbar, which is given in the rightmost of Figure 3(a). Furthermore, we can also find a lot of areas of different color overlap; that is, some parts of light green area are overlapped by the pink area. This phenomenon is due to the coexistence of several attractors, which will be studied in the next section. It should be pointed out that the color of dark black may correspond to four different states, which may be chaotic state, quasiperiodic state, large periodic state, and unfeasible trajectories (or divergent trajectories by some researchers). To distinguish the states in dark black area of Figure 3(a), the 2-D largest Lyapunov exponent diagram is employed (see Figure 3(b)). The gradient colors are shown in Figure 3(b). We can obtain from Figure 3(b) that different colors represent different values of the largest Lyapunov exponent. The dark gray, which represents the relatively small values (negative), is for the stable state. As the color gradually fades, the value of the largest Lyapunov exponent also increases little by little. However, the maximum is set to 2 because of the limitation of the algorithm, and it corresponds to the white area in Figure 3(a). In addition, it is worth pointing out that the largest Lyapunov exponent is zero at the bifurcation point. By comparing the differences between these two figures, we can find the divergent trajectories very clearly, which is shown in the white region of Figure 3(b). Furthermore, we can find that Figures 3(a) and 3(b) have similar fractal structures.Figure 3
(a) The 2-D bifurcation diagram in theα1,α2 parameter space with parameters a=50, c=42, r=0.8, b=1, and β=0.7. (b) The 2-D largest Lyapunov exponent diagram corresponding to Figure 3(a).
(a)(b)Similarly, there is also a symmetric fractal structure in Figure4(a). All the parameters of Figure 4(a) are the same with the parameters in Figure 3(a) except for β. In Figure 4(a), the value of β is set as 0.75, and it is little bigger than that of in Figure 3(a). Figure 4(b) gives the largest Lyapunov exponent diagram against the parameter α1,α2, which corresponds to the 2-D bifurcation diagram of Figure 4(a). Comparing with Figure 3(b), we can find that in addition to the reduction of the stability region of E∗, the divergent area has also increased.Figure 4
(a) The 2-D bifurcation diagram against parametersα1 and α2 with a=50, c=42, r=0.8, b=1, and β=0.75. (b) The corresponding 2-D largest Lyapunov exponent diagram.
(a)(b)Through the above discussion, we can see that when other parameters are fixed, an increase of parametersαi(i=1,2) will lead to the instability of Nash equilibrium through a flip bifurcation or a Neimark–Sacker bifurcation (see Figure 2(b)). For example, if the parameter set is chosen as a=50, c=42, r=0.8, b=1, and β=0.75 and the parameters α1 and α2 are chosen as the bifurcation parameter, we can find that the Nash equilibrium point is stable until a flip bifurcation occurs at which a set of stable period-2 points arise as the parameter α1 grows. Then, with the speed of adjustment α1 further increases, a series of period-doubling bifurcations occurs, and the cycles with higher periods and chaotic state will arise eventually. However, if the value of α2 meets 0.476<α2<0.800, the period-2 cycle can also lose its stability through a Neimark–Sacker bifurcation as the value of α1 increases.Figure5(a) gives a series of period-doubling bifurcations that cause instability, where the set of parameters are chosen as α2=0.15, a=50, c=42, r=0.8, b=1, and β=0.7. We can notice that the Nash equilibrium point is stable for 0≤α1≤0.6105, and the Nash equilibrium point will lose stability through a flip bifurcation at α1≈0.6105. Successively, a series of different state of period-2, period-4, period-8, etc., will appear constantly. At last, the system will enter into chaotic state when α1>0.8012.Figure 5
(a) The 1-D bifurcation diagram ofα1, where the reset of parameters is chosen as α2=0.15, a=50, c=42, r=0.8, b=1, and β=0.7. (b) The largest Lyapunov exponent of α1corresponding to Figure 5(a). (c) The 1-D bifurcation diagram of α1, where the reset of parameters is chosen as α2=0.35, a=50, c=42, r=0.8, b=1, and β=0.7. (d) The largest Lyapunov exponent of α1corresponding to Figure 5(c). (e) The 1-D bifurcation diagram of α1, where the reset of parameters is chosen as α2=0.50, a=50, c=42, r=0.8, b=1, and β=0.7. (f) The largest Lyapunov exponent of α1 corresponding to Figure 5(e). (g) The 1-D bifurcation diagram of α1, where the reset of parameters is chosen as α2=0.76, a=50, c=42, r=0.8, b=1, and β=0.7. (h) The largest Lyapunov exponent of α1 corresponding to Figure 5(g).
(a)(b)(c)(d)(e)(f)(g)(h)Another useful tool to study the complex dynamical behaviors of nonlinear systems is the so-called Lyapunov exponent. The Lyapunov exponent describes the convergence or divergence of adjacent orbits in nonlinear dynamical systems, which is suitable for distinguishing the regular attractors and the strange attractors. If there is no positive Lyapunov exponent, we can make sure that the attractor is a regular attractor. If there is at least one positive Lyapunov exponent, then we can induce that a strange attractor (or chaotic attractor) exists. Zero Lyapunov exponent corresponds to bifurcation point or quasiperiodic orbit, while the especially large Lyapunov exponent always implies that there is an unfeasible trajectory that would evolve to infinity. Therefore, the largest Lyapunov exponent will be employed to distinguish the regular orbits, the quasiperiodic orbits, the chaotic orbit, and the unfeasible orbits (or divergent trajectories). The largest Lyapunov exponent with parameterα1 is calculated and plotted in Figure 5(b), which also corresponds to Figure 5(a). It can be seen from Figure 5(b) that the largest Lyapunov exponent is negative when the bifurcation parameter is chosen as α1<0.6105. When α1≈0.6105, the first period-doubling bifurcation takes place. As the value of parameter α1 further increases, the occurrence of period-doubling bifurcation leads to chaos, eventually.Similarly, the bifurcation diagram of system (10) with α1 under the parameter set α2=0.35, a=50, c=42, r=0.8, b=1, and β=0.7 is shown in Figure 5(c). The largest Lyapunov exponent associated with Figure 5(c) is disposed in Figure 5(d). The bifurcation diagram and the corresponding largest Lyapunov exponents of system (10) with the parameters set α2=0.5, a=50, c=42, r=0.8, b=1, and β=0.7 are shown in Figures 5(e) and 5(f), respectively. From these two figures, we can see that the Nash equilibrium is locally stable when α1<0.4851. While if α1≈0.4581, a period-doubling bifurcation occurs. Then, the period-2 cycle loses stability through a Neimark–Sacker bifurcation at α1=0.715. Subsequently, an attracting invariant circle appears when the parameter α1 exceeds 0.715 and finally the chaotic attractor arises. The Neimark–Sacker bifurcation occurs when the modulus of a pair of complex eigenvalues crosses the unit circle. When a Neimark–Sacker bifurcation appears, the dynamics of market suddenly becomes quasiperiodic, and this is more difficult to deal with for a boundedly rational firm. The time series plot of a quasiperiodic trajectory is sometimes hardly distinguishable from a chaotic state or even a random state.According to the previous method, more complex phenomena can also be revealed. If we chose the parameters asα2=0.76, a=50, c=42, r=0.8, b=1, and β=0.7, the bifurcation diagram of system (10) with respect to α1 and the corresponding largest Lyapunov exponent can be plotted, which are shown in Figures 5(g) and 5(h), respectively. When α1=0.2628, a period-4 cycle appears and then it will lose stability through Neimark–Sacker bifurcation at α1=0.4716. In Figure 5(h), the positive largest Lyapunov exponent corresponds to chaotic state. This means that the dynamics of market becomes unstable and easily accesses to the chaotic state for large value of α1.In the following discussion, we will discuss the impact of other parameters on the complex dynamical behaviors of system (10). It can be seen from Figure 6(a) that the internal equilibrium point E∗ is stable for small value of parameter a. However, as the value of parameter a increases, the internal equilibrium point becomes unstable, suddenly. Gradually, the irregular dynamical behaviors arise, including high-periodic cycles, chaotic states, and so on. In Figure 6(b), the bifurcation diagram about parameter b is plotted. We can see that there is a mixed bifurcation process with period-doubling bifurcation and inverse Neimark–Sacker bifurcation merged together. In this process of bifurcation, there is also a phenomenon called “chaotic bubbles.” The “chaotic bubbles” begin in a multiperiodic state and end in a multiperiodic state, but there exists chaotic state between these two multiperiodic states. As the parameter value of b goes further, there exists a period-halving bifurcation. So, we can conclude that it is an effective way to make system (10) more stable by improving the product differentiation. Or we can say that the greater the degree of product differentiation, the more stable the market is.Figure 6
The bifurcation diagrams and the largest Lyapunov exponent diagram of system (10) with fixed values of α1 and α2, where α1=0.63 and α2=0.76. (a) The bifurcation diagram of parametera, where the rest parameters are b=1, c=42, β=0.7, and γ=0.8. (b) The bifurcation diagram of parameterb, where the rest parameters area=50, c=42, β=0.7, andγ=0.8. (c) The bifurcation diagram of parameter c, where the rest parameters are a=50, b=1, β=0.7, and γ=0.8. (d) The bifurcation diagram of parameter γ, and the rest parameters area=50, b=1, c=42, and β=0.7. (e) The bifurcation diagram of parameter β, where the rest parameters are a=50, b=1, c=42, and γ=0.8. (f) The largest Lyapunov exponents of β corresponding to Figure 6(e).
(a)(b)(c)(d)(e)(f)Similarly, from Figures6(c) and 6(d), we can see that the system is chaotic if the marginal cost c or the efficient measure of R&D investment cost γ is small enough. As c or γ increases, there exist mixed bifurcation processes, too. The bifurcation diagram of investment spillovers parameter β about the R&D investment is even more complicated (see Figure 6(e)). We can see from Figure 6(e) that the stable period-2 cycle loses its stability via a series of period-doubling bifurcations and Neimark–Sacker bifurcation, successively. Irregular states appear ultimately. Moreover, we can also observe “chaotic bubbles” in this bifurcation process. These numerical results fit very well with Proposition 3. The largest Lyapunov exponent plot corresponding to Figure 6(e) is displayed in Figure 6(f). In addition, a lot of irregular points can be found from Figure 6(e). These irregular points correspond to the coexistence of multiple attractors, which is the key issue of our research and will be discussed in the coming discussion.


## 4.2. Basins of Attraction and Multistability

A dramatic change in the long run dynamics occurs if we take different parameters and chooseα1 as the bifurcation parameter. By assuming α2=0.76, a=50, c=42, r=0.8, b=1, and β=0.7, some interesting events take place when the bifurcation parameter α1 belongs to the interval 0.2628,0.6491 (see Figure 5(g)). We have shown in Figure 5(g) that a supercritical Neimark–Sacker bifurcation happens after the period-2 stable cycle loses its stability. Then, the annular attractor, which consists of two closed curves, arises in the phase space, as shown in Figure 7(a). If the parameter α1 further increases, different phase-locking situations will take place and the periodic-7 cycle undergoes a period-doubling sequence that gives rise to a periodic-14 attractor as shown in Figure 7(b). Then, a homoclinic bifurcation of the two period-14 saddle cycle happens, and it leads to the emergence of a “cyclical chaotic attractor,” see Figure 7(c). The enlargement of one piece of the attractor allows us to appreciate that each of the closed curve exhibits loops and self-intersections. It means that the cyclical chaotic attractor is made up of two weakly chaotic rings, which merge into a unique annular chaotic attractor, when α1 is further increased and the homoclinic bifurcation arises again, as shown in Figure 7(d). In addition, the basins of attraction of these attractors are also given in Figure 7. In Figure 7, the white region represents the attracting domain, while the red region indicates the divergent area. We can find that the shape of basin of attraction has scarcely changed in the process of variation of α1. However, the attractors are getting bigger and bigger. In Figure 7(d), the attractor has almost connected with the boundary of the basin of attraction. Then, a global bifurcation called “contact bifurcation” will happen if the value of α1 further increases. The chaotic attractor will be destroyed and the so-called “ghost” will occupy the entire basin of attraction after the “contact bifurcation” happens.Figure 7
The basins of attraction and the attractors with varyingα1, where the rest parameters are chosen as α2=0.76, a=50, c=42, r=0.8, b=1, and β=0.7. In these figures, the red region is the divergent area and the white region is the basin of attractors. (a) α1=0.48, the period-2 limit cycles (two attracting closed orbits) arise after the Neimark–Sacker bifurcation. (b) α1=0.49, the phenomenon of phase-locking appears, where we can find a period-14 cycle. (c) α1=0.50, an attractor made up by 2 weakly chaotic rings comes up. (d) A chaotic attractor and its basin of attraction when α1=0.63.
(a)(b)(c)(d)Another interesting phenomenon may arise in system (10) is the coexistence of multiple attractors. When the parameter α1 is chosen as the bifurcation parameter, and the values of other parameters are fixed as α2=0.8, a=50, c=42, r=0.8, b=1, and β=0.7, the basins of attraction at this circumstances can be plotted numerically, see Figure 8. In these two pictures, we can find the emergence of very complex basins of attraction with fractal boundaries. The value of parameter α1 in Figure 8(a) is chosen as α1=0.4952. From Figure 8(a), we can see that there is a period-2 cycle coexisting with a chaotic attractor, where the purple region is the attracting area of the period-2, while the yellow region is the attracting area of the chaotic attractor, and the red region is the divergent area. However, the purple region and the yellow region are intertwined, and it is difficult to separate them. It means that the initial conditions have a great influence on the final state of system (10). This fact also tells us to be very cautious in choosing the initial state of the system. As the value of α1 further increases to α1=0.4996, the yellow region is replaced by numerous yellow spots, and this is due to the contact between the chaotic attractor and the boundary of its basin of attraction. Hence, only the “ghost” of the chaotic attractor survives in the yellow region, see Figure 8(b). If the value of α1 increases further, we will find that all the yellow spots will eventually disappear, and these yellow spots will also be replaced by purple. This means that the coexistence of multiattractors has also disappeared.Figure 8
The basins of attraction with varyingα1 in the phase plane with the rest parameters are chosen as α2=0.8, a=50, c=42, r=0.8, b=1, and β=0.7. (a) α1=0.4952, a period-2 cycles coexist with a chaotic attractors. (b)α1=0.4996, the basin becomes very complex after the disappearance of the chaotic attractor.
(a)(b)


## 5. Social Welfare

Based on the numerical simulations of R&D investmentxi given above, we can further analyze the dynamic evolution mechanism of the social welfare in built system (10). From equation (5), we can get that the maximum outputs of firm 1 and firm 2, which are given as(24)q1=a−c+2−βx1+2β−1x23b,q2=a−c+2−βx2+2β−1x13b.From equation (6), we can get that the maximum profits of all the two firms, which are(25)π1=19ba−c+2−βx1+2β−1x22−12γx12,π2=19ba−c+2−βx2+2β−1x12−12γx22.The social welfare (denoted by SW) is the total surplus of the market, including the producer’s surplus (for short PS) and the consumer’s surplus (abbreviation: CS) [34]. That is,(26)SW=PS+CS,with PS=π1+π2,CS=UQ−pQ.According to Qiu [4], we know that the utility function of consumers can be represented as UQ=aQ−b/2Q2, where Q=q1+q2 is the total outputs of these two firms. By substituting this relation into (26), we can get the corresponding mathematic expression of the consumer’s surplus, which is given as(27)CS=aQ−b2Q2−a−bQQ=b2Q2.Correspondingly, the social welfare can be rewritten as(28)SW=π1+π2+b2Q2,where the forms of π1 and π2 can be founded in equation (8). On the basis of the above formulas, the final form of the social welfare can be represented as(29)SW=19ba−c+2−βx1+2β−1x22+19ba−c+2−βx2+2β−1x12−12γx12+x22+2a−c+1+βx1+x226b.According to system (10) and equation (29), we can obtain the bifurcation diagram of social welfare, numerically. In Figure 9, both the bifurcation diagram of social welfare and the bifurcation diagram of R&D investment are plotted under the parameter set α2=0.35, a=50, c=42, r=0.8, b=1, and β=0.7, whereas the adjusting speed of firm 1 α1 is chosen as the bifurcation parameter.Figure 9
(a) The bifurcation diagrams of the social welfare with parameterα1, where the parameters are chosen as α2=0.35, a=50, c=42, r=0.8, b=1, and β=0.7. (b) The bifurcation diagram of R&D investments xi(i=1,2) corresponding to Figure 9(a). (c) The average social welfare with parameter α1 and the rest parameters are the same with Figure 9(a). (d) The average profits of firm 1 and firm 2 with parameter α1 and the rest parameters are the same with Figure 9(a). (e) The bifurcation diagrams of the social welfare with parameter α1, where the parameters are chosen as α2=0.55, a=50, c=42, r=0.8, b=1, and β=0.7. (f) The bifurcation diagram of R&D investments xi(i=1,2) corresponding to Figure 9(e). (g) The average social welfare with parameter α1 and the rest parameters are the same with Figure 9(e). (h) The average profits of firm 1 and firm 2 with parameter α1 and the rest parameters are the same with Figure 9(e).
(a)(b)(c)(d)(e)(f)(g)(h)From Figure9(a), we can see that the social welfare is stable when 0≤α1≤0.5715. And then the balanced social welfare loses its stability via a period-doubling bifurcation when the parameter value meets α1≤0.5715. Then, a series of period-doubling bifurcations occur. The period-2 cycle, period-4 cycle, period-8 cycle, and the chaotic attractors appear successively. And then, another period-doubling bifurcation appears when α1>0.8515. As the value of α1 increases further, the chaotic state occurs again. From the above analysis, we can see that both social welfare and R&D investment of these two firms will converge to an equilibrium, when the value of α1 is small enough. By comparing Figures 9(a) and 9(b), we can find that these two figures have same bifurcation structure. Moreover, when the system is in chaotic state, the level of social welfare may be higher than that of periodic state in some certain time periods. However, it does not mean that chaos is conducive to improving the level of social welfare. If we calculate the average social welfare and the average profits (see Figures 9(c) and 9(d)), we will find that both the level of average social welfare and the level of average profit would be reduced, when the system is in a chaotic state.Similarly, if we choose the parameter setα2=0.55, a=50, c=42, r=0.8, b=1, and β=0.7, the bifurcation diagram of social welfare is shown in Figure 9(e). The corresponding bifurcation diagram of R&D investment is given in Figure 9(f). We can observe from these two figures that the equilibrium point will lose stability through flip bifurcation and Neimark–Sacker bifurcation successively, as the value of α1 increases. The varying curves of average social welfare and average profit with parameter α1 are also plotted in Figures 9(g) and 9(h), respectively. However, unlike Figures 9(c) and 9(d), the quasiperiodic state arises at this set of parameters. We can find both the average social welfare and the average profits have increased to a certain extent, when the system is in a quasiperiodic state.


## 6. Conclusion

In this research, a two-stage Cournot duopoly game, where the duopoly firms produce homogeneous goods and conduct R&D to reduce the production costs, is built. In the built model, we also allow the existence of R&D spillovers in order to better simulate the real situation. And then, the existence and the stability of the fixed points are discussed, and the effects of the system parameters on the size of stability region are discussed. We find that the adjusting speeds of the duopoly firms have a negative effect on the stability region. It means that the larger adjusting speed will lead the loss of stability of the Nash equilibrium point. We also find that the system will become more stable, if the R&D efficiency is improved. However, the higher degree of R&D spillover will make the system more unstable. It suggests that the firms should pay attention to the protection of intellectual property rights in the process of R&D.In addition, the complex behaviors of the given model are also investigated numerically. The numerical tools, such as 2-D bifurcation diagram, 2-D largest Lyapunov exponent, and basins of attraction, are employed to study the complex dynamical behaviors. We find that the given system can present very complex dynamical behaviors. The equilibrium may lose stability through a Neimark–Sacker bifurcation or a flip bifurcation. Through analyzing the bifurcation diagram, the phenomenon of “chaotic bubble” is found in our research. The system will fall into a multiperiodic state, or a quasiperiodic state, or a chaotic state, or even divergent state eventually. The four different states can be distinguished using the 2-D bifurcation diagram combined with the 2-D largest Lyapunov exponent diagram. Moreover, the coexistence of multiattractors is also researched in this study using the so-called basin of attraction. We find the basins of the coexisting attractors can be very complex and present fractal structure. The effects of system parameters on the social welfare of the given system are discussed in the end of this paper. We find that chaos is harmful to the economic system, and it can lead to a decline in social welfare.



---

*Source: 2096868-2021-01-18.xml*