Stability , Bifurcation , and Chaos in N-Firm Nonlinear Cournot Games

An N-firm production game known as oligopoly will be examined with isoelastic price function and linear cost under al Cournot competition. After the best responses of the firms are determined, a dynamic system with adaptive expectations is introduced. It is first shown that the local asymptotic behavior of the system is identical with that of the adaptive adjustment process in which the firms cautiously determine their outputs. Dynamic analysis is confined to two special cases, one in which N is divided into two groups and the other in which N is divided into three groups. Then stability conditions will be derived and the global behavior of the equilibria will be illustrated including chaos control. Lastly the twoand three-group models are compared with two-firm duopoly and three-firm triopoly models to shed light on roles of the number of the firms.


Introduction
The state sequence of discrete dynamic systems will be considered as time series, with a deterministic rule to obtain the consecutive state variables.Among the large variety of dynamic economic systems the oligopoly models have a very special place, since the long-term behavior of the state trajectories has many different possibilities including global asymptotic stability, limit cycles with increasing number of periods and even chaos.Following Cournot 1 many researchers worked on developing more realistic models and on examining their properties.The existence and uniqueness of the equilibrium was the main focus in earlier studies, and later the researchers turned their focus to the dynamic extensions of these models.A comprehensive summary of earlier results is given in Okuguchi 2 , and their multiproduct extensions with different model variants are discussed in Okuguchi and Szidarovszky 3 .Most studies considered concave oligopolies with monotonic best responses.Puu 4 replaces a linear demand function the most popular form with an isoelastic demand function the second popular form , which can be derived by assuming a Cobb-Douglas type utility function of the market.In this case the best responses are unimodal, making equilibrium and stability analyses more complicated and more different from those under monotonic best responses.In particular, chaotic dynamics emerges through a Neimark-Sacker bifurcation when inherited nonlinearities becomes stronger.Agiza  This work continues an earlier paper of Matsumoto 10 , where chaos control for nonlinear duopolies was examined.The number of the firms is generalized to N ≥ 2. After determining best responses of the firms and giving detailed equilibrium analysis, it constructs an adjustment process with adaptive expectations in which the expectations are adaptively updated.To simplify dynamic analysis, attentions are confined to two special cases, one in which the N firms are divided into two groups i.e., two-group model and the other in which the N firms are divided into three groups i.e., three-group model .It is shown that these models generate complex dynamics involving chaos.It is further shown that complex dynamics could be stabilized by proper selection of the speeds of adjustment.
The paper develops as follows.In Section 2 the general model will be introduced, the best responses of the firms and the Cournot equilibrium will be determined.In Section 3 we will show that the local stability properties of the adjustment process with adaptive expectations are identical with those of the adaptive adjustment process in which the firms cautiously adjust their outputs in the next period.In Section 4 two special cases, twoand three-group models, will be analyzed both theoretically and numerically in which the dynamics is two-and three-dimensional, respectively.The stability regions, where chaos is controlled, will be shown and their dependence on the number of firms will be illustrated.Section 5 will conclude the paper.

Nonlinear Oligopoly Models
It is assumed that a homogeneous market is supplied by N firms.For the sake of mathematical simplicity only one product is considered.Let x i denote the production output of firm i, then S i j / i x j is the output of the rest of the industry and S N i 1 x i is the total output of the industry.We assume isoelastic price function p 1/S and linear cost functions C i x i c i x i as in the duopoly model of Puu 8 .Since the firms make decisions about their production levels simultaneously, the firms do not know the outputs of the rivals when their decisions are made.Each firm i can have only an expectation prediction of the output of the rest of the industry, S e i .So the expected profit of firm i can be given as

2.1
Notice that this function is strictly concave in x i .
The strict concavity of Π e i implies that with any given value of S e i the profit maximizing output level of firm i can be computed as

2.2
This function is continuous, piecewise differentiable, and in interval 0, 1/c i it is strictly concave in S e i .Then the best response dynamic process is for i 1, 2, . . ., N.

