Local and global analysis of a nonlinear duopoly game with heterogeneous firms

*Correspondence: saskar@ksu.edu.sa; s.e.a.askar@hotmail.co.uk 1College of Science, Department of Statistics and Operations Research, King Saud University, P.O. Box 2455, 11451 Riyadh, Saudi Arabia 2Faculty of Science, Department of Mathematics, Mansoura University, 35516 Mansoura, Egypt Abstract In this paper, we introduce a nonlinear duopoly game whose players are heterogeneous and their inverse demand functions are derived from a more general isoelastic demand. The game is modeled by a discrete time dynamic system whose Nash equilibrium point is unique. The conditions of local stability of Nash point are calculated. It becomes unstable via two types of bifurcations: flip and Neimark–Sacker. Some local and global numerical investigations are performed to show the dynamic behavior of game’s system. We show that the system is noninvertible and belongs to Z2 – Z0 type. We also show some multistability aspects of the system including basins of attraction and regions known as lobes.


Introduction
In economy, literature has reported several studies about the dynamic characteristics of duopoly games. As known, a duopoly game possesses only two competing firms (or players) whose strategies may be quantities or prices. Updating the outputs of those firms depends on the kind of mechanisms they have adopted. There are different types of updating mechanisms that have been reported in the literature. Examples of those mechanisms are gradient based mechanism such as bounded rationality and best reaction approach. There are also naive mechanism and local monopolistic approximation (LMA) approach. Having each firm to update the outputs of the mechanism it adopts, the game's players are able to calculate Nash equilibrium and study its complex dynamic characteristics. Studying these dynamic characteristics shows the routes on where Nash equilibrium may be unstable. These routes are due to different types of bifurcation: flip and Neimark-Sacker bifurcations.
There are several contributions that have reported such dynamic characteristics in the literature. We start here by mentioning the nice work of Puu [1]. In [1], Puu proposed a duopoly game based on an isoelastic demand function which has been derived from Cobb-Douglas utility. He showed that the Cournot equilibrium point of the two competing players can lose its stability through a period doubling bifurcation. This work of Puu opened the route to different kinds of contributions about duopoly and oligopoly compe-titions. For instance, the stability of an economic equilibrium problem was addressed in [2]. Tuinstra in [2] introduced a simple adjustment process competing firms can adopt to set their outputs. This process requires information about the amount sold in the market against prices from the previous period. In [3], an oligopoly game was introduced with competitors having incomplete information about the demand function and adopting the LMA mechanism. It has been proved that competitors adopting this mechanism may have a converging Nash equilibrium point. Naimzada et al. applied the LMA in a monopolistic quantity setting game in [4] and more recently in [5]. In the LMA approach, the competitors do not have any information about the demand function of the market but instead they guess it as linear.
Another well-known mechanism is the gradient-based mechanism such as the bounded rationality. Competitors adopting such mechanism do not need to know the demand and cost functions but instead they make an estimation of the marginal profit in order to update their outputs in the next period of time. The bounded rationality has been intensively adopted in many works in literature. Tramontana et al. [6] analyzed a game of triopoly whose firms are heterogeneous with an isoelastic demand function and adopted bounded rationality and another mechanism to update their quantities in the next period of time. A multi-team Cournot game based on bounded rationality has been introduced and studied by Ahmed et al. in [7]. Heterogeneous oligopolies were considered and studied in many works in the literature such as Askar et al. who analyzed in [8] a Cournot duopoly game based on Cobb-Douglas preferences using different types of decision mechanisms including bounded rationality. More recent works in this direction by Askar have been reported in [9][10][11]. In the sense of heterogeneous framework, literature has reported many works such as Leonard and Nishimura [12], Zhang et al. [13], Ma et al. [14], Peng et al. [15], Ma et al. [16], Elsadany [17], and Naimzada et al. [18]. For extra reading, readers are advised to see [19,20].
The current paper belongs to the above category of heterogeneous oligopoly. It introduces a more general isoelastic utility function than that used in [21] and [5]; however, it has been derived from Cobb-Douglas utility. Even though the decisional mechanisms adopted by both competitors in this work are similar to those in [21], our work generalizes the work studied in [21]. The main results in this paper focus on studying the local stability of the game's equilibrium point. This includes detecting the routes by where the equilibrium point becomes unstable; the global analysis of the game's model. The results show that the region of stability of Nash point becomes larger for both competitors when their marginal costs are more comparable. This region is reduced for both firms when their marginal costs are sufficiently different. This leads to two different types of bifurcations by which the Nash equilibrium point can be destabilized. Those types are flip and Neimark-Sacker bifurcations. Moreover, the nonlinearity and noninvertibility of the system describing the current game give rise to important dynamic characteristics which are not provided by the local analysis. We show that the basins of attraction of periodic cycles can be quite complicated and are characterized by popular regions known in the literature as lobes.
The work in this paper is summarized as follows. In Sect. 2, the system describing the game is modeled. In Sect. 3, the stability conditions of the equilibrium point are discussed. In Sect. 4, numerical simulation experiments that include local and global analysis are carried out. In Sect. 5, we study the noninvertible aspect of the system. Conclusion is finally given in Sect. 6.

