Controlling Chaotic Dynamics in N-firm
Nonlinear Cournot Games∗
Akio Matsumoto†Chuo University

Ferenc Szidarovszky‡
University of Arizona

Abstract
The state trajectories of a special class of discrete dynamic economic
systems will be considered. The control of their long—term behavior is amajor research issue of dynamic markets. An N-firm production game
known as oligopoly will be examined with isoelastic price function andlinear costs. After the reaction functions of the firms are determined, dynamic systems with adaptive expectations and with adaptive adjustments
will be introduced and the equivalence of their local asymptotic behavior
will be verified. Then two special cases will be investigated in details,
with two and three groups of identical firms, in which the dynamics is twoand three dimensional, respectively. Stability conditions will be derived
and the global behavior of the equilibria will be illustrated including chaos
control.

∗ The authors highly appreciate financial supports from the Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research (C), 21530172), Chuo University (Joint
Research Grant 0981) and the Japan Economic Research Foundation. They also want to acknowledge the encouragement and support by Kei Matsumoto for the research leading to thispaper. Part of this paper was done when the first author visited the Department of Systems
and Industrial Engineering of the University of Arizona and the second author visited theDepartment of Economics of Chuo University. The usual disclaimer applies.
† Department of Economics, Chuo University, 742-1, Higashi-Nakano, Hachioji, Tokyo, 1920393, Japan. akiom@tamacc.chuo-u.ac.jp; +81-426-74-3351(tel); +81-426-74-3425 (fax).
‡ Department of Systems and Industrial Engineering, University of Arizona, Tucson, AZ,
85721-0020, USA. szidar@sie.arizona.edu; +01-520-621-6557 (tel); +01-520-621-6555 (fax).

1

1

Introduction

The state sequence of discrete dynamic systems will be considered as time series,
with a deterministic rule to obtain the consecutive state variables. Many important classes of dynamic systems are examined frequently in the literature ofmathematical economics. 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 diﬀerent possibilities including global asymptotic stability, limit cycles with increasing number of nodes and even chaos. In
this paper a special oligopoly model will be investigated and conditions will bederived to avoid chaotic behavior. Following Cournot (1838) 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 ofthese models. A comprehensive summary of earlier results is given in Okuguchi(1976), and their multi-product extensions with diﬀerent model variants are
discussed in Okuguchi and Szidarovszky (1999). Most studies considered concave oligopolies with monotonic reaction functions, and only a little attention
was given to the isoelastic case, which can be derived by assuming a CobbDouglas type utility function of the market. In this case the best response is
not monotonic anymore, making equilibrium and stability analysis more complicated. Agiza (1999) has examined the two-dimensional Kopel map in which
the reaction functions are unimodal, derived stability conditions and studied
bifurcation and chaos by computing the maximum Lyapunov exponents. Agiza and Elsadany (2003) have investigated discrete-time Cournot duopolies withheterogeneous players. Richter and Stolk (2004) have introduced a new methodof controlling coexisting chaotic attractors in Cournot triopolies by means ofsteering the systems dynamics from one attractor to another. See Puu andSushko (2002), Puu (2003) and Bischi et al. (2010) for comprehensive summaryof recent developments in the theory of nonlinear oligopolies.
This work is an extension and generalization of an earlier paper of Matsumoto (2006), where chaos control for nonlinear duopolies was examined. Inline with the economic literature, controlled systems are dynamics with adaptive expectations and uncontrolled systems are those with naive expectations.
We will give detailed equilibrium analysis, show the complexity of the dynamics
of such models and further show that complex dynamics involving chaos could
be stabilized.
The paper develops as follows. In Section 2 the general model will be introduced, the reaction functions and the Cournot equilibrium will be determined.
In Section 3 we will show that the local stability properties of systems withadaptive expectations and with inertia control (i.e., adaptive adjustments) are
identical. In Section 4 two special cases, two and three groups of firms, will beanalyzed both theoretically and numerically in which the dynamics are two andthree 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.

22Nonlinear Oligopoly Models

It is assumed that a homogeneous market is supplied by N firms. For the sakeof mathematical simplicity only one productP is considered. Let xi denote theproduction output of firm i, then yi =j6=i xj is the output of the rest ofPN
the industry and Q = i=1 xi is the total output of the industry. We assume
isoelastic price function p = 1/Q and linear cost functions Ci (xi ) = ci xi as in
the duopoly model of Puu (2003). Since the firms make decisions about theirproduction levels simultaneously, the firms do not know the output of the rivals
when their decisions are made. They can have only an expectation (prediction)of the output of the rest of the industry, yie . So the expected profit of firm i canbe given asxi
Πei =− ci xi , (i = 1, 2, ..., N ).
xi + yie
Notice that this function is strictly concave in xi .

2.1

Reaction Functions

The strict concavity of Πei implies that with any given value of yie the profit
maximizing output level of firm i can be computed as
⎧ r e
1
yi
⎪
⎨
− yie if 0 < yie ≤ ,
e
c
c
i
i
xi = fi (yi ) =
(1)
1
⎪
e
⎩ 0
< yi .
ifciThis function is continuous, and strictly concave in the first segment.
Then the best response dynamic process is
xi (t + 1) = fi (yie (t + 1)) for i = 1, 2, ...N.
Dynamic characteristics are sensitive to the expectation formation. In this study
we first consider naive expectation in which the firms assume that the outputof the rest of the industry remains the same:Xyie (t + 1) =xj (t),
j6=i