2.3
Dynamic characteristics are sensitive to the expectation formation.In this study we first consider naive expectation in which the firms assume that the output of the rest of the industry remains the same as in the previous period: and call it a naive system.It is well known that the naive system is a special case of the best reply dynamics with adaptive expectations.In the appendix, local stability conditions for dynamic systems with adaptive expectations are derived.It is, as will be seen, useful to determine the stability of not only naive systems but also that of the controlled systems.
Without losing generality we may assume that at the equilibrium all firms have positive outputs, otherwise we can ignore the firms with zero equilibrium output values and decrease the value of N. Assuming a positive equilibrium, then, from the definitions of naive expectations and the reaction function of firm i, we get

2.5
Since the total industry output is S x i S e i , from the second equation we have that is,

2.7
Adding this equation for all values of i and denoting the sum of the marginal costs by Therefore there is a trivial equilibrium with S x 1 • • • x N 0, and a nontrivial positive equilibrium with where the Cournot output of firm i becomes The superscript "c" is attached to variables to indicate that they are computed at the Cournot equilibrium.Our concern is on the nontrivial point, x c i i 1,2,...,N , and thus no further considerations will be given to the trivial point.For a positive Cournot output, the following inequality has to be satisfied: This always holds for N 2, and necessarily holds for N > 2 if the marginal costs, c i , are sufficiently close to each other.In the rest of this paper, we assume that this condition is satisfied.We also need to guarantee that in 2.2 the first case applies.Substituting c i S 2 for S e i and arranging terms yield C > c i N − 1 .Thus the condition S e i < 1/c i can be reduced to the nonnegativity condition 2.11 .
Adaptive expectations are generalizations of naive expectation where the expected output of the rest of the industry is computed as x j t .

2.12
Notice that in the case of α i 0, the expectation S e i t of firm i remains constant and therefore the same best response is chosen at all time periods t, which is not equilibrium strategy in general.If α i 1, then this formula reduces to naive expectation.The best response dynamic process with adaptive expectations is a 2N-dimensional system with state variables x i and S e i i 1, 2, . . ., N .It is easy to see that the steady state x i , S i of the system and the Cournot equilibrium of the static N-firm oligopoly coincide: x i x c i and S i S c i .In the appendix the local stability conditions for system with adaptive expectations are derived.Combining inequalities 2.11 and A.15 from the appendix with α i 1 gives and the eigenvalue equation A.16 of the appendix reduces to the following: 1, where γ i ≡ ∂f i S e i ∂S e i at the equilibrium.

2.14
It can be written as the quadratic and cubic equations, 2.15

General Stability Conditions
In this section, we first consider the dynamic process in which firms cautiously adjust their outputs, that is, the output in the next period is a weighted average of the current output and the best reply with naive expectations: or, equivalently, This is commonly known as the adaptive adjustment process.Here we assume that firm i moves into the direction toward its profit maximizing output, and reaches it only for α i 1.
Since this adjustment process describes the best reply dynamics with inertia, we call it the inertia control system.Here α i is the inertia or control parameter of firm i and assumed to be positive and not greater than unity.It is easy to see that the fixed point of the inertia control system is the same as that of the naive system.The Jacobian of the system has the form In the appendix, we show that the nonzero eigenvalues of the Jacobian of the system with adaptive expectations are the eigenvalues of matrix It is easy to see that the characteristic polynomials of matrices H and H c are identical since If any or more γ i 0, then the continuity of the characteristic polynomial coefficients in the matrix elements implies the result.Hence the local stability conditions of dynamics with adaptive expectations and with inertia control are identical.So in the rest of this paper we will consider only inertia control which also contains models with best response dynamics as special case with α i 1.
In the case of concave oligopolies see, e.g., Szidarovszky and Chiarella 11 , Bischi et al. 9 it is proved that −1 < γ i ≤ 0, all eigenvalues are real and they are inside the unit circle if and only if for all i,

3.6
These conditions are clearly satisfied if the speeds α i of adjustments are sufficiently small.However in the case of isoelastic inverse demand functions there is the possibility of complex eigenvalues, so no such simple general stability conditions can be derived.In the next section stability in two special cases will be examined both theoretically and numerically.

Stability Conditions in Special Models
In this section we will focus on two cases: one case where the industry consists of two groups of identical firms and the other case with three groups of identical firms.Our aim is to see whether the unstable naive system can be stabilized by the inertia control method.

