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

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.


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. erefore, the research about R&D has attracted more and more attention of both theoretical and applied researchers.
ere 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. e 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. ey 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][8][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][20][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 twostage 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. e 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. erefore, 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. e rest of this research is organized as follows: the twostage duopoly model will be built in Section 2. 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. e effects of parameters on the social welfare will be analyzed in Section 5. At last, this research is summarized in the final section.

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 by i, (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 quantitysetting 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 where q i ≥ 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. e production cost of firm i can be given by 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). x i ≥ 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. at 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 where c > 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 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. e marginal profit for firm i at any given point can be given by It is obvious that the maximum profit of firm i can be solved by setting the derivative of π i with respected to q i to zero. at is, Equation (6) can be regarded as a linear system of equations about unknowns q 1 and q 2 . rough solving the above two equations simultaneously, we can obtain the optimal output of firm i, which is given by Substituting equation (7) into equation (4), then the expected profit of firm i can be written as It is worth noting that equation (8) can be regarded as a function with variables x i and x j . All the firms will maximize their own profits by choosing the proper R&D investment x i . en, an optimization problem can be obtained. Differentiating the function π (ex) i with respect to the variable x i , then we can get It is necessary to presume that the firms are boundedly rational owing to the limited market information. at is to say, all the firms can only make decisions on the basis of the market information they have. Furthermore, the process of decisionmaking is also dynamic for the boundedly rational firms. Namely, if x i (t) 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., x i (t + 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 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.
is parameter reflects the response speed of firm i to its marginal profit signal. In period t, if the estimated marginal profit firm i Φ i (x i ) 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

Stability Analysis for the Duopoly Model
Obviously, previous duopoly model (10) consists of two nonlinear difference equations. e 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. erefore, if we set x i (t + 1) � x i (t) in equation (10), then a set of nonlinear algebraic equations can be obtained. rough solving the nonlinear algebraic equations, four equilibrium points can be found, which are E 0 � (0, 0), It is clear that the equilibrium points E 0 , E 1 , and E 2 are on the coordinate axes of phase plane (x 1 , x 2 ), 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 (x 1 , x 2 ). We stress that it could be shown that E * corresponds to the Nash equilibrium of the two-stage static game. e 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, c, and β are all positive, so it can be inferred that the boundary equilibrium points E 1 , E 2 and the unique Nash equilibrium point E * are positive if and only if the inequalities given below hold: Otherwise, there will be at least one firm out of the game. In this case, the duopoly market will become a monopoly market. e 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 symbol J(x 1 , x 2 ) to represent the Jacobian matrix of game (10) at any given point (x 1 , x 2 ), which has a specific form as e characteristic equation of J(x 1 , x 2 ) can be given as where Tr and Det are the trace and the determinant of J(x 1 , x 2 ), respectively. According to the theory of nonlinear dynamics, we know that any fixed point E(x 1 , x 2 ) will be locally asymptotically stable if all the eigenvalues λ i (i � 1, 2, . . .) of the Jacobian matrix evaluated at E(x 1 , x 2 ) are located inside the unit circle of the complex plane. at is, all the eigenvalues need to satisfy |λ i | < 1 (i � 1, 2, . . .).

Proposition 1.
e trivial equilibrium point E 0 is a repelling node.
Proof. Obviously, the Jacobian matrix evaluated at the trivial equilibrium point E 0 can be rewritten as Quite evidently, the eigenvalues of J(E 0 ) are Since all the parameters are positive constants and the conditions a > c, b > 0, α i > 0 (i � 1, 2), and 0 < β < 1 always meet according to the economic meaning, then we have λ 1 > 1 and λ 2 > 1. erefore, we can easily draw a conclusion that E 0 is a repelling node.

□
In the same way, the similar results about the boundary equilibrium points E 1 and E 2 can be obtained as follows.

Proposition 2.
e boundary equilibrium points E 1 and E 2 are saddle points or attracting nodes. Proof.
e Jacobian matrix evaluated at the equilibrium point E 1 is It is clear that the Jacobian matrix J(E 1 ) given above is an upper triangular matrix. e eigenvalues of this matrix can be easily gained, which are ). According to inequalities (11) and the nonnegative conditions α 1 > 0, α 2 > 0, we know that then the modules of the eigenvalues will be less than one, i.e., |λ 1 | < 1 and |λ 2 | < 1. erefore, the fixed point E 1 would be an attracting node in such a situation.
Similarly, we can also prove that the boundary equilibrium point E 2 is a saddle point or an attracting node. □

Local Stability Analysis of the Internal Equilibrium Point.
Another topic of this research is the local stability of the internal equilibrium point . Actually, the local stability analysis of the internal equilibrium point is more complicated than that of the boundary equilibrium points. e Jacobian matrix at E * is e specific forms of trace and determinant of J(E * ) are given as where 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. at is, e dual inequalities of these three conditions (18) correspond to three different ways that one of the eigenvalues may cross the unit circle. If P(1) < 0, then one of the real eigenvalues of J(E * ) is larger than 1. While if P(− 1) < 0, then one of the real eigenvalues ofJ(E * ) will cross the unit circle of complex plane from (− 1, 0). Finally, if |Det| > 1, then J(E * ) 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 Mathematical Problems in Engineering e equation (19) can be expressed equivalently as bc > (2/3)(2 − β)(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 x * e above two inequalities together constitute the stability region of system (10) in the parameter space. at 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 and the Neimark-Sacker bifurcation curve is defined by α 1 � α ns 1 , where α ns 1 is shown as where α f 1 and α ns 1 are defined by the system parameters except for α 1 . ey 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.
e 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 parameter b is positively related to the size of the stability region. If we choose the values of the parameters as 6 Mathematical Problems in Engineering b � 1.2, a � 50, c � 42, c � 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, c � 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 c � 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 c 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 c � 0.6, a � 50, c � 42, b � 1, and β � 0.7. If we fix the values of α 1 , α 2 , b, c, 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)). e values of the parameters in Figure 1(f) are given as a � 55, c � 42, b � 1, β � 0.7, and c � 0.8.

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. erefore, 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. is 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 α 2 on 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. is research result has very important economic significance and wide applications.
is 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. e 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.