and call it an 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 arederived. It is, as will be seen, useful to determine the stability of not only naive
systems but also that of the controlled systems.

2.2

Cournot Equilibrium

Without losing generality we may assume that at the equilibrium all firms havepositive outputs, otherwise we can ignore the firms with zero equilibrium output
values and decrease the value of N . Assuming a positive equilibrium, then, fromthe definition of the naive expectation and the reaction function of firm i,
r e
X
yi
e
yi =
xj and xi =− yie .
cij6=i

3

By introducing the notation y = xi + yie , from the second equation we haveyie = ci y 2 ,that is,xi = y − ci y 2 .
Adding this equation
for all values of i and denoting the sum of the marginalPN
costs by C = i=1 ci gives
y = N y − Cy 2 .

Therefore there is a trivial equilibrium with y = x1 = ... = xN = 0, and a
nontrivial positive equilibrium with
yc =

N −1
,
C

where the Cournot output of firm i becomesxci =(N − 1)(C − ci (N − 1))
.
C2

(2)

The superscript ”c” is attached to variables to indicate that they are computedat the Cournot equilibrium. Our concern is on the nontrivial point, (xci )i=1.2,...N ,and thus no further considerations will be given to the trivial point. For nonnegative Cournot outputs, the following inequality has to be satisfied:C − ci (N − 1) ≥ 0 or ci ≤

C
.
N −1

(3)

This always holds for N = 2, and necessarily holds for N > 2 if the marginal
costs, ci , are suﬃciently close to each other. In the rest of this paper, we assume
that this condition is satisfied. Notice that under this condition
µ
¶2
1
N −1
1
e
c 2
≤ ci 2 ≤ ,yi = ci (y ) = ciCcic1
So at the equilibrium the first case of (1) applies.
Adaptive expectations are generalizations of naive expectation where theexpected output of the rest of the industry is computed asXyie (t + 1) = (1 − αi )yie (t) + αixj (t).
j6=i

Notice that by selecting αi = 1 this formula reduces to naive expectation. In the
Appendix the local stability conditions for system with adaptive expectations
are derived. Combining inequalities (3) and (A-1) from the Appendix with
αi = 1 givesCC< ci ≤
for i = 1.2, ...N,4(N − 1)N −1and the eigenvalue equation (A-2) of the Appendix reduces to the following
equation:NX

γi∂fi (yie )at the equilibrium.
= 1, where γ i ≡
γ +λ
∂yie
i=1 i
4

It can be written as the quadratic and cubic equations,
λ2 − γ 1 γ 2 = 0 for N = 2,and−λ3 + (γ 1 γ 2 + γ 1 γ 3 + γ 2 γ 3 )λ + 2γ 1 γ 2 γ 3 = 0 for N = 3.

3

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 averageof the current output and the naively-determined optimal output:
xi (t + 1) = αi fi (yie (t + 1)) + (1 − αi )xi (t) for i = 1, 2, ..., N,or, equivalently,⎛s P

xi (t + 1) = αi ⎝j6=i xj (t)

ci

−

X
j6=i

⎞

xj (t)⎠ + (1 − αi )xi (t).

(4)

This is commonly known as the adaptive adjustment process. Here we assume
that firm i moves into the direction toward its profit maximizing output, andreaches 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⎞⎛
1 − α1 γ 1 α1 · γ 1 α1
⎜ γ 2 α2 1 − α2 · γ 2 α2 ⎟⎟.
Hc = ⎜⎠⎝····
γ N αN γ N αN · 1 − αN

In the Appendix, we show that the nonzero eigenvalues of the Jacobian of
the system with adaptive expectations are the eigenvalues of matrix⎛⎞
1 − α1 γ 2 α1 · γ N α1
⎜ γ 1 α2 1 − α2 · γ N α2 ⎟
⎟.
H=⎜
⎝
⎠
·
·
·
·
γ 1 αNγ 2 αN · 1 − αN

Next we will prove that the characteristic polynomials of matrixes H and Hc
are identical, so as far as local stability is concerned, the asymptotic behaviorof the two systems is identical. Assume first that all γ i 6= 0. Define
R = diag(γ 1 , γ 2 , ..., γ N ).

5

Then clearly
⎛

1
γ1

0

⎜ 0
1
⎜
γ2
R−1 Hc R = ⎜
·
⎝ ·
0
·
⎛ 1−α
1

⎜
⎜
= ⎜
⎝

γ1

α2
·
αN

·
·
·
0

0
0
·

1
γN

α1
1−α2
γ2
·
αN

·
·
·
·

= H.

⎞

⎛
⎞
γ1 0 ·
0⎟0 γ2 ·
0 ⎟
⎟ c⎜
⎟
⎟H ⎜
⎝
·
·
·
· ⎠
⎠
0
· 0 γN⎞⎛⎞α1
γ1 0 ·0α2 ⎟0 γ2 ·
0 ⎟
⎟⎜
⎟
⎟⎜
⎝
·
· ·
· ⎠
⎠
·
1−αN
0
· 0 γNγ
N

