Abstract
In this paper, a dynamic two-stage Cournot duopoly game with R&D efforts is built. Then, the local stability of the equilibrium points are discussed, and the stability condition of the Nash equilibrium point is also deduced through Jury criterion. The complex dynamical behaviors of the built model are investigated by numerical simulations. We found that the unique route to chaos is flip bifurcation, and the increase of adjusting speed will cause the system to lose stability and produce more complex dynamic behavior. In addition, we also found the phenomenon of multistability in the given model. Several kinds of coexistence of attractors are shown. In particular, we found that boundary attractors can coexist with internal attractors, which also aggravates the complexity of the system. At last, the chaotic state in the built system has been successfully controlled.
1. Introduction
In recent years, the market competition mode among firms are no longer subjected to price or quantity, but more towards soft power, including technical level, talent resources, and corporate culture. With more and more systematic and mature development of the organization management for a great quantity of large firms, R&D (the abbreviation of Research and Development) is becoming more and more important for them. The visionary entrepreneurs will spend a lot of money to build their own exclusive R&D team. It may be risky for the capitalists who pursue maximum profits to do such things which may be harmful to their revenues. However, the reasons for these firms that lay much stress on R&D activity can be listed as follows: on the one hand, R&D activities can improve production efficiency, reduce production costs, and then the entrepreneurs can gain surplus value proliferation; on the other hand, R&D activities can also improve quality of products so as to increase advantage competitiveness of firms in the market. Hence, R&D activity gradually attracts a lot of attentions of economists.
In the process of carrying out R&D activities, due to the flow of technical talents, information diffusion, and other factors, R&D spillover effect may emerge among firms [1, 2]. Numerous scholars studied and advanced this problem. Zhou et al. [3] studied the types of bifurcation, intermittent chaos, and coexistences in a dynamical two-stage Cournot duopoly game with R&D spillover effect. Zhou and Wang [4] studied the stability and multistability of a duopoly game model with R&D spillover. While in the case that the firms take highly confidential measures for its R&D achievements, the positive effects of R&D achievements only promote the firms that carry out R&D activities, but not the competitors. In other words, there is no spillover effect. Based on this assumption, D’Aspremont and Jacquemin [5] built the foundation of the performance of R&D cooperation in the product differentiation market. Luckraz [6] established a Cournot duopoly game model with isoelastic demand function and also considered the case without R&D spillovers. The work of Chatterjee [7] showed that the firms in the market undertake R&D activity and studied that R&D incentive may depend on the information structures they had. In addition, Ferreira et al. [8] constructed myopic optimal discrete and continuous R&D dynamics.
The most extensive research on spillover effect is technology spillover. AJ Model [5] and KMZ model [9] are the most representative ones. In both cases, they are both two-stage game models. Also, the two-stage game presents good adaptability in dealing with problems about uncertainty of R&D market and mutuality of competitive strategies. That is to say, on the first stage, both these two firms conduct R&D activities to reduce their production costs. And on the second stage, they compete in market in term of the price or the quantity. The investigations on the two-stage model have gradually become the focus of economists. While at the beginning of the study, many scholars were limited in studying the properties of the static model [10–12]. All the investigations about the static model just illustrated the equilibrium solutions under different R&D strategies. No one can deny that the static model does have some limitations. One limitation is that it can only analyze the single supply and demand balance between firms. Once the supply and demand factors change, the corresponding supply relationship will shift.
In this research, we introduce nonlinear dynamics theory to discuss evolution process among firms’ games. Many scholars had carried out some studies about this issue. Bischi and Lamantia [13, 14] revealed the local and global dynamical properties of the two-stage oligopoly game model with adaptive dynamic mechanism and described its complex evolutionary behaviors. Askar et al. [15] investigated the properties in a dynamic Cournot duopoly game model with nonlinear demand function. A homogeneous Cournot duopoly market possessing two different mechanisms was proposed by Elettreby [16]. Most of the similar works can be seen in [17–19]. Additionally, the dynamical two-stage game has raised much interest in the economic dynamical system. For instance, Zhang et al. [20] proposed a two-stage Cournot duopoly game of semicollusion in production and analyzed its complex dynamics, such as local stabilities and existence. Ma et al. [21, 22] have focused the analysis of complex dynamical behaviors of the two-stage dynamic games.
Furthermore, although in the electronic information world, the corresponding channels for firms to obtain market information have proliferated. However, it cannot be ignored that if the firms want to get complete information, they still need to pay its costs. That is to say, few firms are willing to pay high cost to obtain all the information about market demand and consumer preferences. So, most firms have incomplete information about the market, so they could be regarded to be bounded rational [23] when making decisions. In fact, compared with those having complete information of demand function, it is much easier to obtain the existing information, such as local estimation of marginal profit. That is, gradient adjustment mechanism which has attracted many scholars to carry out deep research studies on this topic (see e.g., [23–26]). Cavalli and Naimzada [27] took the factors such as rationality and information set into consideration and adopted a traditional dynamical adjustment mechanism based on classical best reply function. Elsadany [28] interpreted a Cournot duopoly dynamics model with different strategies including bounded rationality. Fanti and Gori [29] put forward a Bertrand duopoly game model with heterogeneous product expectations and further studied the properties of the model. And Naimzada and Sbragia [30] dealt with an oligopoly game model with nonlinear demand function and bounded rationality.
Additionally, other topics worth studying in nonlinear dynamics are global dynamics, synchronization, and multistability. In recent years, great deals of scholars have deepened the study of these topics. Askar [31] dealt with the dynamical characteristics with log-concave demand of equilibrium points, such as stability, bifurcation, and chaos. Andaluz and Jarne [32] analyzed the dynamics of not only quantity-setting competition but also price-setting competition. Agliari et al. [19] observed that some bifurcations can lead to changes in the structure of attractors and their basins. And a higher degree of production differentiation tends to reduce competition. Tu and Wang [33] built a dynamical R&D two-stage input competition triopoly game model with bounded rationality and further analyzed the dynamics properties of this model. Bischi et al. [34, 35] had studied the phenomena of dynamic properties, such as synchronization, intermittency, coexisting attractors, and complex basins. Fanti et al. [36] provided local and global dynamic characteristics in a duopoly game model with price-setting competition and market share delegation.
Basing on the above analyses, the behaviors including dynamic evolution, intermittency, and multistability are exploited to interpret the complexity of competition. We take use of a Cournot duopoly model with differentiated products to consider the impact of adjustment speeds. And the influence of firm’s adjustment speed on the development of firms is studied through combining a two-stage dynamic game with the method of numerical simulation and the properties of nonlinear dynamics. It is worth to focus on whether the firm will gain more advantages at a lower adjustment speed.
Furthermore, another noticed feature is chaos control. Since chaos could increase the disorder and unpredictability of the established model. In the sense of economics, if actual economic market is in a chaotic state, then the resources of firms will not be controlled by market supply, demand relationship, and the government macrocontrol to some extent. This may disrupt the development of market economy. However, as for firms in the market, they prefer to pursue a long-term and stable development, which can raise their accumulation of capitals to expand themselves industrial scale. Hence, it is necessary to control chaos. A large number of scholars have carried out relevant research studies and achieved some results. For instance, feedback control was given by Elsadany et al. [37, 38]. Ma and Ji [39] studied the dynamical repeated Cournot model by using the chaos control method with feedback control. Parameters adjustment was provided to control chaos by Askar and Al-khedhairi [40].
The main purpose of this research is to investigate the influence of continuous change of system parameters on the eventual dynamical scenarios of the dynamic duopoly model with R&D efforts. The local stability of boundary equilibrium and Nash equilibrium is provided by Jury stability criterion and numerical simulation to obtain the internal complexity of the model. We show that the successive change of system parameters not only causes the system to lose its stability but also causes the high complexity of eventual behaviors of the system, such as the number of coexisting attractors, the critical bifurcation of attractors, and the global bifurcation of attracting domain.
This research is organized as follows. In Section 2, a two-stage dynamical game model without R&D efforts is established, where these two firms’ products are of the same kind but different products. In addition, its stability and the conditions of stability of equilibrium points are discussed. In Section 3, the complex evolution law based on speed of adjustment is discussed in details. Through numerical simulation, the complexity of bifurcations, evolution of strange attractors, intermittent chaos, multistability, and control of chaos are studied. At last, the impacts of adjustment speeds on built model are concluded in Section 4.
2. The model
Assume that there are two technology-based firms operating in the market, labeled by i (i = 1, 2), and they supply homogeneous products to the same market. Research and Development (R&D) is a kind of creative behavior for the firms to increase production efficiency and decrease costs so as to gain more market share. While in the actual economics environment, R&D activity is always carried out in a stepwise manner: research phase and development phase. The process of competition between two R&D firms can be simulated by a two-stage game. On the first stage, these two firms conduct R&D activities to reduce their production costs. And on the second stage, they compete in the market in terms of the price or the quantity. The difference between the R&D levels of the firms can result the difference of goods in quality, function, and externalization. Hence, the inverse demand function can be written aswhere is the price of firm i (i = 1, 2), and represents the outputs of firm i and its rival, respectively, denotes the maximum scale of this market, is defined as the differentiating degree of the products of two firms, and . Moreover, if , it means that the products of firm 1 and firm 2 are mutually independent, and the oligopoly market will be degenerated into two completely independent monopoly markets. If , it indicates that the products of these two firms can replace each other. And if , it means that there is no difference between the products of these two firms. In other words, the products of firm 1 and firm 2 are totally identical.
Consider now the R&D efforts of firm i (i = 1, 2) labeled by . Then, the effective marginal cost of firm i can be represented aswhere represents the marginal cost of two firms without presence of R&D activities and . Equation (2) states that technology is totally confidential and there is no R&D spillover. The investments of firm i (i = 1, 2) can be seen as a quadratic function with respect to its R&D efforts, which can be expressed as (i = 1, 2), where is a parameter about the technical innovation cost of firm i. The innovation ability of firm i is stronger when the technical innovation cost parameter is tinny. Therefore, by following the above assumptions, the profits of firm 1 and firm 2 can be represented as follows:
Substituting equations (1) and (2) into equation (3), we can get the specific expression of profit function, which could be given as
The marginal profits of these two firms captured by taking the derivative of (4) with respect to (i = 1, 2) are given as
Letting (i = 1, 2), the subgame perfect Nash equilibrium (SPNE) of firm i (i = 1, 2) in the second stage can be obtained as
The profit function about R&D efforts captured by taking equation (6) into equation (4) in reverse order can be obtained as
Computing the derivative of equation (7) with respect to (i = 1, 2), the local marginal profits about the R&D efforts of two firms are given as
Actually, it is impossible to make perfect rational decision with incomplete market information, such as incomplete demand information, the investments of their competitors, and the level of technical innovation. These incomplete elements make the firms “myopic” in some aspects. These “myopic” firms can only adjust their strategies through continuous “trial and error.” That is, the output at the next period is decided according to their own experience, so these firms can be seen as bounded rational. The firm shall have a decision on its R&D effort at t + 1 period by using a local estimate of marginal profit at t period. More in detail, if , the firm will increase its R&D investment in the next period. And the firm will keep invariant outputs in next period if the local estimate of marginal profit meets . Then, if , the R&D investment of the firm in the next period will be reduced. On the basis of the above thoughts, a dynamical gradient adjustment mechanism is introduced [20]; then, the dynamic model can be constructed as
Taking equation (8) into equation (9), the two-dimensional nonlinear discrete-time evolution of dynamic system (9) is as follows:where denotes the adjustment speed, which determines the variation ratio of R&D efforts of firm i in the next period.
2.1. Stability Analysis of the Equilibrium Points
System (10) is a two-dimensional noninvertible map whose iteration decides the trajectory of R&D effort. In order to discuss local qualitative behaviors of the solutions of system (10), we need to find equilibrium points of the game between these two firms. Let ; then, we can get the nonlinear equations as follows:
By solving these nonlinear equations, four equilibrium points of system (10) can be obtained. They arewhere . The fixed points , , and are boundary equilibrium points. And the fixed point is an unique Nash equilibrium point. In the sense of economy, the inequalities and should be satisfied. And if the other parameters satisfy the conditionall the four equilibrium points of system (10) will be nonnegative. The corresponding Jacobian matrix of system (10) at any point has the form
In addition, the local stability of equilibrium points of system (10) depends on the eigenvalues of Jacobian matrix at the given equilibrium point.
Proposition 1. The boundary equilibrium point is an unstable node.
Proof. The Jacobian matrix at is given as follows:Obviously, the Jacobian matrix is a diagonal matrix; thus, its corresponding distinct eigenvalues areWith the constricts and , is satisfied, then we can get that the equilibrium point is an unstable node.
Proposition 2. The types of local stability of the equilibrium point can be classified as follows:(1)If the condition meets, the equilibrium point is an unstable node(2)If the conditions , , and meet, then the equilibrium point is an unstable saddle(3)If the conditions , , and meet, then the equilibrium point is a stable node(4)Else if and are satisfied, then the equilibrium point is an unstable point
Proof. The Jacobian matrix at iswhere is a lower triangular matrix; thus, its eigenvalues are main diagonal elements, that is,The type of local stability of this boundary equilibrium point is determined by different values of the parameters:(1)If , then ; thus, the equilibrium point is an unstable node.(2)If and are satisfied, then meets. And if , then meets; hence, the equilibrium point is an unstable saddle.(3)If and , then . And if , then , so the equilibrium point is a stable node.(4)If the parameters meet with and , then ; thus, the equilibrium point is an unstable point.
Proposition 3. The types of equilibrium point are similar to the boundary equilibrium point .
Proof. The Jacobian matrix at can be written as follows:Matrix (19) is an upper triangular matrix and the corresponding eigenvalues areIn addition, the boundary equilibrium points and are symmetrically located around the diagonal line in the plane . And the local stability conditions of are symmetrical to the conditions of listed in Proposition 2, so the conclusions of are no longer repeated here.
In sense of economics, the boundary equilibrium points interpret a situation that one of these two firms has withdrawn from the market. The equilibrium points and can be regarded as the situation that one of these two firms takes the lead in the market. That is to say, the oligopoly market becomes the monopoly market, and the local stability of equilibrium points implies the stability of economic market in a short term. Actually, neither of the above two market situations is expected, only when these two companies restrict each other can they be conducive to the stable development of the market and the country. That is, the game between firms reaches an equilibrium, named “Nash equilibrium.” Nash equilibrium [41] is on the premise of maximizing its own profits. Meanwhile, it also ensures the stable development of the market.
The Jacobian matrix at the Nash equilibrium of system (10) is formulated aswhere . To simplify the process of analyses, the auxiliary variables , , and are introduced. Thus, the matrix can be rewritten as follows:The characteristic polynomial of matrix (22) iswhere and denote the trace and the determinant of matrix (22) at Nash equilibrium point . And the specific formulation of and are given as
Proposition 4. Neimark–Sacker bifurcation will not occur at the Nash equilibrium point .
Proof. The corresponding discriminant of roots of characteristic polynomial (23) is given aswhich means that the dynamical system (10) will not generate complex eigenvalues at the Nash equilibrium point ; hence, the Neimark–Sacker bifurcation will not occur at the Nash equilibrium point .
On the basis of canonical stability analysis, the sufficient condition of local asymptotical stability of Nash equilibrium point is that the eigenvalues of corresponding Jacobian matrix of is within the unit cycle. According to Jury criterion [15], the condition of local asymptotical stability of Nash equilibrium point can be expressed in detail as(i)(ii)(iii)
Proposition 5. When is satisfied, the second condition of Jury criterion holds. That is,
Proof. First, let , and then its roots can be computed as and . We recall that , , and ; then, we can obtainThe numerical relationship between and is that . If the condition holds, then satisfies the condition or . Following the constriction , can be captured. Hence, does not hold. Only when the condition meets, the second condition of Jury criterion holds. The corner case isHence, when and meet, the second condition of Jury criterion holds.
Proposition 6. If the first condition of Jury criterion holds, the third condition of Jury criterion will certainly hold.
Proof. Let of first condition of Jury criterion, and we can get . As for , , , , and , we haveAfter this, we can deduce that .
Let , if , the third condition of Jury criterion will hold. That is,In the same time, it can be easily proved thatwhich is satisfied for any value of . Following the above formulations, we can get . The inequalities andcan be captured through jury criterion and inequality (31). Then, we can deduce the conclusion that if the first condition of Jury criterion holds, the third condition of Jury criterion will certainly hold. In summary, system (10) has an unique Nash equilibrium point . When and , system (10) is local asymptotical stable at Nash equilibrium point . Under the above constrictions of the system parameters, these two firms will reach equilibrium state and get a maximum value for total profit after several games. Then, the market will keep stable development.
3. Analysis of Complex Evolution Law Based on Speed of Adjustment
The local stability of boundary equilibrium points and Nash equilibrium point are discussed by the way of analytics. In this part, in order to investigate the intrinsic laws of evolution of system (10) and to mainly analyze its local and global properties, the method of numerical simulation through two-dimensional bifurcation diagram, one-dimensional bifurcation diagram, two-dimensional largest Lyapunov exponent, and basin of attraction is used.
3.1. The Routes to Chaos
Firstly, the parameters are fixed as , , , and , and then the speeds of adjustment and are chosen as the bifurcation parameters. The global bifurcation route with respect to speed of adjustment is illustrated in Figure 1, where Figure 1(a) is a two-dimensional bifurcation and Figure 1(b) is a two-dimensional largest Lyapunov exponent corresponding to Figure 1(a).
(a)
(b)
In the right of Figure 1(a) is a color bar with different colors, that is to say, these different colors represent different periods. Specifically speaking, the brown area denotes the stable area of Nash equilibrium point, i.e., period-1. The green area denotes period-2, the orange area denotes period-3, the yellow area denotes period-4, and so on; the color of deep black denotes period more than 30. As a result of the constraints of algorithm, the states of quasi-period, chaotic, and escaping of system (10) cannot be classified clearly, so all three states are marked by dark black in Figure 1(a). In addition, when the largest Lyapunov exponent is equal to zero then the corresponding state of system (10) is quasi-period, and when the largest Lyapunov exponent is greater than zero and less than 2, then the state of system (10) is chaotic. The last situation is when the largest Lyapunov exponent is greater than 2, and the corresponding state of system (10) is escaping area. So, the three states mentioned above can be distinguished by this method displayed in Figure 1(b), where the color bar in the right of this figure represents the largest Lyapunov exponent.
The Nash equilibrium point corresponding to Figure 2(a) is . With the increase of , the corresponding eigenvalues at the Nash equilibrium point will increase to and . Then, we can deduce that this Nash equilibrium point is a stable node. And more in detail, when the speed of adjustment , then the state of system (10) is at the Nash equilibrium and will keep stable. While if the speed of adjustment satisfy that , then the stability of this Nash equilibrium will change and a flip bifurcation will occur. Once the speed of adjustment meets the condition that , system (10) will enter into a chaotic state. And in the corresponding largest Lyapunov exponent to Figure 2(a), when the largest Lyapunov exponent is less than zero, the state of system (10) is stable, and system (10) enters into the chaotic state when the largest Lyapunov exponent is greater than zero. And when the largest Lyapunov exponent is equal to zero, a bifurcation takes place in system (10).
(a)
(b)
As it is known to all, the flip bifurcation is a classical trajectory for economic system to change the state from stable into chaotic. A new way of development for firms will generate when the economy develops to a certain stage, and there is only one way for firms to further develop the economy and they will face a dilemma in the choice of the new and old way. As a result, the route of period-doubling is generated in constant choice between firms in dilemma. When the speeds of adjustment of these two firms are relatively small, then the two firms will be in the Nash equilibrium state and maintain its profit maximizing, respectively. In the same time, these two firms mutually hold each other back. Once the speed of adjustment exceeds the threshold value, the stability of Nash equilibrium point will be destroyed and then the period-doubling will occur. Finally, the constant choice of the firms leads to the disorder state of the market.
3.2. The Evolution of Strange Attractor
In the case of a nonlinear dynamic system, its final state inflects the eventual behavior of system because it is a key point to focus on. Furthermore, the final state of the nonlinear system can be indicated by the attractors representing asymptotical behaviors of the nonlinear system if the number of iteration times tends to infinite. In this research, we discuss the impact of adjustment speed of firm 2 on eventual behavior of system (10) by investigating the evolution of attractors, where the parameters are chosen as the same as in Figure 1(a), and the adjustment speed of firm 1 is fixed as . The attractor is a period-4 focus with rough selvedge when the parameter equals to 0.625, which can be seen in Figure 3(a). As shown in Figures 3(b) and 3(c), as the adjustment speed of firm 2 increases, this period-4 focus is transformed to four invariant cycles with selvedge via a Neimark–Sacker bifurcation, and these four invariant cycles become more and more larger and the rough selvedge disappears simultaneously. While when the parameter further grows up to 0.657, these four invariant cycles break to form the attractor shown in Figure 3(d); system (10) enters into “periodic window.” And the broken attractor continues to spread then recombine to form four chaotic attractors, which can be seen in Figure 3(e). When the adjustment speed of firm 2 continues to increase until , the four chaotic attractors mentioned above conduct the dynamic behaviors that spread and merge and finally form a piece of chaotic attractor, see Figure 3(f).
(a)
(b)
(c)
(d)
(e)
(f)
The increase of adjustment speed leads to the complexity of eventual behavior of system (10). This process of evolution presents the properties of “certainty” and “irregular development,” where “certainty” means that the way of bifurcation, that is to say, period-4 Neimark–Sacker bifurcation, and “irregular development” means that the interior structure of the attractors, as well as complexity of interior structure of the chaotic attractor, as shown in Figures 3(e) and 3(f). No matter if from phenomenon development or interior mechanism, system (10) presents a high of certainty and complexity. That means the economic system makes mistakes in predicting, and due to the influence of external condition the development of economy cannot be forecasted anymore.
3.3. Multistability
In the process of studying the nonlinear dynamic system, we can observe that the stability of equilibrium points may vary with the change of parameters and initial conditions by numerical simulation. This variety can cause single-level or multilevel bifurcation so as to lead to the system displaying complex dynamical behaviors, such as coexistence of attractor, fractal, and chaos. The coexistence of the attractor implies the motion of multistability of system. That is to say, the bifurcation of the nonlinear system varies with the change of the parameters and initial conditions, which will lead to the change of number of solution. And the number of bifurcation solution varies raised a phenomenon of coexistence of multiattractors, so there is a close relationship between the motion of multistability and bifurcation. In addition, the phenomenon of coexistence of the attractor also indicates that the global bifurcation behavior of the system and the discussion of global bifurcation behavior are mainly reflected in two aspects. One is about the variation of the structure, shape, and amount of attractors with parameters, and another is about the variation of the structure and shape of basin of attraction with parameters.
Then, we fixed the values of parameters as , , , , and . This set of parameters has a case that the Nash equilibrium point is an unstable node, and the corresponding eigenvalues, respectively, are and . As shown in Figure 4, with the increase of the adjustment speed of system (10) can lead to the phenomenon about the global bifurcation of structure, shape, and amount of attractors.
(a)
(b)
(c)
(d)
When the adjustment speed is equal to 0.4865, the chaotic attractor coexists with a period-2 attractor, see Figure 4(a), where the red area denotes the escaping area. The initial value corresponds to two Lyapunov exponents, and ; the corresponding attractor is period-2 attractor and its attracting area is the light blue area. The initial value corresponds to two Lyapunov exponents, and ; the corresponding attractor is the chaotic attractor and its attracting area is the white area.
When the adjustment speed increases to , the period-4 attractor suddenly adds into the basin of attraction, which can be seen in Figure 4(b). The increase of parameter raises the change of the structure of the chaotic attractor and the number of the attractor, at the same time, also gives rise to the change of corresponding largest Lyapunov exponent of attractor. The initial value corresponds to two Lyapunov exponents, and ; the corresponding attractor is a period-2 attractor. The initial value corresponds to two Lyapunov exponents, and ; the corresponding attractor is period-4 attractor. The initial value corresponds to two Lyapunov exponents, and ; the corresponding chaotic attractor becomes a hyperchaotic attractor. From Figure 4(b), it is easy to know that the attracting area of period-2 attractors is the light blue area, the attracting area of period-4 attractors is the blue area, and the attracting area of super-chaotic attractor is the white area.
And if we further increase the value of parameter to , then we can get an interesting phenomenon that the period-4 attractors become the period-8 attractors via the period-doubling bifurcation, and the chaotic attractors constantly enlarge so the scenario of coexistence, as shown in Figure 4(c). The initial value corresponds to two Lyapunov exponents, and ; the corresponding attractor is period-2 attractor and its attracting area is the light blue area. The initial value corresponds to two Lyapunov exponents, and ; the corresponding attractor is period-4 attractor and its attracting area is the blue area. The initial value corresponds to two Lyapunov exponents, and ; the corresponding attractor is still hyperchaotic attractor and its attracting area is the white area.
And then if we continue to increase the value of parameter , the chaotic attractor will contact with the boundary of the blue area and a “contact bifurcation” will occur, which can be seen in Figure 4(c). As shown in Figure 4(d), this “contact bifurcation” can lead to the chaotic attractor burst and then form a “ghost” existing in the former attracting area. Based on above analysis, it is obvious that, with the increase of the parameter , the largest Lyapunov exponent of period-2 decreases but the corresponding attracting area is the light blue area keeping invariantly, while the largest Lyapunov exponent of the chaotic attractor increases, where the attractor becomes hyperchaotic from chaotic, and its attracting area is the white area. This white area is constantly occupied by the blue attracting area so as to shrink with the presence of period-4 attractors. In addition, the number of attractors has increased from two to three and the structure of the attractor changes. The last bifurcation form of the chaotic attractor is the “critical bifurcation,” when this chaotic attractor contacts the boundary of its basin of attraction.
The parameters are given as , , , , and . A series of special phenomena about the coexistence are discussed in Figure 5, which indicates the one-dimensional bifurcation of different initial values of the adjustment speed ... The initial value corresponds to Figure 5(a) which is , and the initial values corresponds to Figure 5(c) which is , where the blue curve represents the R&D effort of the firm 1 and the red curve represents the R&D effort of the firm 2. The diagrams in Figures 5(b) and 5(d) are the largest Lyapunov exponent corresponding to Figures 5(a) and 5(c), respectively. Let us compare the one-dimensional bifurcation Figure 5(a) with Figure 5(c), and we can draw a result that these two bifurcating diagrams have an obvious difference between them when the adjustment speed satisfies this condition as , which means that due to the different choices of initial values, a huge difference of the eventual state of the system will arise. And if we constantly increase the value of the parameter, tends to a relative large number; then, firm 1 and firm 2 will withdraw from the market in turn. Especially, the way of bifurcation of system (10) is a nonnormal bifurcation in Figure 5(c) when parameter belongs to the interval . As shown in the local enlargement of Figure 5(d), the corresponding largest Lyapunov will present a form of oscillation in this interval and the detailed discussion is shown in Figure 6.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(e)
(f)
In the case of this group of parameters, the boundary equilibrium point is a stable node and the corresponding eigenvalues are and , respectively. When the parameter is , the phenomenon of period-2 attractors coexisting with the stable boundary equilibrium is shown in Figure 6(a), where the red area represents the escaping area. And the attractor corresponding to the initial values is this boundary equilibrium point , and its attracting area is the white area. The attractor corresponding to the initial values is the period-2 attractor, and its attracting area is the light blue area, where the diagram in right upper corner of Figure 6(a) is the enlargement of period-2 attractors.
When the parameter is , the period-2 attractors become the chaotic attractors through a nonnormal bifurcation, see Figure 6(b), where the corresponding largest Lyapunov exponents are and . From the enlargement diagram in the right upper corner of Figure 6(b), we can clearly see these chaotic attractors have an obvious level of structure, by using the colors of green, red, and black separately represent the inner, middle, and outer layers of the attractor. Each layer is surrounded by numerous points, and the outer layer is not connected with the inner layer and the middle layer, but the boundary of the inner layer and the inner layer is nested. In addition, with increase of the adjustment speed , the white attracting area corresponding to the boundary equilibrium point shrinks much more smaller, while the light blue attracting area expands more and more larger. When the value of adjustment speed grows up to 0.472, the three layers of chaotic attractors split from each other and then the outer expands outward continuously to form the attractor, as shown in Figure 5(c). In this case, the attractors are periodic attractors, and the corresponding largest Lyapunov exponent are and . And let the value of parameter further increase, then the inner and middle layers of the periodic attractors vanish. Finally, the outer layer constantly expands outward and mutually merge into the periodic attractors, as shown in Figure 6(c), where the corresponding largest Lyapunov exponent are and .
With the increase of adjustment speed , the periodic attractor will contact its attracting area so as that “critical bifurcation” occurs. And until the value of adjustment speed increases to 1.56, the boundary equilibrium point becomes the boundary chaotic attractor represented by the blue straight line via flip bifurcation, as shown in Figure 6(e). In the same time, the periodic attractor undergoes the “critical bifurcation” and then reforms the boundary chaotic attractor represented by the black straight line, where the attracting area corresponding to the chaotic attractor represented by the blue straight line is the white area, and the attracting area of another boundary chaotic attractor is the light blue area. In addition, these two different attracting areas locate symmetric along the diagonal. If the value of parameter continues to increase, then the basin of attraction will go through a bifurcation, named “global bifurcation,” and the chaotic attractor represented by blue line disappears, forming a “hole,” as shown in Figure 6(f) in the light blue attracting area.
From the study above, we can deduce that there is no relationship between the bifurcation of changing qualitative property of basin of attraction and the bifurcation of changing qualitative property of the attractor. After the local bifurcation, the game between the firms will not converge to Nash equilibrium point representing the global optimal strategy anymore. Even though this game starting near Nash equilibrium point, system (10) will also tend to different attractors, and these attractors may be periodic or aperiodic. Furthermore, these bifurcations may also lead to the lack of predictabilities of system (10). The global bifurcation with respect to the boundary of basin of attraction increases the uncertainty of the game fate of the firms as given in the initial strategy. Hence, a slight variation in initial condition of the firms or a small external disturbance in the process of adjustment, will cause a huge influence on the evolution behaviors of the firms in the long run.
3.4. The Control of Chaos
With the increase of the speed of adjustment, the oligopoly market will enter into the disorder and chaotic state, and the nonlinear dynamic system is at the chaotic state. It is harmful to the development of the firms when this kind of fierce economic fluctuation serves in the market. For the sake of maintaining the longrun and stable development of two firms in this market, we can take control of system (10), aiming at expanding the stable range at Nash equilibrium of system (10) by introducing a control factor. Based on the opinions above, we exploit the method of state feedback and parameter variable control [40] to take control of system (10). And then the form of system (10) after conducting the control factor can be written as follows:where the parameters and are control factors and the restrictive condition is .
Firstly, the parameters are given as , , , , , and , and we choose the control factors and as the parameters of bifurcation. The route of global bifurcation with respect to the control parameters is indicated in Figure 7, where Figure 7(a) is a two-dimensional bifurcation and Figure 7(b) is the largest Lyapunov exponent corresponding to Figure 7(a).
(a)
(b)
It can be observed from Figure 7(a) that there is only one route for system (33) which leads to the Nash equilibrium from chaotic, that is, inverse flip bifurcation. The control parameters and successively pass through the dark black area, yellow area, green area, and eventually enter into brown area, which implies that system (33) changes from chaotic to period-4, period-2, and other period-halving and finally reaches the Nash equilibrium. Figure 7(b) is a diagram of the largest Lyapunov exponent corresponding to Figure 7(a), where the color bar in the right of this figure represents the largest Lyapunov exponent. When the largest Lyapunov exponent is equal to zero then the corresponding state of system (33) is quasi-period, that is, the dark grey area. And when the largest Lyapunov exponent is greater than zero and less than 2, then the state of system (33) is the chaotic state, that is, grey area. The last situation is when the largest Lyapunov exponent is greater than 2, and the corresponding state of system (33) is escaping trajectories, that is, the white area. Fix the control factor as and ; then, the one-dimensional bifurcation and the corresponding largest Lyapunov exponent are displayed in Figures 8(a) and 8(b), respectively, where the red curve represents the R&D effort of firm 1, the blue curve represents the R&D effort of firm 2, and the black curve represents the largest Lyapunov exponent (abbreviated as ).
(a)
(b)
From Figure 8, we can observe that when control factor is fixed as , and if another control factor satisfies the condition that , then system (33) is at the escaping state, that is, the market disappears. If the value of meets , then system (33) is at the chaotic state. If the value of meets , then system (33) conducts an inverse flip bifurcation. And if is met, then system (33) will enter into the Nash equilibrium. And with regards to the corresponding largest Lyapunov exponent, when the value of meets , then . The corresponding largest Lyapunov exponent is less than zero when . When the corresponding largest Lyapunov exponent is equal to zero, then system (33) takes place a bifurcation. Another case is that the control factor increases to 0.65. In this case, when the condition is met, then system (33) is at the escaping state. When the condition is met, then system (33) is at the chaotic state. When the condition is satisfied, then system (33) conducts an inverse flip bifurcation. And if is met, then system (33) will enter into the Nash equilibrium. And with regards to the corresponding largest Lyapunov exponent, when the value of meets then . The corresponding largest Lyapunov exponent is less than zero when . When the corresponding largest Lyapunov exponent is equal to zero, then system (33) takes place at bifurcation. In addition, we can also know that the larger the control factor is, the faster system (33) enters into Nash equilibrium.
For sake of eliminating or slowing down the speed of system (10) entering into chaos, the control factors are introduced and play a stabilizing role in system (33). We will illustrate how they affect the behavior of system (10) entering into chaos through the control factors and the influence on the adjustment speeds and . As shown in Figure 9, a two-dimensional bifurcation with respect to adjustment speeds when the parameters are fixed as , , , and .
(a)
(b)
(c)
(d)
In this group of parameters, the Nash equilibrium point is . The increase of control factors and adjustment speeds have no impacts on Nash equilibrium point . It can be seen from Figure 9 that the stability region with respect to gradually expands at the increase of control factors. However, it does not affect the global shape of two-dimensional bifurcation with respect to , but increases the stability region of the Nash equilibrium. The purpose of introducing the control factor is to expand the region of the Nash equilibrium point attracting area, so as to further delay the speed of system (10) to enter into chaos with increasing of parameters, so as to keep the market stable as long as possible.
Furthermore, we will illustrate the influence of control factors on local bifurcation of system (33). Under the action of control factors, the one-dimensional bifurcation and the corresponding largest Lyapunov exponent about adjustment speed are shown in Figure 10, where the parameters are the same as in Figure 9, and the adjustment speed of firm 1 is fixed as . In Figure 10, the blue curve represents R&D effort of firm 1, the red curve represents R&D effort of firm 2, and the local enlargement is the one-dimensional bifurcation of R&D effort with respect to adjustment speed . When the control factors are chosen as , system (33) is equivalent to system (10). At present, if the adjustment speed of firm 2 is less than 0.465, that is, , then system (33) will enter into the chaotic state and the market enters into a disorder state. And if , then system is at escaping state, and the corresponding largest Lyapunov exponent is greater than zero, that is, . If the largest Lyapunov exponent is less than zero, that is, , then system (33) is at the period state, see Figure 10(a) for more details. When the control factors are chosen as , and the adjustment speed of firm 2 is chosen in the interval that , then system (33) is at the Nash equilibrium state, see Figure 10(b). System (33) will take place a bifurcation, i.e., flip bifurcation, when the adjustment speed of firm 2 is chosen in the interval that . When the condition is satisfied, system (33) enters into chaos. Lastly, when the condition is met, system (33) is at the escaping state. The corresponding largest Lyapunov exponent is negative when . And when , the corresponding largest Lyapunov exponent is positive, that is, . In addition, if we further increase the value of control factors to 0.5 (that is, ), then we can get the results as follows. System (33) is at the Nash equilibrium when . When , system (33) takes place at bifurcation. When , then system (33) will enter into chaos. And when , then system (33) is at the escaping state. If , then the corresponding largest Lyapunov exponent is negative, that is, . If , then the corresponding largest Lyapunov exponent is positive, that is,. The detailed dynamical behaviors under are shown in Figure 10(c). While when and , system (33) is in a stable period-1 state, see Figure 10(d). From Figure 10(d), we find that system (33) takes place at a series of bifurcations when . And when , then system (33) is at the escaping state. If , then the corresponding largest Lyapunov exponent is negative, that is, . The control factors are increased from 0 to 0.95, the stability of system (33) gradually gets stronger and the local bifurcation point is increased from 0.534 to 12.34. It means that the increase of control factors are more beneficial for system (33) to maintain stability, delay the bifurcation point, and keep local stability of the Nash equilibrium point in the long run. However, as for the firms, any slight change of parameters may lead to a huge fluctuation in the market, which makes it difficult for the firms to keep pace with or catch the changing trend of the market. While this fluctuation can be eliminated or slowed down by introducing a control strategy, it is a benefit to a long-term stable development for the market.
(a)
(b)
(c)
(d)
4. Conclusion
In this research, a nonlinear game model of the R&D competition between duopoly firms is studied in detail. The game between these two firms is mainly described by the way of a two-stage game. In the first stage, all of firms determine the level of their R&D efforts to reduce the production cost. And in the second stage, supposing that all firms compete in the market with a Cournot form and determine the outputs to maximize themselves profits, respectively, the local and global properties are discussed by theoretical analysis and numerical simulation. Theoretically, the local stability condition of the equilibrium points and the kind of local bifurcation are analyzed through stability theory and the Jury criterion.
In numerical simulation, we choose the adjusting speeds of these two firms as bifurcation parameters. The evolution law of internal complexity of the nonlinear economic system can be explained by the tools as one-dimensional bifurcation, two-dimensional bifurcation, the one-dimensional largest Lyapunov exponent, the two-dimensional largest Lyapunov exponent, and the basin of attraction. The system will transform from the stable Nash equilibrium state to the chaotic state via a flip bifurcation, which can be verified by the diagrams of one-dimensional bifurcation and the corresponding largest Lyapunov exponent. As it is known to all, the flip bifurcation is a classical route for the economic system to change the state from stable into chaotic. A new way of development for firms will generate when the economy develops to a certain stage, and there is only one way for firms to further develop the economy and they will face a dilemma in the choice of the new and old way. As a result, the route of period-doubling is generated in constant choice between firms in dilemma.
The increase of adjusting speed leads to the complexity of eventual behavior of system (10). This process of evolution presents the properties of “certainty” and “irregular development,” where “certainty” means that the way of bifurcation, that is to say, period-4 Neimark–Sacker bifurcation, and “irregular development” means that the interior structure of strange attractors. No matter from phenomenon development or interior mechanism, system (10) presents a high of certainty and complexity. That means the economic system makes mistakes in predicting, due to the influence of external condition, the development of economy cannot be forecasted anymore. Both the variation of parameters and the choice of initial condition will cause system (10) to generate multistability. With the increase of adjustment speed of firms, the properties of multistability are mainly reflected in the number and structure of coexisting attractors and global bifurcation of the system. There are four distinct coexisting phenomena displayed here, that is, periodic attractors coexist with periodic attractors, boundary equilibrium points coexist with strange attractors, periodic attractors coexist with chaotic attractors, and multiboundary chaotic attractors coexist. In the case of the boundary attractors coexisting with strange attractors with the layered structure, the largest Lyapunov exponent corresponds to strange attractors displaying a fluctuating and oscillating phenomenon as adjusting speed increases. This is a special phenomenon of coexistence and seldom reported by other researchers. As we all know, the boundary attractor corresponds to the fact that the duopoly market has degenerated into a monopoly market. That is, one of these two firms was forced out of the market. Therefore, the coexistence of the boundary attractor and internal attractor implies that not only the parameters can affect the market structure but also the initial state of market can affect the market structure.
The increase of the adjustment speed will lead to change of multistability, and at the same time, it will also cause “critical bifurcation” so as to form “ghost,” as well as the “global bifurcation” of attracting basin so as to form “hole.” The increasing of control factors are more beneficial for system (33) to maintain stability, delay the bifurcation point, and keep local stability of the Nash equilibrium point in the long run. However, as for the firms, any slight change of parameters may lead to a huge fluctuation in the market, which makes it difficult for the firms to keep pace with or catch the changing trend of the market. While this fluctuation can be eliminated or slowed down by introducing a control strategy, it is a benefit to a long-term stable development for the market.
Data Availability
No data were used to support this study.
Conflicts of Interest
The authors declare that there are no conflicts of interest.
Acknowledgments
This investigation was funded by the National Natural Science Foundation of China (Project no. 61463027), Young Scholars Science Foundation of Lanzhou Jiaotong University (Project no. 2015029), and Foundation of Humanities and Social Sciences from the Ministry of Education of China (Project no. 15YJC82007).