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.
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 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. 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 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 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 aswhere is the production output of firm i (i = 1, 2) and are positive constants. Parameter 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 bywhere 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). is the R&D investment of firm i, and measures the R&D spillovers between firm i and its rival. Here, 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 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 aswhere 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. The 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 with respected to to zero. That is,
Equation (6) can be regarded as a linear system of equations about unknowns and . Through 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 and . All the firms will maximize their own profits by choosing the proper R&D investment . Then, an optimization problem can be obtained. Differentiating the function with respect to the variable , then we can get
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, if 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., ) 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 , .
Where is marginal relative profit of firm i, the specific form of it has been given in equation (9). And () 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 is positive/negative, then firm i will increase/decrease its respective R&D investment at period t + 1 with a rate of . According the above analysis, the final dynamic game can be expressed as
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 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 , , , and , where .
It is clear that the equilibrium points , , and are on the coordinate axes of phase plane , so they are often called boundary equilibrium points. In addition, is usually called internal equilibrium point as it is located in the interior of phase plane . We stress that it could be shown that 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 , , , , and are all positive, so it can be inferred that the boundary equilibrium points , and the unique Nash equilibrium point 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. 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 symbol to represent the Jacobian matrix of game (10) at any given point , which has a specific form as
The characteristic equation of can be given aswhere and are the trace and the determinant of , respectively.
According to the theory of nonlinear dynamics, we know that any fixed point will be locally asymptotically stable if all the eigenvalues of the Jacobian matrix evaluated at are located inside the unit circle of the complex plane. That is, all the eigenvalues need to satisfy .
Proposition 1. The trivial equilibrium point is a repelling node.
Proof. Obviously, the Jacobian matrix evaluated at the trivial equilibrium point can be rewritten asQuite evidently, the eigenvalues of are and . Since all the parameters are positive constants and the conditions , , , and always meet according to the economic meaning, then we have and . Therefore, we can easily draw a conclusion that is a repelling node.
In the same way, the similar results about the boundary equilibrium points and can be obtained as follows.
Proposition 2. The boundary equilibrium points and are saddle points or attracting nodes.
Proof. The Jacobian matrix evaluated at the equilibrium point isIt is clear that the Jacobian matrix given above is an upper triangular matrix. The eigenvalues of this matrix can be easily gained, which are and . According to inequalities (11) and the nonnegative conditions , , we know that and if and only if the inequalities and meet. While if the condition and meet, then the modules of the eigenvalues will be less than one, i.e., and . Therefore, the fixed point would be an attracting node in such a situation.
Similarly, we can also prove that the boundary equilibrium point 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 point , (). Actually, the local stability analysis of the internal equilibrium point is more complicated than that of the boundary equilibrium points. The Jacobian matrix at is
The specific forms of trace and determinant of are given aswhere and .
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,
The dual inequalities of these three conditions (18) correspond to three different ways that one of the eigenvalues may cross the unit circle. If , then one of the real eigenvalues of is larger than . While if , then one of the real eigenvalues of will cross the unit circle of complex plane from . Finally, if , then 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
The equation (19) can be expressed equivalently as . 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 . In other words, the pitchfork bifurcation will never happen. Nevertheless, the other two conditions of (18) are equivalent to
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 on the system’s dynamical behavior in this research, the bifurcation curve of flip bifurcation is defined byand the Neimark–Sacker bifurcation curve is defined by , where is shown aswhere and are defined by the system parameters except for . 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 point is a(i)sink if (ii)source if and (iii)saddle if and or and (iv)nonhyperbolic point if either or 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 for the internal equilibrium point is shown in Figure 1(a), where the parameters are fixed as , , , , and , 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 . It means that the position of parameters and is equivalent, which can also be discovered from system (10). Furthermore, we can also find that an increasing of and/or 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 is positively related to the size of the stability region. If we choose the values of the parameters as , , , , and , where only the parameter 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 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 , , , , and . However, the value of the parameter is negatively correlated with the stability region. If the parameters set are taken as , , , , and , 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 , , , , and . If we fix the values of , , , , , and and increase the value of parameter , 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 , , , , and .
(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 , , , , and . 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 and are chosen in the green area of Figure 2(a). However, an increasement of and/or will bring the point 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 and on the local stability of the Nash equilibrium point. In particular, a further increasement of parameter (), 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 , 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.
(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 Figure 3(a) are chosen as the same with that of in Figure 2(a). Firstly, we analyze the dynamical behavior around the line , 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 passes through the line or , 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.
(a)
(b)
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 , 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 , the divergent area has also increased.
(a)
(b)
Through the above discussion, we can see that when other parameters are fixed, an increase of parameters () 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 , , , , and and the parameters and 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 grows. Then, with the speed of adjustment 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 meets , the period-2 cycle can also lose its stability through a Neimark–Sacker bifurcation as the value of increases.
Figure 5(a) gives a series of period-doubling bifurcations that cause instability, where the set of parameters are chosen as , , , , , and . We can notice that the Nash equilibrium point is stable for , and the Nash equilibrium point will lose stability through a flip bifurcation at . 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 .
(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 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 . When , the first period-doubling bifurcation takes place. As the value of parameter further increases, the occurrence of period-doubling bifurcation leads to chaos, eventually.
Similarly, the bifurcation diagram of system (10) with under the parameter set , , , , , and 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 , , , , , and 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 . While if , a period-doubling bifurcation occurs. Then, the period-2 cycle loses stability through a Neimark–Sacker bifurcation at . Subsequently, an attracting invariant circle appears when the parameter exceeds 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 , , , , , and , the bifurcation diagram of system (10) with respect to and the corresponding largest Lyapunov exponent can be plotted, which are shown in Figures 5(g) and 5(h), respectively. When , a period-4 cycle appears and then it will lose stability through Neimark–Sacker bifurcation at . 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 .
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 is stable for small value of parameter . However, as the value of parameter 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 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 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.
(a)
(b)
(c)
(d)
(e)
(f)
Similarly, from Figures 6(c) and 6(d), we can see that the system is chaotic if the marginal cost or the efficient measure of R&D investment cost is small enough. As 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 as the bifurcation parameter. By assuming , , , , , and , some interesting events take place when the bifurcation parameter belongs to the interval (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 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 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 . 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 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.
(a)
(b)
(c)
(d)
Another interesting phenomenon may arise in system (10) is the coexistence of multiple attractors. When the parameter is chosen as the bifurcation parameter, and the values of other parameters are fixed as , , , , , and , 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 in Figure 8(a) is chosen as . 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 further increases to , 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 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.
(a)
(b)
5. Social Welfare
Based on the numerical simulations of R&D investment 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
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,
According to Qiu [4], we know that the utility function of consumers can be represented as , where 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 aswhere the forms of and 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 , , , , , and , whereas the adjusting speed of firm 1 is chosen as the bifurcation parameter.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
From Figure 9(a), we can see that the social welfare is stable when . And then the balanced social welfare loses its stability via a period-doubling bifurcation when the parameter value meets . 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 . As the value of 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 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 , , , , , and , 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 increases. The varying curves of average social welfare and average profit with parameter 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.
Data Availability
No data were used to support the study.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research was supported by the National Natural Science Foundation of China (grant no. 61863022), the Young Scholars Science Foundation of Lanzhou Jiaotong University (grant no. 2015029), and the Foundation of Humanities and Social Sciences from the Ministry of Education of China (grant no. 15YJC820007).