Since Hc and H are similar, their characteristic polynomials are identical. If any
or more γ i = 0, then the continuity of the characteristic polynomial coeﬃcientsin the matrix elements implies the result. Hence the local stability conditionsof dynamics with adaptive expectations and with inertia control are identical.
So in the rest of this paper we will consider only inertia control.
In the case of concave oligopolies (see for example, Szidarovszky and Chiarella(2001), Bischi et al. (2010)) it is proved that −1 < γ i ≤ 0, all eigenvalues arereal and they are inside the unit circle if and only if for all i,
αi <
and

2
1 + γiNP

αi γ i> −1.
i=1 2 − αi (1 + γ i )These conditions are clearly satisfied if the speeds αi of adjustments are suﬃciently small. However in the case of isoelastic inverse demand functions thereis the possibility of complex eigenvalues, so no such simple general stability conditions can be derived. In the next section two special cases will be examinedboth theoretically and numerically.

4

Stability Conditions in Special Models

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

4.1

Two Groups of Firms

Assume there are two groups of firms in a sense that firms of the same groupproduce the same output with the same marginal cost and have the same speedof adjustment. So the N firms are divided into two groups. Without loss ofgenerality, we can assume that the first Na firms are in the first group and theremaining Nb firms in the second group, where Na + Nb = N and Na = w1 N

6with w > 1. We denote the outputs produced by the firms of the two groups by
x and y, so
x1 = · · · = xNa = x and xNa +1 = · · · = xN = y,
the two kinds of marginal costs are denoted by a and b, so
c1 = · · · = cNa = a and cNa +1 = · · · = cN = b,
and the two types of speeds of adjustment by α and β, that is,α1 = · · · = αNa = α and αNa +1 = · · · = αN = β.
We further assume, without any losses of generality, that the firms in the firstgroup are more eﬃcient than the ones in the second in a sense thata < b.
Accordingly, the marginal cost ratio of b over a, which we denote by h, is greaterthan unity,bh = > 1.
a
The sum of the marginal costs and the derivatives of the reaction functions
evaluated at the Cournot equilibrium are
C = Na a + Nb b
and
γa =

C
C
− 1 and γ b =− 1.
2a(N − 1)2b(N − 1)

From (2), the equilibrium outputs of the firms arexc =(N − 1)(C − a(N − 1))
,
C2

yc =

(N − 1)(C − b(N − 1))
,
C2

where
C − a(N − 1) = a(Nb (h − 1) + 1) > 0,
C − b(N − 1) = a(Na − 1)

µ

¶
N
−h .
N −w

The first inequality is always true because h > 1, so xc is always positive. The
second equation indicates that y c is nonnegative if Na > 1 and the marginal
cost ratio is bounded from above:h ≤ hN ≡

N
.
N −w

(5)

Notice that Na > 1 is equivalent to N > w implying the positive denominator in
the upper bound hN of the marginal cost ratio. It can be seen that hN decreases7

in N and is approaching unity as N converges to infinity with fixed values of
w. Since the ratio h is assumed to be greater than unity, this implies that it
becomes more diﬃcult to have a positive Cournot equilibrium as the numberof the firms in the industry increases. Similarly, condition (A-1) for each groupcan be written as⎧4(N − 1) − αNa⎪⎪h<for the first group,
⎪
⎨
αN
b

⎪
⎪
⎪
⎩ h>

βNa
for the second group.
4(N − 1) − βNb

However neither inequality is eﬀective, since for N ≥ 3 and 0 < (α, β) ≤ 1,
N
4(N − 1) − αNa
βNa
<1≤
,
<
4(N − 1) − βNb
N −w
αNbwhich indicates that if h fulfills (5), then it also satisfies (A-1).
Notice that in this case, equation (A-2) assumes the form
1 + Na

αγ a
βγ b
+ Nb
= 0.
1 − α(1 + γ a ) − λ
1 − β(1 + γ b ) − λSince Na and Nb firms have identical parameters, in which case αi = α, γ i =γ a (1 ≤ i ≤ Na ) and αi = β, γ i = γ b (Na + 1 ≤ i ≤ Na + Nb ). Notice that withnotation Ra = 1 + γ a and Rb = 1 + γ b , the above equation can be written as(1 − αRa − λ)(1 − βRb − λ) + Na αγ a (1 − βRb − λ) + Nb βγ b (1 − αRa − λ) = 0.
This is a quadratic equation in λ,λ2 + pλ + q = 0(6)

with coeﬃcientsp = α(Ra (1 − Na ) + Na ) + β(Rb (1 − Nb ) + Nb ) − 2,andq

= −α(Ra (1 − Na ) + Na ) − β(Rb (1 − Nb ) + Nb ) +αβ(Ra Rb (1 − N ) + Ra Nb + Rb Na ) + 1.

It is well-known (see for example, Bischi et al. (2010)) that the roots of thequadratic equation are inside the unit circle if and only if
q < 1,
q + p + 1 > 0,q − p + 1 > 0.
The left hand side of the second condition can be rewritten as
q+p+1=(Na + hNb )2 αβ,4h(N − 1)8which is always positive. This indicates that the characteristic equation (6) does
not have a root equal to unity. To simplify the other two stability conditions,
we make a specializing assumption that the speeds of adjustment are the samefor the two groups.
Assumption α = β.
Solving q − 1 = 0 for h yields the Hopf boundaryhH =√w2 (N − 1) + (w − 1)(2 − α)N + DH
(w − 1)(2(N − w) + (w − 1)αN )(7)

