# Dynamic Models of Pollution Penalties and Rewards with Time Delays

Authors: Akio Matsumoto; Ferenc Szidarovszky
Journal: Abstract and Applied Analysis
 (2020)

---

## Abstract

In cases of nonpoint pollution sources, the regulator can observe the total emission but unable to distinguish between the firms. The regulator then selects an environmental standard. If the total emission level is higher than the standard, then the firms are uniformly punished, and if lower, then uniformly awarded. This environmental regulation is added ton-firm dynamic oligopolies, and the asymptotical behavior of the corresponding dynamic systems is examined. Two particular models are considered with linear and hyperbolic price functions. Without delays, the equilibrium is always (locally) asymptotically stable. It is shown how the stability can be lost if time delays are introduced in the output quantities of the competitors as well as in the firms’ own output levels. Complete stability analysis is presented for the resulting one- and two-delay models including the derivations of stability thresholds, stability switching curves, and directions of the stability switches.


---

## Body

## 1. Introduction

Oligopolies are among the most frequently examined models in mathematical economics. The early results up to the mid-70s are summarized in Okuguchi [1], and their multiproduct extensions are discussed in Okuguchi and Szidarovszky [2]. These models and the corresponding dynamic systems are linear, the asymptotical behavior of which is simple, since local asymptotical stability implies global asymptotical stability. From the early 90s, increasing attention was given to nonlinear oligopolies. The asymptotical stability of these models was examined by a variety of concepts including linearization, Lyapunov functions, and the critical curve methods among others. Bischi et al. [3] offer a comprehensive summary of these developments. In all previous models, instantaneous information was assumed about the actions of the competitors as well as about the own output selections of the firms. However, data collection, determining appropriate actions, and their implementations need time; therefore, delayed models describe reality more accurately. More recently, Matsumoto and Szidarovszky [4] offer a collection of delayed dynamic oligopolies with a brief summary of the used mathematical methodology as well as with discussions on different types of oligopoly models.A large variety of oligopoly models consider environmental issues. The effects of different environmental regulation policies are examined by many authors including Downing and White [5], Segerson [6], Jung et al. [7], Montero [8], and Okuguchi and Szidarovszky [9, 10]. In the case of nonpoint source pollution in which neither the source nor the size of specific emissions can be observed, the government can define a cut-off value for the total emission level of the entire industry but cannot distinguish among the firms due to the information asymmetry. If the total emission level exceeds the cut-off value, then the firms are punished, otherwise rewarded. In the early stages, the existence of the Nash equilibrium was the main focus and how the governmental policy affects the total pollution level of the industry. Depending on the selected model, increased ambient charge in duopolies can lead to higher pollution levels [11], and in other cases, to lower pollution [12]. This result is generalized for n-firm Cournot oligopolies by Matsumoto et al. [13]. Dynamic models are introduced, and their asymptotical behavior was examined without and with time delays. The corresponding Bertrand models are considered and investigated in Ishikawa et al. [14]. Hyperbolic duopolies and triopolies are studied in Matsumoto et al. [15] with static and dynamic analysis. Three-stage optimum models are introduced in Matsumoto et al. [16] for Cournot duopolies without product differentiation.This paper extends and further generalizes the earlier methodology and stability results for two particular models. Linear and hyperbolic oligopolies are discussed with ambient charges or rewards, with selected cut-off pollution levels for the entire industry. As the dynamic models are very similar, we will present the complete analysis in detail for the general case including both particular models. The main results can be summarized as follows:(i)
In the case of one delay, the model with symmetric firms has a critical value of the delayτ0∗ and stability is preserved for τ<τ0∗, is lost for τ=τ0∗ via Hopf bifurcation, and cannot be regained for τ>τ0∗(ii)
In the case of two delays, stability switching curves are constructed and the directions of the stability switches are analytically shownTwo notices of our methodology are given. The first one concerns whether dynamic consideration should be global or local. The asymptotic behavior of trajectories of differential equations without and with delays is usually examined by two different methods. One is the construction and investigation of the properties of an appropriate Lyapunov function. Unfortunately, there is no general tool to create Lyapunov functions in general, which makes the application of this methodology very difficult. In our future study, we will elaborate on the possibility of applying it to different extended versions of oligopolies. The other approach is based on local linearization. This is the method we use in this paper. The other concerns whether the model should be deterministic or stochastic. In our study, the dynamic model is given and the model parameters are considered deterministic. However, in many cases, there is uncertainty in the mathematical formulation and parameter values. This uncertainty can be modeled either by stochastic (e.g., [17], or [18]) or by fuzzy [19, 20] approach. We plan to extend our study into this direction in the future.This paper is developed as follows. In Section2, the mathematical models are introduced. In Sections 3 and 4, one-delay and two-delay models are examined. In both cases, two special cases are discussed in detail: symmetric firms and general duopolies. Conclusions and further research directions are outlined in Section 5.


## 2. The Mathematical Models

