Abstract
In this paper, a dynamical two-stage game with R&D competition and joint profit maximization is built. The stability of all the equilibrium points is discussed through Jury condition, and the stability region of the Nash equilibrium point is then given. The influence of the parameters on the system is discussed, and we find that the firm can even benefit from chaos, when it has higher innovation efficiency and higher adjusting speed. And then the coexistence of multiple attractors is studied using basin of attraction. Our research result shows that the coexisting attractors can be observed in the two-parameter bifurcation diagram. At last, the boundary of feasible region, global bifurcations, and formation mechanism of fractal structure of attracting basin are analyzed through critical curves and noninvertible map theory.
1. Introduction
With the rapid development of economics, the competition among firms has become more and more fierce. The firms, which are the most important part of the market, must form a core competitiveness that is different from other firms, if they want to get a good survival and remain invincible in the market. Research and development (or R&D for short) has become the main driving force for the development of enterprises and it is also an important way for enterprises to obtain the core competitiveness. In fact, the R&D activities of enterprises can reduce production cost, improve the quality of products, expand the market share, and further improve its market competitiveness. However, in the real economic environment, R&D spillovers often arise when the R&D activities are carried out. In fact, the R&D spillovers are inevitable due to the information exchange of R&D among firms and the flow of human resources. In recent years, the problem of competition and cooperation in R&D investment has attracted more and more attention of entrepreneurs and economists. The classic research begins in the 1980s; a famous duopoly model with technology spillover was built by D’Aspremont and Jacquemin [1]. In this model, the competition between firms was divided into two stages. At the first stage, both firms decided their own R&D level, while, at the second stage, all the firms determined their outputs to the market. This famous model is also called AJ model by followers. Then, Kamien et al. [2] extended the model of D’Aspremont and Jacquemin to firms and got the KMZ model. In KMZ model, the forms of R&D are also extended to total competition, R&D cartel, research joint ventures (RJVs), and RJVs cartel.
Both AJ model and KMZ model are two-stage game models. Actually, the two-stage game is better than one-stage game for solving the uncertainty of the market and the interaction of competitors’ strategies when the R&D activity is considered. Following the AJ model and the KMZ model, Ziss [3] built a two-stage duopoly game model with R&D spillovers. Amir [4] studied the effects of R&D spillover and pointed out that the R&D cost of firms would increase with an increasing amount of technology outflow. Lambertinti et al. [5] extended the research results of Amir and proposed that the cooperative innovation will increase information sharing of market. Katsoulacos and Ulph [6] have led to a study of endogenous spillovers in R&D cooperation. Yin [7] studied the asymmetric R&D cooperation structure and pointed out that one of the important motivations for providing cooperative innovation is to increase the market share. Alallah [8] analyzed the influence of spillover rate on cooperative motivation and joint profit under asymmetrical R&D spillover. Bischi et al. [9, 10] established a two-stage dynamic game according to knowledge share and R&D competition in the market and analyzed the stability of the built model. Zhang et al. [11] considered a two-stage duopoly model with semi-collusion in production, and the complex chaotic behavior of the given model was researched in their works.
The influence of chaos theory on modern science is limited not only to nature science, but also to other fields, such as philosophy [12], social science [13], and economics [14]. The economic system is essentially nonlinear, so the chaotic phenomenon must exist in the economic system. In recent years, it has attracted the attention of a growing number of scholars to explore the evolution of the economic system and explain the complicated economic phenomena by using chaos theory. Gangopadhyay [15] built a dynamic model of enterprise merger, and the bifurcation behavior and the coexistence of multiple attractors in the built model were discussed. Canovas et al. [16] set up a piecewise smooth model for market competition, and the asymptotic dynamic behaviors of this model were discussed. Bischi et al. [17] proposed an inertial model of market share based on the optimal response and naïve expectation, the influence of the heterogeneity of the enterprise was also analyzed, and it was found that there are very rich local and global bifurcation phenomena in the built system. Li and Ma [18] considered a bounded rational dual-channel game, and the complex dynamical behavior of their model was simulated in their works. Cavalli et al. [19] built a nonlinear dynamical duopoly game by using bounded rationality and “LMA” (the abbreviation of “local monopolistic approximation”) adjustment strategies, and then they analyzed the coexistence of periodic points and the “lobes” basin of attraction was formed by periodic points. Based on the research results of [19], Cavalli and Naimzada [20] built a heterogeneous duopoly Cournot game with reduced rationality, isoelastic demand function, and linear total cost. And they found that the Nash equilibrium point may lose its stability through a flip bifurcation or a Neimark-Sacker bifurcation. In addition, a lot of researchers have discussed the complex dynamical behaviors of nonlinear oligopolies from different aspects, such as differentiated goods [21–25], heterogeneous firms [26–31], and delayed decisions [32, 33].
Another important issue for dynamical economical model is the coexistence of attractors. During the long run of dynamical economical system, the initial conditions play a very important role in determining the asymptotic behavior of the system. And thus, it is especially important to discuss the dependence of system on the initial conditions using basin of attraction and to analyze the evolution of attracting basin. In recent years, this topic has attracted a lot of attention of many scholars. Bischi et al. [34] discussed the coexistence of multi-attractors and cyclic attractors in a dynamic Cournot game. Agiza et al. [35] analyzed the multi-stability of a three-dimensional dynamical oligopoly model. Both the routes to complex attractors and the creation of basins with complex structures were studied in their paper. Cavalli and Naimzada [19] built an oligopoly model with different rational players, and then they explored the complex dynamical behaviors of the given model induced by coexistence of multiple attractors. Gori and Sodini [36] presented a nonlinear duopoly game with price competition, and the stability and coexistence of attractors were studied using the basin of attraction.
Through analyzing the above related research results, it is not difficult to find out that the R&D competition between duopoly firms that produce complementary products is seldom studied. However, the R&D competition between firms that produce complementary products is very common in the real market. A typical example is the cooperation and competition between “small program” and “APP” software. One firm developing “small program” and another firm developing “APP” always compete in R&D level as the development of them has a lot of commonalities. But they collaborate in the production level because they are complementary in functionality. They choose to cooperate in producing and selling their products, so as to maximize their joint profits. Therefore, a two-stage duopoly game under the background of R&D spillover, complementary products, and joint profits maximization is proposed in this present work.
In addition, the stability and the complex dynamical behavior in the built two-stage duopoly game is another topic in this research. As the stability of equilibrium point is often influenced by various parameters, the stability condition and stability region with various sets of parameters will be discussed in this article. And the complicated dynamical phenomena maybe arise when the equilibrium points lose their stability. The coexistence of many different attractors is the most important dynamical behaviors, which can be found in our model. The coexistence of attractors (i.e., multistability) indicates the long-run dependence of market dynamics on the initial conditions. Different initial conditions not only will result in a vast difference for firms after several games, but also may cause the coexistence of different situations. So by modeling the practical problem and analyzing the stability and multistability of the built model, it is beneficial for firms to make their decisions and predictions.
The plan of this article is given as follows. In the second section, a duopoly two-stage game based on R&D spillover, complementary goods, and maximization of joint profits is built. At the first stage of this game, two firms compete in the R&D level, but the R&D spillovers are permitted in this model. However, at the second stage, these two firms produce complementary goods and choose cooperation at the production level. In Section 3, the equilibrium points of this model are solved and the stability conditions are discussed. In Section 4, the stability of the Nash equilibrium point and the bifurcation of this system are analyzed using numerical simulation. And then, in Section 5, the coexistence among many attractors is investigated. The evolution and formulation mechanism of attracting basins are studied. Results are summarized in the last section.
2. The Model
Suppose there are two firms, labeled by (=1,2), in a market, which conduct R&D and produce complementary goods. In order to better analyze the R&D investment behavior of the two firms, we consider the following two-stage game. In the first stage of the game, both firms compete in the R&D level and choose the R&D efforts as their decision variables to reduce the production cost and improve the quality of their products, while, at the second stage of the game, both firms choose the market prices as their decision variables. And they choose cooperation, share the market information, and pursue the maximization of their joint profit at the production level. Obviously, the R&D investments of the two firms in the first stage will affect their profits in the second stage.
Let and be the price and the quantity of product of firm (), respectively. And the output of firm is determined by the prices of the goods produced by firm and its rival through an inverse demand function, which has the form [27]where represent the market sizes and , are the price sensitivity of consumers.
As we know, the spillover of knowledge is always inevitable in the process of R&D. That is, the R&D results of firm not only can reduce the production cost of its own, but also can reduce the production cost of its rival. So if we suppose and are the R&D efforts of firms and , respectively, then the effective marginal cost of firm with R&D efforts can be given aswhere represents the unit cost of produced goods without R&D efforts. If we assume that the two firms have equal knowledge levels when they are not engaged in R&D, that is, the productive efficiency of these two firms is the same, then we can let . And is related to the R&D spillover. If , it means that the technology between firms is completely confidential and the R&D spillover is ignored. And if , it means that R&D technology is fully shared among competing firms. In this paper, we only consider the intermediate condition .
However, the firms have to pay certain R&D expenditures when they conduct R&D activities. And the R&D expenditure can be assumed to be a quadratic function of R&D effort , which reflects the decreasing returns to R&D effort. And the R&D expenditure of firm can be given as where is the cost parameter of firm’s technological innovation, which indicates the efficiency of using or producing the unique technology or knowledge resources for an enterprise. The smaller the parameter , the stronger the innovation ability of firm .
The individual profit of firm can be given as . Through substituting (2) and (3) into , we can get
For convenience, the parameters and are set as the same in the following discussion. That is, . And the similar hypothesis is also imposed on parameters and ; it means that . The marginal profit of firm can be obtained by taking the first derivative of profits with respect to prices . That is,
Through solving the first order conditions , we can get the optimal prices of products produced by firm 1 and firm 2, respectively. And they are given as
By substituting (8) and (9) into (1), the optimal quantity of firm can be obtained,
Then the profit of firm can be obtained by substituting (8), (9), and (10) into (4) and (5), respectively.
Since, as we have assumed, the two firms choose collaboration in the production level, we consider the joint profit, namely, . By combining equation (11) and (12), the specific forms of the joint profit of these two firms can be given.
The marginal joint marginal profits of each firm can be obtained by taking the first derivative of joint profit with respect to the R&D efforts . That is,
As we have presumed above, the two firms choose to cooperate in the production level. Therefore, they will share the market information in order to improve their joint profit. However, the two firms will not share the information in the R&D level, as they are competitive in this process. And the decision-making of them cannot be completely rational, due to the lack of market information. So, the firms are not able to optimize their joint profit with respect to the R&D effort in one shot. They can only adjust their R&D strategies over time towards an optimal strategy along a direction of local estimation of the expected profits gradient, according to gradient adjustment mechanism. For a detailed description of the economic rationale of the gradient adjustment mechanism we refer to the description in [37].
According to the gradient adjustment mechanism, we know that if firm has positive marginal profit in period , it will increase its R&D effort in period t+1. On the contrary, if its marginal profit is negative, it will reduce its R&D effort. In the actual market, the firm dynamically adjusts its strategy at each discrete time, and the resulting dynamical system can be described by the following couple of equations:where is the R&D effort of firm at discrete time periods and is the speed of adjustment for firm . The dynamic game model can be obtained by substituting equations (14) and (15) into (16).
3. Stability Analysis of Equilibrium Point
In this section, the local stability around the equilibrium point will be analyzed through Jury condition [38]. The equilibrium points of system (17) can be solved by setting , that is,
Through solving the above nonlinear algebraic equations, four equilibrium points of system (17) can be obtained, where , . The equilibrium points , , and are called boundary equilibrium points as they locate at the coordinate axes. And the equilibrium point is called interior equilibrium point, which is also the unique Nash equilibrium. Obviously, if we assume the parameters meet the conditions as and , then the economic meaningfulness of these equilibria can be guaranteed.
It is obvious that all the coordinate axes are invariant sub-manifold of system (17), since implies . It means that any initial conditions starting from coordinate axis will stay on the coordinate axis forever after any time of iteration. And this situation can be regarded as the market competitors have degenerated from duopoly to monopoly. Or we can say that one of the firms has gone out of the market. The eventuality in which one firm exits the market and how the dynamics behave have been discussed by Gori et al. [38]. At this moment, the system (17) can be rewritten as a one-dimensional map,where represents the evolution operator of . The above equation is obtained by letting and in system (17), where and . Obviously, the above equation is conjugate to the famous Logistic map topologically through a proper linear transformation, which is given asand the corresponding value of is
As all the boundary equilibrium points locate on the coordinate axes of , the stability of all the boundary equilibrium points of system (17) is equivalent to the corresponding points of the Logistic map. According to the stability theory, the local stability of an equilibrium point can be determined by analyzing the eigenvalues of Jacobian matrix of system (17) evaluated at the equilibrium point on the complex plane. The Jacobian matrix of system (17) at any point can be given as
Proposition 1. The trivial equilibrium point is an unstable node.
Proof. At the trivial equilibrium point , the Jacobian matrix (23) becomes a diagonal matrix, which can be given asTherefore the eigenvalues of this diagonal matrix are the diagonal entries, which are and . According to the nonnegativity condition , we know that . And thus, both of the eigenvalues are larger than the unity, that is, and . And hence the trivial equilibrium point is an unstable node along the direction of the coordinate axes.
In the actual market, is on behalf of the zero R&D efforts of these two firms. It is obvious that this boundary equilibrium means that there is no investment in R&D for both firms. And if the initial condition is chosen from the coordinate axes of the plane (except the origin ), the system will eventually evolve to the equilibrium point or as is a repelling node. However, the system will stay at forever, if the initial condition is chosen as , even though the trivial equilibrium is unstable.
Proposition 2. The boundary equilibriums point is a saddle point when . But is an unstable node when .
Proof. At the boundary equilibrium point , the Jacobian matrix (23) becomes a lower triangular matrix, which is given asThe eigenvalues of this lower triangular matrix are given by the diagonal entries and . Through the assumption given above, we know that , so . Again with the nonnegative condition we know that , so will also be satisfied. However, when , the eigenvalue will be less than -1, so the boundary equilibrium point is an unstable node. But when , the eigenvalue will be located in the unit circle, so the boundary equilibrium point is an unstable saddle.
It is clear that the boundary equilibrium point presents that only firm 2 conducts R&D, while firm 1 is eliminated or left out of the competitive market and therefore it can be viewed as a single monopoly market. As the equilibrium point is symmetric with in the rectangular coordinate system, the stability analysis of is very similar to that of and is not repeated. In order to analyze the stability of the Nash equilibrium point, the following proposition is firstly given.
Proposition 3. The Neimark-Sacker bifurcation cannot occur at the Nash equilibrium point .
Proof. At the Nash equilibrium point , the Jacobian matrix becomes the following form:where and . Let and be the trace and the determinant of the Jacobian matrix , respectively. The following equation gives the characteristic equation of matrix :where the trace and the determinant of the Jacobian matrix are given as According to Eq. (28) and Eq. (29), the discriminant of the characteristic polynomial is given as As the parameters and meet and , the discriminant of the characteristic polynomial will be always nonnegative. It means that there is no complex root for this quadratic equation. And hence, the Neimark-Sacker bifurcation cannot occur.
Proposition 4. The unique Nash equilibrium is stable when the system parameters meet the following conditions:
Proof. According to the Jury condition [38], the Nash equilibrium point is local asymptotic stable if and only if all the following conditions are satisfied.The condition apparently holds as all the parameters , , , , and are all positive. And the condition is unnecessary as a consequence of Proposition 3. In order to solve the inequality , we let . And then the following threshold of can be gained through solving the given equation. And the inequality gives . Through the discussion given above, the stability region of Nash equilibrium point can be given as
In system (17), the equilibrium is stable when the parameters satisfy the stability condition, while, when one or more parameters violate the stability condition, the Nash equilibrium will lose its stability through some types of bifurcation. In the next section, the stability of the equilibrium point is studied using numerical simulation. Then the bifurcation with varying parameters and the complex dynamical behaviors of this system when the equilibrium point loses its stability are also analyzed.
4. Bifurcation Analysis Using Numerical Simulation
In this section, some numerical simulations will be conducted in order to analyze the stability of the equilibrium point and the complex dynamical behaviors of system (17).
Firstly, the system parameters are fixed as , , , , and , and the initial condition is chosen as . The stability region of the Nash equilibrium point is given in Figure 1. In Figure 1, the area with dark color represents the stability region of , while the white area represents the instability region of . We can also see from Figure 1 that the stability region is situated on the lower left corner of this figure. The shape of the stability region just likes a curved edge rectangle in the parameter plane of . And the boundary between stability region and instability region corresponds to a bifurcation curve. As there is no possibility for Neimark-Sacker bifurcation, this bifurcation curve refers to a period doubling bifurcation curve. It means that if the parameter or/and cross this curve into the white area from the dark area, then a period doubling bifurcation occurs. The endpoints of this bifurcation curve also can be solved through the equation , if we let and in this equation, respectively. And the coordinate axis of these two endpoints can be given as
In order to further analyze the influence of the system parameters on the stability region of the Nash equilibrium point on the plane , the stability curves with different values of are plotted in the plane . For example, if we fix all the other parameters as , , , and , the value of parameter is chosen as 0.1, 0.44, and 1.0, respectively. The stability curves in the parameter plane with different values of are shown in Figure 2. In Figure 2, the stability curves with black color, blue color, and red color represent , , and , respectively. And the area formed by these stability curves and coordinate axes are the stability regions with different values of . From Figure 2 we can see that the stability region increases very quickly as the value of increases from 0.1 to 0.44. However, the change of stability region is very small when increases from 0.44 to 1.0. In fact, there is hardly any change in the stability region as continues to increase further. And the shape of stability curve is also changing as varies; see Figure 2 for details.
The above analysis is mainly about the stability region of the Nash equilibrium point in the parameter plane . However, the complex dynamical behaviors outside the stability region have not been revealed. In order to analyze the complex dynamical behaviors of system (17) when the Nash equilibrium point loses its stability, the bifurcation diagram is employed in this article. However, the tool we used here is not only the one-parameter bifurcation diagram but also two-parameter bifurcation diagram, in order to analyze the complex dynamical behavior of system (17) in depth.
Figures 3(a) and 3(b) show the two-parameter bifurcation diagram with fixed parameter as , , , , and and initial conditions as in the parameter plane of . The different colors in these figures represent different number of points of the numerically evaluated attractor, which can be seen from the color bar in the right of these figures. In Figure 3(a) deep yellow represents the period-1 attractor, which is the stability region of the Nash equilibrium point. Furthermore, the deep yellow region in Figure 3(a) also corresponds to the stability region in Figure 1. The dark green color represents the region of period-2 attractor. And the yellow color represents the region of period-4 attractor, and so on. From Figures 3(a) and 3(b), we can see that the system will lose stability through period doubling bifurcation as one/both of the parameters , increase. In particular, in the black region the trajectories either converge toward a large period attractor, a chaotic or quasi-periodic one or they diverge. It is worth pointing out that the quasi-periodic attractor is a special attractor, which will arise after the Neimark-Sacker bifurcation. And the largest Lyapunov exponent equals zeros, if the system is in a quasi-periodic state.
(a)
(b)
(c)
(d)
Figure 3(b) is a partial enlargement of Figure 3(a). From Figure 3(b), we can find some fractal structures. And symmetry of these figures can also be detected. It also means that the symmetry of parameters and also exists in system (17). Therefore, only the changes of dynamical behavior with parameter are studied in this article. There is no need for us to study the effects of parameter on system (17) due to the symmetry.
If we fix the value of as in Figure 3(a), the detailed dynamical behavior of system (17) can be found in Figures 4(a) and 4(b). In Figures 4(a) and 4(b), the one-parameter bifurcation diagram and the corresponding largest Lyapunov exponent are given, where the bifurcation parameter is chosen as the speed of adjustment of firm 1. And the parameters chosen in Figure 4 are same with those of Figure 3(a). From Figure 4(a), we can see that the Nash equilibrium point is stable when . At about , the Nash equilibrium point loses its stability through a period doubling bifurcation, and a stable state of period-2 appears. As the value of further increases, the stable period-2 orbit loses its stability again, and a period-4 cycle arises, and so on. At last, the system enters into chaos. However, some periodic windows can also be found in the chaotic region, which can be seen from Figure 4(a). Through the above analysis, a common conclusion can be gained that the speed of adjustment plays a vital role in system (17). In other words, smaller speed of adjustment means more stable R&D investments and more stable market environment. In fact, the R&D investments for firms can exhibit a more stable behavior if the reaction speeds for them to the market information are slow enough. Another useful tool to analyze the complex dynamical behavior is the largest Lyapunov exponent. The largest Lyapunov exponent diagram can be used to show the degree of separation of nonlinear system variables running over time, and it can also be used to prove whether the system is chaotic. When the system is in a stable state, the largest Lyapunov exponent is less than zero. While the system enters chaos, its largest Lyapunov exponent is greater than zero. In Figure 4(b), the largest Lyapunov exponent diagram corresponding to Figure 4(a) is plotted. The evolution process of system (17) is confirmed again.
(a)
(b)
(c)
(d)
Moreover, the role of parameter is also very important for the system (17). If we fix the parameter sets as , , , , and , the two-parameter bifurcation diagrams with parameters and are shown in Figures 3(c) and 3(d). In Figures 3(c) and 3(d), the initial conditions are chosen as and . From Figure 3(c), we can see that the dynamical behaviors of system (17) are very complex, when the value of parameter meets (to be honest, this value is not very accurate), while, if , we can see that the system is in a chaotic state when the value of parameter is relatively small. And then the stable periodic states will arise. However, the system will lose stability again, when the value of is large enough. And the escaping trajectories will appear. Figure 3(d) is the partial enlargement of Figure 3(c). Although Figure 3(c) looks ugly, Figure 3(d) does show a good fractal structure. Furthermore, if we fix the value of as in Figures 3(c) and 3(d), the one-parameter bifurcation diagram is given in Figure 4(c). From Figure 4(c), we can see that the system is in a chaotic state when . And then a multi-periodic window appears. The chaotic window reappears, subsequently. But when the value of meets , the system enters a multi-periodic state again. And then a secondary Neimark-Sacker bifurcation occurs at . The Nash equilibrium point becomes stable at . Finally, the Nash equilibrium point loses its stability again at , and the escaping trajectories arise. The escaping trajectory cannot be detected, as the value(s) of or/and will go to infinity. Such changes can also be seen from the corresponding largest Lyapunov exponent diagram, which is given in Figure 4(d).
In order to further explain the influence of parameter on system (17), the two-dimensional bifurcation diagram with different parameter sets is given in Figure 5. The parameters are chosen as , , , and , which are same with those of Figure 3(a). But the value of parameter is different from Figure 3(a). In Figures 5(a) and 5(b), the value of is chosen as , and Figure 5(b) is the partial enlargement of Figure 5(a). However, the value of is chosen as in Figures 5(c) and 5(d), and Figure 5(d) is the corresponding partial enlargement of Figure 5(c). By comparing Figures 3(a) and 5, we can see that the Nash equilibrium point always loses its stability through a period doubling bifurcation, and then a period-2 cycle appears. After that the system enters to chaotic state through a series of period doubling bifurcation. The detailed bifurcation process can be found in Figure 6(a), where the parameters are same with Figure 5(a) except that the value of is fixed as . From Figure 6(a), a clear secondary Neimark-Sacker bifurcation can be found. In addition, we can also find the existence of periodic windows on the way to chaos.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
(e)
(f)
From Figures 3(a) and 5, we can also see that the parameter window into quasi-periodic state becomes very narrow when increases to some certain value. When , the parameter window of Neimark-Sacker bifurcation is even hard to observe; see Figure 5(c). But this phenomenon can be revealed from the partial enlargement of Figure 5(c); see Figure 5(d).
The symmetry structure and fractal structure can also be observed from Figures 3(a) and 5. As the value of increases, the fractal structure of these figures has also changed. But the symmetry structure of all the figures will never change. This phenomenon also directly reflects the symmetry of parameters and . In addition, we can also find that there are a lot of sparse black points, which have covered the light blue area and the yellow area in the upper right corner of Figure 5(d). Here, these sparse black points are called noisy points, as they have affected the beauty of Figure 5(d). However, the presence of the noisy points often implies the coexistence of multiple attractors. And the problem of coexistence of attractors will be discussed in the following section.
The parameter plays a very important role in the analysis of high-tech enterprises and it can reflect the R&D capability of enterprises. The smaller value of parameter means the stronger innovation ability of enterprises. The following analysis is mainly focused on the influence of parameter on the dynamical behaviors and profits of firms in the given model. To be specific, the equilibrium profits and the average profits of firm 1 and firm 2 will be discussed in order to enclose the relation between parameter and the behaviors of firms. The equilibrium profits of firm can be gained by substituting the value of into (11) and (12), respectively. The specific mathematical expressions of equilibrium profits will be given as
And the average profits of firm 1 and firm 2 can be calculated through the mean value of the profit of firm at different period. The detailed expressions of average profits of firm 1 and firm 2 can be written as where is the number of iteration. In the following numerical simulation, the value of is chosen as . Figure 6 gives the one-parameter bifurcation diagrams, the corresponding equilibrium profits, and average profits of firm 1 and firm 2, where the blue color corresponds to firm 1, while the red color represents firm 2. In Figure 6, the fixed parameters are given as , , , , and . However, the value of is fixed as , , and , respectively. From the value of , we can see that the R&D efforts of firm 1 and firm 2 are equal at the Nash equilibrium point. And thus the equilibrium profits are same as well for firm 1 and firm 2. This phenomenon can also be seen from Figure 6. However, the R&D efforts and average profits of firm 1 and firm 2 are no longer equal any more, once the Nash equilibrium point loses its stability, as increases. Through comparing Figures 6(b), 6(d), and 6(f), we can find that the average profits will fall when the value of increases. It is easy to understand that larger value of represents a lower R&D capability, so the profits will fall with increasing of . And the amplitude of fluctuation of R&D efforts is smaller for firm with slower speed of adjustment; see Figures 6(a), 6(c), and 6(e).
From Figure 6, we can see that the average profits would be smaller than the equilibrium profits when the Nash equilibrium point loses stability. However, the average profits are close to the equilibrium profits, if the firm’s speed of adjustment is large enough; see Figures 6(d) and 6(f). It is worth pointing out that the average profit is even higher than the equilibrium profit for the firm with large speed of adjustment, once the system enters into chaos; see the blue curves in Figures 6(d) and 6(f), while the average profit of another firm fluctuates very fierce and is far below its equilibrium profit; see the red curves in Figures 6(d) and 6(f). This phenomenon can be concluded as chaos is usually harmful to the entire market, but firms with large speed of adjustment sometimes can even benefit from chaos.
5. Multistability
As we have noticed, there are some noisy points in the upper right corner of Figure 5(d). These noisy points are probably caused by the coexistence of multiple attractors as we have mentioned above. To give a more explicit picture about this problem, the one-parameter bifurcation diagrams with different initial conditions are plotted if we fix the parameters as , , , , and and choose the parameter as the bifurcation parameter. The one-parameter bifurcation diagrams with initial conditions and are given in Figures 7(a) and 7(b), respectively. The differences between these two figures are very clear. In order to highlight these differences, the bifurcation diagrams are plotted with different colors. In Figure 7(a), a period doubling bifurcation firstly happens at , and then a Neimark-Sacker bifurcation arises at . Later, a period-2 point appears again. After that the system enters into chaos through a period doubling bifurcation. In Figure 7(b), the dynamical behavior of system (17) is relatively simple, but a transcritical bifurcation can be observed.
(a)
(b)
In Figure 7, different colors represent different meanings. From left to right, there is no coexistence attractor in the red region. Blue corresponds to the coexistence of attractors. Three different forms of coexistence of multiple attractors, which are the coexistence of two period-2 attractors, the coexistence of a period-2 attractor with a 2-cyclic attractor, and the coexistence of a period-2 attractor with a chaotic attractor, can be found in the blue region. However, the coexistence of multiple attractors will disappear again in the green area.
The multi-stability in economic system has a very important research meaning. In fact, the choices of initial conditions are very crucial for firms in the real market. As the system we built here can be regarded as a repeated process, the state of firms at time can be considered as the initial conditions for firms at time t+1. Furthermore, once the initial conditions are given for system (17), the long and asymptotic behaviors of this system after many iterations can be predicted. All the time, a recognized tool to analyze the coexistence of attractors is the basin of attraction. And thus, the basin of attraction is employed in our article. Firstly, the definition of basin of attraction is given here.
Suppose there is an attractor of a nonlinear system, which can be expressed by . Then all the basins of attractor can be written as , where is the immediate basin of . Obviously, , and the points mapped to after a finite number of iterations also belong to . represents the set of rank-n preimages of . The basin of attraction can be divided into the attractor, domain of attraction, and domain of escape. And the domain of attraction refers to the set of initial conditions that converge to the same attractor after a series of iterations. The domain of attraction can be divided into connected sets and nonconnected sets. The domain of escape refers to the set of initial conditions that diverge to infinity after a series of iterations. If we fix the parameters as , , , , , , and , a typical basin of attraction of system (17) is given in Figure 8. In Figure 8, the domain of attraction and the escaping domain are represented by yellow area and dark blue area, respectively.
As it can be seen from Figure 8, the shape of the basin of Nash equilibrium point is a quadrilateral , where is the origin of coordinate axis and the rest vertices of this quadrilateral , , and are the rank-1 preimages of the origin . It is clear that the points of and are located on the horizontal axis and the vertical axis , respectively. And their coordinates can be gained through solving the algebraic equation in formula (20), which are given as
For the convenience of analysis, the line-segments and are labeled as and , i.e., . And the rank-1 preimages of , are labeled as and , respectively. The rank-1 preimages of , can be understood as a set of all the preimages of the points on . In order to get the coordinates of points located on , it is supposed that the coordinate of any point on can be represented as ; then the preimages of should satisfy the following equations:
From the second equation of (39) we can find that the rank-1 preimage of should meet or the linear function given below:The function represents that is one of the rank-1 preimages of itself. However, the linear function (40) gives another rank-1 preimage of , that is, . Through similar method, the coordinates of points located on also can be solved. And the following linear function can determine the straight line :
It is obvious that the line determined by (40) and the line determined by (41) have a unique intersection , which is one of the rank-1 preimages of the origin . The coordinate of point can be gained through solving the combination of (40) and (41), which can be represented as
A more intuitive explanation of the points , , , and these straight lines , , and can be found in Figure 8. It is clear that the boundaries of the attracting domain of the Nash equilibrium point are constituted by these lines. And the points , , , and are just the intersections of these straight lines. From the detailed expressions of (40) and (41), we find that the structure of domain of attraction is related to the system parameters. Any change in these parameters will lead to changes of structure of basins.
For example, we fix all the parameters as , , , , and , which are same with the parameters given in Figure 7. And the parameter varies. Firstly, the parameter is fixed as 0.53325, and the corresponding basins can be found in Figure 9(a). At this very moment, a chaotic attractor coexists with a period-2 cycle. The attracting domain of the chaotic attractor is represented by yellow color, and the attracting domain of the period-2 cycle is represented by light green, while the dark blue area means the domain of escaping. And the shape of Figure 9(a) is just like a quadrilateral “sieve” with fractal structure. As parameter further increases to 0.55675, the chaotic attractor enlarges, but the fractal structure of the attracting basin has hardly changed. The chaotic attractor has almost contacted with its attracting domain; see Figure 9(b). It means that a global bifurcation called “contact bifurcation” will arise. In Figure 9(c), the parameter is fixed as 0.557, we can find that both the chaotic attractor and its corresponding domain of attraction have been destroyed, and the destroyed attractor will fill the whole domain of attraction. This phenomenon can also be observed from the time series plot, which is given in Figure 10. In Figure 10, the values of parameters are same with those in Figure 9(c) and the initial conditions are chosen as . We can find that the evolution is disorderly in Figure 10. Comparing Figures 9(b) and 9(c), it can be found that the truth of the matter is that the attracting basin of the “dead” chaotic attractor is also occupied by the period-2 attractor. We refer to the long-run transition as the “ghost” of the destroyed attractor; see Bischi, etc. [17, 39] for more detail. Furthermore, if we choose as 0.58117, it can be found that the “ghost” of the destroyed attractor still exists, and the coexistence of a period-4 point (red points in Figure 9(d)) and an 8-piece quasi-periodic attractor (black cycle in Figure 9(d)) can be found. But the attracting domain of the period-4 points is confined to the diagonal. As further increases to 0.58225, the 8-piece quasi-periodic attractor has evolved into a 4-piece quasi-periodic attractor. The coexisting period-4 attractor is hard to distinguish from the quasi-periodic attractor; see Figure 9(e). As continues to increase, the 4-piece quasi-periodic attractor becomes a 4-piece chaotic attractor. And some disconnected sets arise in the attracting basins of these two attractors; see Figure 9(f). The 4-piece chaotic attractor starts to contact with its attracting domain. This means that global bifurcation occurs, and the structure of attracting domain will change once again; see Figure 9(g). In Figure 9(g), we fix the value of as 0.6025, and the coexistence of a 4-piece chaotic attractor and a multi-periodic cycle in the diagonal can be found. The attracting basin of the multi-periodic cycle is only the diagonal line. With a further increase of the parameter , the 4-piece chaotic attractor starts to merge into a one-piece chaotic attractor at ; see Figure 9(h). However, a weak chaotic attractor is beginning to appear on the diagonal, too. And its attracting domain is just the diagonal line.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
In Figure 9, the multistability of system (17) with has been discussed. Then, the cases with will be discussed. Firstly, the parameter sets are fixed as , , , , and , and the parameters and are chosen as the bifurcation parameters. In Figure 11(a), a two-parameter bifurcation diagram is plotted. And the local enlargement of Figure 11(a) is given in Figure 11(b). The noisy points can be found in the lower right corner of Figure 11(b). In the previous section, we have mentioned that these noisy points are related with the coexistence of attractors. To throw further light on this problem, the values of and are chosen as 0.524 and 0.2037, respectively. This set of parameter corresponds to the overlapping region of fluorescein and turquoise in Figure 11(b). The basins of attraction at this set of parameters are given in Figure 11(c). From Figure 11(c), we can find that there are two attractors coexisting: one is a period-2 cycle (labeled by ), and the other is a period-6 cycle (labeled by ). The immediate basin of is the area where locates, while the total basin of , which can be written as , is the entire region with light blue. In fact, the light blue region occupies the vast majority of the feasible region; see Figure 11(c). Accordingly, the total basin of is very sparse, which can be represented as and is expressed with yellow color in Figure 11(c). The feasible region can be represented as .
(a)
(b)
(c)
(d)
If the values of and are changed to 0.5665 and 0.2037, respectively, other coexisting attractors can be found, see Figure 11(d), where a period-4 cycle (labeled by ) coexists with a chaotic attractor constituted by 7-pieces (labeled by ). In Figure 11(d), the total basin of is the area of light blue, and the total basin of is the area of yellow. In Figure 11(d), the yellow region is further increased, while the light blue region is reduced compared with Figure 11(c). The feasible regions are the total basins of these two attractors. The shape of Figure 11(d) is similar to that of Figure 11(c), and it is just like a knife coin in ancient China.
It is generally known that the critical set is the primary tool to analyze the topological structure of attracting basin in noninvertible maps. If we let , in system (17), and solve the algebraic equations about and , which are given as Then, the rank-1 preimages of point can be obtained. Actually, the preimages of point may be more than one, as the nonlinear equations (43) are constituted by two quadratic equations. And therefore, the system can be regarded as a two-dimensional noninvertible map, which can be abbreviated as . The noninvertible map represents many-to-one. It means that the image of two distinct points on the two-dimensional plane may be same. That is, even if , may be equal to , too. In fact, if we solve the nonlinear algebraic equations (43), there are maybe four, maybe two, or even no solutions for equations (43) in the real field. And the critical curves can be used to divide the phase space into the regions called , , and . The symbol represents the region in phase space, where the points have preimages. The boundary of these regions, also called critical curves, can be denoted as . Suppose that are the preimages sets of points on , that is, . And can be defined as a locus of the points where the determination of Jacobian matrix vanishes, that is,
As the Jacobian matrix of system (17) has given in (23), thus the points on can be determined by the following equation:where , , , , , and .
It can be found that there are two branches in , i.e., and . Therefore, there are also two branches in , that is, and . A more intuitive image of and can be found in Figure 12. The sketch of has been given in Figure 12(a), while the sketch of has been given in Figure 12(b). In Figure 12(b), we find that the branch separates from , while the other branch separates from . For example, the origin point has four preimages as and itself, so the origin .
(a)
(b)
In the above discussion, the coexistence of attractors in system (17) has been shown. However, another interesting phenomenon in the basin of attraction called “hole” by Bischi, etc., [40] can be found in system (17). The formation mechanism of these holes and their evolutions will be discussed in the following through the critical curves, which we have given above.
First of all, the two-parameter bifurcation diagram is given in Figure 13, where the parameters are fixed as , , , , and . The specific meaning of the color in Figure 13 can be referred to the color bar in the right side of these figures. Figure 13(b) is the local enlargement of Figure 13(a). The obvious fractal structure can be found in these two figures. But there are also a lot of noisy points in these figures, which have already affected the beauty of these figures. Similarly, these noisy points also give us the information of the coexistence of multiple attractors.
(a)
(b)
As a matter of fact, if we fix and as 1.572 and 1.4893, the basin of attraction is given in Figure 14(a). The rough boundaries of feasible region can be determined using Eqs. (40) and (41), but there is also a small gap in the feasible region; see the topmost of Figure 14(a) for detail. And the boundary of basin is tangent to the critical curve at the gap. No coexisting of attractors can be found at such a situation. Only one chaotic attractor exists in the feasible region. That is, the feasible region is just the attracting basin of that chaotic attractor. The critical curve is also the boundary of that chaotic attractor. It means that the chaotic attractor has contacted with the boundary of its attracting domain; then a contact bifurcation will happen soon.
(a)
(b)
(c)
(d)
(e)
(f)
When increases to 1.49, the contact bifurcation happens and the chaotic attractor disappears. And the feasible region is full of the “ghost” of the disappeared attractor; see Figure 14(b). At the same time, the critical curve starts intersecting with the gap, and a very tiny hole in the infeasible region arises. As the critical curve separates the phase space into and , so the closed hole should have two preimages. And the two preimages of , which are represented as and , form another two holes in the feasible region; see Figure 14(b) for detail. The holes and lie in the region, which is different from the research result of other scholars [40]. Therefore, there should be four preimages for the hole and , respectively. One of the preimages of hole and is given as and , respectively. And other preimages are too small to see.
As further increases to 1.52, the main hole becomes larger than that of Figure 14(b), and more holes arise. These holes can be regarded as the rank-k (k=0,1,2,…) preimages of ; see Figure 14(c). By comparison with Figure 14(b), the “ghost” of chaotic attractor has disappeared and a new period-2 cycle emerges in Figure 14(c). With a further increase of to 1.62, the main hole is getting larger, and the preimages of are going to get bigger, too. In the meanwhile, the attractor has evolved from a period-2 cycle to a period-4 cycle; see Figure 14(d). It is very clear that the holes in the feasible region will grow bigger and bigger and burst into larger holes, as the main hole is no longer closed; see Figure 14(e). In Figure 14(e), we can see that the number of holes has been reduced, but the scale of these holes has increased. Furthermore, a chaotic attractor arises when . However, the newborn chaotic attractor is going to die, as there is no place for it to live. And then a contact bifurcation happens at , the newborn chaotic attractor has already died, and only its ghost still remains in the feasible region; see Figure 14(f).
6. Conclusion
In this paper, a two-stage dynamical duopoly game of R&D competition between two high-tech enterprises is studied. At the first stage, both firms that conduct R&D so as to reduce their production cost and improve the quality of products compete in the R&D level. And at the second stage, both firms choose cooperation and pursue the maximization of joint profit. And then, the existence and stability of all the equilibrium points are discussed. It has been proved that the boundary equilibrium points are always unstable, and the Nash equilibrium point is stable only when the parameters meet some conditions. The stability region of Nash equilibrium point is derived through Jury condition. We find that the Neimark-Sacker bifurcation cannot occur as both eigenvalues of the Jacobian matrix at Nash equilibrium point are always real numbers. And a numerical simulation is employed in order to check the theoretical results. The stability region is plotted numerically with different values of parameter . We find that the stability region is easy to be influenced by the parameter . The stability region is bigger when the value of is larger. But the stability region will not change when the parameter is large enough. We also find that the parameter window into quasi-periodic state becomes very narrow when increases to some certain value. In addition, the influence of adjusting speed on the system is also discussed in this paper.
Another topic of this paper is the multi-stability. The multistability always implies path-dependence, which means the long run behavior of firms is deeply influenced by historical accidents. That is, a small perturbation of initial conditions will affect the evolution of the system. And hence, the choice of initial conditions plays a very important role in the long run of the built system. For the sake of analyzing the effect of initial conditions on the long run evolution of the system, the basin of attraction is employed. In this research, the two-parameter bifurcation diagram is used to find the coexistence of multiple attractors. We found that there exist noisy points in the two-parameter bifurcation diagram, when the phenomenon of coexistence of multiple attractors arises. And then the coexistence of multi-attractors is studied with two different cases: and , respectively. Our research result shows that the topological structures of these two cases are totally different. The global bifurcation, such as “contact” bifurcation, and the phenomenon of “holes” in attracting basins have been found through numerical simulation. And at last, the formation mechanisms of these phenomena have been analyzed through critical curves and noninvertible maps.
Data Availability
All the authors solemnly declare that there is no data used to support the findings of this study.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This project is supported by the Talents’ Innovation and Entrepreneurship Project of Lanzhou City (Project 2015-RC-3), the Young Scholars Science Foundation of Lanzhou Jiaotong University (Project No. 2015029), and the Foundation of Humanities and Social Sciences from the Ministry of Education of China (Project No. 15YJC820007).