with the discriminant©
ª
DH = w2 (N − 1) (N − 1)w2 + 4(w − 1)((1 − α)N − 1) .

On the Hopf boundary the eigenvalues are complex and their absolute values areequal to unity. Thus this boundary is valid for p2 < 4q. It is rather diﬃcult tosimplify this condition, however it can be checked numerically in all particularsituations. Notice that DH is linear in α. At α = 0,ª©
DH = w2 (N − 1)2 w2 + 4(w − 1) > 0,
and at α = 1,DH©
ª
= w2 (N − 1) (N − 1)w2 − 4(w − 1)©
ª= w2 (N − 1) (N − 2)w2 + (w − 2)2 ≥ 0 for N ≥ 2.

Therefore DH ≥ 0 for all α ∈ (0, 1].
By gradually varying the parameter values α and h, our simulation study indicates that if pair (α, h) crosses the Hopf boundary, then the Cournot point bifurcates to periodic cycles, quasi-periodic cycles and then to chaotic fluctuations.
The same phenomenon has been observed earlier by Puu (2003) in his nonlinear
duopoly model. As shown below, the Hopf bifurcation of the n-firm oligopoly,
however, occurs only in the infeasible region in which the non-negativity condition (4) is violated. In consequence we do not pay much attention to the Hopf
bifurcation.
Solving 1 − p + q = 0 for h yields the flip boundary√(4 − α)(w − 1)αN 2 + 2w2 (N − 1)(N α − 4) + 2 DFhF =(8)(w − 1)(4(N − w) + (w − 1)αN )αN
with©
ª
DF = w2 (N − 1)(αN − 4) 2(2 − α)(w − 1)αN 2 + w2 (N − 1)(αN − 4) .

On the flip boundary, at least one of the eigenvalues is equal to −1. Crossing
this boundary the equilibrium point goes through a period doubling cascade tochaos. Note that the sign of this discriminant depends on the particular value

9of α. It can be confirmed that DF = 0 has three real, positive roots:√−w2 + N (w2 + 4(w − 1)) − ∆αS =,4N (w − 1)√
−w2 + N (w2 + 4(w − 1)) + ∆αL =,4N (w − 1)
αM

=

4
,
N

£
¡
¢¤2
with ∆ = w2 − N w2 + 4(w − 1) −32(N −1)(w−1)w2 , since it can be shownthat ∆ ≥ 0. Here S, L and M stand for "smaller", "larger" and "middle," asαS < αM < αL holds. Therefore, DF > 0 for α < αS and αM < α < αL , andDF ≤ 0 otherwise. The adaptive coeﬃcient α is assumed to be positive and not
greater than unity, and the number of firms N must be greater than 2. Since
αL ≥ 2 for N ≥ 2, αL is outside the unit domain of α. When N ≥ 4,
DF ≥ 0 for either αM ≤ α ≤ 1 or 0 < α ≤ αS ,andDF < 0 for αS < α < αM .
When 2 ≤ N < 4, αM becomes greater than unity and thus it is outside thedomain. The sign of DF is positive or negative according to α is less or greaterthan αS .
Judging from the above considerations, the stability domain for non-negative
Cournot equilibrium is surrounded by the non-negativity boundary hN , the Hopf
boundary hH , and the flip boundary hF . The last two boundaries, however, still
have complicated expressions. So, instead of analytic study we will consider an
important special case and perform numerical simulation. For further simplification, we set w = 2 (that is, the industry consists of two groups with equalsize) and increase the number of the firms by two from N = 4 to N = 10 tosee how the shape of the stability domain changes as the number of firms in
the industry increases. We note that the present model with N = 2 reduces to
the Puu√ (duopoly) model in which the Hopf bifurcation is shown to occurs for3 + 2 2 < h ≤ 25/4.
The feasible region in terms of the adjustment speeds and the marginal cost
ratios for which the Cournot equilibrium is non-negative is
PN = {(α, h) | 0 < α ≤ 1 and 1 < h ≤ hN } for N = 4, 6, 8, 10,
where hN is the upper bound of the marginal cost ratio defined in (5). Since
the upper bound is determined by N and w, the feasible region is presented bya rectangular, PN = (0, 1] × (1, hN ] which decreases in N (i.e., PN +1 ⊂ PN ).
Then the stability region is given as
SN = {(α, h) | 0 < α ≤ 1 and M ax{1, hFN } < h ≤ hHN },where hFN and hHN denote the flip boundary (7) and the Hopf boundary (6).
The feasible and stable region is the intersection of these two regions, PN ∩ SN .

10