An oligopoly ofn firm is considered. Let qk denote the output of firm k, Qk=∑i≠kqi the output of the rest of the industry, and Q=qk+Qk the industry output. Assume that the price function of the product by firm k is given as Pkqk,Qk. Firm k emits pollutions ekqk in connection with its production, so the total amount of pollutions is ∑k=1nekqk. The government can measure this total amount and unable to distinguish between the firms. An exogenously determined environmental standard E is selected by the government, and a θ>0 is chosen to determine the emission penalty or award for the firms. If ck is the production unit cost of firm k, then its payoff is given as
(1)πk=Pkqk,Qk−ckqk−θ∑i=1neiqi−E.Assume that at timet, each firm k has only delayed information about the outputs of the competitors, so the payoff of firm k is the following:
(2)πkt=Pkqkt,Qkt−τk−ckqkt−θekqkt+∑i≠kneiqit−τk−E.The gradient adjustment process of firmk is driven by the delay differential equation
(3)q̇kt=Kk∂πkt∂qkt,where Kk>0 is the speed of adjustment of firm k and the marginal profit is
(4)∂πkt∂qkt=∂Pk∂qktqkt+Pk−ck−θek.Let(5)gkqkt,Qkt−τk=∂πkt∂qkt.Then, equation (3) can be rewritten as
(6)q̇kt=Kkgkqkt,Qkt−τk.To linearize this equation, let(7)Uk=∂gk∂qkt,Vk=∂gk∂Qkt−τkand notice that for all i≠k,
(8)∂gk∂qit−τk=Vk,since by differentiation, the last term of gkqkt,Qkt−τk cancels out. Then, the linearized equation has the form
(9)q̇kεt=KkUkqkε+KkVk∑i≠kqiεt−τk,where qiε is the difference between qi and its equilibrium level.Before proceeding to the stability analysis of the system, two important cases are introduced.Case 1.
Assume differentiated products, when the price of the product of firmk is as follows:
(10)Pk=αk−qk−γk∑i≠kqi,where αk is the maximum price and γk represents the substitutability of the products, 0≤γk≤1. In this case,
(11)πkt=αk−qkt−γk∑i≠kqit−τk−ckqkt−θekqkt+∑i≠kneiqit−τk−E,therefore,
(12)gkqkt,Qkt−τk=αk−2qkt−γkQkt−τk−ck−θek,and so
(13)Uk=−2,Vh=−γk.
Notice that for allk,
(14)Uk<Vk≤0.Case 2.
Assume a hyperbolic oligopoly without product differentiation and common price function(15)P=αqk+Qk,where α is a positive constant. In this case, the profit function is rewritten as
(16)πkt=αqkt+Qkt−τk−ckqkt−θekqkt+∑i≠kneiqit−τk−E,consequently,
(17)gkqkt,Qkt−τk=αQkt−τkqkt+Qkt−τk2−ck−θek.
Then at the equilibrium,(18)Uk=−2αQkt−τkqkt+Qkt−τkqkt+Qkt−τk4=−2αQ¯kQ¯3,where Q¯k and Q¯ are the equilibrium levels of Qk and Q. Similarly,
(19)Vk=αqkt+Qkt−τk2−2αQkt−τkqkt+Qkt−τkqkt+Qkt−τk4=αq¯k−Q¯kQ¯3,where q¯k is the equilibrium level of qk.
Notice thatUk<0,Uk<Vk, and if there is no dominant firm, then Vk≤0. In the rest of this paper, we will assume the absence of a dominant firm. So we will assume that
(20)Uk<Vk≤0.


## 3. Single Delay Stability

To determine the characteristic equation of model (5), assume exponential solution forms, qiεt=eλtui, then substituting them into (5) gives
(21)λeλtuk=KkUkeλtuk+KkVkeλt−τk∑i≠kui.So the characteristic equation can be written as(22)φλ=detA1−λB1e−λτ1⋯B1e−λτ1B2e−λτ2A2−λ⋯B2e−λτ2··········Bne−λτnBne−λτn⋯An−λ=0,where the simplifying notations Ak=KkUk and Bk=KkVk are used. Let
(23)a=B1e−λτ1B2e−λτ2···Bne−λτn,b=11···1,D=diagA1−λ−B1e−λτ1,⋯,An−λ−Bn−λτnto have
(24)φλ=detD+abT=detDdetI+D−1abT=detD1+bTD−1a,where the identity discussed in Bischi et al. [3] is used. So
(25)φλ=∏k=1nAk−λ−Bke−λτk1+∑k=1nBke−λτkAk−λ−Bke−λτk.We have now two possibilities. Consider the first equation(26)Ak−λ−Bke−λτk=0.Without delay τk=0 and from (26), λ=Ak−Bk<0. Stability switch might occur if λ=iω with some ω>0. Then, (26) implies that
(27)Ak−iω−Bkcosωτk−isinωτk=0.Separation of the real and imaginary parts shows that(28)Bkcosωτk=Ak,Bksinωτk=ω.Adding the squares of these equations gives that(29)ω2=Bk2−Ak2<0,since Ak<Bk≤0. There is no solution for ω.Consider next the equation(30)1+∑k=1nBke−λτkAk−λ−Bke−λτk=0.Proposition 1.
System (5) is locally asymptotically stable if all roots of equation (30) are negative real numbers or complex with negative real parts.In Case1, system (5) is linear, so the asymptotical stability is global. Equation (30) is very complicated in general, so two special cases are examined.


