Discussion Paper No.240
Nonlinear Cournot Duopoly with Implementation Delays
Akio Matsumoto
Chuo University
Ferenc Szidarovszky
University of Pécs
January

2015

INSTITUTE OF ECONOMIC RESEARCH
Chuo University
Tokyo, Japan

Nonlinear Cournot Duopoly
with Implementation Delays∗
Akio Matsumoto†

Ferenc Szidarovszky‡

Abstract
We study the effects of two delays on the local as well as on global
stability of nonlinear Cournot duopoly dynamics. The two major findings
are an analytical construction of the stability switching curve on which
stability is lost and the numerical confirmation of the birth of aperiodic
global dynamics when the stationary state is locally unstable. The delays matters and can generate various dynamics ranging from simple tocomplicated dynamics.
Keywords: Implementation delays, Gradient dynamics, Stability switches,
Cournot duopoly, Hopf bifurcation∗ The authors highly appreciate the financial supports from the MEXT-Supported Program
for the Strategic Research Foundation at Private Universities 2013-207, the Japan Society for
the Promotion of Science (Grant-in-Aid for Scientific Research (C) 24530202, 25380238 and26380316). The usual disclaimers apply.
† Professor, Department of Economics, Senior Researcher, International Center for further
Development of Dynamic Economic Research, Chuo University, 742-1, Higashi-Nakano, Hachioji, Tokyo, 192-0393, Japan; akiom@tamacc.chuo-u.ac.jp‡ Professor, Department of Applied Mathematics, University of Pécs, Ifjúság, u. 6, H-7624,
Pécs, Hungary; szidarka@gmail.com

1

1

Introduction

Oligopoly theory has a long history since the pioneering work of Cournot (1838).
It has played a central role in mathematical economics and developed in various
directions. The existence and uniqueness of the equilibrium were the researchissues in early stage and then dynamic extensions become the main topic of researchers. Linear models were first examined, where local asymptotical stability
implies global stability. Each model is based on a particular output adjustment
scheme. In applying best response dynamics, global information is needed about
the profit function while in the case of gradient adjustments only local information is needed to assess the marginal profit. Okuguchi (1976) summarizes theearly results on static and dynamic oligopolies and Okuguchi and Szidarovszky(1999) discuss their multiproduct generalizations. During the last two decades
an increasing attention has been given to nonlinear dynamics. Bischi et al.
(2010) give a comprehensive summary of the newer developments. Concerning oligopoly dynamics, it is now well-known that oligopoly models formulated
in discrete-time framework can generate various dynamics ranging from simple dynamics to complex dynamics including chaos if endogenous nonlinearities
become stronger. It is also well-known that continuous-time oligopoly models
behave better than discrete-time models in a sense that the former models havelarger stability regions of the parameters. In the existing literature, however,not much has been revealed with the delay oligopoly model that is a hybridof the discrete-time and the continuous-time models.1 In particular, it has not
yet been known what dynamic behavior arises when more than one delays are
involved in economic variables.
The main purpose of this paper is to provide rigorous theoretical analysis
on multiple-delay dynamic behavior and to confirm the analytical results byperforming numerical simulations. For this end, we construct an á la Cournot
duopoly model in which each firm has a delay in implementing information
about its own output. This is a continuity of Matsumoto et al. (2011) in whicha linear duopoly model with two fixed delay is considered. It is a special case of
Howroyd and Russell (1984) in a sense that the number of the firms is reduced to
two from a general N and, at the same time, considered to be an extended casein a sense that the number of delays is increased to two from one. It is also anapplication of the results obtained in delay dynamic monopolies of Matsumoto
and Szidarovszky (2014a, b). We apply an analytical method developed by Lin
and Wang (2012) and construct a stability switching curve with respect to twodelays. We provide some examples in which periodic as well as aperiodic delay
dynamics may appear in a 2D system of differential equations.
The paper is organized as follows. In Section 2, the traditional continuoustime duopoly model is presented and then delays are implemented. The stability
switching curves are analytically derived. In Section 3, we give numerical ex1 The delay equations describing the population of the species are frequently discussed intheoretical biology. Even in early 80s, Shibata and Saito (1980) investigate the time delayed
saturation effect in a logistic model of two competing species and numerically confirm thebirth of chaotic solutions.

2amples to illustrate the theoretical findings and to find how the delays affectdynamic behavior. In the final section, concluding remarks are given.