Two Groups of Firms
Assume there are two groups of firms in a sense that firms of the same group produce the same output with the same marginal cost and have the same speed of adjustment.So the N firms are divided into two groups.Without loss of generality, we can assume that the first N a firms are in the first group and the remaining N b firms in the second group, where N a N b N and N a 1/w N with w > 1.By assuming N/N − 1 < w < N, we will not consider the extreme division in which N is divided into 1 and N − 1.We denote the outputs produced by the firms of the two groups by x and y, so the two kinds of marginal costs are denoted by a and b, so It is also assumed that the firms use adaptive output adjustments 3.2 with naive expectation with the two values of speeds of adjustment We further assume, without any losses of generality, that the firms in the first group are more efficient than the ones in the second in a sense that a < b.

4.4
Accordingly, the marginal cost ratio of b over a, which we denote by h, is greater than unity, The sum of the marginal costs and the derivatives of the reaction functions evaluated at the Cournot equilibrium are, respectively,

4.6
From 2.10 , the equilibrium outputs of the firms are where

4.8
The first inequality is always true because h > 1, so x c is always positive.The second equation indicates that y c is positive if N a > 1 and the marginal cost ratio is bounded from above: Notice that N a > 1 implies that the denominator of the marginal cost ratio, h N , is positive.It can be seen that h N decreases in N and is approaching unity as N converges to infinity with fixed values of w or N a converges to infinity regardless of the value of w.Since the ratio h is assumed to be greater than unity, this implies that it becomes more difficult to have a positive Cournot equilibrium as the number of the firms in the industry increases.Similarly, condition A.15 for each group can be written as for the second group.

4.10
However neither inequality is effective, since for N ≥ 3 and α, β ∈ 0, 1 , which indicates that if h fulfills 4.9 , then it also satisfies A.15 .In this case, A.16 assumes the form since N a and N b firms have identical parameters, in which case This is a quadratic equation in λ, λ 2 pλ q 0 4.14

4.15
It is well-known that the roots of the quadratic equation 4.14 are inside the unit circle if and only if

4.16
The left hand side of the second stability condition can be rewritten as which is always positive.This indicates that the characteristic equation 4.14 does not have a root equal to unity.The other two conditions are hard to be simplified and explained in general.Therefore we make a specializing assumption that the speeds of adjustment are the same for the two groups.
Assumtion 4.1 α β .The feasible region in terms of the adjustment speeds and the marginal cost ratios for which the Cournot equilibrium is positive is where h N is the upper bound of the marginal cost ratio defined in 4.9 .Since it is determined by N and w, the feasible region is presented by a rectangular, P N 0, 1 × 1, h N which decreases in N i.e., P N 1 ⊂ P N .
The locus of q − 1 0 is the Neimark-Sacker boundary on which the Cournot equilibrium changes stability through a pair of complex conjugates and a closed invariant curve is born.We first examine the possibility of the Neimark-Sacker bifurcation.Under α β and h h N , the left hand side of the first stability condition, q − 1, has the form It is natural to confine our analysis to the case where N, N a and N b are integers.In consequence, the relation N ≥ 2w must hold to have that N a N/w is an integer, since N > w has already been assumed.Then substituting N/2 for w yields

4.20
It then follows from the assumption α ≤ 1 that q − 1 < 0. In other words, the q − 1 0 locus does not cross the h h N locus.We may go from this to the conclusion that the first stability condition q − 1 < 0 is always satisfied for any α, h ∈ P N , given w and N. Our first result on Cournot dynamics is summarized as follows.We then examine the possibility of flip bifurcation.The locus of 1 − p q 0 is the flip boundary on which at least one of the eigenvalues is equal to −1.Crossing this boundary the equilibrium point goes through a period-doubling cascade to chaos.Solving 1 − p q 0 for α under α β and h h N yields

4.21
Given w and N, the above equation determines a threshold value of α at which the flip boundary crosses the h h N locus.It, however, still has a complicated expression.So, instead of analytic study, we will consider an important special case and perform numerical simulations.For further simplification, we set w 2 i.e., the industry consists of two groups with equal size and increase the number of the firms by two from N 4 to N 10 to find the threshold value of α when the number of the firms in the industry increases:

4.22
For N 4, the threshold value α 4 is greater than unity.This implies that the flip boundary and the h h N locus do not intersect.Thus we see that the Cournot equilibrium is stable in the feasible region P N when N 4. On the other hand, the threshold value α i for i 6, 8, 10 is less than unity.This implies that the flip boundary with N ≥ 6 intersects the h h N locus and then divides the corresponding feasible region P N into two parts, stability region and instability region.Figure 1 illustrates three intersections denoted by three dots of the flip boundaries and the h h i locus for i 6, 8, 10 and the separation of the feasible regions.When N 6, the corresponding feasible region P 6 is a rectangular having the h h 6 locus as its upper boundary.The locus denoted by F 6 is the flip boundary with N 6 and divides the feasible region P 6 into two regions: the stability region is the largest rectangular with lightgray in the left of the flip boundary and the instability region is the white region in the right of the flip boundary and below the h h 6 locus.In the same way, it is shown in Figure 1 that the flip boundaries with N 8 and N 10 denoted by F 8 and F 10 divide the corresponding feasible regions into two parts.The smallest rectangular with dark-gray is the stability region with N 10 and the remaining middle rectangular is the one with N 8. Since the h h N locus and the flip boundary shift downward and leftward, respectively, the feasible and stable region becomes smaller as the number of firms N increases.
Our extension from the two-firm i.e., duopoly model to the two-group model reveals interesting features of the nonlinear oligopoly.See Puu 8 for a theoretical and numerical analysis of the two-firm duopoly and three-firm triopoly models.To see them, we perform numerical simulations in two models and compare the results.To this end, in the two-group model, we take the set of parameters N 6, w 2 and choose the two bifurcation parameters: one is the common adjustment coefficient α and the other is the cost ratio h b/a between the two groups where a 1 is taken for simplicity.The dynamic system 3.2 with the two groups has now the special form: x

4.23
On the other hand, in the two-firm model, replacing N 6 with N 2 reduces the above dynamic system to

4.24
Changing α from 0.6 to 1.0 and h actually b from 1 to 1.5 h 6 generates the twoparameter bifurcation diagram shown as the left part of Figure 2 while changing α from 0.9 to 1.0 and h from 3 2 √ 2 to 6.5 yields the two-parameter bifurcation diagram, the right part of Figure 2. Different colors in the α, h plane indicate different periods of cycles up to 16. Periodic points with a period larger than 16 and aperiodic points are colored in gray.

Three Groups of Firms
Assume that the industry consists of three groups with N a , N b and N c firms where N a N b N c N and N a 1/w a N and N b 1/w b N with w a > and w b > 1.The common outputs of the three groups of firms are denoted by x, y and z, the marginal costs are a, b and c with and the speeds of adjustment are denoted by α, β and γ , As in the previous section, we assume that the firms use adaptive output adjustments with naive expectations, and the firms in the first group are the most efficient in a sense that their marginal cost is the smallest, a < min{b, c}.

4.28
The marginal cost ratios are denoted by h and k,

4.29
The Cournot outputs are obtained from 2.10 , Then it can be seen that left hand sides of the positivity conditions 2.11 have the forms

4.31
The first inequality is always true.In order to have y c , z c > 0, we have to assume that

4.32
Let the right hand sides of the first and the second inequalities be denoted by f b h and f c h .There is no guarantee that f b h > 1, but f c h > 1 always for h > 1.Then the set of the marginal cost ratios that generate positive Cournot equilibrium can be defined by

4.33
Characteristic equation A.16 now can be written in the form where the coefficients are The Cournot equilibrium is locally asymptotically stable if all eigenvalues are less than unity in absolute value.The most simplified form of the necessary and sufficient conditions for the cubic equation to have roots only inside the unit circle is as it has been proved in Farebrother 12 and Okuguchi and Irie 13 .It is easy to show that which implies that unity is not a root of the cubic equation.It seems tedious to examine the remaining two conditions in general, instead we numerically confirm the stability region in the special case with N 6 and w 1 w 2 3.The qualitatively same results can be obtained for any N > 4 , w 1 and w 2 .In this case we have three groups with two firms in each of them.The region for positive Cournot equilibrium is surrounded by the k f a h and k f b h loci.Substituting N a N b N c 2 into these functions determines the nonnegativity region