### 3.1. Case of Symmetric Firms

AssumeAk=A,Bk=B,τk=τ, then (30) is specialized as
(31)A−λ+n−1Be−λτ=0.Without delay,λ=A+n−1B<0. Stability switching might occur with λ=iωω>0, and substituting it into (31) yields
(32)A−iω+n−1Bcosωτ−isinωτ=0,which implies that
(33)n−1Bcosωτ=−A,n−1Bsinωτ=−ω.By adding the squares of these equations, we have(34)ω2=n−12B2−A2.Ifn−1B≥A, then the right-hand side is nonpositive with no solution and without stability switch. Otherwise,
(35)ω∗=n−12B2−A2.From (33), it is clear that cosωτ<0 and sinωτ>0; furthermore, the critical values of the delays are
(36)τm∗=1ω∗cos−1−An−1B+2mπ,m≥0.The direction of stability switches can be determined by using Hopf bifurcation. Letτ be selected as the bifurcation parameter and consider λ as a function of τ:λ=λτ. By implicitly, differentiating equation (31) with respect to τ shows that
(37)λ′=−λ−Aλ1+λ−Aτ,where equation (31) is used again. At the critical value λ=iω,(38)λ′=ω2+iAω1−Aτ+iωτ·1−Aτ−iωτ1−Aτ−iωτ,the real part of which has the same sign as
(39)ω21−Aτ+ωAωτ=ω2>0.Proposition 2.
In the symmetric case, the equilibrium is locally asymptotically stable ifn−1B≥A, otherwise if τ<τ0∗. At τ=τ0∗ stability is lost via Hopf bifurcation and stability cannot be regained with larger values of τ.In the linear case,(40)Uk=−2,Vk=−γksoγk=γ.Furthermore withKk=K,(41)A=−2K,B=−Kγ,so
(42)n−1B−A=K−n−1γ+2,which is nonnegative if
(43)γ≤2n−1.One numerical example is given to confirm Proposition2. Only for analytical simplicity, γk=γ=1 (that is, no product differentiation) and Kk=K=1 are assumed. From (36), the critical value of the delay at which stability is lost is given by
(44)τ0∗K=1n−3n+1Kcos−1−2n−1.The down-slopingτ0∗Kcurve is illustrated in Figure 1(a) with n=4. It divides the parameter space of K and τ into two subregions. Stability is held for K and τ in the shaded region below the curve. For K=1, the critical value of the delay is τ0∗1≃1.029 denoted as the red point on the curve, and stability is lost. A resultant limit cycle is depicted in Figure 1(b). When n increases, the τ0∗Kcurve shifts downward, leading to a smaller stability (shaded) region.Figure 1
One delay destabilizing effect.


### 3.2. General Duopolies

In the case ofn=2, equation (22) shows that
(45)φλ=A1−λA2−λ−B1B2e−λτ1+τ2=0,which is a single-delay equation with τ=τ1+τ2:
(46)λ2−A1+A2λ+A1A2−B1B2e−λτ=0.Without delay,τ=0, and (46) becomes
(47)λ2−A1+A2λ+A1A2−B1B2=0.SinceAk>Bk for k=1,2, both the linear coefficient and the constant term are positive implying that the roots are negative real values or complex with negative real parts. Stability switch might occur if λ=iω with ω>0 when from (46) we have
(48)−ω2−iA1+A2ω+A1A2−B1B2cosωτ−isinωτ=0,implying that
(49)B1B2cosωτ=−ω2+A1A2,B1B2sinωτ=A1+A2ω.By adding the squares of these equations and arranging the terms, we have(50)ω4+A12+A22ω2+A12A22−B12B22=0.SinceAk>Bk for k=1,2, all coefficients are positive showing that no stability switch can occur.Proposition 3.
The equilibrium in a duopoly is always locally asymptotically stable with allτ1,τ2≥0.


## 3.1. Case of Symmetric Firms