2

Delay Model

We consider an industry of two firms producing a homogeneous good. Let x
and y be the output quantities produced by firms x and y, respectively. The
price function is assumed to be linear,p = a − b(x + y) with a > 0 and b > 0.
Production costs are also assumed to be linear and the marginal costs are denoted by cx and cy . The profit function of firm z (z = x, y) is given by
π z = (a − b(x + y))z − cz z.
Firm z determines its output to maximize profit, taking quantity supplied by
the competition as given. Assuming an interior optimal solution and solving thefirst-order condition for the profit maximizing problem yield the best replies (orreaction functions),
Rx (y) =11(a − cx − by) and Ry (x) = (a − cy − bx) .
2
2

A Cournot equilibrium is a pair of outputs (xe , y e ) satisfying xe = Rx (ye ) andye = Ry (xe ), that is,xe =a − 2cx + cy
a − 2cy + cx
and y e =3b3bwhere the following is assumed to guarantee positive equilibrium,a > max[2cx − cy , 2cy − cx ].
Concerning the adjustment process of output, we make two assumptions.
One is that the firms are boundedly rational in a sense that they modify theirproduction according to the sign of the marginal profits. That is, the firm
increases its production if the marginal profit is positive, decreases it if negativeand does not change it if zero. Such an adjustment process is called gradient
dynamics and satisfies the following condition,
dπ z
sign [ż(t)] = sign.
dz
Specifically, two formulations are often used,
ż(t)
dπz
dπz
= kz
or ż(t) = fz (z)with fz′ (z) > 0,
z(t)
dz
dz
3

where the growth rate of the variable is portional to the marginal profit in thefirst formulation and the adjustment rate is determined in the second formulation by a product of the marginal profit and the degree of the adjustment that
depends on the level of the variable. Notice that the two are mathematicallyequivalent if fz (z) = kz z although their economic interpretations are different.
The other assumtption is that the firms obtain only delayed information on production activities. There are two ways to introduce a delay in our framework.
The first is a delay in implementing information about firm’s own output andthe second is a delay in obtaining information about the competitor’s output.
However, it has been shown by Howroyd and Russell (1984) that the second
delay is harmless in a Cournot oligopoly model.2 Thus we focus on only the
first delay in this study to simplify the analysis. The two assumptions can begiven in the following way.
Assumption 1. The firms continuously adjust their output growth rates
proportional to changes in their profits.
Assumption 2. The firms experience delays in implementing information
about their own outputs.
With positive adjustment coefficients αx and β y of firms x and y, this gradient dynamic system is described by a 2D system of delay differential equations,ẋ(t) = αx x(t) [a − cx − 2bx(t − τ x ) − by(t)] ,(1)ẏ(t) = β y y(t) [a − cy − bx(t) − 2by(t − τ y )] ,where τ x ≥ 0 and τ y ≥ 0 are time delays. The positive steady state of system(1) is identical with Cournot equilibrium (xe , ye ).3 The linear approximation ina neighborhood of the equilibrium isẋδ (t) = α [−2bxδ (t − τ x ) − byδ (t)] ,(2)ẏδ (t) = β [−bxδ (t) − 2byδ (t − τ y )] ,wherexδ (t) = x(t) − xe , yδ (t) = y(t) − ye , α = αx xe and β = β y y e .
Substituting exponential solutions
x(t) = eλt u and y(t) = eλt v,
2 A more general result obtained in Hofbauer and So (2000).
3 This model has multiple equilibria. Other than the positive equilibrium, there are one
zero solution and two corner solutions that correspond to monopoly equilibrium. Our concern
is on the positive equilibrium in this study.

4into the linearized equations in (2) and arranging terms yield a 2D simultaneous
equation system with respect to u and v,


  
λ + 2αbe−λτ x
αb
u
0


 =  .
v
0βb
λ + 2βbe−λτ yThus the corresponding characteristic equation is

P (λ, τ x , τ y ) = P0 (λ) + P1 (λ)e−λτ x + P2 (λ)e−λτ y + P3 (λ)e−λ(τ x +τ y ) = 0 (3)where

P0 (λ) = λ2 − αβb2 ,P1 (λ) = 2αbλ,P2 (λ) = 2βbλ,P3 (λ) = 4αβb2 .