4.39
The positivity and stability regions under naive expectation with α 1 and those under adaptive adjustment with α 0.7 are given in Figures 3 a and 3 b .In each figure, the outer curve is the Neimark-Sacker boundary, the inner curve is the flip boundary and the two straight lines are the k f b h and the k f c h loci where the former is steeper than the latter.The light gray area bounded by the two straight lines is the positivity region and the dark gray area illustrates the stability region.Their intersection is the feasible and stability region.It is seen in Figure 3 a that the positivity region is completely outside the stability region when the expectations are naively formed.This implies that the positive Cournot equilibrium is always unstable under naive expectation.This result reminds us the classical result of Theocharis 14 that the stability of a nondifferentiated Cournot equilibrium can be confirmed only in the duopoly framework if the expectations are naively formed and the price and cost functions are linear.It is, in turn, seen in Figure 3 b that about half of the positivity region is included in the stability region for adaptive adjustment, which implies that an unstable Cournot equilibrium under naive expectations could be stabilized by the adaptive adjustment system.
Figure 4 is an enlargement of the south-west corner of Figure 3 b .It shows how the flip boundary shifts if the value of α changes.In the current example, solving 1 p −q r 0 for α yields 0.8, given k h 2. This implies that the flip boundary with α 0.8 passes through the point 2, 2 , the vertex of the triangular part of the stable region.Thus if α ≥ 0.8, then the positivity region is completely outside the stability region for adaptive adjustment, so no positive equilibrium becomes stable.On the other hand, solving 1 p −q r 0 for α yields 2/3, given k h 1.This implies that the flip boundary with α 2/3 passes through the point 1, 1 .If α ≤ 2/3, then the positivity region is entirely inside the stability region, in which case all nonnegative equilibria become stable.If 2/3 < α < 0.8, then only a certain part of the nonnegativity region belongs to the stability region, and this part becomes larger if the value of α decreases.In particular, for α 0.7, the stability region is horizontally-striped and is the triangle with a base of the flip boundary with α 0.7, the most outer circle-wise curve.In the remaining light-gray area of the positivity region, the Cournot output is locally unstable.When the adaptive parameter α decreases, the flip boundary shifts inside accordingly.As a consequence, the stability region enlarges and the unstable region shrinks.Our third result on Cournot dynamics is numerically confirmed and then summarized as follows.
Theorem 4.4.In the special case of the three-group model where α β γ, N 6 and w 1 w 2 2, the Cournot equilibrium is locally unstable if α ≥ 4/5, locally asymptotically stable if α ≤ 2/3 while the feasible region is divided into the stability and instability regions by the flip boundary if 4/5 > α > 2/3.
As the comparison of the three-group and the three-firm triopoly models, we, again, perform numerical simulations for these two models.The dynamic system in the three-group model becomes

4.40
where for simplicity, we assume that α β γ , h k, and a 1 implying b c.The dynamic system in the three-firm model i.e., triopoly can be constructed similarly.Selecting α and h in particular b as the bifurcation parameters, Figure 5 illustrates the bifurcation diagrams in the α, h plane.The numerical investigations clearly reveal qualitatively the same issues as the ones we saw in Figure 2. Namely, first of all, the destabilizing process goes to chaos through a flip bifurcation with a lower production ratio in the three-group model and a Neimark-Sacker bifurcation with a higher ratio in the three-firm model.Second, the adjustment speed can be used to control unstable trajectories in the two-and three-group models.Comparing the bifurcation diagram of the two-group model and that of the threegroup model shows the similarity of the destabilizing process in which period-doubling bifurcation takes place.
In Figure 6, we present the bifurcation diagram in the h, k plane to draw attention to the feasibility of the solutions of the three-group model.The value α 0.8 is fixed.As already explored in Figure 4, h and k range from 1 to 2 and the stationary state is unstable for any combination of h and k from this region.Color has the same meanings as before.The feasible region of the three-group model is defined by two upward sloping curves, k f b h and k f c h .Notice that although the area outside the feasible region is colored in the same way as the feasible region, the stationary point defined in that area becomes negative and thus economically meaningless.Bifurcation makes sense only in the feasible region from an economic point of view.Notice further that a trajectory that oscillates around the stationary point periodically or aperiodically may take negative values and thus become economically meaningless.One way is to confine the parameter choice in such a way that the resultant dynamics does not become infeasible.The other way is to reconstruct the dynamic system by taking into account the nonnegativity constraint explicitly.However, the former has the difficulty of deriving the confinement conditions in the currently considered model as many   parameters are involved, and the latter makes the asymptotic behavior of the dynamic system significantly different and more difficult to analyze.Since our main concern is to control the unstable trajectories and our main conclusion is that the adjustment speed is an effective control parameter, which is supposed to hold in those models, we used the model without such modifications at the expense of some economic precision.