AssumeAk=A,Bk=B,τk=τ, then (30) is specialized as
(31)A−λ+n−1Be−λτ=0.Without delay,λ=A+n−1B<0. Stability switching might occur with λ=iωω>0, and substituting it into (31) yields
(32)A−iω+n−1Bcosωτ−isinωτ=0,which implies that
(33)n−1Bcosωτ=−A,n−1Bsinωτ=−ω.By adding the squares of these equations, we have(34)ω2=n−12B2−A2.Ifn−1B≥A, then the right-hand side is nonpositive with no solution and without stability switch. Otherwise,
(35)ω∗=n−12B2−A2.From (33), it is clear that cosωτ<0 and sinωτ>0; furthermore, the critical values of the delays are
(36)τm∗=1ω∗cos−1−An−1B+2mπ,m≥0.The direction of stability switches can be determined by using Hopf bifurcation. Letτ be selected as the bifurcation parameter and consider λ as a function of τ:λ=λτ. By implicitly, differentiating equation (31) with respect to τ shows that
(37)λ′=−λ−Aλ1+λ−Aτ,where equation (31) is used again. At the critical value λ=iω,(38)λ′=ω2+iAω1−Aτ+iωτ·1−Aτ−iωτ1−Aτ−iωτ,the real part of which has the same sign as
(39)ω21−Aτ+ωAωτ=ω2>0.Proposition 2.
In the symmetric case, the equilibrium is locally asymptotically stable ifn−1B≥A, otherwise if τ<τ0∗. At τ=τ0∗ stability is lost via Hopf bifurcation and stability cannot be regained with larger values of τ.In the linear case,(40)Uk=−2,Vk=−γksoγk=γ.Furthermore withKk=K,(41)A=−2K,B=−Kγ,so
(42)n−1B−A=K−n−1γ+2,which is nonnegative if
(43)γ≤2n−1.One numerical example is given to confirm Proposition2. Only for analytical simplicity, γk=γ=1 (that is, no product differentiation) and Kk=K=1 are assumed. From (36), the critical value of the delay at which stability is lost is given by
(44)τ0∗K=1n−3n+1Kcos−1−2n−1.The down-slopingτ0∗Kcurve is illustrated in Figure 1(a) with n=4. It divides the parameter space of K and τ into two subregions. Stability is held for K and τ in the shaded region below the curve. For K=1, the critical value of the delay is τ0∗1≃1.029 denoted as the red point on the curve, and stability is lost. A resultant limit cycle is depicted in Figure 1(b). When n increases, the τ0∗Kcurve shifts downward, leading to a smaller stability (shaded) region.Figure 1
One delay destabilizing effect.


## 3.2. General Duopolies

In the case ofn=2, equation (22) shows that
(45)φλ=A1−λA2−λ−B1B2e−λτ1+τ2=0,which is a single-delay equation with τ=τ1+τ2:
(46)λ2−A1+A2λ+A1A2−B1B2e−λτ=0.Without delay,τ=0, and (46) becomes
(47)λ2−A1+A2λ+A1A2−B1B2=0.SinceAk>Bk for k=1,2, both the linear coefficient and the constant term are positive implying that the roots are negative real values or complex with negative real parts. Stability switch might occur if λ=iω with ω>0 when from (46) we have
(48)−ω2−iA1+A2ω+A1A2−B1B2cosωτ−isinωτ=0,implying that
(49)B1B2cosωτ=−ω2+A1A2,B1B2sinωτ=A1+A2ω.By adding the squares of these equations and arranging the terms, we have(50)ω4+A12+A22ω2+A12A22−B12B22=0.SinceAk>Bk for k=1,2, all coefficients are positive showing that no stability switch can occur.Proposition 3.
The equilibrium in a duopoly is always locally asymptotically stable with allτ1,τ2≥0.


## 4. Two-Delay Stability

Assume next that firmk has a delay τ1k in its output and delay τ2k in the outputs of its competitors. Then, the dynamic equation (9) modifies as follows:
(51)q̇kεt=Akqkεt−τ1k+Bk∑i≠kqiε1−τ2k.Similarly to equations (22) and (25), the characteristic equation can be derived as
(52)φλ=∏k=1nAke−λτ1k−λ−Bke−λτ2k1+∑k=1nBke−λτ2kAke−λτ1k−λ−Bke−λτ2k=0.Proposition 4.
The equilibrium of system (51) is locally asymptotically stable if all roots of (52) are negative reals or complex with negative real parts.Similarly to the previous model, two special cases will be reexamined.


### 4.1. Case of Symmetric Firms