Figure 1(A) illustrates these boundaries and the upper bound of the marginal
cost ratio, hN , for each value of N (= 4, 6, 8, 10). The Hopf boundary is negativesloping and shifts downward as N increases. In addition it is bounded frombelow and its lower limit is defined by√
6−α+4 2−αlim hHN =
.
N →∞
2+α
It can be checked that the condition p2 < 4q holds if h > 2.19 for N = 4, h >2.62 for N = 6, h > 2.94 for N = 8 and h > 3.17 for N = 10. It is also observedin Figure 1(A) that the Hopf boundaries hHN are located far above the h = hN
locus for each N. Thus the Hopf bifurcation does not occur in the feasible region
PN . The Flip boundary is bounded from above with its upper limitp12 − α + 4 2(4 − α)lim hFN =.
N→∞
4+α
It is positive-sloping for the appropriate values of α and shifts leftward as N
increases. Although the speed of adjustment α is positive and not greater thanunity, it is measured to the right up to the point 1.2 on the horizontal line in
Figure 1(A) to confirm that the flip boundary with N = 4, the most inner flip
boundary, takes unity for α = 1. Thus P4 ∩ S4 = P4 , so all positive Cournot
points are stable for N = 4. However, the flip boundary with N ≥ 6 divides
the corresponding feasible region PN into two sub regions, stability region withhFN < h and instability region with hFN > h. Figure 1(B) enlarges the lowerleft part of Figure 1(A) and depicts three stability regions: the largest rectangular with light-gray is the stability region with N = 6, the smallest rectangularwith dark-gray is the one with N = 10 and the remaining middle rectangular isthe one with N = 8. Since the h = hN locus and the flip boundary shift downward and leftward, respectively, the feasible and stable region becomes smaller
as the number of N increases.

(A) Various boundaries(B) Enlargment of Figure 1(A)Figure 1. Stability and feasible regions for N = 4, 6, 8, 10
Our extension from the two-firm (i.e., duopoly) model to the two-groupmodel reveals interesting features of the nonlinear oligopoly.1 To see them, we
1 See Puu (2003) for the theoretical and numerical analyses of the two-firm (duopoly) andthe three-firm (triopoly) models.

11

perform numerical simulations in two models and compare the results. To thisend, 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
coeﬃcient α = β and the other is the cost ratio h = b/a between the two groups
where a = 1 is assumed for simplicity. The dynamic system has now the from(4)Ãr!
2x(t) + 3y(t)x(t + 1) = (1 − α)x(t) + α− (2x(t) + 3y(t)) ,ay(t + 1) = (1 − α)y(t) + α

Ãr

!
3x(t) + 2y(t)− (3x(t) + 2y(t)) .
b

On the other hand, in the two-firm model, replacing N = 6 with N = 2 reduces
the this dynamic system to
Ãr
!
y(t)
x(t + 1) = (1 − α)x(t) + α− y(t) ,ay(t + 1) = (1 − α)y(t) + α

Ãr

!
x(t)
− x(t) .
b

Changing α from 0.6 to 1.0 and h (actually b) from 1 to 1.5 generates thetwo-parameter bifurcation diagram√ shown in Figure 2 (A) while changing α
from 0.9 to 1.0 and h from 3 + 2 2 to 6.5 yields the two-parameter bifurcation
diagram presented in Figure 2 (B). Diﬀerent colors in the (α, h) plane indicate
diﬀerent periods of the cycles up to 16. Periodic points with a period larger than
16 and aperiodic points are colored in gray. The solution becomes imaginary ifparameters are selected from the white region. The upward sloping black curvein Figure 2(A) and the downward sloping black curve in Figure 2(B) are the flip
boundary of the two-group model and the Hopf boundary of the two-firm (i.e.,
duopoly) model, respectively. As the destabilizing scenario is concerned, thecomparison shows two issues. First the stationary state is destabilized througha flip bifurcation in the two-group model and with Hopf bifurcation in the twofirm model. Even though both models can generate complex dynamics involving
chaos, the way to chaos is diﬀerent. Second, the loss of stability occurs for a
relatively lower production cost ratio in the two-group model and for a higher
ratio in the two-firm model. Furthermore, Figure 2(A) indicates that the valueof the adjustment speed seems to be a main source of the flip bifurcation scenario
since similar bifurcation takes place as α increases regardless of the value of h.
This implies that the adjustment speed can be an eﬀective control variable in
the two-group model.

12(A) The two-group model(B) The two-firm model
Figure 2. Two-parameter bifurcation diagrams in the (α, h) plane

4.2

Three Groups of Firms

Assume that the industry is now divided into three groups with Na , Nb and Nc
firms where Na + Nb + Nc = N and Na = w1a N and Nb = w1b N with wa > 1
and wb > 1. The outputs of the three groups of firms are denoted by x, y andz,x1 = · · · = xNa = x, xNa +1 = · · · = xNa +Nb = y and xNa +Nb +1 = · · · = xN = z,
the marginal costs are by a, b and c,
c1 = · · · = cNa = a, cNa +1 = · · · = cNa +Nb = b and cNa +Nb +1 = · · · = cN = c,and the speeds of adjustment are by α, β and γ,
α1 = ··· = αNa = α, αNa +1 = ··· = αNa +Nb = β and αNa +Nb +1 = ··· = αN = γ.
As in the previous section, we assume that the firms in the first group are themost eﬃcient in a sense that their marginal cost is the smallest,a < min {b, c} .
The marginal cost ratios are denoted by h and k,h=cb> 1 and k = > 1.
aa

The Cournot outputs are obtained from equation (2) as
xc =

N −1
(C − a(N − 1)),
C2

yc =

N −1
(C − b(N − 1)),
C2

zc =

N −1
(C − c(N − 1))C2
13