Before proceeding, we examine the non-delay case with τ x = τ y = 0 in whichcase the characteristic equation is simplified as
λ2 + 2b(α + β)λ + 3b2 αβ = 0.
Since all coefficients of this equation and the discriminant are positive, thecharacteristic roots are real and negative, implying that the steady state withno delays is locally stable. We can mention that it is still stable as long as
positive delays are sufficiently small by continuous dependency of λ on thevalues of delays. To see what extent the stationary point can preserve stability,
we determine the threshold values of the delays for which stability is just lost.
The characteristic equation (3) is now investigated with a constructive method
developed by Lin and Wang (2012) to solve a class of two delay differential equations. We look for a pair of the delays for which the characteristic equation has
purely imaginary roots. The set of such pairs is called stability switching curve.
Since λ = 0 is not a solution and roots of a real function come in conjugate
pairs, we assume, without loss of generality, that λ = iω with ω > 0 on thiscurve. Substituting this solution into (3), we have two different forms,

 

(4)
P0 (iω) + P1 (iω)e−iωτ x + P2 (iω) + P3 (iω)e−iωτ x e−iωτ y = 0

and


 

P0 (iω) + P2 (iω)e−iωτ y + P1 (iω) + P3 (iω)e−iωτ y e−iωτ x = 0(5)Using equations (4) and (5), we can derive the following theorem that yields
explicit expressions for the stability switching curves.

5

Theorem 1 Suppose that P2 P̄3 − P0 P̄1 = 0 and P1 P̄3 − P0 P̄2 = 0. Then thestability switching curves are the loci of LR (m, n) and LB (m, n) defined by−
LR (m, n) = (τ +x (ω, m), τ y (ω, n))ω ∈ Ω1 , m, n ∈ Z+
LB (m, n) = (τ −x (ω, m), τ y (ω, n))ω ∈ Ω2 , m, n ∈ Z

andwhere Ω1 and Ω2 are the sets of ω > 0 for which relations (A-5) and (A-6)/ of
the Appendix are satisfied, furthermoreτ±
x (ω, m) =±ψx (ω) − ϕx (ω) + 2mπ
,
ω

and

±ψy (ω) − ϕy (ω) + 2nπ,ωwith ϕx (ω) = arg(P2 P̄3 − P0 P̄1 ), ϕy (ω) = arg(P1 P̄3 − P0 P̄2 ),τ±y (ω, n) =ψx (ω) = cos−1|P0 |2 + |P1 |2 − |P2 |2 − |P3 |2, ψx ∈ [0, π],
2 |B1 (ω)|ψy (ω) = cos−1|P0 |2 − |P1 |2 + |P2 |2 − |P3 |2, ψy ∈ [0, π].
2 |B2 (ω)|

and

Proof. Proof is given in the Appendix.

To check the existence and domain of ω for which inequalities (A-5) and(A-6) (given in the Appendix) hold, we introduce two functions,2
F (ω) = |P0 |2 + |P1 |2 − |P2 |2 − |P3 |2 − 4B12 ≤ 0
and

2

2
2
2
2
G(ω) = |P0 | − |P1 | + |P2 | − |P3 |− 4B22 ≤ 0.

It can be shown that the right hand side expressions of F (ω) and G(ω) areequivalent and, with the notation of z = ω/b, are reduced toz 8 + 4a3 z 6 + 2a2 z 4 + 4a1 z 2 + a0
wherea3 = −2α2 + αβ − 2β 2 ,a2 = 8α4 + 8α3 β − 19α2 β 2 − 8β 4 + 16α2 β(4β − α),a1 = −30α4 β 2 − 15α3 β 3 + 30α2 β 4 − 4α2 β 2 (4β − α)2 ,a0 = 225α4 β 4 .
6Let Z = z 2 andf (Z) = Z 4 + 4a3 Z 3 + 2a2 Z 2 + 4a1 Z + a0 .
Solving f(Z) = 0 for Z yields four solutions,Z1 = 2α2 − αβ + 2β 2 + 2 |α − β|Z2 = 2α2 − αβ + 2β 2 − 2 |α − β|Z3 = 2α2 − αβ + 2β 2 + 2 |α − β|Z4 = 2α2 − αβ + 2β 2 − 2 |α − β|α2 + αβ + β 2 ,
α2 + αβ + β 2 ,
α2 − 3αβ + β 2 ,
α2 − 3αβ + β 2 .