Assume now thatAk=A,Bk=B,τ1k=τ1 and τ2k=τ2. From (52), we have to consider two cases. In the first case, we examine the equation
(53)Ae−λτ1−λ−Be−λτ2=0which can be rewritten as
(54)1+a1λe−λτ1+a2λe−λτ2=0with
(55)a1λ=−Aλ,a2λ=Bλ.We will apply the method introduced by Gu et al. (2005) and discussed in detail in Matsumoto and Szidarovszky [4]. Notice that
(56)a1iω=iAω,a2iω=−iBω,a1iω=−Aω,a2iω=−Bω,arga1iω=3π2,arga2iω=π2.The range ofω is determined by conditions,
(57)a1iω+a2iω≥1,−1≤a1iω−a2iω≤1which simplify in our case as
(58)−Aω−Bω≥1orω≤−A+B,−1≤−Aω+Bω≤1orω≥B−A.Soω runs through interval B−A,−A+B. Moreover, by the law of cosine,
(59)θ1ω=cos−11+a1iω2−a2iω22a1iω=cos−1ω2+A2−B2−2Aω,θ2ω=cos−11+a2iω2−a1iω22a2iω=cos−1ω2+B2−A2−2Bω.The stability switching curves are given by pairsτ1±k,τ2∓m with
(60)τ1±kω=1ω3π2+2k−1π±θ1ω,τ2∓mω=1ωπ2+2m−1π∓θ2ω.In the second case, we consider the other equation(61)Ae−λτ1−λ+n−1Be−λτ2=0,which can be rewritten as
(62)1+a1λe−λτ1+a2λe−λτ2=0with
(63)a1λ=−Aλ,a2λ=−n−1Bλ.Then,(64)a1iω=iAω,a2iω=in−1Bω,a1iω=−Aω,a2iω=−n−1Bω,arga1iω=arga2iω=3π2.The range ofω is determined again based on conditions (57), which are the following in this case:
(65)−Aω−n−1Bω≥1orω≤−A+n−1B,−1≤−Aω+n−1Bω≤1orω≥A−n−1B,so a range of ω is in the interval,
(66)A−n−1B,−A+n−1B.By the law of cosines,(67)θ~1ω=cos−1ω2+A2−n−12B2−2Aω,θ~2ω=cos−1ω2+n−12B2−A2−2n−1Bω.The stability switching curves are given by pairsτ1±k,τ2∓m with
(68)τ~1±kω=1ω3π2+2k−1π±θ~1ω,τ~2∓mω=1ω3π2+2m−1π∓θ~2ω.Proposition 5.
The stability switching curves are formed by pointsτ1±kω,τ2∓mω with k=m=0 when ω runs through interval B−A,−A+B and points τ~1±kω,τ~2∓mω with k=m=0 when ω runs through A−n−1B,−A+n−1B.We numerically confirm Proposition5. To this end, we take K=1 and n=4 in the second case. The stability switching curve in the first case consists of the pairs of τ1±0ω,τ2∓0ω and illustrated as the dotted red-blue curve in Figure 2(a). It can be verified by equation (53) that the model with symmetric firms is stable without delays (i.e., τ1=τ2=0). Hence, stability is preserved for τ1,τ2 to the left of the curve. On the other hand, the stability switching curve in the second case, the locus of τ~1±0ω,τ~2∓0ω, is illustrated as the solid red-blue curve. Due to equation (61), the model in the second case is also stable without delays, implying that stability is preserved for τ1,τ2 below the curve. The lower-left corner surrounded by the dotted line is enlarged in Figure 2(b). The model is stable in the region including the origin and bounded by the solid red-blue curve and the dotted blue curve.Figure 2
Two-delay destabilizing effect.The directions of stability switching in the first case can be assessed by computing the following expressions:(69)a1ωe−iωτ1=iAωcosωτ1−isinωτ1,a2ωe−iωτ2=−iBωcosωτ2−isinωτ2with real and imaginary parts,
(70)R1=Rea1ωe−iωτ1=Aωsinωτ1,I1=Ima1ωe−iωτ1=Aωcosωτ1,R2=Rea2ωe−iωτ2=−Bωsinωτ2,I2=Ima2ωe−iωτ2=−Bωcosωτ2,and finally,
(71)S1=R2I1−R1I2=ABω2sinωτ1cosωτ2−cosωτ1sinωτ2,which has the same sign as sinωτ1−τ2.The directions of stability switches in the second case can be determined similarly to the first case. Notice that(72)a1ωe−iωτ1=iAωcosωτ1−isinωτ1,a2ωe−iωτ2=in−1Bωcosωτ2−isinωτ2with real and imaginary parts,
(73)R1=Rea1ωe−iωτ1=Aωsinωτ1,I1=Ima1ωe−iωτ1=Aωcosωτ1,R2=Rea2ωe−iωτ2=n−1Bωsinωτ2,I2=Ima2ωe−iωτ2=n−1Bωcosωτ2,and therefore,
(74)S2=R2I1−R1I2=n−1ABω2sinωτ2cosωτ1−cosωτ2sinωτ1,which has the same sign as sinωτ2−τ1. Then, we have the following:Proposition 6.
Letτ1,τ2 be a point on the stability switching curve and assume that the curve is crossed at this point from right to left when we are looking forward to increasing values of ω on the curve. If S1orS2 is positive, then at least one pair of eigenvalues changes the sign of the real part from negative to positive. If S1orS2 is negative, then the sign change is in the opposite direction.


### 4.2. General Duopolies