Global Dynamical Behavior Analysis.
e 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 quantitysetting game, the 2-D bifurcation diagram is mainly employed. To clarify the problem more clearly, the parameters in Figure 3(a) are chosen as the same with that of in Figure 2(a). Firstly, we analyze the dynamical behavior around the line I 3 , 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 I 1 or I 2 , 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. is 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)). e 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. e 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.
Similarly, there is also a symmetric fractal structure in Figure 4(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.
rough 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. en, 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. Figure 5(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.
Another useful tool to study the complex dynamical behaviors of nonlinear systems is the so-called Lyapunov exponent.
e 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. erefore, 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). e 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.
e 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. en, 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.
e 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. e 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  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. is 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." e "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.
Similarly, from Figures 6(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 c is small enough. As c or c increases, there exist mixed bifurcation processes, too. e 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. ese numerical results fit very well with Proposition 3. e 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). ese 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.

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. en, 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). en, a homoclinic bifurcation of the two period-14 saddle cycle happens, and it leads to the   e 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. en, a global bifurcation called "contact bifurcation" will happen if the value of α 1 further increases. e chaotic attractor will be destroyed and the so-called "ghost" will occupy the entire basin of attraction after the "contact bifurcation" happens.
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. e 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). is 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. is means that the coexistence of multiattractors has also disappeared.

Social Welfare
Based on the numerical simulations of R&D investment x i 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 From equation (6), we can get that the maximum profits of all the two firms, which are e 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]. at is, According to Qiu [4], we know that the utility function of consumers can be represented as U(Q) � aQ − (b/2)Q 2 , where Q � q 1 + q 2 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 Correspondingly, the social welfare can be rewritten as 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 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.
From Figure 9(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. en, a series of period-doubling bifurcations occur. e 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). e 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. e 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.

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. e 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. e equilibrium may lose stability through a Neimark-Sacker bifurcation or a flip bifurcation. rough analyzing the bifurcation diagram, the phenomenon of "chaotic bubble" is found in our research. e system will fall into a multiperiodic state, or a quasiperiodic state, or a chaotic state, or even divergent state eventually.
e 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. e 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.

Data Availability
No data were used to support the study.

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