- Email: [email protected]
Delay dynamics of a Cournot game with heterogeneous duopolies
Applied Mathematics and Computation 269 (2015) 699–713
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Delay dynamics of a Cournot game with heterogeneous duopolies✩ Akio Matsumoto a,∗, Ferenc Szidarovszky b a 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 b Department of Applied Mathematics, University of Pécs, Ifjúság u. 6, Pécs, H-7624, Hungary
a r t i c l e
i n f o
Keywords: Delay duopoly dynamics Bounded rationality Continuous-time system Stability switch
a b s t r a c t We construct a Cournot duopoly model in a continuous-time framework and consider its dynamic behavior when the firms are heterogeneous in determining their output decision: one firm has an information delay in the competitor’s output as well as an implementation delay in its own output while the other firm does not have any delays. Two main results are obtained. One is that the information delay does not affect dynamics of the commodities and the other is that the implementation delay can destabilize the otherwise stable equilibrium. Furthermore, it is numerically confirmed that the two delay duopoly model can generate a wide spectrum of dynamics ranging from cyclic dynamics to chaos. © 2015 Elsevier Inc. All rights reserved.
1. Introduction Oligopoly theory is one of the most frequently studied research areas in mathematical economics. Based on the pioneering work [6] many scientists devoted their effort to this interesting and important model family. The earlier results up to the mid 70s are summarized in [21] and their multiproduct generalizations are presented with additional applications in [22]. Recently nonlinear models become the main focus, and a summary of more recent developments is presented in [5]. There is a wide variety of oligopoly models: single product models with and without product differentiation, multiproduct models, labor managed oligopolies among others. In most studies instantaneous information is assumed to the firms about the price and the output levels of the competitors, however in real economies there is always a time delay due to information collection, decision making and implementation. Time delays with continuous time scales can be modeled as fixed or continuously distributed delays. Differential-difference equations are used to model fixed delays [2] while continuously distributed delays are described by Volterra-type integro-differential equations [7]. The central issue in examining dynamic models with delays is the way how the delays affect the asymptotic behavior of the systems. If discrete time scales are assumed, then the lengths of the delays are assumed to be integer multiples of the time unit and higher-order difference equations model the dynamism of the systems. This paper presents a new characterization of price formation of a single commodity by which regular and irregular price oscillations can be constructed in a continuous-time framework with fixed time delays. In the existing literature of macroeconomic dynamics, it has been well-known since the 1930s that a delay in investment or production is one of the key factors to ✩ The authors are indebted to two anonymous referees for their comments that are most helpful in preparing the final version. They highly appreciate the financial supports from the MEXT-Supported Program for the Strategic Research Foundation at Private Universities 2013–2017, the Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research (C) 24530202, 25380238 and 26380316) and Chuo University (Joint Research Grant). The usual disclaimers apply. ∗ Corresponding author. Tel.: +81 42 674 3351; fax: +81 42 674 3425. E-mail addresses: [email protected] (A. Matsumoto), [email protected] (F. Szidarovszky).
http://dx.doi.org/10.1016/j.amc.2015.07.051 0096-3003/© 2015 Elsevier Inc. All rights reserved.
700
A. Matsumoto, F. Szidarovszky / Applied Mathematics and Computation 269 (2015) 699–713
generate cyclic dynamics in national income, capital accumulation, employment rate, etc. See the pioneering papers, [10] and [15]. On the other hand, in the existing literature of microeconomic dynamics such as price and commodity oscillations, it has been less-known that production delays are also able to generate economic fluctuations of a continuous-time unstable economic system.1 As early as in the 1930s, it has been already shown in [12] that cyclic behavior can arise in a simple price dynamics model with production delay, coaxing an idea from theoretical biology.2 A continuous-time cobweb model is constructed in [16] for the pork market in which pork is inelastically supplied with 12 month delay and gives rise to cyclic behavior. [14] is the first work to consider delay dynamics in a N-firm oligopoly model. On the other hand, [17] formulates the price dynamics of a single commodity market as a nonlinear delay differential equation and rigorously derives the stability conditions and the birth of a cyclic oscillation via Hopf bifurcation. Recently [18,19] reconsider delay dynamics in various monopoly models and emphasize that delay may explain various dynamic behavior of micro economic variables. In the existing literature, however, analysis is mostly confined to cases of a single dimensional model with one delay or multiple delay or a multiple dimensional model with one delay. In this paper we construct a multiple dimensional model with multiple delays and examine cyclic behavior in the model. This paper revisits the delay duopoly game investigated in [8] in order to reconsider the delay effect on stability on the Nash equilibrium in a continuous-time framework. The game in [8] is built in a discrete-time framework and its distinguished feature is the heterogeneity of the duopolists. In predicting competitor’s output, one duopolist uses a combination of the current and delayed information while the other adopts naive expectation. The main result is that the delay has a stabilizing effect: it enlarges the stability region of the Nash equilibrium. Generally speaking, a continuous-time model has a larger stability region than a discrete-time model (see, for example, [20]). Furthermore, the delay in a continuous-time model often exhibits stability switching from stability to instability if it becomes larger than some threshold value. However, it is not known yet whether such stability switch occurs in the heterogenous duopoly game. Hence in this paper, we retain the heterogeneity of the discrete-time dynamic system and consider its effect in a delay continuous-time framework . The main results can be summarized as follows: (I) When only one firm has information delays in the competitor’s output, then this heterogeneous expectation formation has no destabilizing effect; (II) When the firm has a delay in implementing information about its own output, the otherwise stable equilibrium can be destabilized and give rise to simple as well as to complex fluctuations when the delay is sufficiently long. This paper is organized as follows. In Section 2, the basic elements of the model of [8] are recapitulated. In Section 3 the corresponding continuous-time duopoly model is constructed and the delay effects due to heterogenous expectation formation are examined. In Section 4, the destabilizing effect caused by the implementation delay is analytically and numerically examined. In the final section some concluding remarks are given and further research directions are outlined. 2. Discrete-time duopoly model A duopoly game with heterogeneous bounded rationality is formulated in discrete-time framework. Two firms produce a homogenous good. Firm 1 produces quantity x with marginal cost cx , while firm 2 produces quantity y at marginal cost cy . The market demand is linear and depends on the total output of the industry,
p = a − b(x + y)
(1)
where a and b are positive constants. The profit of firm k is
πk = [a − b(x + y)]z − cz z.
(2)
where k = 1, 2 and z = x, y. [4] assumes bounded rational firms and introduce gradient dynamics in which the firms adjust production levels according to their marginal profits in such a way that a firm increases output if the marginal profit is positive, decreases it if it is negative and does not change it if it is zero. The gradient dynamics with boundedly rational firm k can be described by
z(t + 1) = z(t ) + αk z(t )
∂πk ∂z
(3)
where α k > 0 is an adjustment coefficient. Hence in the duopoly game, the output adjustment process with gradient method is presented by a two dimensional system of difference equations,
x(t + 1) = x(t ) + αx x(t )[a − cx − 2bx(t ) − bye (t + 1)], y(t + 1) = y(t ) + αy y(t )[a − cy − bxe (t + 1) − 2by(t )], xe (t
+ 1) and + 1) are expected output levels at period t + 1. The positive stationary outputs, where x(t + 1) and y∗ = y(t ) = ye (t ) = y(t + 1), are given by
x∗ =
ye (t
a − 2cx + cy a − 2cy + cx and y∗ = 3b 3b
(4) x∗
=
xe (t )
= x(t ) = (5)
1 Since Cobweb dynamics examined by [9], the discrete-time economic models have a long history in establishing cyclic behavior of the price and the commodity. 2 Surprisingly, it was mentioned in the postscript of the paper in which his theory was completed in 1924.
A. Matsumoto, F. Szidarovszky / Applied Mathematics and Computation 269 (2015) 699–713
701
with the assumption
a > max[2cx − cy , 2cy − cx ]
(6)
to ensure the positivity of the outputs. In the literature, traditional expectation formations such as naive expectation and adaptive expectation are adopted in dynamic games with homogenous firms. [8] constructs a dynamic game with heterogeneous expectation formation, one firm makes output decision based on delayed information on the competitor’s output while the other firm makes its output decision on current information. 3 In particular, it is assumed that
ye (t + 1) = ηy(t ) + (1 − η)y(t − 1) and xe (t + 1) = x(t )
(7)
with the weight η being positive and less than unity. It is demonstrated that dynamic system (4) with (7) enlarges the stability region and generates complex dynamics via period doubling cascade when stability is lost. 3. Continuous-time duopoly models We modify three issues of the discrete-time model (4) in order to consider “delay dynamics” in a continuous-time framework. The first one is to replace z(t + 1) − z(t ) with z˙ (t ), the second is to replace ye (t + 1) with ye (t) and the last one is to substitute one unit discrete time delay with τ continuous time delay. Then the discrete-time model can be converted to a continuous-time model,
x˙ (t ) = αx x(t )[a − cx − 2bx(t ) − bye (t )], y˙ (t ) = αy y(t )[a − cy − bx(t ) − 2by(t )],
(8)
where ye (t) is the expected output formed at time t. We introduce the following four different expectation formations and then consider how the different formations affect dynamics in the continuous-time framework. Forming the expectation on the competitor’s output, (E1) uses the realized output at time t − τ , which is similar to the naive expectation in a discrete-time model, (E2) uses the weighted average of two past outputs at times t − τ1 and t − τ2 , (E3) is an extension of (E2) and employs three past outputs at times t − τ1 , t − τ2 and t − τ3 to obtain the weighted average and finally (E4) generalizes the weighted average using all past outputs from time 0 to t and the weight is exponentially declining with the most weight given to the most current output:
(E1) ye (t ) = y(t − τ ); (E2) ye (t ) =
2
2
ηk y(t − τk ),
k=1
(E3) ye (t ) =
3
3
ηk y(t − τk ),
k=1
(E4) ye (t ) =
0
t
ηk = 1;
k=1
ηk = 1;
k=0
1 −τ e T y(t − τ )dτ = T
0
t
(9)
1 − t−s e T y(s)ds. T
3.1. One fixed delay The point (x∗ , y∗ ) is also the stationary point of the continuous-time model. Assuming formation (E1) in (9) and introducing new notation, xδ (t ) = x(t ) − x∗ and yδ (t ) = y(t ) − y∗ , we obtain a two dimensional system of linearized equations
x˙ δ (t ) = α [−2bxδ (t ) − byδ (t − τ )], y˙ δ (t ) = β [−bxδ (t ) − 2byδ (t )]
(10)
α = αx x∗ and β = αy y∗ .
(11)
with
Substituting exponential solutions
xδ (t ) = eλt u and yδ (t ) = eλt v
(12)
into the linearized system (10) and arranging terms yield an alternative linear system,
λ + 2bα bβ
bα e−λτ λ + 2bβ
u
v
=
0 . 0
(13)
Excluding the trivial solutions (i.e., u = v = 0), we obtain non-trivial solutions by solving the characteristic equation
(λ + 2bα)(λ + 2bβ) − b2 αβ e−λτ = 0. 3
[23] and [13] assume homogenous expectation formations in which both firms use past production data to determine their current outputs.
(14)
702
A. Matsumoto, F. Szidarovszky / Applied Mathematics and Computation 269 (2015) 699–713
For τ = 0, the characteristic equation is reduced to
λ2 + 2b(α + β)λ + 3αβ b2 = 0.
(15) 3αβ b2
> 0, the non-delay characteristic equation has roots with negative real parts, Since, by assumption, 2b(α + β) > 0 and implying that the stationary point is locally stable if there is no delay. It is true that the stationary state preserves stability as far as the delay is positive but sufficiently small. We are concerned with a threshold value (if exists) of the delay for which the stationary point is just destabilized. It is well known that if the stability of the stationary point switches at τ = τ¯ , then the characteristic equation must have a pair of pure conjugate imaginary roots there. To examine the stability switches of the dynamic system (10), we determine this threshold value τ¯ . For this purpose, we substitute a purely imaginary solution λ = iω with ω > 0 into Eq. (14) and separate the real and imaginary parts of the resulted equation,
−ω2 + 4b2 αβ = b2 αβ cos τ ω, −2b(α + β)ω = b2 αβ sin τ ω.
(16)
Squaring both sides of these equations and adding the resulted expressions yield, after arranging terms, a quartic equation in ω,
ω4 + 4b2 (α 2 + β 2 )ω2 + 15(b2 αβ)2 = 0.
(17)
Since the left hand side is positive for all real values of ω, stability switch cannot occur and therefore the steady state is always stable for any τ > 0. In other words, time delay is harmless. Theorem 1. Under the expectation formation (E1), the continuous-time dynamic system (8) is locally asymptotically stable for any value of τ > 0. 3.2. Two fixed delays Assuming (E2) as the output expectation formation and substituting the exponential forms, x(t ) = eλt u and y(t ) = eλt v into the linearized Eq. (10) yields
λ + 2bα bβ
bα(η1 e−λτ1 + η2 e−λτ2 ) λ + 2bβ
u
v
0 . 0
=
(18)
Assuming nontrivial solution gives the corresponding characteristic equation,
(λ + 2bα)(λ + 2bβ) − b2 αβ(η1 e−λτ1 + η2 e−λτ2 ) = 0.
(19)
This is equivalent to equation
a(λ, τ1 , τ2 ) = 1 + a1 (λ)e−λτ1 + a2 (λ)e−λτ2 = 0
(20)
−b2 αβ η for k = 1, 2. (λ + 2bα)(λ + 2bβ) k
(21)
where
ak (λ) =
Suppose λ = iω with ω > 0, then
ak (iω) =
−αβ b2 [4b2 αβ − ω2 − i2b(α + β)ω] ηk for k = 1, 2 [4b2 αβ − ω2 ]2 + [2b(α + β)ω]2
(22)
and their absolute values are
|ak (iω)| =
b2 αβ [4b2 αβ − ω2 ]2 + [2b(α + β)ω]2
ηk for k = 1, 2.
(23)
[11] geometrically investigates the two delay equation based on the idea that solving a(iω, τ1 , τ2 ) = 0 is equivalent to forming a triangle by the three terms in (20) as three vectors in the complex plane placed head to tail. Since the triangle consists of three line segments, it is a necessary and sufficient condition for the existence of a positive solution ω > 0 of Eq. (20) that the sum of the lengths of any two adjacent line segments is not shorter than the length of the remaining line segment. Theorem 2. A purely imaginary root λ = iω with ω > 0 is a solution of a(iω, τ1 , τ2 ) = 0 if and only if the following three inequalities hold:
|a1 (iω)| + |a2 (iω)| ≥ 1, |a1 (iω)| + 1 ≥ |a2 (iω)|, |a2 (iω)| + 1 ≥ |a1 (iω)|.
(24)
Proof. Let a, b and c be three line segments. The end points of a and b with a common starting point can be connected with a segment of length c if and only if
|a − b| ≤ c ≤ a + b.
(25)
A. Matsumoto, F. Szidarovszky / Applied Mathematics and Computation 269 (2015) 699–713
703
The second inequality gives one condition. The first inequality can be rewritten as
−c ≤ a − b ≤ c,
(26)
which means that
b≤a+c
(27)
a ≤ b + c.
(28)
and
These inequality conditions complete the proof. We check whether the first condition of Theorem 2 holds. Substituting (23) into it presents the inequality condition in parametric terms,
ω4 + 4b2 (α 2 + β 2 )ω2 + 15(b2 αβ)2 ≤ 0.
(29)
Notice that the left hand side is positive with all real values of ω implying that stability switch cannot occur for τ 1 > 0 and τ 2 > 0. Theorem 3. Under the expectation formation (E2), the continuous-time dynamic system (8) is locally asymptotically stable for any values of τ 1 > 0 and τ 2 > 0. 3.3. Three fixed delays In the similar way, we can show that the expectation formation (E3) cannot destabilize the stationary state as well. As before, substituting the exponential forms, x(t ) = eλt u and y(t ) = eλt v into the linearized system under (E3) yields
λ + 2bα bβ
bα
−λτ u 0 η1 e 1 + η2 e−λτ2 + η3 e−λτ3 = . v 0 λ + 2bβ
(30)
Nontrivial solution exists if and only if
(λ + 2bα)(λ + 2bβ) − b2 αβ(η1 e−λτ1 + η2 e−λτ2 + η3 e−λτ3 ) = 0.
(31)
As before, this equation can be rewritten as
a(λ, τ1 , τ2 , τ3 ) = 1 + a1 (λ)e−λτ1 + a2 (λ)e−λτ2 + a3 (λ)e−λτ3 = 0
(32)
where
ak (λ) =
−b2 αβ η for k = 1, 2, 3, (λ + 2bα)(λ + 2bβ) k
(33)
Suppose λ = iω with ω > 0, then
ak (iω) =
−b2 αβ [4b2 αβ − ω2 − i2b(α + β)ω] ηk for k = 1, 2, 3, [4b2 αβ − ω2 ]2 + [2b(α + β)ω]2
(34)
and their absolute values are
|ak (iω)| =
b2 αβ [4b2 αβ − ω2 ]2 + [2b(α + β)ω]2
ηk for k = 1, 2, 3.
(35)
[1] gives a straightforward extension of the two delay case considered in [11] to the more general case of three delays and obtains the following result: Theorem 4. Pure imaginary root λ = iω with ω > 0 is a solution of a(iω, τ1 , τ2 , τ3 ) = 0 if and only if the following four inequalities hold:
|a1 (iω)| + |a2 (iω)| + |a3 (iω)| ≥ 1, |a2 (iω)| + |a3 (iω)| + 1 ≥ |a1 (iω)|, |a1 (iω)| + |a3 (iω)| + 1 ≥ |a2 (iω)|, |a1 (iω)| + |a2 (iω)| + 1 ≥ |a3 (iω)|.
(36)
We again check the first condition, which, after substituting (35), can be written as
ω4 + 4b2 (α 2 + β 2 )2 ω2 + 15(b2 αβ)2 ≤ 0.
(37)
Similarly to the two delay case, the left hand side is positive with any real value of ω, implying that stability switch cannot occur for any positive delays. Theorem 5. Under the expectation formation (E3), the continuous-time dynamic system (8) is locally asymptotically stable for any values of τ 1 > 0, τ 2 > 0 and τ 3 > 0.
704
A. Matsumoto, F. Szidarovszky / Applied Mathematics and Computation 269 (2015) 699–713
Fig. 1. Stable trajectories converging to x∗ .
3.4. Continuously distributed delays The expectation formation (E4) assumes that the expected output is a weighted average of all past outputs from time 0 to time t. It is a generalization of (E2) and (E3) with infinitely many continuously distributed delays. Adding the time derivative of the expectation (E4) to linearized system (10) yields a three dimensional system of differential equations,
x˙ (t ) = α [−2bx(t ) − bye (t )], y˙ (t ) = β [−bx(t ) − 2by(t )], 1 y˙ e (t ) = [y(t ) − ye (t )]. T
(38)
The corresponding characteristic equation is given by
λ3 + a1 λ2 + a2 λ + a3 = 0
(39)
where
a1 =
1 + 2(α + β)bT > 0, T
(40)
a2 =
2(α + β)b + 4αβ b2 T >0 T
(41)
a3 =
3αβ b2 > 0. T
(42)
and
A set of necessary and sufficient conditions for all roots of the cubic characteristic equation to have negative real parts is
ai > 0 for i = 1, 2, 3 and a1 a2 − a3 > 0,
(43)
which is a special case of the Routh–Hurwitz criterion. Since all coefficients are positive, we need to check whether the last condition is satisfied or not. Clearly
a1 a2 − a3 =
b [8b2 αβ(α + β)T 2 + b(4α 2 + 9αβ + 4β 2 )T + 2(α + β)] T2
(44)
is always positive for any T > 0. Hence, the stationary state of the three dimensional system is locally stable. Theorem 6. Under the expectation formation (E4), the continuous-time dynamic system (8) is locally asymptotically stable for any values of T > 0. Analytical results obtained so far are numerically confirmed in Fig. 1 in which time trajectories of x(t) in various cases are plotted in different colors.4 Constant initial functions are selected, ϕx (t ) = 5 and ϕy (t ) = 1 for t ≤ 0, the parameters are fixed at αx = αy = 1, a = 4, b = 1, cx = cy = 1, the delays are selected to be τ = 5 in one delay case, τ1 = 5 and τ2 = 6 in two delay case, τ1 = 5, τ2 = 6 and τ3 = 7 in three delay case and T = 5 in the continuously distributed case. It can be seen that all trajectories converge to the steady state of x as analytically shown above. Two issues are also observed: one is that the black real and black dotted trajectories monotonically converge and the latter approaches the steady state more rapidly; the other is that three trajectories in the fixed delay cases start to fluctuate at t = 5 and convergence is oscillatory. 4
“C delay” in the legend of Fig. 1 means continuously distributed delay.
A. Matsumoto, F. Szidarovszky / Applied Mathematics and Computation 269 (2015) 699–713
705
4. Own and competitor’s delays In this section we examine the effects caused by the two delays on dynamics: one delay in the competitor’s output and the other delay in the firm’s own output. Thus the dynamic system is modified in the following way:
x˙ (t ) = αx x(t )[a − cx − 2bx(t − τx ) − by(t − τy )], y˙ (t ) = αy y(t )[a − cy − bx(t ) − 2by(t )]
(45)
where τ x > 0 is the implementation delay and τ y > 0 is the competitor’s information delay. The stationary point is the same as given in (5). The linearized system is
x˙ δ (t ) = α [−2bxδ (t − τx ) − byδ (t − τy )], y˙ δ (t ) = β [−bxδ (t ) − 2byδ (t )].
(46)
Substituting exponential solutions xδ (t ) = eλt u and yδ (t ) = eλt v reduces the linear system to an alternative form,
λ + 2α be−λτx bβ
bα e−λτy λ + 2bβ
u
=
v
0 . 0
(47)
Nontrivial solution exists if and only if
λ2 + 2bβλ + 2bα(λ + 2bβ)e−λτx − b2 αβ e−λτy = 0.
(48)
Before turning to a closer examination of Eq. (48), a few remarks should be made concerning the effect caused by the implementation delay. Assuming that τy = 0 and λ = iω with ω > 0, we can separate Eq. (48) into real and imaginary parts,
2bα(2bβ cos τx ω + ω sin τx ω) = ω2 + b2 αβ , (49)
2bα(ω cos τx ω − 2bβ sin τx ω) = −2bβω. Adding the squares of these equations yields a fourth-order polynomial equation in ω,
ω4 + 2b2 [αβ + 2(β 2 − α 2 )]ω2 − 15b4 α 2 β 2 = 0 that can be solved for
ω2
(50)
and the positive solution is
ω+2 = b2 {−[αβ + 2(β 2 − α 2 )] +
4β 2 + 4αβ 3 + 7α 2 β 2 + (2α 2 − αβ)2 } > 0
(51)
Solving (49) for cos τ x ω and sin τ x ω , substituting ω+ into cos τ x ω and then solving it for τ give the threshold value of the delay,
τxn =
1
cos−1
ω+
b2 β 2 2 ω+ + 4b2 β 2
+ 2nπ , for n = 0, 1, 2, . . .
(52)
A routine bifurcation analysis shows that stability is lost at the smallest threshold value,
τx0 =
1
ω+
cos−1
b2 β 2 2 ω+ + 4b2 β 2
.
(53)
and stability cannot be regained later. In the bifurcation analysis we consider the eigenvalues as functions of the bifurcation parameter; λ = λ(τx ), which is substituted into the characteristic Eq. (48). By taking τy = 0 and implicitly differentiating the resulted equation with respect to the bifurcation parameter, we can have,
{2λ + 2bβ + 2bα [1 − (λ + 2bβ)τx ]e−λτx }
dλ − 2bα(λ + 2bβ)λe−λτx = 0 d τx
(54)
which is solved for (dλ/dτx )−1 for analytical simplicity,
dλ d τx
−1
=
2λ + 2bβ + 2bα [1 − (λ + 2bβ)τx ]e−λτx 2bα(λ + 2bβ)λe−λτx
(55)
Since from(48) with τy = 0,
λ2 + 2bβλ − b2 αβ = −2bα(λ + 2bβ)e−λτx ,
(56)
the right hand side of (55) can be written as
dλ d τx
−1
=
2(λ + bβ) + −λ(λ2 + 2bβλ − b2 αβ)
1
λ(λ + 2bβ)
−
τx . λ
(57)
706
A. Matsumoto, F. Szidarovszky / Applied Mathematics and Computation 269 (2015) 699–713
Its real part at λ = iω is, after arranging terms,
Re
dλ d τx
−1
= Re λ=iω
2(iω + bβ) −iω( − ω2 − b2 αβ + i2bβω)
+ Re
1 iω(iω + 2bβ)
1 ω2 + (2bβ)2 + b2 αβ − (2bβω)2 + (ω2 + b2 αβ)2 ω2 + (2bβ)2 b2 β [(4β − α)(ω2 + b2 αβ) + (4bβ)2 ] = [(2bβω)2 + (ω2 + b2 αβ)2 ][ω2 + (2bβ)2 ] =
(58)
The last expression is positive if 4β > α or the adjustment coefficients do not differ so much, which is assumed to be henceforth. Thus we obtain the derivatives λ (τxn ) at the critical values τxn . Since these derivatives are positive with 4β > α , at all critical points at least one eigenvalue pair changes the sign of its real part from negative to positive. So stability is lost at the smallest critical point τx0 and the stability cannot be regained later. Theorem 7. If 4β > α holds, then the stationary point of (45) with τy = 0 is locally asymptotically stable for τx < τx0 , loses stability at τx = τx0 and bifurcates to a limit cycle via Hopf bifurcation for τx > τx0 . We now investigate the characteristic equation with two positive delays by applying Gu’s method [11]. Dividing both sides by
λ2 + 2bβλ and introducing new functions a1 (λ) =
2bα
λ
and a2 (λ) =
−b2 αβ λ(λ + 2bβ)
(59)
simplify Eq. (48) as
a(λ, τx , τy ) = 1 + a1 (λ)e−λτx + a2 (λ)e−λτy = 0.
(60)
Suppose that λ = iω with ω > 0. Substituting it into ai (λ) yields
a1 (iω) = −i
2bα
(61)
ω
and
a2 (iω) =
b2 αβω 2b3 αβ 2 + i . ω(ω2 + 4b2 β 2 ) ω(ω2 + 4b2 β 2 )
(62)
Their absolute values are
|a1 (iω)| =
2bα
(63)
ω
and
|a2 (iω)| =
b2 αβ
ω ω2 + 4b2 β 2
.
(64)
We check whether the three conditions of Theorem 2 are satisfied. Substituting the absolute values reduces the three inequalities to the following two conditions,
f (ω) = ω4 − 4bαω3 + 4b2 (α 2 + β 2 )ω2 − 16b3 αβ 2 ω + 15b4 α 2 β 2 ≤ 0
(65)
g(ω) = ω4 + 4bαω3 + 4b2 (α 2 + β 2 )ω2 + 16b3 αβ 2 ω + 15b4 α 2 β 2 ≥ 0.
(66)
and
Notice that g(ω) ≥ 0 is always true for ω ≥ 0. However, it is ambiguous whether f(ω) ≤ 0 holds or not. Using the variable transformation
4bα ω =x− − = x + bα ,
(67)
4
we rewrite f(ω) as
F (x) = x4 − 2b2 (α 2 − 2β 2 )x2 − 8b3 αβ 2 x + b4 α 2 (α 2 + 3β 2 ).
(68)
α2
> (+, +, −, +) if < and (+, −, +) if The sequence of the signs of the polynomial coefficients are (+, −, −, +) if α 2 = 2β 2 . In any sequence, the number of sign changes is two. According to Descartes’ rule of sign, this polynomial has either two or zero positive roots. It is clear that F(0) > 0 and F (∞) = ∞. Substituting x = α b yields
F (α b) = α 2 b4 ( − β 2 ) < 0,
2β 2 ,
α2
2β 2 ,
(69)
A. Matsumoto, F. Szidarovszky / Applied Mathematics and Computation 269 (2015) 699–713
707
which implies that F (x) = 0 has two positive solutions, x1 and x2 where 0 < x1 < α b and x2 > α b. Hence f (ω) = 0 has two positive solutions
ω1 = x1 + α b > 0 and ω2 = x2 + α b > 0
(70)
and f(ω) ≤ 0 for ω ∈ [ω1 , ω2 ], and a routine investigation can show that f(ω) > 0 otherwise.5 We will next find all pairs of (τ 1 , τ 2 ) satisfying a(iω, τ1 , τ2 ) = 0. The three terms in (60) are three vectors in the complex plane that construct a triangle. Let us suppose that 1 is its base and let us denote the angle between 1 and |a1 (iω)e−iωτx | by θ 1 and an angle between 1 and |a2 (iω)e−iωτy | by θ 2 . Since, |e−iωτx | = |e−iωτx | = 1, we have, by the law of cosine,
θ1 (ω) = cos−1 = cos−1 and
θ2 (ω) = cos
−1
= cos−1
2
2 2 · 1 · a1 (iω)
12 + a1 (iω) − a2 (iω)
ω4 + 4b2 (α 2 + β 2 )ω2 + 15b4 α 2 β 2 4bαω(ω2 + 4b2 β 2 )
(71)
2
2 2 · 1 · a2 (iω)
12 + a2 (iω) − a1 (iω)
ω4 − 4b2 (α 2 − β 2 )ω2 − 15b4 α 2 β 2 . 2b2 αω ω2 + 4b2 β 2
(72)
Since the triangle may be located above and also under the horizontal axis, we get two equations as
{arg[a1 (iω)e−iτx ω ] + 2mπ } ± θ1 (ω) = π {arg[a2 (iω)e−iτy ω ] + 2nπ } ∓ θ2 (ω) = π which yield the threshold values of the delays:
τx± (ω, n) =
τ (ω, n) =
1 3π + (2m − 1)π ± θ1 (ω) ω 2
and ∓ y
(73)
1
ω
tan
−1
2bβ
ω
+ (2n − 1)π ∓ θ2 (ω) .
(74)
(75)
Let be the set of all ω values for which f(ω) ≤ 0 holds. Then we can find all pairs of (τ 1 , τ 2 ) constructing the stability switching curves for ω ∈ which consists of two sets of parametric segments,
L1 (m, n) = {τx+ (ω, m), τy− (ω, n)}
(76)
L2 (m, n) = {τx− (ω, m), τy+ (ω, n)}.
(77)
and
To illustrate the stability switching curves we first specify the parameter values that make the mathematical derivations simple without changing the main conclusions of the analysis: Assumption: α = β = b = 1. Under this assumption, we have
f (ω) = ω4 − 4ω3 + 8ω2 − 16ω + 15.
(78)
It is not difficult to show that f (ω) = 0 has two real and positive roots,
1
ωs = 1 +
2 3/K 1
ωe = 1 +
2 3/K
with
5
−
1 2
+
1 2
− 13 (12 + K ) + 16 3/K 1.611
− 13
(79)
(12 + K ) + 16 3/K 2.326
√ √ 1 1 K = −4 + (584 − 48 87) 3 + 2(73 + 6 87) 3 .
A formal proof is given in the Appendix.
(80)
708
A. Matsumoto, F. Szidarovszky / Applied Mathematics and Computation 269 (2015) 699–713
The domain of ω is the interval [ωs , ωe ]( = ). So the segment (76) starts at point (τx+ (ωs , m), τy− (ωs , n)) and terminates at point (τx+ (ωe , m), τy− (ωe , n)) and so does the segment (77) at points (τx− (ωs , m), τy+ (ωs , n)) and (τx− (ωe , m), τy+ (ωe , n)). Furthermore the angles of the triangle are
ω4 + 8ω2 + 15 θ1 (ω) = cos 4ω(ω2 + 4) ω4 − 15 . θ2 (ω) = cos−1 √ 2ω ω 2 + 4
−1
(81)
Since
ωs4 + 8ωs2 + 15 = 1, 4ωs (ωs2 + 4)
ωe4 + 8ωe2 + 15 =1 4ωe (ωe2 + 4)
(82)
ωs4 − 15 = −1, 2ωs ωs2 + 4
ωe4 − 15 = 1, 2ωe ωe2 + 4
(83)
and
we have
θ1 (ωs ) = 0, θ1 (ωe ) = 0
(84)
θ2 (ωs ) = π , θ2 (ωe ) = 0.
(85)
and
We can now formulate the following results concerning the locations of the curve segments: Theorem 8. Given m and n, the segments L1 (m, n + 1) and L2 (m, n) have the same starting point whereas the segments L1 (m, n) and L2 (m, n) have the same end point. Proof. For analytical simplicity we assume m = 0. The starting points of L1 (0, n) and L2 (0, n) are
Ls1 (0, n) = {τx+ (ωs , 0), τy− (ωs , n)}
(86)
with
τx+ (ωs , 0) =
1
π
ωs 2
and τy− (ωs , n) =
1
tan−1
ωs
2 + 2(n − 1)π , ωs
(87)
and
Ls2 (0, n) = {τx− (ωs , 0), τy+ (ωs , n)}
(88)
where
τx− (ωs , 0) =
1
π
ωs 2
and τy+ (ωs , n) =
1
ωs
tan−1
2 + 2nπ . ωs
(89)
Hence we have
Ls1 (0, n + 1) = Ls2 (0, n).
(90)
In the same way, the end points of L1 (0, n) and L2 (0, n) are
Le1 (0, n) = {τx+ (ωe , 0), τy− (ωe , n)}
(91)
with
τx+ (ωe , 0) =
1
π
ωe 2
and τy− (ωe , n) =
1
ωe
tan−1
2 + (2n − 1)π ωe
(92)
and
Le2 (0, n) = {τx− (ωe , 0), τy+ (ωe , n)}
(93)
where
τx− (ωe , 0) =
1
π
ωe 2
and τy+ (ωe , n) =
1
ωe
tan−1
2 + (2n − 1)π . ωe
(94)
Hence we have
Le1 (0, n) = Le2 (0, n). The same result is obtained for any integer m > 0.
(95)
A. Matsumoto, F. Szidarovszky / Applied Mathematics and Computation 269 (2015) 699–713
709
Fig. 2. Stability switching curve.
This result is numerically confirmed in Fig. 2 where m = 0. The stability switching curve consists of the red and blue segments that correspond to L1 (0, n) and L2 (0, n) for n = 0, 1, 2. The red and blue segments shift upward when n increases and to the right when m increases. The upward-sloping blue segment crossing the horizontal axis at τx = τx0 is the L2 (0, 0) segment.6 It connects to the red L1 (0, 1) segment at a point where Ls1 (0, 1) = Ls2 (0, 0) holds. As before, the upper script “s” means the starting (i.e., initial) point of the segment and
Ls1 (0, 1) = (τx+ (ωs , 1), τy− (ωs , 1)) and Ls2 (0, 0) = (τx− (ωs , 0), τy+ (ωs , 0)). This L1 (0, 1) segment then connects to the blue L2 (0, 1) segment at a point where “e” means again the end (i.e., terminal) point of the segment and
Le1 (0, 1) = (τx+ (ωe , 1), τy− (ωe , 1)) and Le2 (0, 1) = (τx− (ωe , 1), τy+ (ωe , 1)).
(96) Le1 (0, 1)
=
Le2 (0, 1)
holds. The upper script (97)
As n increases, the two segments are connected in the same way to construct the continuous stability switching curve. The eigenvalues are purely imaginary on this curve. The stability switching curves L1 (0, n) and L2 (0, n) divide the first quadrant of the (τ 1 , τ 2 ) plane into two parts. One contains the origin and its every point can be reached from the origin via continuous curve not crossing the stability switching curve. At the points in this region the real parts of the eigenvalues are negative, so the system is locally asymptotically stable. The complement of this region except the curves gives points when the system is unstable. We again notice that the blue L2 (0, 0) segment crosses the horizontal axis at τx = τx0 ( 0.824). This indicates that stability is preserved for τx < τx0 and lost for τx > τx0 when τy = 0. In other words, this is the threshold value of τ x obtained in the case of one delay (Theorem 7), which is equal to 0.824 under Assumption. So we could confirm this threshold value in two different ways. The direction of change (that is, stability loss or regain) can be checked based on the Hopf bifurcation theory and at each point (τ x , τ y ) the number of stability losses or stability regains determine if this point is stable or not as in [18]. Instead of following this lengthy and complicated method, we perform numerical simulations to investigate how stability changes in case of two positive delays. We now examine the effects caused by changing the value of τ x , keeping the value of τ y at some positive fixed value. In Fig. 3(A) we increase the value of τ x from 0 to 1.3 in 0.001 increments along the horizontal dotted line with fixed value τy = τya ( = 2). For each value of τ x , the dynamic system runs for 0 ≤ t ≤ 1000 and only the data for 950 ≤ t ≤ 1000 are used to get rid of the transients. The local maxima and minima are plotted against each value of τ x . If the bifurcation diagram has one point against τ x , then the maximum point is identical with the minimum point, implying that the stationary state is locally stable. If it has two points, then a limit cycle having one maximum and one minimum emerges and if many points, then a limit cycle exhibits many ups and downs. The horizontal dotted line crosses the blue L2 (0, 1) segment at τxa 0.652 denoted by the green dot in Fig. 2. It is seen that the stationary state is locally asymptotically stable for τx < τxa and becomes unstable for τx > τxa as shown in Fig. 3(A). When stability is lost at τx = τxa , then the stationary point bifurcates to a limit cycle and does not regain stability for larger values of τ x . We move to the second simulation in which the fixed value of τ y is increased to τyb = 4.4 and the same procedure is repeated. As is seen in Fig. 2, the horizontal dotted line at τyb crosses the stability switching curve three times as denoted by three green dots. As described in Fig. 3(B), both stability regain and stability loss occur in this example. In particular, stability is first lost at the first crossing point, the left most green dot denoted by τx1 ( 0.699). Then limit cycle emerges for τ x in the interval [τx1 , τx2 ] where τx2 ( 0.892) corresponds to the middle green point. If the τ x -value of the most right green point is denoted by τx3 ( 0.964), then stability is regained for τx ∈ [τx2 , τx3 ]. The stationary state again loses its stability for τx > τx3 . The bifurcation diagram shows that a trajectory gradually exhibits more complex dynamics through a quasi period-doubling process as τ x become larger and then returns to the simple limit cycle through a quasi period-halving process as τ x further increases. 6
The L1 (0, 0) segment is located in the fourth quadrant so it is not illustrated.
710
A. Matsumoto, F. Szidarovszky / Applied Mathematics and Computation 269 (2015) 699–713
Fig. 3. Bifurcation diagram with respect to τ x , given τ y .
Fig. 4. Bifurcation diagram with respect to τ y , given τ x .
We now turn attention to the delay effect caused by changing values of τ y , keeping the value of τ x at some positive value. We also perform two simulations. In Fig. 2, τ x is fixed at τxA = 0.75 and τ y increases along the vertical dotted line crossing the stability switching curve three times denoted by three black dotes. Although they are not labelled on the vertical axis to avoid confusion in the figure, we denote their τ y -values by τy1 ( 1.410), τy2 ( 3.138) and τy3 ( 4.215) with ascending order as shown in Fig. 4(A). The bifurcation diagram indicates that stability is lost at τy = τy1 , a limit cycle emerges for larger values than τy1 and then regains stability at τy = τy2 . Stability is preserved for τy < τy3 and is lost again for τy > τy3 . The diagram implies a similar dynamic cascade in which the process of stability loss, birth of limit cycle and regain of stability is repeated as the value of τ y increases. In Fig. 4(B) the fixed value of τ x is increased to τxB = 1.2 from τxA = 0.75 and the value of τ y increases along the vertical dotted line at τxB = 1.2 which is located to the right of the stability switching curve in Fig. 2. Even for τy = 0, the dynamic system is unstable as τxB > τx0 and generates a limit cycle with one maximum and one minimum. As the value of τ y increases, the bifurcation diagram implies that the diameter of the limit cycle varies and then more complex dynamics appears through a period-doubling like process. 5. Concluding remarks Delay dynamics of heterogeneous duopolies in continuous-time framework was investigated. This is an extension of the earlier study [8] on heterogeneous duopolies in a discrete-time framework where the authors investigated a heterogenous duopoly market in which one firm forms its expectation on the competitor’s output only with delayed information and the other firm uses the current information. Following their spirit, we assumed that one firm has a delay in obtaining information about the competitor’s output (i.e., information delay) in four different ways in which one, two or three delayed information is used to form expectation on the output of the competitor, or all past output values are used as a continuously distributed delay. Applying a recently developed method to examine stability of delay differential equations, we analytically demonstrated that the information delay is harmless to stability in the continuous-time duopoly model. In other words, no stability switch occurs regardless
A. Matsumoto, F. Szidarovszky / Applied Mathematics and Computation 269 (2015) 699–713
711
of the expectation formation and the lengths of the delays, meaning that the system is stable with any values of the delays, so delays are harmless in all examined cases. This shows a sharp difference from the main result of the discrete case that a delay in the competitor’s production can have a stabilizing effect. Then we assumed that one firm had implementation and information delays and the other firm could make decisions without any delay. We first constructed the stability switching curve which consists of a sequence of connected segments forming a continuous curve as function of the two delays. We confirmed that stability can be switched to instability when the length of the implementation delay takes over some threshold value which depends on the length of the information delay. It is numerically verified that various dynamics ranging from simple to complex can emerge according to the values of the delays. If one firm had only implementation delay, then stability was lost when the delay exceeded a certain threshold value. Therefore it was demonstrated that the implementation delay matters. We also examined the dynamics of the system with changing the lengths of the delays. With fixed small information delay, the system lost stability at a threshold of the implementation delay and the stability could not be regained later. However with larger fixed information delay the stability is lost at a threshold of the implementation delay and the stability was gained for a while with the further increase of the implementation delay. Similar situation occurs with small fixed implementation delay, and if this delay is larger, then the system is unstable for all lengths of the information delay. Hopf bifurcation occurs at the critical values giving the possibility of the birth of limit cycles. There was a variety of dynamic behavior in the case of instability including limit cycles and period-doubling like processes in losing stability and stability can be regained by a period-halving process. The methodology used in this paper to determine stability switching curves and the stability regions can be used in many other delay models arising in several fields of mathematical economics. There are many directions for future research. In this study, we assumed heterogenous duopolies in which only one firm had delayed information formulated in (45). One possible extension is given below in which duopolies are heterogeneous and one firm has an implementation delay and the other has an information delay,
x˙ (t ) = αx x(t )[a − cx − 2bx(t − τx ) − by(t )], y˙ (t ) = αy y(t )[a − cy − bx(t − τy ) − 2by(t )].
(98)
However, the characteristic equation derived from this dynamic system is identical with Eq. (48). It may be interesting to consider duopoly models with different economic backgrounds generating exactly the same dynamics. The second extension concerns the way of the transformation to a continuous-time model from a discrete-time model. We make two assumptions: one is replacing the discrete time difference z(t + 1) − z(t ) with time derivative z˙ (t ) and the other is introducing continuous-time expectation formations concerning ye (t) instead of the discrete-time formation. The delay discrete-time formation was given as
ye (t + 1) = ηy(t ) + (1 − η)y(t − 1)
(99)
which can be written as
ye (t + 1) = =
η[y(t ) − y(t − 1)] + y(t − 1) ηy˙ (t ) + y(t − 1).
(100)
If the unit time difference is replace with delay τ , then we have the following form of a continuous-time system:
x˙ (t ) = αx x(t )[a − cx − 2bx(t ) − b{ηy˙ (t ) + y(t − τ )}], y˙ (t ) = αy y(t )[a − cy − bx(t ) − 2by(t )].
(101)
This can be solved for x˙ (t ) and y˙ (t ) to derive a dynamic system of explicit delay differential equations. Further it is possible to introduce an implementation delay on the firm’s own output although analysis will be more complicated. The third extension is to adopt the modelling method proposed in [3] and applied in [20] when considering delay monopoly dynamics. The continuoustime system is described by the following delay equations,
σx x˙ (t ) + x(t ) = x(t − τx ) + αx x(t − τx )[a − cx − 2bx(t − τx ) − bye (t )], σy y˙ (t ) + y(t ) = y(t − τx ) + βy y(t − τy )[a − cy − bxe (t ) − 2by(t − τy )],
(102)
in which σ x ≥ 0 and σ y ≥ 0 denote the inertias (or frictions) in the production process. It should be noticed that this system can be reduced to system (4) if σx = σy = 0 and τx = τy = 1. So when σ x and σ y take smaller values, the generated dynamics is similar to the discrete-time case (4). On the other hand, when σ x and σ y take larger values, it may generate different dynamics. Appendix In this Appendix, we will prove that equation
x4 − 2b2 (α 2 − 2β 2 )x2 − 8b3 αβ 2 x + b4 α 2 (α 2 + 3β 2 ) = 0
(103)
has two distinct positive roots and no negative root exists. Let x = α bz, then
α 4 b4 z4 − 2b2 (α 2 − 2β 2 )α 2 b2 z2 − 8b3 αβ 2 α bz + b4 α 2 (α 2 + 3β 2 ) = 0
(104)
which can be simplified as
h(z) = α 2 z4 − 2(α 2 − 2β 2 )z2 − 8β 2 z + (α 2 + 3β 2 ) = 0.
(105)
712
A. Matsumoto, F. Szidarovszky / Applied Mathematics and Computation 269 (2015) 699–713
Notice that
h(0) = α 2 + 3β 2 > 0,
lim h(z) = ∞ and h(1) = −β 2 < 0.
z→±∞
(106)
The derivative of h(z) has the form
h (z) = 4(z − 1)(α 2 z2 + α 2 z + 2β 2 ).
(107)
The sign of this derivative depends on the sign of the quadratic function
g(z) = α 2 z2 + α 2 z + 2β 2
(108)
with the discriminant
D = α 2 (α 2 − 8β 2 ).
(109)
If α 2 = 8β 2 , then it is positive except its negative vertex, so h (z) < 0 as z < 1 and differs from the vertex, h (z) > 0 as z > 1 and h (z) = 0 as z = 1 or equals the vertex. Therefore h(z) strictly decreases as z < 1 and strictly increases as z > 1. If α 2 < 8β 2 , then h (z) > 0 as z > 1 and h (z) < 0 as z < 1. Since h(0) > 0 and both limits at ± ∞ are positive there are two positive roots of h(z) and no negative root exists. Assume next that α 2 > 8β 2 . Then the quadratic function g(z) has two negative roots:
z1,2 = −
1 ± 2
α 2 − 8β 2 (z1 < z2 ). 2α
(110)
,
(111)
So h (z) < 0 if z < z1 or z2 < z < 1, and h (z) > 0 if z1 < z < z2 or z > 1. At these roots
z2 = −z − so
2β 2
α2
h(z) = =
α
2
−z −
2β 2
2
α2
(α 2 − 8β 2 )z +
2β 2 − 2(α − 2β ) −z − 2 − 8β 2 z + (α 2 + 3β 2 ) α 2
1
α2
2
(α 4 + 5α 2 β 2 − 4β 4 ).
(112)
Since z1 and z2 are negative and their average value is −1/2, both are larger than −1. Since α 2 > 8β 2 , at the roots
(α 2 − 8β 2 )z > 8β 2 − α 2 .
(113)
Hence
h(z) > 8β 2 − α 2 +
1
α2
(α 4 + 5α 2 β 2 − 4β 4 )
β2 (5α 2 − 4β 2 ) α2 β2 > 8β 2 + 2 (40β 2 − 4β 2 ) > 0. α = 8β 2 +
(114)
So h(z) decreases from the ∞ limit as z → −∞ to h(z1 ) > 0, then increases to h(z2 ) > 0 and decreases again until h(1) < 0 and then increases and tends to ∞ as z → ∞ . Since h(0) > 0, there are no negative roots, only two distinct positive roots. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17]
E. Almoderesi, M. Bozorg, Stability crossing surfaces for linear time delay systems with three delays, Int. J. Control 82 (2009) 2304–2310. R. Bellman, K. Cooke, Differential-Difference Equations, Academic Press, New York, 1956. M. Berezowski, Effect of time delay on the generation of chaos in continuous systems, Chaos, Solitons Fractals 12 (2001) 83–89. G.-I. Bischi, A. Namzada, Global analysis of a dynamic duopoly game with bounded rationality, in: Advances in Dynamic Games and Application, vol. 5, Birkhauser, Boston, Massachussetts, USA, 2000, pp. 361–385. (Chapter 20). G.-I. Bischi, C. Chiarella, M. Kopel, F. Szidarovszky, Nonlinear Oligopolies: Stability and Bifurcations, Springer, Berlin/Heidelberg/New York, 2010. A. Cournot, Recherches sur les Principes Mathématiques de la Théorie des Richessess, Hachette, Paris, 1838.(English translation (1960): Researches into the Mathematical Principles of the Theory of Wealth, Kelley, New York) J. Cushing, Integro-Differential Equations and Delay Models in Population Dynamics, Springer-Verlag, Berlin/Heidelberg/New York, 1977. A. Elsadany, A. Matouk, Dynamic Cournot duopoly game with delay, J. Complex Syst. (2014). http://dx.doi.org/10.1155/2014/384843 M. Ezekiel, The cobweb theorem, Q. J. Econ. 52 (1933) 255–280. R. Goodwin, The nonlinear accelerator and the persistence of business cycles, Econometrica 19 (1951) 1–17. K. Gu, S.-I. Nicolescu, J. Chen, On stability crossing curves for general systems with two delays, J. Math. Anal. Appl. 311 (2005) 231–253. J. Haldane, A contribution to the theory of price fluctuations, Rev. Econ. Stud. 1 (1933) 186–195. S.Z. Hassan, On delayed dynamical duopoly, App. Math. Comput. 151 (2004) 275–286. T. Howroyd, A. Russel, Cournot oligopoly models with time delays, J. Math. Econ. 13 (1984) 97–103. M. Kalecki, A macrodynamic theory of business cycles, Econometrica 3 (1935) 327–344. A. Larson, The hog cycle as harmonic motion, J. Farm Econ. 46 (1964) 375–386. M. Mackey, Commodity price fluctuations: Price dependent delays and nonlinearities as explanatory factors, J. Econ. Theory 48 (1989) 497–509.
A. Matsumoto, F. Szidarovszky / Applied Mathematics and Computation 269 (2015) 699–713
713
[18] A. Matsumoto, F. Szidarovszky, Learning monopolies with delayed feedback on price expectations, Commun. Nonlinear Sci. Numer. Simul. 28 (2015) 151– 165. [19] A. Matsumoto, F. Szidarovszky, Discrete-time delay dynamics of boundedly rational monopoly, Dec. Econ. Finance 37 (2014) 53–70. [20] A. Matsumoto, F. Szidarovszky, On the Comparison of Discrete and Continuous Dynamic Systems, IERCU DP #248, Chuo University (http://www.chuou.ac.jp/research/institutes/economic/publication/discussion/pdf/discussno248.pdf?1429781458990). [21] K. Okuguchi, Expectations and Stability in Oligopoly Models, Springer-Verlag, Berlin/Heidelberg/New York, 1976. [22] K. Okuguchi, F. Szidarovszky, The Theory of Oligopoly with Multi-Product Firms, Springer-Verlag, Berlin/Heidelberg/New York, 1999. [23] M.T. Yassen, H.N. Agiza, Analysis of a duopoly game with delayed bounded rationality, Appl. Math. Comput. 138 (2003) 387–402.
Contents lists available at ScienceDirect
Applied Mathematics and Computation journal homepage: www.elsevier.com/locate/amc
Delay dynamics of a Cournot game with heterogeneous duopolies✩ Akio Matsumoto a,∗, Ferenc Szidarovszky b a 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 b Department of Applied Mathematics, University of Pécs, Ifjúság u. 6, Pécs, H-7624, Hungary
a r t i c l e
i n f o
Keywords: Delay duopoly dynamics Bounded rationality Continuous-time system Stability switch
a b s t r a c t We construct a Cournot duopoly model in a continuous-time framework and consider its dynamic behavior when the firms are heterogeneous in determining their output decision: one firm has an information delay in the competitor’s output as well as an implementation delay in its own output while the other firm does not have any delays. Two main results are obtained. One is that the information delay does not affect dynamics of the commodities and the other is that the implementation delay can destabilize the otherwise stable equilibrium. Furthermore, it is numerically confirmed that the two delay duopoly model can generate a wide spectrum of dynamics ranging from cyclic dynamics to chaos. © 2015 Elsevier Inc. All rights reserved.
1. Introduction Oligopoly theory is one of the most frequently studied research areas in mathematical economics. Based on the pioneering work [6] many scientists devoted their effort to this interesting and important model family. The earlier results up to the mid 70s are summarized in [21] and their multiproduct generalizations are presented with additional applications in [22]. Recently nonlinear models become the main focus, and a summary of more recent developments is presented in [5]. There is a wide variety of oligopoly models: single product models with and without product differentiation, multiproduct models, labor managed oligopolies among others. In most studies instantaneous information is assumed to the firms about the price and the output levels of the competitors, however in real economies there is always a time delay due to information collection, decision making and implementation. Time delays with continuous time scales can be modeled as fixed or continuously distributed delays. Differential-difference equations are used to model fixed delays [2] while continuously distributed delays are described by Volterra-type integro-differential equations [7]. The central issue in examining dynamic models with delays is the way how the delays affect the asymptotic behavior of the systems. If discrete time scales are assumed, then the lengths of the delays are assumed to be integer multiples of the time unit and higher-order difference equations model the dynamism of the systems. This paper presents a new characterization of price formation of a single commodity by which regular and irregular price oscillations can be constructed in a continuous-time framework with fixed time delays. In the existing literature of macroeconomic dynamics, it has been well-known since the 1930s that a delay in investment or production is one of the key factors to ✩ The authors are indebted to two anonymous referees for their comments that are most helpful in preparing the final version. They highly appreciate the financial supports from the MEXT-Supported Program for the Strategic Research Foundation at Private Universities 2013–2017, the Japan Society for the Promotion of Science (Grant-in-Aid for Scientific Research (C) 24530202, 25380238 and 26380316) and Chuo University (Joint Research Grant). The usual disclaimers apply. ∗ Corresponding author. Tel.: +81 42 674 3351; fax: +81 42 674 3425. E-mail addresses: [email protected] (A. Matsumoto), [email protected] (F. Szidarovszky).
http://dx.doi.org/10.1016/j.amc.2015.07.051 0096-3003/© 2015 Elsevier Inc. All rights reserved.
700
A. Matsumoto, F. Szidarovszky / Applied Mathematics and Computation 269 (2015) 699–713
generate cyclic dynamics in national income, capital accumulation, employment rate, etc. See the pioneering papers, [10] and [15]. On the other hand, in the existing literature of microeconomic dynamics such as price and commodity oscillations, it has been less-known that production delays are also able to generate economic fluctuations of a continuous-time unstable economic system.1 As early as in the 1930s, it has been already shown in [12] that cyclic behavior can arise in a simple price dynamics model with production delay, coaxing an idea from theoretical biology.2 A continuous-time cobweb model is constructed in [16] for the pork market in which pork is inelastically supplied with 12 month delay and gives rise to cyclic behavior. [14] is the first work to consider delay dynamics in a N-firm oligopoly model. On the other hand, [17] formulates the price dynamics of a single commodity market as a nonlinear delay differential equation and rigorously derives the stability conditions and the birth of a cyclic oscillation via Hopf bifurcation. Recently [18,19] reconsider delay dynamics in various monopoly models and emphasize that delay may explain various dynamic behavior of micro economic variables. In the existing literature, however, analysis is mostly confined to cases of a single dimensional model with one delay or multiple delay or a multiple dimensional model with one delay. In this paper we construct a multiple dimensional model with multiple delays and examine cyclic behavior in the model. This paper revisits the delay duopoly game investigated in [8] in order to reconsider the delay effect on stability on the Nash equilibrium in a continuous-time framework. The game in [8] is built in a discrete-time framework and its distinguished feature is the heterogeneity of the duopolists. In predicting competitor’s output, one duopolist uses a combination of the current and delayed information while the other adopts naive expectation. The main result is that the delay has a stabilizing effect: it enlarges the stability region of the Nash equilibrium. Generally speaking, a continuous-time model has a larger stability region than a discrete-time model (see, for example, [20]). Furthermore, the delay in a continuous-time model often exhibits stability switching from stability to instability if it becomes larger than some threshold value. However, it is not known yet whether such stability switch occurs in the heterogenous duopoly game. Hence in this paper, we retain the heterogeneity of the discrete-time dynamic system and consider its effect in a delay continuous-time framework . The main results can be summarized as follows: (I) When only one firm has information delays in the competitor’s output, then this heterogeneous expectation formation has no destabilizing effect; (II) When the firm has a delay in implementing information about its own output, the otherwise stable equilibrium can be destabilized and give rise to simple as well as to complex fluctuations when the delay is sufficiently long. This paper is organized as follows. In Section 2, the basic elements of the model of [8] are recapitulated. In Section 3 the corresponding continuous-time duopoly model is constructed and the delay effects due to heterogenous expectation formation are examined. In Section 4, the destabilizing effect caused by the implementation delay is analytically and numerically examined. In the final section some concluding remarks are given and further research directions are outlined. 2. Discrete-time duopoly model A duopoly game with heterogeneous bounded rationality is formulated in discrete-time framework. Two firms produce a homogenous good. Firm 1 produces quantity x with marginal cost cx , while firm 2 produces quantity y at marginal cost cy . The market demand is linear and depends on the total output of the industry,
p = a − b(x + y)
(1)
where a and b are positive constants. The profit of firm k is
πk = [a − b(x + y)]z − cz z.
(2)
where k = 1, 2 and z = x, y. [4] assumes bounded rational firms and introduce gradient dynamics in which the firms adjust production levels according to their marginal profits in such a way that a firm increases output if the marginal profit is positive, decreases it if it is negative and does not change it if it is zero. The gradient dynamics with boundedly rational firm k can be described by
z(t + 1) = z(t ) + αk z(t )
∂πk ∂z
(3)
where α k > 0 is an adjustment coefficient. Hence in the duopoly game, the output adjustment process with gradient method is presented by a two dimensional system of difference equations,
x(t + 1) = x(t ) + αx x(t )[a − cx − 2bx(t ) − bye (t + 1)], y(t + 1) = y(t ) + αy y(t )[a − cy − bxe (t + 1) − 2by(t )], xe (t
+ 1) and + 1) are expected output levels at period t + 1. The positive stationary outputs, where x(t + 1) and y∗ = y(t ) = ye (t ) = y(t + 1), are given by
x∗ =
ye (t
a − 2cx + cy a − 2cy + cx and y∗ = 3b 3b
(4) x∗
=
xe (t )
= x(t ) = (5)
1 Since Cobweb dynamics examined by [9], the discrete-time economic models have a long history in establishing cyclic behavior of the price and the commodity. 2 Surprisingly, it was mentioned in the postscript of the paper in which his theory was completed in 1924.
A. Matsumoto, F. Szidarovszky / Applied Mathematics and Computation 269 (2015) 699–713
701
with the assumption
a > max[2cx − cy , 2cy − cx ]
(6)
to ensure the positivity of the outputs. In the literature, traditional expectation formations such as naive expectation and adaptive expectation are adopted in dynamic games with homogenous firms. [8] constructs a dynamic game with heterogeneous expectation formation, one firm makes output decision based on delayed information on the competitor’s output while the other firm makes its output decision on current information. 3 In particular, it is assumed that
ye (t + 1) = ηy(t ) + (1 − η)y(t − 1) and xe (t + 1) = x(t )
(7)
with the weight η being positive and less than unity. It is demonstrated that dynamic system (4) with (7) enlarges the stability region and generates complex dynamics via period doubling cascade when stability is lost. 3. Continuous-time duopoly models We modify three issues of the discrete-time model (4) in order to consider “delay dynamics” in a continuous-time framework. The first one is to replace z(t + 1) − z(t ) with z˙ (t ), the second is to replace ye (t + 1) with ye (t) and the last one is to substitute one unit discrete time delay with τ continuous time delay. Then the discrete-time model can be converted to a continuous-time model,
x˙ (t ) = αx x(t )[a − cx − 2bx(t ) − bye (t )], y˙ (t ) = αy y(t )[a − cy − bx(t ) − 2by(t )],
(8)
where ye (t) is the expected output formed at time t. We introduce the following four different expectation formations and then consider how the different formations affect dynamics in the continuous-time framework. Forming the expectation on the competitor’s output, (E1) uses the realized output at time t − τ , which is similar to the naive expectation in a discrete-time model, (E2) uses the weighted average of two past outputs at times t − τ1 and t − τ2 , (E3) is an extension of (E2) and employs three past outputs at times t − τ1 , t − τ2 and t − τ3 to obtain the weighted average and finally (E4) generalizes the weighted average using all past outputs from time 0 to t and the weight is exponentially declining with the most weight given to the most current output:
(E1) ye (t ) = y(t − τ ); (E2) ye (t ) =
2
2
ηk y(t − τk ),
k=1
(E3) ye (t ) =
3
3
ηk y(t − τk ),
k=1
(E4) ye (t ) =
0
t
ηk = 1;
k=1
ηk = 1;
k=0
1 −τ e T y(t − τ )dτ = T
0
t
(9)
1 − t−s e T y(s)ds. T
3.1. One fixed delay The point (x∗ , y∗ ) is also the stationary point of the continuous-time model. Assuming formation (E1) in (9) and introducing new notation, xδ (t ) = x(t ) − x∗ and yδ (t ) = y(t ) − y∗ , we obtain a two dimensional system of linearized equations
x˙ δ (t ) = α [−2bxδ (t ) − byδ (t − τ )], y˙ δ (t ) = β [−bxδ (t ) − 2byδ (t )]
(10)
α = αx x∗ and β = αy y∗ .
(11)
with
Substituting exponential solutions
xδ (t ) = eλt u and yδ (t ) = eλt v
(12)
into the linearized system (10) and arranging terms yield an alternative linear system,
λ + 2bα bβ
bα e−λτ λ + 2bβ
u
v
=
0 . 0
(13)
Excluding the trivial solutions (i.e., u = v = 0), we obtain non-trivial solutions by solving the characteristic equation
(λ + 2bα)(λ + 2bβ) − b2 αβ e−λτ = 0. 3
[23] and [13] assume homogenous expectation formations in which both firms use past production data to determine their current outputs.
(14)
702
A. Matsumoto, F. Szidarovszky / Applied Mathematics and Computation 269 (2015) 699–713
For τ = 0, the characteristic equation is reduced to
λ2 + 2b(α + β)λ + 3αβ b2 = 0.
(15) 3αβ b2
> 0, the non-delay characteristic equation has roots with negative real parts, Since, by assumption, 2b(α + β) > 0 and implying that the stationary point is locally stable if there is no delay. It is true that the stationary state preserves stability as far as the delay is positive but sufficiently small. We are concerned with a threshold value (if exists) of the delay for which the stationary point is just destabilized. It is well known that if the stability of the stationary point switches at τ = τ¯ , then the characteristic equation must have a pair of pure conjugate imaginary roots there. To examine the stability switches of the dynamic system (10), we determine this threshold value τ¯ . For this purpose, we substitute a purely imaginary solution λ = iω with ω > 0 into Eq. (14) and separate the real and imaginary parts of the resulted equation,
−ω2 + 4b2 αβ = b2 αβ cos τ ω, −2b(α + β)ω = b2 αβ sin τ ω.
(16)
Squaring both sides of these equations and adding the resulted expressions yield, after arranging terms, a quartic equation in ω,
ω4 + 4b2 (α 2 + β 2 )ω2 + 15(b2 αβ)2 = 0.
(17)
Since the left hand side is positive for all real values of ω, stability switch cannot occur and therefore the steady state is always stable for any τ > 0. In other words, time delay is harmless. Theorem 1. Under the expectation formation (E1), the continuous-time dynamic system (8) is locally asymptotically stable for any value of τ > 0. 3.2. Two fixed delays Assuming (E2) as the output expectation formation and substituting the exponential forms, x(t ) = eλt u and y(t ) = eλt v into the linearized Eq. (10) yields
λ + 2bα bβ
bα(η1 e−λτ1 + η2 e−λτ2 ) λ + 2bβ
u
v
0 . 0
=
(18)
Assuming nontrivial solution gives the corresponding characteristic equation,
(λ + 2bα)(λ + 2bβ) − b2 αβ(η1 e−λτ1 + η2 e−λτ2 ) = 0.
(19)
This is equivalent to equation
a(λ, τ1 , τ2 ) = 1 + a1 (λ)e−λτ1 + a2 (λ)e−λτ2 = 0
(20)
−b2 αβ η for k = 1, 2. (λ + 2bα)(λ + 2bβ) k
(21)
where
ak (λ) =
Suppose λ = iω with ω > 0, then
ak (iω) =
−αβ b2 [4b2 αβ − ω2 − i2b(α + β)ω] ηk for k = 1, 2 [4b2 αβ − ω2 ]2 + [2b(α + β)ω]2
(22)
and their absolute values are
|ak (iω)| =
b2 αβ [4b2 αβ − ω2 ]2 + [2b(α + β)ω]2
ηk for k = 1, 2.
(23)
[11] geometrically investigates the two delay equation based on the idea that solving a(iω, τ1 , τ2 ) = 0 is equivalent to forming a triangle by the three terms in (20) as three vectors in the complex plane placed head to tail. Since the triangle consists of three line segments, it is a necessary and sufficient condition for the existence of a positive solution ω > 0 of Eq. (20) that the sum of the lengths of any two adjacent line segments is not shorter than the length of the remaining line segment. Theorem 2. A purely imaginary root λ = iω with ω > 0 is a solution of a(iω, τ1 , τ2 ) = 0 if and only if the following three inequalities hold:
|a1 (iω)| + |a2 (iω)| ≥ 1, |a1 (iω)| + 1 ≥ |a2 (iω)|, |a2 (iω)| + 1 ≥ |a1 (iω)|.
(24)
Proof. Let a, b and c be three line segments. The end points of a and b with a common starting point can be connected with a segment of length c if and only if
|a − b| ≤ c ≤ a + b.
(25)
A. Matsumoto, F. Szidarovszky / Applied Mathematics and Computation 269 (2015) 699–713
703
The second inequality gives one condition. The first inequality can be rewritten as
−c ≤ a − b ≤ c,
(26)
which means that
b≤a+c
(27)
a ≤ b + c.
(28)
and
These inequality conditions complete the proof. We check whether the first condition of Theorem 2 holds. Substituting (23) into it presents the inequality condition in parametric terms,
ω4 + 4b2 (α 2 + β 2 )ω2 + 15(b2 αβ)2 ≤ 0.
(29)
Notice that the left hand side is positive with all real values of ω implying that stability switch cannot occur for τ 1 > 0 and τ 2 > 0. Theorem 3. Under the expectation formation (E2), the continuous-time dynamic system (8) is locally asymptotically stable for any values of τ 1 > 0 and τ 2 > 0. 3.3. Three fixed delays In the similar way, we can show that the expectation formation (E3) cannot destabilize the stationary state as well. As before, substituting the exponential forms, x(t ) = eλt u and y(t ) = eλt v into the linearized system under (E3) yields
λ + 2bα bβ
bα
−λτ u 0 η1 e 1 + η2 e−λτ2 + η3 e−λτ3 = . v 0 λ + 2bβ
(30)
Nontrivial solution exists if and only if
(λ + 2bα)(λ + 2bβ) − b2 αβ(η1 e−λτ1 + η2 e−λτ2 + η3 e−λτ3 ) = 0.
(31)
As before, this equation can be rewritten as
a(λ, τ1 , τ2 , τ3 ) = 1 + a1 (λ)e−λτ1 + a2 (λ)e−λτ2 + a3 (λ)e−λτ3 = 0
(32)
where
ak (λ) =
−b2 αβ η for k = 1, 2, 3, (λ + 2bα)(λ + 2bβ) k
(33)
Suppose λ = iω with ω > 0, then
ak (iω) =
−b2 αβ [4b2 αβ − ω2 − i2b(α + β)ω] ηk for k = 1, 2, 3, [4b2 αβ − ω2 ]2 + [2b(α + β)ω]2
(34)
and their absolute values are
|ak (iω)| =
b2 αβ [4b2 αβ − ω2 ]2 + [2b(α + β)ω]2
ηk for k = 1, 2, 3.
(35)
[1] gives a straightforward extension of the two delay case considered in [11] to the more general case of three delays and obtains the following result: Theorem 4. Pure imaginary root λ = iω with ω > 0 is a solution of a(iω, τ1 , τ2 , τ3 ) = 0 if and only if the following four inequalities hold:
|a1 (iω)| + |a2 (iω)| + |a3 (iω)| ≥ 1, |a2 (iω)| + |a3 (iω)| + 1 ≥ |a1 (iω)|, |a1 (iω)| + |a3 (iω)| + 1 ≥ |a2 (iω)|, |a1 (iω)| + |a2 (iω)| + 1 ≥ |a3 (iω)|.
(36)
We again check the first condition, which, after substituting (35), can be written as
ω4 + 4b2 (α 2 + β 2 )2 ω2 + 15(b2 αβ)2 ≤ 0.
(37)
Similarly to the two delay case, the left hand side is positive with any real value of ω, implying that stability switch cannot occur for any positive delays. Theorem 5. Under the expectation formation (E3), the continuous-time dynamic system (8) is locally asymptotically stable for any values of τ 1 > 0, τ 2 > 0 and τ 3 > 0.
704
A. Matsumoto, F. Szidarovszky / Applied Mathematics and Computation 269 (2015) 699–713
Fig. 1. Stable trajectories converging to x∗ .
3.4. Continuously distributed delays The expectation formation (E4) assumes that the expected output is a weighted average of all past outputs from time 0 to time t. It is a generalization of (E2) and (E3) with infinitely many continuously distributed delays. Adding the time derivative of the expectation (E4) to linearized system (10) yields a three dimensional system of differential equations,
x˙ (t ) = α [−2bx(t ) − bye (t )], y˙ (t ) = β [−bx(t ) − 2by(t )], 1 y˙ e (t ) = [y(t ) − ye (t )]. T
(38)
The corresponding characteristic equation is given by
λ3 + a1 λ2 + a2 λ + a3 = 0
(39)
where
a1 =
1 + 2(α + β)bT > 0, T
(40)
a2 =
2(α + β)b + 4αβ b2 T >0 T
(41)
a3 =
3αβ b2 > 0. T
(42)
and
A set of necessary and sufficient conditions for all roots of the cubic characteristic equation to have negative real parts is
ai > 0 for i = 1, 2, 3 and a1 a2 − a3 > 0,
(43)
which is a special case of the Routh–Hurwitz criterion. Since all coefficients are positive, we need to check whether the last condition is satisfied or not. Clearly
a1 a2 − a3 =
b [8b2 αβ(α + β)T 2 + b(4α 2 + 9αβ + 4β 2 )T + 2(α + β)] T2
(44)
is always positive for any T > 0. Hence, the stationary state of the three dimensional system is locally stable. Theorem 6. Under the expectation formation (E4), the continuous-time dynamic system (8) is locally asymptotically stable for any values of T > 0. Analytical results obtained so far are numerically confirmed in Fig. 1 in which time trajectories of x(t) in various cases are plotted in different colors.4 Constant initial functions are selected, ϕx (t ) = 5 and ϕy (t ) = 1 for t ≤ 0, the parameters are fixed at αx = αy = 1, a = 4, b = 1, cx = cy = 1, the delays are selected to be τ = 5 in one delay case, τ1 = 5 and τ2 = 6 in two delay case, τ1 = 5, τ2 = 6 and τ3 = 7 in three delay case and T = 5 in the continuously distributed case. It can be seen that all trajectories converge to the steady state of x as analytically shown above. Two issues are also observed: one is that the black real and black dotted trajectories monotonically converge and the latter approaches the steady state more rapidly; the other is that three trajectories in the fixed delay cases start to fluctuate at t = 5 and convergence is oscillatory. 4
“C delay” in the legend of Fig. 1 means continuously distributed delay.
A. Matsumoto, F. Szidarovszky / Applied Mathematics and Computation 269 (2015) 699–713
705
4. Own and competitor’s delays In this section we examine the effects caused by the two delays on dynamics: one delay in the competitor’s output and the other delay in the firm’s own output. Thus the dynamic system is modified in the following way:
x˙ (t ) = αx x(t )[a − cx − 2bx(t − τx ) − by(t − τy )], y˙ (t ) = αy y(t )[a − cy − bx(t ) − 2by(t )]
(45)
where τ x > 0 is the implementation delay and τ y > 0 is the competitor’s information delay. The stationary point is the same as given in (5). The linearized system is
x˙ δ (t ) = α [−2bxδ (t − τx ) − byδ (t − τy )], y˙ δ (t ) = β [−bxδ (t ) − 2byδ (t )].
(46)
Substituting exponential solutions xδ (t ) = eλt u and yδ (t ) = eλt v reduces the linear system to an alternative form,
λ + 2α be−λτx bβ
bα e−λτy λ + 2bβ
u
=
v
0 . 0
(47)
Nontrivial solution exists if and only if
λ2 + 2bβλ + 2bα(λ + 2bβ)e−λτx − b2 αβ e−λτy = 0.
(48)
Before turning to a closer examination of Eq. (48), a few remarks should be made concerning the effect caused by the implementation delay. Assuming that τy = 0 and λ = iω with ω > 0, we can separate Eq. (48) into real and imaginary parts,
2bα(2bβ cos τx ω + ω sin τx ω) = ω2 + b2 αβ , (49)
2bα(ω cos τx ω − 2bβ sin τx ω) = −2bβω. Adding the squares of these equations yields a fourth-order polynomial equation in ω,
ω4 + 2b2 [αβ + 2(β 2 − α 2 )]ω2 − 15b4 α 2 β 2 = 0 that can be solved for
ω2
(50)
and the positive solution is
ω+2 = b2 {−[αβ + 2(β 2 − α 2 )] +
4β 2 + 4αβ 3 + 7α 2 β 2 + (2α 2 − αβ)2 } > 0
(51)
Solving (49) for cos τ x ω and sin τ x ω , substituting ω+ into cos τ x ω and then solving it for τ give the threshold value of the delay,
τxn =
1
cos−1
ω+
b2 β 2 2 ω+ + 4b2 β 2
+ 2nπ , for n = 0, 1, 2, . . .
(52)
A routine bifurcation analysis shows that stability is lost at the smallest threshold value,
τx0 =
1
ω+
cos−1
b2 β 2 2 ω+ + 4b2 β 2
.
(53)
and stability cannot be regained later. In the bifurcation analysis we consider the eigenvalues as functions of the bifurcation parameter; λ = λ(τx ), which is substituted into the characteristic Eq. (48). By taking τy = 0 and implicitly differentiating the resulted equation with respect to the bifurcation parameter, we can have,
{2λ + 2bβ + 2bα [1 − (λ + 2bβ)τx ]e−λτx }
dλ − 2bα(λ + 2bβ)λe−λτx = 0 d τx
(54)
which is solved for (dλ/dτx )−1 for analytical simplicity,
dλ d τx
−1
=
2λ + 2bβ + 2bα [1 − (λ + 2bβ)τx ]e−λτx 2bα(λ + 2bβ)λe−λτx
(55)
Since from(48) with τy = 0,
λ2 + 2bβλ − b2 αβ = −2bα(λ + 2bβ)e−λτx ,
(56)
the right hand side of (55) can be written as
dλ d τx
−1
=
2(λ + bβ) + −λ(λ2 + 2bβλ − b2 αβ)
1
λ(λ + 2bβ)
−
τx . λ
(57)
706
A. Matsumoto, F. Szidarovszky / Applied Mathematics and Computation 269 (2015) 699–713
Its real part at λ = iω is, after arranging terms,
Re
dλ d τx
−1
= Re λ=iω
2(iω + bβ) −iω( − ω2 − b2 αβ + i2bβω)
+ Re
1 iω(iω + 2bβ)
1 ω2 + (2bβ)2 + b2 αβ − (2bβω)2 + (ω2 + b2 αβ)2 ω2 + (2bβ)2 b2 β [(4β − α)(ω2 + b2 αβ) + (4bβ)2 ] = [(2bβω)2 + (ω2 + b2 αβ)2 ][ω2 + (2bβ)2 ] =
(58)
The last expression is positive if 4β > α or the adjustment coefficients do not differ so much, which is assumed to be henceforth. Thus we obtain the derivatives λ (τxn ) at the critical values τxn . Since these derivatives are positive with 4β > α , at all critical points at least one eigenvalue pair changes the sign of its real part from negative to positive. So stability is lost at the smallest critical point τx0 and the stability cannot be regained later. Theorem 7. If 4β > α holds, then the stationary point of (45) with τy = 0 is locally asymptotically stable for τx < τx0 , loses stability at τx = τx0 and bifurcates to a limit cycle via Hopf bifurcation for τx > τx0 . We now investigate the characteristic equation with two positive delays by applying Gu’s method [11]. Dividing both sides by
λ2 + 2bβλ and introducing new functions a1 (λ) =
2bα
λ
and a2 (λ) =
−b2 αβ λ(λ + 2bβ)
(59)
simplify Eq. (48) as
a(λ, τx , τy ) = 1 + a1 (λ)e−λτx + a2 (λ)e−λτy = 0.
(60)
Suppose that λ = iω with ω > 0. Substituting it into ai (λ) yields
a1 (iω) = −i
2bα
(61)
ω
and
a2 (iω) =
b2 αβω 2b3 αβ 2 + i . ω(ω2 + 4b2 β 2 ) ω(ω2 + 4b2 β 2 )
(62)
Their absolute values are
|a1 (iω)| =
2bα
(63)
ω
and
|a2 (iω)| =
b2 αβ
ω ω2 + 4b2 β 2
.
(64)
We check whether the three conditions of Theorem 2 are satisfied. Substituting the absolute values reduces the three inequalities to the following two conditions,
f (ω) = ω4 − 4bαω3 + 4b2 (α 2 + β 2 )ω2 − 16b3 αβ 2 ω + 15b4 α 2 β 2 ≤ 0
(65)
g(ω) = ω4 + 4bαω3 + 4b2 (α 2 + β 2 )ω2 + 16b3 αβ 2 ω + 15b4 α 2 β 2 ≥ 0.
(66)
and
Notice that g(ω) ≥ 0 is always true for ω ≥ 0. However, it is ambiguous whether f(ω) ≤ 0 holds or not. Using the variable transformation
4bα ω =x− − = x + bα ,
(67)
4
we rewrite f(ω) as
F (x) = x4 − 2b2 (α 2 − 2β 2 )x2 − 8b3 αβ 2 x + b4 α 2 (α 2 + 3β 2 ).
(68)
α2
> (+, +, −, +) if < and (+, −, +) if The sequence of the signs of the polynomial coefficients are (+, −, −, +) if α 2 = 2β 2 . In any sequence, the number of sign changes is two. According to Descartes’ rule of sign, this polynomial has either two or zero positive roots. It is clear that F(0) > 0 and F (∞) = ∞. Substituting x = α b yields
F (α b) = α 2 b4 ( − β 2 ) < 0,
2β 2 ,
α2
2β 2 ,
(69)
A. Matsumoto, F. Szidarovszky / Applied Mathematics and Computation 269 (2015) 699–713
707
which implies that F (x) = 0 has two positive solutions, x1 and x2 where 0 < x1 < α b and x2 > α b. Hence f (ω) = 0 has two positive solutions
ω1 = x1 + α b > 0 and ω2 = x2 + α b > 0
(70)
and f(ω) ≤ 0 for ω ∈ [ω1 , ω2 ], and a routine investigation can show that f(ω) > 0 otherwise.5 We will next find all pairs of (τ 1 , τ 2 ) satisfying a(iω, τ1 , τ2 ) = 0. The three terms in (60) are three vectors in the complex plane that construct a triangle. Let us suppose that 1 is its base and let us denote the angle between 1 and |a1 (iω)e−iωτx | by θ 1 and an angle between 1 and |a2 (iω)e−iωτy | by θ 2 . Since, |e−iωτx | = |e−iωτx | = 1, we have, by the law of cosine,
θ1 (ω) = cos−1 = cos−1 and
θ2 (ω) = cos
−1
= cos−1
2
2 2 · 1 · a1 (iω)
12 + a1 (iω) − a2 (iω)
ω4 + 4b2 (α 2 + β 2 )ω2 + 15b4 α 2 β 2 4bαω(ω2 + 4b2 β 2 )
(71)
2
2 2 · 1 · a2 (iω)
12 + a2 (iω) − a1 (iω)
ω4 − 4b2 (α 2 − β 2 )ω2 − 15b4 α 2 β 2 . 2b2 αω ω2 + 4b2 β 2
(72)
Since the triangle may be located above and also under the horizontal axis, we get two equations as
{arg[a1 (iω)e−iτx ω ] + 2mπ } ± θ1 (ω) = π {arg[a2 (iω)e−iτy ω ] + 2nπ } ∓ θ2 (ω) = π which yield the threshold values of the delays:
τx± (ω, n) =
τ (ω, n) =
1 3π + (2m − 1)π ± θ1 (ω) ω 2
and ∓ y
(73)
1
ω
tan
−1
2bβ
ω
+ (2n − 1)π ∓ θ2 (ω) .
(74)
(75)
Let be the set of all ω values for which f(ω) ≤ 0 holds. Then we can find all pairs of (τ 1 , τ 2 ) constructing the stability switching curves for ω ∈ which consists of two sets of parametric segments,
L1 (m, n) = {τx+ (ω, m), τy− (ω, n)}
(76)
L2 (m, n) = {τx− (ω, m), τy+ (ω, n)}.
(77)
and
To illustrate the stability switching curves we first specify the parameter values that make the mathematical derivations simple without changing the main conclusions of the analysis: Assumption: α = β = b = 1. Under this assumption, we have
f (ω) = ω4 − 4ω3 + 8ω2 − 16ω + 15.
(78)
It is not difficult to show that f (ω) = 0 has two real and positive roots,
1
ωs = 1 +
2 3/K 1
ωe = 1 +
2 3/K
with
5
−
1 2
+
1 2
− 13 (12 + K ) + 16 3/K 1.611
− 13
(79)
(12 + K ) + 16 3/K 2.326
√ √ 1 1 K = −4 + (584 − 48 87) 3 + 2(73 + 6 87) 3 .
A formal proof is given in the Appendix.
(80)
708
A. Matsumoto, F. Szidarovszky / Applied Mathematics and Computation 269 (2015) 699–713
The domain of ω is the interval [ωs , ωe ]( = ). So the segment (76) starts at point (τx+ (ωs , m), τy− (ωs , n)) and terminates at point (τx+ (ωe , m), τy− (ωe , n)) and so does the segment (77) at points (τx− (ωs , m), τy+ (ωs , n)) and (τx− (ωe , m), τy+ (ωe , n)). Furthermore the angles of the triangle are
ω4 + 8ω2 + 15 θ1 (ω) = cos 4ω(ω2 + 4) ω4 − 15 . θ2 (ω) = cos−1 √ 2ω ω 2 + 4
−1
(81)
Since
ωs4 + 8ωs2 + 15 = 1, 4ωs (ωs2 + 4)
ωe4 + 8ωe2 + 15 =1 4ωe (ωe2 + 4)
(82)
ωs4 − 15 = −1, 2ωs ωs2 + 4
ωe4 − 15 = 1, 2ωe ωe2 + 4
(83)
and
we have
θ1 (ωs ) = 0, θ1 (ωe ) = 0
(84)
θ2 (ωs ) = π , θ2 (ωe ) = 0.
(85)
and
We can now formulate the following results concerning the locations of the curve segments: Theorem 8. Given m and n, the segments L1 (m, n + 1) and L2 (m, n) have the same starting point whereas the segments L1 (m, n) and L2 (m, n) have the same end point. Proof. For analytical simplicity we assume m = 0. The starting points of L1 (0, n) and L2 (0, n) are
Ls1 (0, n) = {τx+ (ωs , 0), τy− (ωs , n)}
(86)
with
τx+ (ωs , 0) =
1
π
ωs 2
and τy− (ωs , n) =
1
tan−1
ωs
2 + 2(n − 1)π , ωs
(87)
and
Ls2 (0, n) = {τx− (ωs , 0), τy+ (ωs , n)}
(88)
where
τx− (ωs , 0) =
1
π
ωs 2
and τy+ (ωs , n) =
1
ωs
tan−1
2 + 2nπ . ωs
(89)
Hence we have
Ls1 (0, n + 1) = Ls2 (0, n).
(90)
In the same way, the end points of L1 (0, n) and L2 (0, n) are
Le1 (0, n) = {τx+ (ωe , 0), τy− (ωe , n)}
(91)
with
τx+ (ωe , 0) =
1
π
ωe 2
and τy− (ωe , n) =
1
ωe
tan−1
2 + (2n − 1)π ωe
(92)
and
Le2 (0, n) = {τx− (ωe , 0), τy+ (ωe , n)}
(93)
where
τx− (ωe , 0) =
1
π
ωe 2
and τy+ (ωe , n) =
1
ωe
tan−1
2 + (2n − 1)π . ωe
(94)
Hence we have
Le1 (0, n) = Le2 (0, n). The same result is obtained for any integer m > 0.
(95)
A. Matsumoto, F. Szidarovszky / Applied Mathematics and Computation 269 (2015) 699–713
709
Fig. 2. Stability switching curve.
This result is numerically confirmed in Fig. 2 where m = 0. The stability switching curve consists of the red and blue segments that correspond to L1 (0, n) and L2 (0, n) for n = 0, 1, 2. The red and blue segments shift upward when n increases and to the right when m increases. The upward-sloping blue segment crossing the horizontal axis at τx = τx0 is the L2 (0, 0) segment.6 It connects to the red L1 (0, 1) segment at a point where Ls1 (0, 1) = Ls2 (0, 0) holds. As before, the upper script “s” means the starting (i.e., initial) point of the segment and
Ls1 (0, 1) = (τx+ (ωs , 1), τy− (ωs , 1)) and Ls2 (0, 0) = (τx− (ωs , 0), τy+ (ωs , 0)). This L1 (0, 1) segment then connects to the blue L2 (0, 1) segment at a point where “e” means again the end (i.e., terminal) point of the segment and
Le1 (0, 1) = (τx+ (ωe , 1), τy− (ωe , 1)) and Le2 (0, 1) = (τx− (ωe , 1), τy+ (ωe , 1)).
(96) Le1 (0, 1)
=
Le2 (0, 1)
holds. The upper script (97)
As n increases, the two segments are connected in the same way to construct the continuous stability switching curve. The eigenvalues are purely imaginary on this curve. The stability switching curves L1 (0, n) and L2 (0, n) divide the first quadrant of the (τ 1 , τ 2 ) plane into two parts. One contains the origin and its every point can be reached from the origin via continuous curve not crossing the stability switching curve. At the points in this region the real parts of the eigenvalues are negative, so the system is locally asymptotically stable. The complement of this region except the curves gives points when the system is unstable. We again notice that the blue L2 (0, 0) segment crosses the horizontal axis at τx = τx0 ( 0.824). This indicates that stability is preserved for τx < τx0 and lost for τx > τx0 when τy = 0. In other words, this is the threshold value of τ x obtained in the case of one delay (Theorem 7), which is equal to 0.824 under Assumption. So we could confirm this threshold value in two different ways. The direction of change (that is, stability loss or regain) can be checked based on the Hopf bifurcation theory and at each point (τ x , τ y ) the number of stability losses or stability regains determine if this point is stable or not as in [18]. Instead of following this lengthy and complicated method, we perform numerical simulations to investigate how stability changes in case of two positive delays. We now examine the effects caused by changing the value of τ x , keeping the value of τ y at some positive fixed value. In Fig. 3(A) we increase the value of τ x from 0 to 1.3 in 0.001 increments along the horizontal dotted line with fixed value τy = τya ( = 2). For each value of τ x , the dynamic system runs for 0 ≤ t ≤ 1000 and only the data for 950 ≤ t ≤ 1000 are used to get rid of the transients. The local maxima and minima are plotted against each value of τ x . If the bifurcation diagram has one point against τ x , then the maximum point is identical with the minimum point, implying that the stationary state is locally stable. If it has two points, then a limit cycle having one maximum and one minimum emerges and if many points, then a limit cycle exhibits many ups and downs. The horizontal dotted line crosses the blue L2 (0, 1) segment at τxa 0.652 denoted by the green dot in Fig. 2. It is seen that the stationary state is locally asymptotically stable for τx < τxa and becomes unstable for τx > τxa as shown in Fig. 3(A). When stability is lost at τx = τxa , then the stationary point bifurcates to a limit cycle and does not regain stability for larger values of τ x . We move to the second simulation in which the fixed value of τ y is increased to τyb = 4.4 and the same procedure is repeated. As is seen in Fig. 2, the horizontal dotted line at τyb crosses the stability switching curve three times as denoted by three green dots. As described in Fig. 3(B), both stability regain and stability loss occur in this example. In particular, stability is first lost at the first crossing point, the left most green dot denoted by τx1 ( 0.699). Then limit cycle emerges for τ x in the interval [τx1 , τx2 ] where τx2 ( 0.892) corresponds to the middle green point. If the τ x -value of the most right green point is denoted by τx3 ( 0.964), then stability is regained for τx ∈ [τx2 , τx3 ]. The stationary state again loses its stability for τx > τx3 . The bifurcation diagram shows that a trajectory gradually exhibits more complex dynamics through a quasi period-doubling process as τ x become larger and then returns to the simple limit cycle through a quasi period-halving process as τ x further increases. 6
The L1 (0, 0) segment is located in the fourth quadrant so it is not illustrated.
710
A. Matsumoto, F. Szidarovszky / Applied Mathematics and Computation 269 (2015) 699–713
Fig. 3. Bifurcation diagram with respect to τ x , given τ y .
Fig. 4. Bifurcation diagram with respect to τ y , given τ x .
We now turn attention to the delay effect caused by changing values of τ y , keeping the value of τ x at some positive value. We also perform two simulations. In Fig. 2, τ x is fixed at τxA = 0.75 and τ y increases along the vertical dotted line crossing the stability switching curve three times denoted by three black dotes. Although they are not labelled on the vertical axis to avoid confusion in the figure, we denote their τ y -values by τy1 ( 1.410), τy2 ( 3.138) and τy3 ( 4.215) with ascending order as shown in Fig. 4(A). The bifurcation diagram indicates that stability is lost at τy = τy1 , a limit cycle emerges for larger values than τy1 and then regains stability at τy = τy2 . Stability is preserved for τy < τy3 and is lost again for τy > τy3 . The diagram implies a similar dynamic cascade in which the process of stability loss, birth of limit cycle and regain of stability is repeated as the value of τ y increases. In Fig. 4(B) the fixed value of τ x is increased to τxB = 1.2 from τxA = 0.75 and the value of τ y increases along the vertical dotted line at τxB = 1.2 which is located to the right of the stability switching curve in Fig. 2. Even for τy = 0, the dynamic system is unstable as τxB > τx0 and generates a limit cycle with one maximum and one minimum. As the value of τ y increases, the bifurcation diagram implies that the diameter of the limit cycle varies and then more complex dynamics appears through a period-doubling like process. 5. Concluding remarks Delay dynamics of heterogeneous duopolies in continuous-time framework was investigated. This is an extension of the earlier study [8] on heterogeneous duopolies in a discrete-time framework where the authors investigated a heterogenous duopoly market in which one firm forms its expectation on the competitor’s output only with delayed information and the other firm uses the current information. Following their spirit, we assumed that one firm has a delay in obtaining information about the competitor’s output (i.e., information delay) in four different ways in which one, two or three delayed information is used to form expectation on the output of the competitor, or all past output values are used as a continuously distributed delay. Applying a recently developed method to examine stability of delay differential equations, we analytically demonstrated that the information delay is harmless to stability in the continuous-time duopoly model. In other words, no stability switch occurs regardless
A. Matsumoto, F. Szidarovszky / Applied Mathematics and Computation 269 (2015) 699–713
711
of the expectation formation and the lengths of the delays, meaning that the system is stable with any values of the delays, so delays are harmless in all examined cases. This shows a sharp difference from the main result of the discrete case that a delay in the competitor’s production can have a stabilizing effect. Then we assumed that one firm had implementation and information delays and the other firm could make decisions without any delay. We first constructed the stability switching curve which consists of a sequence of connected segments forming a continuous curve as function of the two delays. We confirmed that stability can be switched to instability when the length of the implementation delay takes over some threshold value which depends on the length of the information delay. It is numerically verified that various dynamics ranging from simple to complex can emerge according to the values of the delays. If one firm had only implementation delay, then stability was lost when the delay exceeded a certain threshold value. Therefore it was demonstrated that the implementation delay matters. We also examined the dynamics of the system with changing the lengths of the delays. With fixed small information delay, the system lost stability at a threshold of the implementation delay and the stability could not be regained later. However with larger fixed information delay the stability is lost at a threshold of the implementation delay and the stability was gained for a while with the further increase of the implementation delay. Similar situation occurs with small fixed implementation delay, and if this delay is larger, then the system is unstable for all lengths of the information delay. Hopf bifurcation occurs at the critical values giving the possibility of the birth of limit cycles. There was a variety of dynamic behavior in the case of instability including limit cycles and period-doubling like processes in losing stability and stability can be regained by a period-halving process. The methodology used in this paper to determine stability switching curves and the stability regions can be used in many other delay models arising in several fields of mathematical economics. There are many directions for future research. In this study, we assumed heterogenous duopolies in which only one firm had delayed information formulated in (45). One possible extension is given below in which duopolies are heterogeneous and one firm has an implementation delay and the other has an information delay,
x˙ (t ) = αx x(t )[a − cx − 2bx(t − τx ) − by(t )], y˙ (t ) = αy y(t )[a − cy − bx(t − τy ) − 2by(t )].
(98)
However, the characteristic equation derived from this dynamic system is identical with Eq. (48). It may be interesting to consider duopoly models with different economic backgrounds generating exactly the same dynamics. The second extension concerns the way of the transformation to a continuous-time model from a discrete-time model. We make two assumptions: one is replacing the discrete time difference z(t + 1) − z(t ) with time derivative z˙ (t ) and the other is introducing continuous-time expectation formations concerning ye (t) instead of the discrete-time formation. The delay discrete-time formation was given as
ye (t + 1) = ηy(t ) + (1 − η)y(t − 1)
(99)
which can be written as
ye (t + 1) = =
η[y(t ) − y(t − 1)] + y(t − 1) ηy˙ (t ) + y(t − 1).
(100)
If the unit time difference is replace with delay τ , then we have the following form of a continuous-time system:
x˙ (t ) = αx x(t )[a − cx − 2bx(t ) − b{ηy˙ (t ) + y(t − τ )}], y˙ (t ) = αy y(t )[a − cy − bx(t ) − 2by(t )].
(101)
This can be solved for x˙ (t ) and y˙ (t ) to derive a dynamic system of explicit delay differential equations. Further it is possible to introduce an implementation delay on the firm’s own output although analysis will be more complicated. The third extension is to adopt the modelling method proposed in [3] and applied in [20] when considering delay monopoly dynamics. The continuoustime system is described by the following delay equations,
σx x˙ (t ) + x(t ) = x(t − τx ) + αx x(t − τx )[a − cx − 2bx(t − τx ) − bye (t )], σy y˙ (t ) + y(t ) = y(t − τx ) + βy y(t − τy )[a − cy − bxe (t ) − 2by(t − τy )],
(102)
in which σ x ≥ 0 and σ y ≥ 0 denote the inertias (or frictions) in the production process. It should be noticed that this system can be reduced to system (4) if σx = σy = 0 and τx = τy = 1. So when σ x and σ y take smaller values, the generated dynamics is similar to the discrete-time case (4). On the other hand, when σ x and σ y take larger values, it may generate different dynamics. Appendix In this Appendix, we will prove that equation
x4 − 2b2 (α 2 − 2β 2 )x2 − 8b3 αβ 2 x + b4 α 2 (α 2 + 3β 2 ) = 0
(103)
has two distinct positive roots and no negative root exists. Let x = α bz, then
α 4 b4 z4 − 2b2 (α 2 − 2β 2 )α 2 b2 z2 − 8b3 αβ 2 α bz + b4 α 2 (α 2 + 3β 2 ) = 0
(104)
which can be simplified as
h(z) = α 2 z4 − 2(α 2 − 2β 2 )z2 − 8β 2 z + (α 2 + 3β 2 ) = 0.
(105)
712
A. Matsumoto, F. Szidarovszky / Applied Mathematics and Computation 269 (2015) 699–713
Notice that
h(0) = α 2 + 3β 2 > 0,
lim h(z) = ∞ and h(1) = −β 2 < 0.
z→±∞
(106)
The derivative of h(z) has the form
h (z) = 4(z − 1)(α 2 z2 + α 2 z + 2β 2 ).
(107)
The sign of this derivative depends on the sign of the quadratic function
g(z) = α 2 z2 + α 2 z + 2β 2
(108)
with the discriminant
D = α 2 (α 2 − 8β 2 ).
(109)
If α 2 = 8β 2 , then it is positive except its negative vertex, so h (z) < 0 as z < 1 and differs from the vertex, h (z) > 0 as z > 1 and h (z) = 0 as z = 1 or equals the vertex. Therefore h(z) strictly decreases as z < 1 and strictly increases as z > 1. If α 2 < 8β 2 , then h (z) > 0 as z > 1 and h (z) < 0 as z < 1. Since h(0) > 0 and both limits at ± ∞ are positive there are two positive roots of h(z) and no negative root exists. Assume next that α 2 > 8β 2 . Then the quadratic function g(z) has two negative roots:
z1,2 = −
1 ± 2
α 2 − 8β 2 (z1 < z2 ). 2α
(110)
,
(111)
So h (z) < 0 if z < z1 or z2 < z < 1, and h (z) > 0 if z1 < z < z2 or z > 1. At these roots
z2 = −z − so
2β 2
α2
h(z) = =
α
2
−z −
2β 2
2
α2
(α 2 − 8β 2 )z +
2β 2 − 2(α − 2β ) −z − 2 − 8β 2 z + (α 2 + 3β 2 ) α 2
1
α2
2
(α 4 + 5α 2 β 2 − 4β 4 ).
(112)
Since z1 and z2 are negative and their average value is −1/2, both are larger than −1. Since α 2 > 8β 2 , at the roots
(α 2 − 8β 2 )z > 8β 2 − α 2 .
(113)
Hence
h(z) > 8β 2 − α 2 +
1
α2
(α 4 + 5α 2 β 2 − 4β 4 )
β2 (5α 2 − 4β 2 ) α2 β2 > 8β 2 + 2 (40β 2 − 4β 2 ) > 0. α = 8β 2 +
(114)
So h(z) decreases from the ∞ limit as z → −∞ to h(z1 ) > 0, then increases to h(z2 ) > 0 and decreases again until h(1) < 0 and then increases and tends to ∞ as z → ∞ . Since h(0) > 0, there are no negative roots, only two distinct positive roots. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17]
E. Almoderesi, M. Bozorg, Stability crossing surfaces for linear time delay systems with three delays, Int. J. Control 82 (2009) 2304–2310. R. Bellman, K. Cooke, Differential-Difference Equations, Academic Press, New York, 1956. M. Berezowski, Effect of time delay on the generation of chaos in continuous systems, Chaos, Solitons Fractals 12 (2001) 83–89. G.-I. Bischi, A. Namzada, Global analysis of a dynamic duopoly game with bounded rationality, in: Advances in Dynamic Games and Application, vol. 5, Birkhauser, Boston, Massachussetts, USA, 2000, pp. 361–385. (Chapter 20). G.-I. Bischi, C. Chiarella, M. Kopel, F. Szidarovszky, Nonlinear Oligopolies: Stability and Bifurcations, Springer, Berlin/Heidelberg/New York, 2010. A. Cournot, Recherches sur les Principes Mathématiques de la Théorie des Richessess, Hachette, Paris, 1838.(English translation (1960): Researches into the Mathematical Principles of the Theory of Wealth, Kelley, New York) J. Cushing, Integro-Differential Equations and Delay Models in Population Dynamics, Springer-Verlag, Berlin/Heidelberg/New York, 1977. A. Elsadany, A. Matouk, Dynamic Cournot duopoly game with delay, J. Complex Syst. (2014). http://dx.doi.org/10.1155/2014/384843 M. Ezekiel, The cobweb theorem, Q. J. Econ. 52 (1933) 255–280. R. Goodwin, The nonlinear accelerator and the persistence of business cycles, Econometrica 19 (1951) 1–17. K. Gu, S.-I. Nicolescu, J. Chen, On stability crossing curves for general systems with two delays, J. Math. Anal. Appl. 311 (2005) 231–253. J. Haldane, A contribution to the theory of price fluctuations, Rev. Econ. Stud. 1 (1933) 186–195. S.Z. Hassan, On delayed dynamical duopoly, App. Math. Comput. 151 (2004) 275–286. T. Howroyd, A. Russel, Cournot oligopoly models with time delays, J. Math. Econ. 13 (1984) 97–103. M. Kalecki, A macrodynamic theory of business cycles, Econometrica 3 (1935) 327–344. A. Larson, The hog cycle as harmonic motion, J. Farm Econ. 46 (1964) 375–386. M. Mackey, Commodity price fluctuations: Price dependent delays and nonlinearities as explanatory factors, J. Econ. Theory 48 (1989) 497–509.
A. Matsumoto, F. Szidarovszky / Applied Mathematics and Computation 269 (2015) 699–713
713
[18] A. Matsumoto, F. Szidarovszky, Learning monopolies with delayed feedback on price expectations, Commun. Nonlinear Sci. Numer. Simul. 28 (2015) 151– 165. [19] A. Matsumoto, F. Szidarovszky, Discrete-time delay dynamics of boundedly rational monopoly, Dec. Econ. Finance 37 (2014) 53–70. [20] A. Matsumoto, F. Szidarovszky, On the Comparison of Discrete and Continuous Dynamic Systems, IERCU DP #248, Chuo University (http://www.chuou.ac.jp/research/institutes/economic/publication/discussion/pdf/discussno248.pdf?1429781458990). [21] K. Okuguchi, Expectations and Stability in Oligopoly Models, Springer-Verlag, Berlin/Heidelberg/New York, 1976. [22] K. Okuguchi, F. Szidarovszky, The Theory of Oligopoly with Multi-Product Firms, Springer-Verlag, Berlin/Heidelberg/New York, 1999. [23] M.T. Yassen, H.N. Agiza, Analysis of a duopoly game with delayed bounded rationality, Appl. Math. Comput. 138 (2003) 387–402.