Let us denote the discriminant of the last two equations by D = α2 − 3αβ +β 2 . Notice that Z3 and Z4 are real if D ≥ 0 and complex if D < 0, and thefollowing result is the simple consequence of the above expressions:Lemma 1 Given α > 0 and β > 0,
If α = β,then 0 < Z1 = Z2 = Z3 = Z4 ;If α = β and D ≥ 0, then 0 < Z2 < Z4 ≤ Z3 < Z1 ;If α = β and D < 0, then 0 < Z2 < Z1 .
This lemma immediately implies the following result:
Theorem 2 Given α > 0, β > 0 and b > 0,
If α = β,then F (ω) ≥ 0 for all ω ≥ 0;

If α = β and D ≥ 0, then F (ω) ≤ 0 for ω ∈ [ω 2 , ω 4 ] and ω ∈ [ω 3 , ω1 ];If α = β and D < 0, then F (ω) ≤ 0 for ω ∈ [ω 2 , ω 1 ]where for i = 1, 2, 3, 4,
ωi = zi b with zi =Zi .

In examining the sign of the discriminant D, we solve α2 − 3αβ + β 2 = 0 for
α giving two solutions describing two lines with positive slopes,
√
√
3− 5
3+ 5
α=
β and α =β.
22The β = α line divides the positive quadrant of the (α, β) plane into two parts
as shown in Figure 1. Each part is further divided into two parts by one ofthese two lines since the first line is steeper than the diagonal and the second7

is flatter. Since the situation is symmetric with respect to the diagonal line, we
focus on the left part to it by making the following assumption.
Assumption 3: α < β√It is clear that in Figure 1, D > 0 in the left to the α = 3−2 5 β line andD < 0 in the right.

Figure 1. Division of the (α, β) plane
We reveal the effects caused by the changing values of the adjustment speedson the shapes of the stability switching curves. Figures 2(A), (B) and (C)illustrate the F (ω) curves with β = 1,
√
√
√
3− 5
3− 5
3− 5
(A) α =
β − 0.01, (B) α =β and (C) α =β + 0.01.
2
2
2
As is already mentioned in Theorem 2, F (ω) has a finite number of roots, inparticular, there are four intersections with the axis of abscissa in Figure 2(A)with D > 0, two intersections and one tangency point in Figure 2(B) with
D = 0 and two intersections in Figure 2(C) with D < 0.

(A) Four roots(B) Three roots
Figure 2. Shapes of F (ω)8(C) Two roots

Accordingly, the corresponding stability switching curves are determined.
For simplicity we confine attention to the case of m = 0, 1 and n = 0 andillustrate LR (0, 0) and LB (1, 0) as red and blue curves, respectively, in Figures
3(A), (B) and (C). We start with Figure 3(B), which is a boundary case between
Figures 3(A) and 3(C). Since D = 0, F (ω) ≤ 0 for ω ∈ [ω2 , ω1 ] with ω 2 < ω 4 =ω3 < ω 1 . The continuous red curve LR (0, 0) consists of the two segments, one
defined over the interval [ω 2 , ω4 ] is unimodal and the other defined over theinterval [ω3 , ω1 ] is positive-sloping. It has a kink at point (τ 0x , τ 0y ) with0−
τ 0x = τ +x (ω m , 0) and τ y = τ y (ω m , 0)

where ω m = ω4 = ω 3 . In the same way, the blue curve LB (1, 0) has two segments, one over [ω 2 , ω 4 ] has a steep positive slope and the other over [ω 3 , ω 1 ]has a flat negative slope. These segments are connected at the same kinked
point,
0+τ 0x = τ −x (ωm , 1) and τ y = τ y (ωm , 1).
√
When the value of α increases from (3 − 5)β/2, the discriminant becomes
negative so the kinked point disappears. In consequence, LR (0, 0) as well as
LB (1, 0) defined over interval [ω 2 , ω 3 ] become√ smoother as depicted in Figure
3(C). When the value of α decreases from (3− 5)β/2, the discriminant becomes
positive and the interval is separated into two distinct intervals, [ω 2 , ω 4 ] and[ω3 , ω1 ] with ω 3 > ω4 . Further, the kinked point (τ 0x , τ 0y ) is broken into two
different points, (τ 1x , τ 1y ) and (τ 2x , τ 2y ). In particular, the stability switching curve
defined over [ω 2 , ω4 ] takes an inequality-wise shaped curve and has the blue
curve LB (1, 0) connected with the red curve LR (0, 0) at the point (τ 1x , τ 1y ) where−1−+τ 1x = τ +x (ω 4 , 0) = τ x (ω 4 , 1) and τ y = τ y (ω 4 , 0) = τ y (ω4 , 1).