Model
Let us consider the following inverse demand function: where a is a constant and Q = 2 i=1 q i , q i ≥ 0; i = 1, 2. If a = 0, we get a demand function that is derived from a Cobb-Douglas utility. If Q = 0, which means market is not provided by quantities, the price becomes constant and equal to 1 a . We suppose that the industrial market consists of two firms (or players) that produce goods denoted by q 1 and q 2 . Both firms adopt linear costs with constant marginal c i , i = 1, 2. So the profits of firms are and the marginal profits become where q -i refers to the output of the other competitor. Since both players want to maximize their profits, the Nash equilibrium E = (q 1 ,q 2 ) for this game can be calculated by vanishing the marginal profits as follows: Now, we consider two heterogeneous players. The first player adopts a gradient based mechanism to update its production outputs. This mechanism is known with bounded rationality and depends on the change occurring in the marginal profit where it increases or decreases. It is governed by the following equation: The parameter k is known as the speed of adjustment. While the second player is considered as a naive player and he chooses to compete in the market using the best reaction function given by [21] Now the proposed game can be described by the following discrete-time dynamical system: It is clear that system (7) is undefined at the line, q 2 = -aq 1 , and this may give some interesting characteristics that will be discussed later. Furthermore, the system may give negative trajectories due to the second equation of (7) as will be shown in the numerical simulation. For this reason we focus on parameter values, for which q 2 stays positive for all iterations, otherwise max{q 2 (t + 1), 0} is taken. Moreover, if the initial condition q 10 > 1-ac 2 c 2 , the first iteration of the second equation of (7) gives negative values, while k > 1 c 1 gives negative values for the first equation as long as q 20 goes to infinity. We should highlight that the model discussed in [21] becomes a special case of our proposed system describing the current game by setting a = 0. It is also more general than the model studied in [5]. Studying the system's qualitative behaviors requires recalling the Jacobian matrix in order to find out the characteristics of its eigenvalues and hence discover at which types of bifurcations the Nash equilibrium point E becomes unstable. The next section discusses analytically the stability conditions of Nash point.

Local stability analysis
The Jacobian at E takes the form whose trace and determinant are given by τ = 1 + k (a +q 2 )(a -q 1 +q 2 ) (a +q 1 +q 2 ) 3 c 1 , δ = kq 1 (a -q 1 +q 2 )(1 -2 c 2 (a +q 1 )) 2 c 2 (a +q 1 )(a +q 1 +q 2 ) 3 (9) and then the characteristic polynomial becomes It is known that the local stability of Nash point is governed by the eigenvalues λ 1 and λ 2 that are the solutions of (10). These solutions must be in the unit circle |λ 1,2 | < 1. In our case it is hard to get these eigenvalues, and then we recall Jury conditions [22] that take the following form: Relations (11) have the following properties: • The Nash point is locally asymptotically stable provided that g(1) > 0, g(-1) > 0 and > 0. • The Nash point becomes unstable due to flip bifurcation if g(1) > 0, > 0, and g(-1) = 0. It means we have two real eigenvalues passing through -1. • The Nash point becomes unstable due to transcritical or fold bifurcation if g(-1) > 0, > 0, and g(1) = 0. It means we have two real eigenvalues passing through 1. • The Nash point becomes unstable due to Neimark-Sacker bifurcation if g(1) > 0, g(-1) > 0, and < 0. It means we have two complex and conjugate eigenvalues whose modulus is passing through 1.
Proof The proof is straightforward. Substituting (9) in (11) completes the proof.