with C = Na a + Nb b + Nc c. It can also be seen that in the nonnegativity
conditions (3),C − a(N − 1) = a {Nb (h − 1) + Nc (k − 1) + 1} > 0,
C − b(N − 1) = a {Nc k − (Na + Nc − 1)h + Na }andC − c(N − 1) = a {−(Na + Nb − 1)k + Nb h + Na } .

The first inequality is always true. In order to have positive values of y c and z c ,we have to assume that
k>
and
k≤

Nc + Na − 1
Na
h−
Nc
Nc

Nb
Na
h+
.
Na + Nb − 1
Na + Nb − 1

Let the right hand sides of the first and the second inequalities be denoted byfb (h) and fc (h). There is no guarantee that fb (h) > 1, but fc (h) > 1 always forh > 1. Therefore the set of the marginal cost ratios that generate non-negative
Cournot equilibrium can be defined by
NN = {(h, k) | max{1, fb (h)} ≤ k ≤ fc (h) and h > 1} .
The characteristic equation (A-2) can be written in the form
αγ a
βγ b
γγ c
1 + Na 1−α(1+γ)−λ + Nb 1−β(1+γ )−λ + Nc 1−γ(1+γ )−λ = 0.
abcIntroducing the notation Ra = 1 + γ a , Rb = 1 + γ b and Rc = 1 + γ c , it can be
reduced to a cubic equation,−λ3 + pλ2 + qλ + r = 0,
where the coeﬃcients are
p = (−Na + (Na − 1)Ra )α + (−Nb + (Nb − 1)Rb )β + (−Nc + (Nc − 1)Rc )γ + 3,
q

= −2 {(−Na + (Na − 1)Ra )α + (−Nb + (Nb − 1)Rb )β + (−Nc + (Nc − 1)Rc )γ}+(−Nb Ra − Na Rb + (Na + Nb − 1)Ra Rb )αβ+(−Nc Ra − Na Rc + (Na + Nc − 1)Ra Rc )αγ+(−Nc Rb − Nb Rc + (Nb + Nc − 1)Rb Rc )βγ − 3andr= (−Na + (Na − 1)Ra )α + (−Nb + (Nb − 1)Rb )β + (−Nc + (Nc − 1)Rc )γ−(−Nb Ra − Na Rb + (Na + Nb − 1)Ra Rb )αβ−(−Nc Ra − Na Rc + (Na + Nc − 1)Ra Rc )αγ−(−Nc Rb − Nb Rc + (Nb + Nc − 1)Rb Rc )βγ+(−Nc Ra Rb − Nb Ra Rc − Na Rb Rc − (Na + Nb + Nc − 1))αβγ.
14The Cournot equilibrium is locally asymptotically stable if all eigenvalues areless than unity in absolute value. Okuguchi and Irie (1990) have shown thatthe most simplified form of the necessary and suﬃcient conditions for the cubic
equation to have roots only inside the unit circle is1 − p + (−q) − r > 0,
1 + p + (−q) + r > 0,
1 + q + pr − r2 > 0.
It is easy to show that1 − p + (−q) − r =(n − 1)(Na + hNb + kNc )3 αβγ> 0,8hk(n − 1)3

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 w1 = w2 = 3. The same qualitative conclusions can be obtained for any
N (> 4), w1 and w2 . In this case we have three groups with two firms in eachof them. The region for the non-negative Cournot equilibrium is surrounded by
the k = fa (h) and k = fb (h) loci. Substituting Na = Nb = Nc = 2 into thesefunctions determines the non-negativity region
½
¾
2
2
3
P6 = (h, k) | h + ≥ k ≥ h − 1, h > 1 and k > 1 .
3
3
2
The nonnegativity and stability regions under naive expectation (with α = 1)
and those under adaptive adjustment (with α = 0.7) are shown in Figures 3A
and 3B, in each of which the outer curve is the Hopf boundary, the inner curve isthe flip boundary and the two straight lines are the k = fa (h) and the k = fb (h)loci. The light gray area is the nonnegativity region and the dark gray area
illustrates the stability region. Their intersection is the feasible and stability
region. It is seen in Figure 3A that the non-negativity region is completely
outside the stability region for naive expectation. This implies that the nonnegative Cournot equilibrium is always unstable under naive expectation.2 It
is, in turn, seen in Figure 3B that about half of the non-negativity region is
included in the stability region for adaptive adjustment, which implies that an
unstable Cournot equilibrium under naive expectations could be stabilized byadaptive adjustments.
2 This result reminds us the classical result of Theocaris (1960) that the a non-diﬀerentiated
Cournot equilibrium is stable only in the duopoly case if the expectatitions are naively formed
and the price and the cost functions are linear.

15

(A) α = 1(B) α = 0.7
Figure 3. Stability and feasible regions of the naive and adaptive systems
Figure 4 is an enlargement of the south-west corner of Figure 3B. 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 thatthe flip boundary with α = 0.8 passes through the point (2, 2), the vertex of thetriangular part of the stablility region. Thus if α ≥ 0.8, then the nonnegativity
region is completely outside the stability region for adaptive adjustment, so
no nonnegative equilibrium becomes stable. On the other hand, solving 1 +
p + (−q) + r = 0 for α yields 23 , given k = h = 1. This implies that theflip boundary with α = 23 passes through the point (1, 1). If α ≤ 23 , then thenonnegativity region is entirely inside the stability region, in which case allnonnegative equilibria become stable. If 23 < α < 0.8, then only a certain partof the nonnegativity region belongs to the stability region, and this part becomeslarger 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 circular curve. In the remaining light-gray area of
the nonnegativity region, the Cournot output is locally unstable. When theadaptive parameter α decreases, the flip boundary shifts inside accordingly. As