Conclusion
Discrete dynamic systems always generate time series.The asymptotic behavior of them is a fundamental research issue.In this paper a special class of economic models was examined.
For the sake of mathematical simplicity we selected N-firm Cournot oligopolies without product differentiation, and with isoelastic price function.The reaction functions and the equilibrium were determined first, and then the asymptotic behavior of the equilibrium was illustrated in two special cases with two or three groups of identical firms.Stability conditions could be derived analytically in the first case, and the dependence of the asymptotic properties of the equilibrium on the number of firms was illustrated by computer simulation in the second case.The results of the nonlinear duopoly and triopoly models show that the Cournot equilibrium can be destabilized through a Neimark-Sacker bifurcation.We also found the following new dynamic phenomenon.In the multigroup models, the stationary state is destabilized through the Feigenbaum period doubling sequence, furthermore a Neimark-Sacker bifurcation can occur only in the infeasible regions in which the stationary state is negative.For N > 4, the multigroup models are unstable if α is close to one and become stable if α is below a certain threshold, regardless of the production cost ratios.This implies that the main source of instability is the speeds of adjustment and thus the stationary state could be stabilized by selecting sufficiently small speeds of adjustment.That is, the multigroup models with N firms are unstable under naive expectations but are controllable under adaptive adjustment process in which the speed of adjustment is the control parameter.

Appendix
In this appendix, we derive the stability conditions of the adaptive system in which the expectation is adaptively formulated, y e i t 1 1 − α i y e i t α i j / i x j t .A.1 Here α i ∈ 0, 1 is the speed of adjustment of firm i in assuming the expected aggregate output of the other firms.The stability conditions are also useful to determine the local dynamic behavior of the naive system as well as that of the inertia system, that is, the system where the firms with naive expectations adaptively adjust toward the best reply solution.
We consider the adjustment process with adaptive expectations first, since the one with naive expectations can be obtained by selecting the speeds of adjustment equal to unity i.e., α i 1, i 1, 2, . . ., N .For i 1, 2, . . ., N, The Jacobian at the equilibrium has the form A.6 Subtracting the γ i -multiple of the second equation from the first one gives The value λ 0 cannot destroy stability, so we may assume λ / 0. Then u i γ i v i , and by substituting it into the second equation, we have This is the usual eigenvalue problem of the N × N matrix Notice that

A.11
The characteristic polynomial of H can be decomposed by using the simple fact that if x, y ∈ R N , then

A.13
The roots of the first factor are 1 − α i 1 γ i which are inside the unit circle if and only if which occurs if and only if The other eigenvalues are the roots of and Elsadany 5 have investigated discrete-time Cournot duopolies with heterogeneous players.Richter and Stolk 6 have introduced a new method of controlling coexisting chaotic attractors in Cournot triopolies by means of steering the systems dynamics from one attractor to another.See Puu and Sushko 7 , Puu 8 and Bischi et al. 9 for comprehensive summary of recent developments in the theory of nonlinear oligopolies.

Theorem 4 . 2 .
If α β, N/ N − 1 < w < N and all of N, N a and N b are integers, then the Neimark-Sacker bifurcation does not occur in the feasible region P N .

Figure 2 :Theorem 4 . 3 .α β, w 2 and N 6 ,
Figure 2: Bifurcation diagrams in the α, h plane for the two-group model a and for the two-firm model b .

7 Figure 3 :Figure 4 :
Figure 3: Stability and feasible regions of the naive and adaptive systems.

Figure 5 :
Figure 5: Bifurcation diagrams in the α, h plane for the three-group model a and for the three-firm model b .

det I xy T 1 y
T x, A.12 where I is the N × N identity matrix.So we have det H − λI det D ab T − λI det D − λI det I D − λI −1 ab T .