In the case ofn=2 from (52), we have
(75)φλ=A1e−λτ11−λ−B1e−λτ21A2e−λτ12−λ−B2e−λτ22+B1e−λτ21A2e−λτ12−λ−B2e−λτ22+B2e−λτ22A1e−λτ11−λ−B1e−λτ21=0,which can be simplified as
(76)A1e−λτ11−λA2e−λτ12−λ−B1B2e−λτ21+τ22=0.This equation is analytically intractable since it has four delays,τ11,τ12,τ11+τ12,τ21+τ22. Therefore, we make the following simplifying assumption:
(77)τ11=τ12=τ21=τ22=τ,when (76) becomes
(78)λ2−A1+A2λe−λτ+A1A2−B1B2e−2λτ=0.By multiplying both sides byeλτ, we have
(79)λ2eλτ−A1+A2λ+A1A2−B1B2e−λτ=0.Without delay, a quadratic equation is obtained byλ(80)λ2−A1+A2λ+A1A2−B1B2=0.The linear coefficient and constant term are positive; the roots are negative real values since the discriminant is positive. Stability switch might occur ifλ=iω with ω>0, then from (79),
(81)−ω2cosωτ+isinωτ−iωA1+A2+A1A2−B1B2cosωτ−isinωτ=0.By separating the real and imaginary parts, we have(82)−ω2+A1A2−B1B2cosωτ=0,−ω2−A1A2−B1B2sinωτ=ωA1+A2.We have to consider now two possibilities from the first equation of (82).
(83)icosωτ≠0,then,
(84)ω2=A1A2−B1B2,and from the second equation of (82),
(85)−2A1A2−B1B2sinωτ=ωA1+A2.However,(86)ω2A1+A22−4A1A2−B1B22=A1A2−B1B2A12+A22+2A1A2−4A1A2−B1B22=A1A2−B1B2A1−A22+4B1B2,where the first term is positive and the second is being nonnegative. If it is positive, (85) has no solution. If it is zero, then sinωτ=+1 implying that cosωτ=0, which is a contradiction. The other possibility is
(87)iicosωτ=0,then, from the second equation of (82), sinωτ=1 and
(88)ω2+A1+A2ω+A1A2−B1B2=0.The discriminant is(89)Δ=A1−A22+4B1B2≥0,furthermore,
(90)A1+A22−Δ=4A1A2−B1B2>0,so both roots are positive,
(91)ω±=−A1+A2±Δ2>0.Notice thatΔ=0 if A1=A2 and one of B1 and B2 equals to zero. Then, ω±=−A1=−A2. From the first equation of (82), the critical value of the delay is
(92)τm±=1ω±π2+2mπform=0,1,2,⋯The direction of the stability switches are assessed by consideringλ as the function of the bifurcation parameter, λ=λτ and implicitly differentiating equation (79):
(93)2λλ′eλτ+λ2eλτλ′τ+λ−λ′A1+A2+A1A2−B1B2e−λτ−λ′τ−λ=0.The multiplier ofλ′ and the constant term are
(94)2λeλτ+λ2τeλτ−A1−A2−A1A2−B1B2e−λττ,λ3eλτ−A1A2−B1B2e−λτλ,respectively, so
(95)λ′−1=−2λeλτ+A1+A2λ3eλτ−A1A2−B1B2λe−λτ−τλ.Notice that at the critical values ofτm±,cosωτ=0 and sinωτ=1 implying that
(96)eiωτ=cosωτ+isinωτ=i,e−iωτ=cosωτ−isinωτ=−i.Atλ=iω, we are interested in the real part of (95), where the second term is pure complex. Therefore,
(97)Reλ′−1=Re−2λeλτ+A1+A2λ3eλτ−A1A2−B1B2λe−λτ=2ω+A1+A2ωω2−A1A2−B1B2.Substituting (91) with numerator ±Δ into the denominator presents
(98)ω42A1−A22+8B1B2∓2A1+A2Δ=ω2ΔΔ∓A1+A2.We know however thatA1+A2>Δ, so this expression is positive at ω+ and negative at ω−. Hence, Reλ′−1 is always positive.Proposition 7.
The equilibrium is locally asymptotically stable forτ>τ0+, stability is lost at τ=τ0+ via Hopf bifurcation, and stability cannot be regained with larger values of τ.


## 4.1. Case of Symmetric Firms