On the other hand, the stability switching curve defined over [ω 3 , ω 1 ] takes aflat-roof shaped curve and has the blue curve LB (1, 0) connected with the red
curve LR (0, 0) at the point (τ 2x , τ 2y ) where−2−+τ 2x = τ +x (ω 3 , 0) = τ x (ω 3 , 1) and τ y = τ y (ω 3 , 0) = τ y (ω3 , 1).

(A) D > 0

(B) D = 0

(C) D < 0

Figure 3. Stability crossing curves for m = 0, 1 and n = 0
93Simulations

We now perform simulations to illustrate the theoretical findings. For this purpose, we specify some of the parameters in (1) as follows:
Assumption 4: a = 4, b = 1, cx = 1, cy = 1 and β y = 1
Under this assumption, we havexe = y e = 1, β = 1 and α = αx .
√In Section 3.1, we take αx = 0.5 for which D < 0 as (3 − 5)/2 ≃ 0.382, selectthree values of τ y and examine how changing the value of τ x affects dynamics
generated by system (1). Then in Section 3.2, we decrease the value of αx to
0.35 for which D > 0 and repeat the same procedure.

3.1

Simulation Study I

The adjustment coefficient of firm x is taken to be αx = 0.5 in this subsection.
As in Figure 3(C), we have the smoother red curve LR (0, 0) and the blue curve
LB (1, 0).4 A larger value of α shifts the blue curve in Figure 3(C) rightward
and the red curve leftward, resulting in the location of the curves as shown in
Figure 4. Since it is shown that the stationary point is locally stable in theyellow region below the red curve, the red curve is the stability switching curveon which stability is lost. We select three different values of τ y and perform
three simulations with the increasing value of τ x from zero to detect the effects
caused by the delays.

Figure 4. Division of the time delay plane
In the first example we take τ y = 0.5 and increase the value of τ x along
the dotted line that crosses the red curve at τ x = τ 0.5 (≃ 1.348) as shown in
4 The segments of L (0, 0) and L (1, 0) do not appear in this region of (τ , τ ).
x y
B
R

10

Figure 4. A bifurcation diagram with respect to τ x is given in Figure 5(A).
It is constructed in the following way. With fixed value of τ x , we run the
delay system (1) with the specified parameter values for 0 ≤ t ≤ 1000. To takeaway the initial disturbances, we discard the numerical data of x(t) and y(t)for t ≤ 950 and plot the local maximum and local minimum of y(t) obtained
from the data for 950 ≤ t ≤ 1000 against this value of τ x . Then we increase thevalue of τ x and repeat the same procedure until τ x arrives at 2.8. It is shown
there that until τ x reaches τ 0.5 , the stationary point is asymptotically stable,loses stability at τ x = τ 0.5 and a limit cycle of y(t) emerges for τ x > τ 0.5 . To
observe how cyclic fluctuations emerge, we gradually increase further the valueof τ x . We find first a unique minimum and maximum, then two minima andmaxima, even further increase of τ x results in three minima and maxima, andso on. Accordingly, we have a regular cycle, which is distorted more and more
to be a cycle having several ups and downs by the increase of τ x . For example,
the vertical dotted line at τ x = 2 crosses the bifurcation diagrams eight times,
four time with the upper branches of the diagram and four times with the lower
branches. This implies that the motion of y(t) is oscillatory and the trajectoryhas four increasing segments and four decreasing segments. In Figure 5(B) aperiodic solution for 900 ≤ t ≤ 1000 is illustrated in the (x, y) plane in red while
the two best reply lines are shown in black.

(A) Bifurcation diagram(B) Distorted limit cycle