Global analysis via numerical simulation
Due to the complicated expressions of Nash point and Jury conditions obtained in (12), we perform some numerical experiments in order to get more information about the stability of Nash point to predict the future evolution of the system. We start our experiments by assuming the parameter set, c 1 = 0.2, c 2 = 0.3, a = 0.9, k = 3, with the initial datum (q 0,1 , q 0,2 ) = (0.11, 0.12). At this set, the Nash point is positive and becomes (1.2440, 0.52933), and then the Jacobian matrix (8) is given by which has two real eigenvalues λ 1 ≈ 0.47057 and λ 2 ≈ -0.02897; that means the Nash point is stable as |λ 1,2 | < 1. As one can see, conditions (12) depend on three parameters, which are c 1 , c 2 , and k. Now, we analyze the influences of the parameter k on the stability of Nash point and on the dynamics of system (7) as we take it as the bifurcation parameter. This makes us to start with the 1D bifurcation diagram. As in Fig. 1a, the Nash point is stable for all k until k ≈ 10.49, it gets unstable due to the coexistence of flip bifurcation. Substituting c 1 = 0.2, c 2 = 0.3, a = 0.9, k = 10.49 in (11), we get g(1) > 0, > 0 and g(-1) < 0, which confirm the appearance of flip bifurcation, and hence E gets unstable. As the  Figure 1b displays this period which has quite complicated basins of attraction with lobes. This cycle is positive but its basin contains negative points because the map is undefined at q 2 = -aq 1 as mentioned before. The gray color denotes the divergent points, while the other colors are for the basins of period-2 cycle. Increasing k to 12.29 gives rise to a period-4 cycle represented by square around the Nash point. It is a positive cycle with quite complicated basins as shown in Fig. 1c. The gray color in Fig. 1c refers to the divergent and infeasible points in the phase plane. The basins of attraction become more complicated as k increases. This is plotted in Fig. 1d where period-8 cycle takes place at k = 10.64. Further increasing in the parameter k gives rise to higher periodic cycles with complicated basins of attraction until k = 12.747 on where the dynamic of map (7) changes to four unconnected chaotic areas. Figure 1e shows the phase plane of these chaotic areas. They continue to appear until k approaches 12.757, where a period-12 cycle emerges. After that these four areas gather together to form a one-piece chaotic attractor at k = 13.8 as depicted in Fig. 1f. In addition, the numerical experiments show that the dynamic behavior of the map is affected when we increase the parameter a. Negative trajectories arise when we take higher values for the parameter a. We should notice here that the marginal costs are close to each other. Now, we study another parameter set, that is, c 1 = 0.2, c 2 = 1.2, and a = 0.1. At this set, the Nash point is positive and becomes (0.67418, 0.029031). The Jacobian matrix (8) at this set with k = 2.84 is given by which has two complex conjugate eigenvalues λ 1,2 ≈ 0.13245 ± 0.42747i with |λ 1,2 | < 1; that means the Nash point is stable. It gets unstable due to the coexistence of a mixed situation, that is, a Neimark-Sacker followed by a flip bifurcation as k increases. Substituting c 1 = 0.2, c 2 = 1.2, a = 0.1, k = 2.84 in (11), we get g(1) > 0, g(-1) > 0, and < 0, which confirm the appearance of Neimark-Sacker bifurcation, and hence E gets unstable. This is depicted in Fig. 2a. Keeping the parameter set fixed and increasing k to 2.914, a spiral is formed around the Nash point as given in Fig. 2b. The dynamic behavior given in Fig. 2b contains negative values for the second player, which means it may leave the market. At k = 2.925 and the other parameter values are fixed, a period-4 cycle is born with negative values. It is plotted in Fig. 2c with its basins of attraction that contain lobes. The basins of attraction get quite complicated as k increases more. This is clear from the basins given in Fig. 2d for the period-8 cycle. This cycle emerges at k = 3.25. Further increase in the parameter k gives rise to different chaotic attractors of the system that change as k approaches 3.488 to one chaotic attractor as shown in Fig. 2e and Fig. 2f. Any other increasing in the bifurcation parameter k at this set of parameter values gives only one chaotic attractor. In addition, numerical experiments show that negative trajectories may coexist when a takes higher values, and it must be small should we require to get positive trajectories with respect to k. Now, we study the influence of the cost parameter c 1 on the system dynamics when the other parameter values including k are fixed. Assuming the following parameter set, c 2 = 0.3, k = 3, and a = 0.9, it is clear in Fig. 3a that the Nash point loses its stability due to flip bifurcation. Since the Nash point depends on c 1 , c 2 , and a, we assume c 1 = 0.5, and then E = (0.1335402095, 0.8225670156). It is a positive point and the Jacobian becomes J E = 0.78416 -0.099553 -0.10205 0 which has two real eigenvalues λ 1 ≈ 0.79691 and λ 2 ≈ -0.01275; that means the Nash point is stable as |λ 1,2 | < 1. As previously, increasing c 1 and keeping the other set of parameter values fixed give rise to different periodic cycles and chaotic attractors. For instance, at c 1 = 1.31 a period-2 cycle arises as shown in Fig. 3b. It has complicated basins of attraction which contain lobes. It is also clear that this period cycle contains negative values for decision variable of the first firm. Figure 3c and Fig. 3d show other dynamic behaviors of the system at c 1 = 1.38 (the basins of period-4 cycle with gray color refer to divergent

Critical curves and noninvertibility
From the above discussion and analysis, we find that the structure of basins of some periodic cycles of system (7) is quite complicated. For that, we set q 1,t+1 =q 1 and q 2,t+1 =q 2 in system (7), which means the time evolution is obtained by the iteration T : (q 1 , q 2 ) → (q 1 ,q 2 ) given by T : System (13) is a noninvertible system because for a given point (q 1 ,q 2 ) ∈ R 2 it may have up to four real rank-1 preimages; these preimages can be obtained by solving algebraically system (13) with respect to the variables q 1 and q 2 . Furthermore, the phase plane contains these points in different regions. Therefore, the phase plane may be divided into several regions represented by Z i , where i refers to the rank-1 preimages number in each region. This requires to calculate both LC and LC -1 . The later is determined by vanishing the determinant of the Jacobian of system (7) as follows: kq 1 (aq 1 + q 2 )(1 -2 c 2 (a + q 1 )) (a + q 1 + q 2 ) 3 c 2 (a + q 1 ) = 0.
It is clear in Fig. 4 that at this set of parameter values the critical curve LC divides the phase plane of period-2 cycle into two different regions: Z 0 and Z 2 . Therefore, system (7) is noninvertible. The region Z 0 does not have rank-1 preimages, while the region Z 2 possesses two real rank-1 preimages. Because of the complexity, to get an analytical expressions for those preimages, we can obtain real preimages for the points (q 1 , 0) ∈ R 2 and (0,q 2 ) ∈ R 2 as follows.
Case 1: Substituting q 1 = -a in the first equation of (19) with simple calculations gives Solving algebraically (20) with respect to q 2 , we get the first part of Proposition 3 proved. Case 2: Substituting q 1 = 1-ac 2 c 2 in the first equation of (19) and solving it as in case 1 complete the proof.
Another important aspect of system (13) we should mention here is that it is undefined at the line q 2 = -aq 1 . This makes it belonging to the family of maps cited in [23][24][25] and defined in the form (q 1 ,q 2 ) = ( N(q 1 ,q 2 ) D(q 1 ,q 2 ) , G(q 1 , q 2 )). Our system (13) is defined for all (q 1 , q 2 ) ∈ R 2 in the phase plane except at the point (q 1 , -aq 1 ) that makes D(q 1 , q 2 ) = 0 and its preimages too. Furthermore, system (13) is in the form 0/0 because N(q 1 , q 2 ) = 0 and D(q 1 , q 2 ) = 0 at the point (0, -a). The point (0, -a) is called a focal point. But the line q 2 = -aq 1 makes D(q 1 , q 2 ) = 0 only and forms a set of non-definition points denoted by L S = {(q 1 , q 2 ) ∈ R 2 : q 2 = -aa 1 } and its set of prefocal points gets L Q : q 2 = -a. More information about the properties of this point and its lobes is well discussed in [5].

Conclusion
The current paper has introduced and studied a nonlinear duopoly game whose players are heterogeneous. The players have adopted a gradient based mechanism known as bounded rationality and a naive mechanism based on the best reaction. The unique equilibrium point of the game, which is a Nash equilibrium, has been calculated and its stability conditions have been obtained. We have investigated the two routes by which Nash equilibrium point can be destabilized. Those routes have given a destabilization behavior of the game's system via flip and Neimark-Sacker bifurcations. In order to detect these routes, we have used numerical simulation that includes 1D bifurcation diagram, largest Lyapunov exponent, phase diagram, and basins of attraction. All the numerical simulations carried out have shown quite complicated basins of attraction for the periodic cycles constructed by the game's dynamics. Furthermore, we have shown the appearance of lobes regions in the basins of periodic cycles. In addition, we have proved that the game's system is noninvertible and belongs to Z 2 -Z 0 type. Finally, our game's system has generalized some works in literature.