16

a consequence, the stability region enlarges and the unstable region shrinks.

Figure 4. Dependency of the stability region on α
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 isÃr!
x(t) + 2y(t) + 2z(t)x(t + 1) = (1 − α)x(t) + α− (x(t) + 2y(t) + 2z(t))
ay(t + 1) = (1 − β)y(t) + β

Ãr

!
2x(t) + y(t) + 2z(t)− (2x(t) + y(t) + 2z(t))
bz(t + 1) = (1 − γ)z(t) + γ

Ãr!

2x(t) + 2y(t) + z(t)− (2x(t) + 2y(t) + z(t))
cwhere for simplicity, α = β = γ, h = k and a = 1 are selected implying b = c.
The dynamic system in the three-firm model (i.e., triopoly) can be constructedsimilarly. Selecting α and h (in particular b) as the bifurcation parameters,
Figure 5 illustrates the bifurcation diagrams in the (α, h) plane. The numericalinvestigations clearly reveal qualitatively the same issues as those we saw inFigure 2. Namely, first of all, the destabilizing process goes to chaos througha flip bifurcation with a lower production ratio in the three-group model andthrough a Hopf bifurcation with a higher ratio in the three-firm model. Second,
the adjustment speed can be used to control unstable trajectories. Comparing
the bifurcation diagrams of the two-group and the three-group models shows thesimilarity of the destabilizing process of these models in which period-doubling17

bifurcation takes place.

(A) The three-group model(B) The three-firm model
Figure 5. Two-parameter bifurcation diagrams in the (α, h = k) plane
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 andthe 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-groupmodel is defined by two upward sloping curves, k = fb (h) and k = fc (h). Noticethat 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 takenegative values and thus becomes economically meaningless. One way is toconfine 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 non-negativity constraint explicitly. However, theformer has the diﬃculty of deriving the confinement conditions in our model as
many parameters are involved, and the latter makes the asymptotic behavior ofthe dynamic system significantly diﬀerent and more diﬃcult to analyze. Since
our main concern is to control the unstable trajectories and our main conclusion
is that the adjustment speed is an eﬀective control parameter, which is supposed
to hold in those models, we used the model without such modifications at the18expense of some economic precision.

Figure 6. Two-parameter bifurcation diagram in the (h, k) plane

5

Conclusion

The asymptotic behavior of discrete dynamic systems is a fundamental researchissue. In this paper a special class of economic models was examined. For thesake of mathematical simplicity we selected N -firm Cournot oligopolies without product diﬀerentiation, 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 threegroups of identical firms. Stability conditions could be derived analytically inthe first case, and the dependence of the asymptotic properties of the equilibrium on the number of firms was illustrated by computer simulation in thesecond case. The results of the nonlinear duopoly and triopoly models show that
the Cournot equilibrium can be destabilized through a Hopf bifurcation. We
also found the following new dynamic phenomenon. In the multi-group models,
the stationary state is destabilized through the Feigenbaum period doubling sequence, and a Hopf bifurcation can occur only in the infeasible regions in which
the stationary state is negative. For N > 4, the multi-group models are unstable
if α is close to one and become stable if α is below a certain threshold, regardlessof the production cost ratios. This implies that the main source of instability
is the speed of adjustment and thus the stationary state could be stabilized by
selecting suﬃciently small speed of adjustment. That is, the multi-group models
with N firms are unstable under naive expectations but are controllable withthe adaptive adjustment process in which the speed of adjustment is the control
parameter.

19

Appendix
In this appendix, we derive the stability conditions of the adaptive system
in which the expectation is adaptively formulated,X
xj (t).
yie (t + 1) = (1 − αi )yie (t) + αi
j6=iHere αi ∈ (0, 1] is the speed of adjustment of firm i. The stability conditions
are also useful to determine the local dynamic behavior of the naive system aswell as that of the inertia system.
We consider the adjustment process with adaptive expectations first, since
the one with naive expectations can be obtained by selecting the speed of adjustment equal to unity (i.e., αi = 1). For i = 1, 2, ..., N,
⎧
s P
´
⎪
αi j6=i xj (t) + (1 − αi )yie (t) ³ P⎪⎪⎨ xi (t + 1) =− αi j6=i xj (t) + (1 − αi )yie (t) ,
ci
⎪
⎪
⎪
P
⎩ e
yi (t + 1) = αi j6=i xj (t) + (1 − αi )yie (t).

The Jacobian at the equilibrium has the form⎛0
⎜ γ 2 α2⎜⎜·⎜⎜ γ N αN
J=⎜
⎜
0
⎜
⎜ α2
⎜
⎝
·
αN

γ 1 α1
0
·
·
α1
0
·
·

·
·
·
γ N αN
·
·
·
αN

where

γi ≡

γ 1 α1γ 2 α2
·
0
α1
α2
·
0