Figure 5. Delay effects of τ x , given τ y = 0.5
In the second example, we increase the value of τ y to 0.9. The horizontal
dotted line at τ y = 0.9 in Figure 4 crosses the red curve two times at τ x = τ a0.9 (≃0.593) and τ x = τ b0.9 (≃ 1.058) and the blue curve one time at τ x = τ c0.9 (≃1.885). A bifurcation diagram of y(t) with respect to τ x is given in Figure 6(A)and shows that increasing the value of τ x along the dotted line gives rise to a wide
variety of dynamics of y(t) ranging from a stable trajectory to chaotic motions.
There, a limit cycle emerges for τ x < τ a0.9 , stability is gained in the interval(τ a0.9 , τ b0.9 ) for which the dotted line is located in the yellow region of Figure
4. Stability is lost again at τ x = τ b0.9 and bifurcates to a periodic cycle having
11

several local maxima and minima for τ x < τ c0.9 and finally more complicateddynamics arises for larger values of τ x . Figure 6(B) illustrates a phase diagramfor 9500 ≤ t ≤ 10000 in the (log[x(t)],log[y(t)]) plane for τ x = τ̄ (= 2) in which
output exhibits many ups and downs.

(A) Bifurcation diagram(B) Chaotic trajectory

Figure 6. Delay effects of τ x , given τ y = 0.9
In the third example, the value of τ y is further increased to 1.1 at which thedotted horizontal line starts and crosses the blue curve at τ x = τ 1.1 (≃ 1.409) in
Figure 4. A bifurcation diagram shown in Figure 7(A) indicates that the steadystate is locally unstable for τ x = 0 and the dynamic system (1) generates alimit cycle for τ x > 0. The limit cycle then bifurcates to complicated cycles
through a period-doubling-like cascade and then comes back to a simple limit
cycle through period-halving-like cascade for τ x ≤ τ 1.1 . Further increasing
τ x indicates transformation to aperiodic oscillations from periodic oscillations.
Figure 7(B) illustrates a phase diagrams of log[x(t)] and log[y(t)] for 9500 ≤ t ≤10000 when τ x = τ̄ (= 2.25).

12

(A) Bifurcation diagram(B) Chaotic movement

Figure 7. Delay effects of τ x , given τ y = 1.1

3.2

Simulation Study II√
In this section we present numerical simulations when α = 0.35 < (3 − 5)/2.
The last inequality implies D > 0 and as in Figure 3(A), there are two independent intervals, [ω 2 , ω 4 ] with ω 2 ≃ 0.563 and ω4 ≃ 1.081 and [ω 3 , ω 1 ]with ω3 ≃ 1.619 and ω 1 ≃ 1.863. For ω ∈ [ω2 , ω 4 ], the red curve LR (0, 0)
is connected to the blue curve LB (1, 0) and thus both curves take distorted
inequality-shaped profiles while only a small portion of the blue curve LB (1, 0)
is illustrated in Figure 8. On the other hand, for ω ∈ [ω 3 , ω 1 ], the red curve
LR (0, 0) and the blue curve LB (0, 0) together construct a flatterly roof-shaped
profile. In the yellow region surrounded by these curves, the stationary point islocally stable and its stability is lost on these curves. To examine what kinds ofdynamics emerge for (τ x , τ y ) in the white region, we will perform three numerical simulations as before.

13

Figure 8. Dividion of the (τ x , τ y ) plane
In the first simulation, we take τ y = 0.5 and increase the value of τ x along thedotted line that crosses the downward sloping red curve at τ x = τ 0.5 (≃ 1.934)
in Figure 8. A bifurcation diagram is presented in Figure 9(A) in which the
dynamics seems to be similar to the one in Figure 5(A), that is, stability is lostat τ x = τ 0.5 and a limit cycle is born for τ x > τ 0.5 . It is seen that the numberof increasing and decreasing segments increases as the value of τ x increases.
In Figure 9(B) a time trajectory of y(t) for 955 ≤ t ≤ 1000 is depicted when
τ x = 3. Within one cycle, the trajectory hits maximum and minimum seventimes each.

(A) Bifurcation diagram(B) Time trajectory

Figure 9. Delay effect of τ x , given τ y = 0.5
In the second example, τ y is increased to 0.9. As can be seen in Figure
8, the dotted line at τ y = 0.9 crosses the positive-sloping red curve at τ x =BτA
0.9 (≃ 0.954), the negative sloping blue curve at τ x = τ 0.9 (≃ 1.467) and the14negative sloping red curve at τ x = τ C
0.9 (≃ 1.743). So as the value of τ x increases
Bfrom zero, the stationary state gains stability for τ A
0.9 < τ x < τ 0.9 and loses
it otherwise. The bifurcation diagram presented in Figure 10(A) reveals whatkind of dynamics emerges when stability is lost. A simple limit cycle arisesfor τ x < τ A
0.9 while more complex dynamics emerges via a period-increasing
cascade for τ x > τ B
0.9 . With further increasing the value of τ x , a very long
periodic cycle or chaotic behavior might emerge and its oscillations are aperiodic
but not erratic. Figure 10(B) gives time trajectories for τ x = τ̄ (= 2.8), thered one for log[y(t)] and the blue one for log[x(t)]. A phase diagram in the(log[x(t)], log[y(t)]) plane is very similar to the one shown in Figure 9(B).