Assume now thatAk=A,Bk=B,τ1k=τ1 and τ2k=τ2. From (52), we have to consider two cases. In the first case, we examine the equation
(53)Ae−λτ1−λ−Be−λτ2=0which can be rewritten as
(54)1+a1λe−λτ1+a2λe−λτ2=0with
(55)a1λ=−Aλ,a2λ=Bλ.We will apply the method introduced by Gu et al. (2005) and discussed in detail in Matsumoto and Szidarovszky [4]. Notice that
(56)a1iω=iAω,a2iω=−iBω,a1iω=−Aω,a2iω=−Bω,arga1iω=3π2,arga2iω=π2.The range ofω is determined by conditions,
(57)a1iω+a2iω≥1,−1≤a1iω−a2iω≤1which simplify in our case as
(58)−Aω−Bω≥1orω≤−A+B,−1≤−Aω+Bω≤1orω≥B−A.Soω runs through interval B−A,−A+B. Moreover, by the law of cosine,
(59)θ1ω=cos−11+a1iω2−a2iω22a1iω=cos−1ω2+A2−B2−2Aω,θ2ω=cos−11+a2iω2−a1iω22a2iω=cos−1ω2+B2−A2−2Bω.The stability switching curves are given by pairsτ1±k,τ2∓m with
(60)τ1±kω=1ω3π2+2k−1π±θ1ω,τ2∓mω=1ωπ2+2m−1π∓θ2ω.In the second case, we consider the other equation(61)Ae−λτ1−λ+n−1Be−λτ2=0,which can be rewritten as
(62)1+a1λe−λτ1+a2λe−λτ2=0with
(63)a1λ=−Aλ,a2λ=−n−1Bλ.Then,(64)a1iω=iAω,a2iω=in−1Bω,a1iω=−Aω,a2iω=−n−1Bω,arga1iω=arga2iω=3π2.The range ofω is determined again based on conditions (57), which are the following in this case:
(65)−Aω−n−1Bω≥1orω≤−A+n−1B,−1≤−Aω+n−1Bω≤1orω≥A−n−1B,so a range of ω is in the interval,
(66)A−n−1B,−A+n−1B.By the law of cosines,(67)θ~1ω=cos−1ω2+A2−n−12B2−2Aω,θ~2ω=cos−1ω2+n−12B2−A2−2n−1Bω.The stability switching curves are given by pairsτ1±k,τ2∓m with
(68)τ~1±kω=1ω3π2+2k−1π±θ~1ω,τ~2∓mω=1ω3π2+2m−1π∓θ~2ω.Proposition 5.
The stability switching curves are formed by pointsτ1±kω,τ2∓mω with k=m=0 when ω runs through interval B−A,−A+B and points τ~1±kω,τ~2∓mω with k=m=0 when ω runs through A−n−1B,−A+n−1B.We numerically confirm Proposition5. To this end, we take K=1 and n=4 in the second case. The stability switching curve in the first case consists of the pairs of τ1±0ω,τ2∓0ω and illustrated as the dotted red-blue curve in Figure 2(a). It can be verified by equation (53) that the model with symmetric firms is stable without delays (i.e., τ1=τ2=0). Hence, stability is preserved for τ1,τ2 to the left of the curve. On the other hand, the stability switching curve in the second case, the locus of τ~1±0ω,τ~2∓0ω, is illustrated as the solid red-blue curve. Due to equation (61), the model in the second case is also stable without delays, implying that stability is preserved for τ1,τ2 below the curve. The lower-left corner surrounded by the dotted line is enlarged in Figure 2(b). The model is stable in the region including the origin and bounded by the solid red-blue curve and the dotted blue curve.Figure 2
Two-delay destabilizing effect.The directions of stability switching in the first case can be assessed by computing the following expressions:(69)a1ωe−iωτ1=iAωcosωτ1−isinωτ1,a2ωe−iωτ2=−iBωcosωτ2−isinωτ2with real and imaginary parts,
(70)R1=Rea1ωe−iωτ1=Aωsinωτ1,I1=Ima1ωe−iωτ1=Aωcosωτ1,R2=Rea2ωe−iωτ2=−Bωsinωτ2,I2=Ima2ωe−iωτ2=−Bωcosωτ2,and finally,
(71)S1=R2I1−R1I2=ABω2sinωτ1cosωτ2−cosωτ1sinωτ2,which has the same sign as sinωτ1−τ2.The directions of stability switches in the second case can be determined similarly to the first case. Notice that(72)a1ωe−iωτ1=iAωcosωτ1−isinωτ1,a2ωe−iωτ2=in−1Bωcosωτ2−isinωτ2with real and imaginary parts,
(73)R1=Rea1ωe−iωτ1=Aωsinωτ1,I1=Ima1ωe−iωτ1=Aωcosωτ1,R2=Rea2ωe−iωτ2=n−1Bωsinωτ2,I2=Ima2ωe−iωτ2=n−1Bωcosωτ2,and therefore,
(74)S2=R2I1−R1I2=n−1ABω2sinωτ2cosωτ1−cosωτ2sinωτ1,which has the same sign as sinωτ2−τ1. Then, we have the following:Proposition 6.
Letτ1,τ2 be a point on the stability switching curve and assume that the curve is crossed at this point from right to left when we are looking forward to increasing values of ω on the curve. If S1orS2 is positive, then at least one pair of eigenvalues changes the sign of the real part from negative to positive. If S1orS2 is negative, then the sign change is in the opposite direction.


## 4.2. General Duopolies