∂fi (yie )∂yie⎞
γ 1 (1 − α1 )
0
·
0
⎟
0
γ 2 (1 − α2 ) ·
0
⎟
⎟
·
·
·
·
⎟
0
·
0 γ N (1 − αN ) ⎟⎟,⎟1 − α10·0⎟⎟0
1 − α2
·
0
⎟
⎠
·
·
·
·
0
·
0
1 − αN=C− 1.
2ci (N − 1)

Notice that condition (3) implies that γ i ≥ − 12 for all i.
The eigenvalue equation has the formT

Jx =λx with x = (u1 , ...uN , v1 , ...vN ) ,or equivalently,⎧P⎨ γ i αi j6=i uj + γ i (1 − αi )vi = λui , 1 ≤ i ≤ N,
⎩

αi

P

j6=i uj + (1 − αi )vi = λvi ,1 ≤ i ≤ N.

Subtracting the γ i -multiple of the second equation from the first one gives
λ(ui − γ i vi ) = 0.
The value λ = 0 cannot destroy stability, so we may assume λ 6= 0. Thenui = γ i vi , and by substituting it into the second equation, we haveXγ j vj + (1 − αi )vi = λvi , 1 ≤ i ≤ N.
αi
j6=i

20

This is the usual eigenvalue problem of the N × N matrix⎞⎛
1 − α1 γ 2 α1 · γ N α1
⎜ γ 1 α2 1 − α2 · γ N α2 ⎟
⎟.
H=⎜
⎠
⎝
·
·
·
·
γ 1 αNγ 2 αN · 1 − αN

Notice thatH = D + abT ,with
aT = (α1 , α2 , ..., αN ), bT = (γ 1 , γ 2 , ..., γ N ),and
D = diag(1 − α1 (1 + γ 1 ), 1 − α2 (1 + γ 2 ), ..., 1 − αN (1 + γ N )).
The characteristic polynomial of H can be decomposed by using the simple factthat if x, y ∈ RN , then¡¢det I + xyT = 1 + yT x,

where I is the N × N identity matrix. So we havedet(H − λI) = det(D + abT − λI)= det(D − λI) det(I + (D − λI)−1abT ).

The roots of the first factor are 1 − αi (1 + γ i ) which are inside the unit circle if
and only if−1 < 1 − αi (1 + γ i ) < 1,which occurs if and only if
αi C.
4(N − 1)(A-1)αi γ i= 0.
1 − αi (1 + γ i ) − λ
i=1

(A-2)

ci >

The other eigenvalues are the roots of equation
1+N
X

The Cournot equilibrium is locally asymptotically stable if all eigenvalues areless than unity in absolute value and is unstable if at least one eigenvalue isoutside the unit circle. However, since parameters α1 , α2 , ..., aN can be selected
arbitrarily in interval (0, 1], there is a large flexibility in the location of theeigenvalues.

21

References[1] Agiza, H. N. (1999), "On the Analysis of Stability, Bifurcation, Chaos andChaos Control of Kopel Map," Chaos, Solitons and Fractals, vol. 10, no.
11, 1909-1916.
[2] Agiza, H. N. and A. A. Elsadany (2003), "Nonlinear Dynamics in the
Cournot Duopoly Game with Heterogeneous Players," Physica A, vol. 320,
512-524.
[3] Bischi, G-I, C. Chiarella, M. Kopel and F. Szidarovszky (2010), Nonlinear
Oligopolies: Stability and Bifurcations, Springer, Berlin/Heidelberg/New
York.
[4] Cournot, A. (1838), Recherches surles Principes Mathematiques de la Theorie de Richessess. Hachette, Paris. (English translation (1960): Researches
into the Mathematical Principle of the Theory of Wealth, Kelly, New York)[5] Matsumoto, A. (2006), "Controlling the Cournot-Nash Chaos," Journal of
Optimization Theory and Applications, vol.128, 379-392.
[6] Okuguchi, K. (1976), Expectations and Stability in Oligopoly Models,
Springer-Verlag, Berlin/Heidelberg/New York.
[7] Okuguchi, K. and K. Irie (1990), "The Schur and Samuelson Conditions fora Cubic Equation," The Manchester School of Economic & Social Studies,vol. 58, no.4, 414-418.
[8] Okuguchi, K. and F. Szidarovszky (1999), The Theory of Oligopoly with
Multi-product Firms, 2nd edn., Springer-Verlag, Berlin/Heidelberg/New
York.
[9] Puu, T. (2003), Attractors, Bifurcations and Chaos, 2nd edn., SpringerVerlag, Berlin/Heidelberg/New York.
[10] Puu, T. and I. Sushko (eds) (2002), Oligopoly Dynamics, Springer-Verlag,
Berlin/Heidelberg/New York.
[11] Richter, H. and A. Stolk (2004), "Control of the Triple Chaotic Attractor
in Cournot Triopoly Model," Chaos, Solitons and Fractals, vol. 20, no. 2,
409-413.
[12] Szidarovszky, F. and C. Chiarella (2001), "Dynamic Oligopolies, Stability
and Bifurcation," Cubo Matematica Educacional, vol. 3, no. 2, 267-284.
[13] Theocharis, R. D., (1960) "On the Stability of the Cournot Solution in the
Oligopoly Problem," Review of Economic Studies, vol. 27, 133-134.

22