(A) Bifurcation diagram(B) Time trajectories

Figure 10. Delay effect of τ x , given τ y = 0.9
In the third example, τ y is further increased to 1.1 at which the dottedhorizontal line crosses the positive-sloping segment of the red curve at τ x =τ 1.1 (≃ 1.624). Since the line is located in the white region of Figure 8, thestationary state is locally unstable for any τ x ≥ 0 and a bifurcation diagramin Figure 11(A) shows that periodic and aperiodic motions of y(t) alternate
as the value of τ x increases. The time trajectories of log[x(t)] and log[y(t)] atτ x = τ̄ (= 2.5) are illustrated again as the blue and the red curves in Figure15

11(B). They are aperiodic but not erratic.

(A) Bifurcation diagrams(B) Time trajectories

Figure 11. Delay effect of τ x , given τ y = 1.1

4

Concluding Remarks

In this paper we have focused on dynamics of the Cournot duopoly model inwhich each firm has a delay in implementing information about its own output.
First we analytically derived the stability switching curve that divides the delay
region into two regions, one in which stability is preserved and the other inwhich stability is lost. Then we numerically examined dynamic behavior of
output when the stationary state loses stability. It is found that the shape ofthe stability switching curve is parameter-dependent. It is demonstrated thatthe two delay model can generate rich dynamics involving chaotic oscillations.
This finding implies that the delay can be a source of periodic as well as aperiodic
dynamics in the continuous-time framework. It may be possible to extend thenumerical results in the cases of higher dimension. However the construction of
the stability switching surfaces in higher dimesions is still in progress.

16

Appendix
Following Lin and Wang (2012), we now present the proof of Theorem 1.
Since e−iωτ x = e−iωτ y = 1, equations (4) and (5) can be written as
P0 (iω) + P1 (iω)e−iωτ x = P2 (iω) + P3 (iω)e−iωτ x
andP0 (iω) + P2 (iω)e−iωτ y = P1 (iω) + P3 (iω)e−iωτ y .
Using the relationa + be−iωτ2

we can rewrite (4) as
2

2

2




= a + be−iωτ ā + b̄eiωτ ,2|P0 | + |P1 | − |P2 | − |P3 | = 2A1 (ω) cos(ωτ x ) − 2B1 (ω) sin(ωτ x )(A-1)

and (5) as|P0 |2 − |P1 |2 + |P2 |2 − |P3 |2 = 2A2 (ω) cos(ωτ y ) − 2B2 (ω) sin(ωτ y )

(A-2)

where
A1 (ω) = Re[P2 P̄3 − P0 P̄1 ] and B1 (ω) = Im[P2 P̄3 − P0 P̄1 ]andA2 (ω) = Re[P1 P̄3 − P0 P̄2 ] and B2 (ω) = Im[P1 P̄3 − P0 P̄2 ].
It is clear that both P2 P̄3 − P0 P̄1 and P1 P̄3 − P0 P̄2 are pure complex with zeroreal parts. So A1 (ω) = A2 (ω) = 0. Assume B1 (ω) = 0 and B2 (ω) = 0 and thencompute the arguments of these complex values as π 2 if B1 (ω) > 0,ϕx (ω) := arg(P2 P̄3 − P0 P̄1 ) = 3π if B1 (ω) < 02andϕy (ω) := arg(P1 P̄3 − P0 P̄2 ) = π 2 if B2 (ω) > 0,
 3π if B2 (ω) < 0.
2
With A1 (ω) = 0 and cos ϕx (ω) = 0, the right-hand side of equation (A-1) canbe rewritten as−2B1 (ω) sin ωτ x