In the case ofn=2 from (52), we have
(75)φλ=A1e−λτ11−λ−B1e−λτ21A2e−λτ12−λ−B2e−λτ22+B1e−λτ21A2e−λτ12−λ−B2e−λτ22+B2e−λτ22A1e−λτ11−λ−B1e−λτ21=0,which can be simplified as
(76)A1e−λτ11−λA2e−λτ12−λ−B1B2e−λτ21+τ22=0.This equation is analytically intractable since it has four delays,τ11,τ12,τ11+τ12,τ21+τ22. Therefore, we make the following simplifying assumption:
(77)τ11=τ12=τ21=τ22=τ,when (76) becomes
(78)λ2−A1+A2λe−λτ+A1A2−B1B2e−2λτ=0.By multiplying both sides byeλτ, we have
(79)λ2eλτ−A1+A2λ+A1A2−B1B2e−λτ=0.Without delay, a quadratic equation is obtained byλ(80)λ2−A1+A2λ+A1A2−B1B2=0.The linear coefficient and constant term are positive; the roots are negative real values since the discriminant is positive. Stability switch might occur ifλ=iω with ω>0, then from (79),
(81)−ω2cosωτ+isinωτ−iωA1+A2+A1A2−B1B2cosωτ−isinωτ=0.By separating the real and imaginary parts, we have(82)−ω2+A1A2−B1B2cosωτ=0,−ω2−A1A2−B1B2sinωτ=ωA1+A2.We have to consider now two possibilities from the first equation of (82).
(83)icosωτ≠0,then,
(84)ω2=A1A2−B1B2,and from the second equation of (82),
(85)−2A1A2−B1B2sinωτ=ωA1+A2.However,(86)ω2A1+A22−4A1A2−B1B22=A1A2−B1B2A12+A22+2A1A2−4A1A2−B1B22=A1A2−B1B2A1−A22+4B1B2,where the first term is positive and the second is being nonnegative. If it is positive, (85) has no solution. If it is zero, then sinωτ=+1 implying that cosωτ=0, which is a contradiction. The other possibility is
(87)iicosωτ=0,then, from the second equation of (82), sinωτ=1 and
(88)ω2+A1+A2ω+A1A2−B1B2=0.The discriminant is(89)Δ=A1−A22+4B1B2≥0,furthermore,
(90)A1+A22−Δ=4A1A2−B1B2>0,so both roots are positive,
(91)ω±=−A1+A2±Δ2>0.Notice thatΔ=0 if A1=A2 and one of B1 and B2 equals to zero. Then, ω±=−A1=−A2. From the first equation of (82), the critical value of the delay is
(92)τm±=1ω±π2+2mπform=0,1,2,⋯The direction of the stability switches are assessed by consideringλ as the function of the bifurcation parameter, λ=λτ and implicitly differentiating equation (79):
(93)2λλ′eλτ+λ2eλτλ′τ+λ−λ′A1+A2+A1A2−B1B2e−λτ−λ′τ−λ=0.The multiplier ofλ′ and the constant term are
(94)2λeλτ+λ2τeλτ−A1−A2−A1A2−B1B2e−λττ,λ3eλτ−A1A2−B1B2e−λτλ,respectively, so
(95)λ′−1=−2λeλτ+A1+A2λ3eλτ−A1A2−B1B2λe−λτ−τλ.Notice that at the critical values ofτm±,cosωτ=0 and sinωτ=1 implying that
(96)eiωτ=cosωτ+isinωτ=i,e−iωτ=cosωτ−isinωτ=−i.Atλ=iω, we are interested in the real part of (95), where the second term is pure complex. Therefore,
(97)Reλ′−1=Re−2λeλτ+A1+A2λ3eλτ−A1A2−B1B2λe−λτ=2ω+A1+A2ωω2−A1A2−B1B2.Substituting (91) with numerator ±Δ into the denominator presents
(98)ω42A1−A22+8B1B2∓2A1+A2Δ=ω2ΔΔ∓A1+A2.We know however thatA1+A2>Δ, so this expression is positive at ω+ and negative at ω−. Hence, Reλ′−1 is always positive.Proposition 7.
The equilibrium is locally asymptotically stable forτ>τ0+, stability is lost at τ=τ0+ via Hopf bifurcation, and stability cannot be regained with larger values of τ.


## 5. Conclusions

Environmental regulations were added to the classicaln-firm Cournot model. In cases of nonpoint source pollution, the regulator can measure only the total emission level without knowing the individual emissions of the firms. Therefore in the regulation, the firms are uniformly punished if the total emission is higher than a regulator selected standard, and awarded otherwise. In the dynamic extensions, we considered three cases. First, no delays were introduced about the information on the outputs of all firms; second, delayed data were assumed about the output levels of the competitors; and third, additional delays were added in the firms’ own output levels. The stability analysis was conducted under general conditions that are satisfied in cases of linear and hyperbolic price functions. Models without delays, with a single delay, and with two delays were analyzed in detail. It was demonstrated how the stability of the no-delay models can be lost by introducing delays. The stability thresholds, stability switching curves, and directions of stability switches were analytically derived.It will be an interesting project to extend the results of this paper to more general cases including nonlinear Cournot models, multiproduct, labor-managed oligopolies, and rent-seeking games among others. This will be the subject of our continued research project. In addition, the possible uncertainty of model parameters will be considered in the extended oligopoly models with a fuzzy and stochastic approach.



---

*Source: 3162634-2020-07-06.xml*