= −2 |B1 (ω)| sin ϕx (ω) sin ωτ x
= 2 |B1 (ω)| [cos ϕx (ω) cos ωτ x − sin ϕx (ω) sin ωτ x ]= 2 |B1 (ω)| cos[ϕx (ω) + ωτ x ].
17

The right-hand side of equation (A-2) can be rewritten in the same way. Therefore, (A-1) and (A-2) becomes
2

2

2

2

(A-3)

|P0 |2 − |P1 |2 + |P2 |2 − |P3 |2 = 2 |B1 (ω)| cos(ϕy (ω) + ωτ y ).

(A-4)

|P0 | + |P1 | − |P2 | − |P3 | = 2 |B1 (ω)| cos(ϕx (ω) + ωτ x )andA sufficient and necessary condition for the existence of τ x > 0 and τ y > 0satisfying the above equations areand|P0 |2 + |P1 |2 − |P2 |2 − |P3 |2 ≤ 2 |B1 (ω)|(A-5)

|P0 |2 − |P1 |2 + |P2 |2 − |P3 |2 ≤ 2 |B1 (ω)| .

(A-6)

Denote by Ω1 and Ω2 the set of ω > 0 values satisfying (A-5) and (A-6), respectively. Let us define new functions bycos[ψx (ω)] =|P0 |2 + |P1 |2 − |P2 |2 − |P3 |2, ψx ∈ [0, π]
2 |B1 (ω)|(A-7)

cos[ψy (ω)] =|P0 |2 − |P1 |2 + |P2 |2 − |P3 |2, ψy ∈ [0, π],
2 |B1 (ω)|

(A-8)

and

then from (A-3) and (A-7), (A-4) and (A-8), we havecos[ϕx (ω) + ωτ x ] = cos[ψx (ω)]
andcos[ϕy (ω) + ωτ y ] = cos[ψy (ω)]
both of which yield, for arbitrary integers m and n,τ±
x (ω, m) =±ψx (ω) − ϕx (ω) + 2mπωand±ψy (ω) − ϕy (ω) + 2nπ.
ωThese are the values of the delays that solve equations (4) and (5). Therefore,−1−+given m and n, (τ +x (ω, m), τ y (ω, n)) for ω ∈ Ω and (τ x (ω, m), τ y (ω, n)) for2ω ∈ Ω construct the stability switching curves in the (τ x , τ y ) plane. Q.E.D.
τ±y (ω, n) =18

References
Bischi, G-I, Chiarella, C., Kopel, M. and Szidarovszky, F. Nonlinear Oligopolies:
Stability and Bifurcations, 2010, Heidelberg Dordrecht, London, New York,
Springer.
Cournot, A. Recherches sur les Principes Mathématiques de la Théorie desRichessess, 1838, Hachette, Paris (English translation (1960): Researches
into the Mathematical Principles of the Theory of Wealth, Kelley, NewYork).
Hofbauer, J. and So, J. Diagonal Dominance and Harmless Off-diagonal Delays,
Proceedings of the American Mathematical Society, 128, 2000, 2675-2682.
Howroyd, T. and Russell, A. Cournot Oligopoly Models with Time Delays,
Journal of Mathematical Economics, 13, 1984, 7-103.
Lin, X. and Wang, H., Stability Analysis of Delay Differential Equations with
Two Discrete Delays, Canadian Applied Mathematics Quarterly, 20, 2012,
519-533.
Matsumoto, A. and Szidarovszky, F. Learning Monopolies with Delayed Feedback on Price Expectations, IERCU Discussion Paper #223 (http://www.chuou.ac.jp/research/institutes/economic/publication/discussion/pdf/discussno223.pdf),
Institute of Economic Research, Chuo University, 2014a.
Matsumoto, A. and Szidarovszky, F., Discrete and Continuous Dynamics in
Nonlinear Monopolies, Applied Mathematics and Computation, 232, 632642, 2014b.
Matsumoto, A., Szidarovszky, F. and Yoshida, H. Dynamics in Linear Cournot
Duopolies with Two Time Delays, Computational Economics, 38, 2011,311-327.
Okuguchi, K. Expectations and Stability in Oligopoly Models, 1976, Berlin,
Springer.
Okuguchi, K. and Szidarovszky, F. The Theory of Oligopoly with Multi-productFirms (2nd ed), 1999, Berlin, Springer.
Shibata, A. and Saito, N. Time Delays and Chaos in Two Competitng Species,
Mathematical Biosciences, 51, 1980, 199-211

19