Dynamic Models of Pollution Penalties and Rewards with Time Delays

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 to n-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 oneand two-delay models including the derivations of stability thresholds, stability switching curves, and directions of the stability switches.


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 shown Two 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 Section 2, 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.

The Mathematical Models
An oligopoly of n firm is considered. Let q k denote the output of firm k, Q k = ∑ i≠k q i the output of the rest of the industry, and Q = q k + Q k the industry output. Assume that the price function of the product by firm k is given as P k ðq k , Q k Þ. Firm k emits pollutions e k q k in connection with its production, so the total amount of pollutions is ∑ n k=1 e k q k . 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 c k is the production unit cost of firm k, then its payoff is given as Assume that at time t, each firm k has only delayed information about the outputs of the competitors, so the payoff of firm k is the following: The gradient adjustment process of firm k is driven by the delay differential equation where K k > 0 is the speed of adjustment of firm k and the marginal profit is Let Then, equation (3) can be rewritten as To linearize this equation, let and notice that for all i ≠ k, Abstract and Applied Analysis since by differentiation, the last term of g k ðq k ðtÞ, Q k ðt − τ k ÞÞ cancels out. Then, the linearized equation has the form where q iε is the difference between q i 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 firm k is as follows: where α k is the maximum price and γ k represents the substitutability of the products, 0 ≤ γ k ≤ 1. In this case, therefore, and so Notice that for all k, Case 2. Assume a hyperbolic oligopoly without product differentiation and common price function where α is a positive constant. In this case, the profit function is rewritten as consequently, Then at the equilibrium, where Q k and Q are the equilibrium levels of Q k and Q. Similarly, where q k is the equilibrium level of q k .
Notice that U k < 0, U k < V k , and if there is no dominant firm, then V k ≤ 0. In the rest of this paper, we will assume the absence of a dominant firm. So we will assume that

Single Delay Stability
To determine the characteristic equation of model (5), assume exponential solution forms, q iε ðtÞ = e λt u i , then substituting them into (5) gives So the characteristic equation can be written as where the simplifying notations A k = K k U k and B k = K k V k are used. Let 3 Abstract and Applied Analysis where the identity discussed in Bischi et al. [3] is used. So We have now two possibilities. Consider the first equation Without delay τ k = 0 and from (26), λ = A k − B k < 0. Stability switch might occur if λ = iω with some ω > 0. Then, (26) implies that Separation of the real and imaginary parts shows that Adding the squares of these equations gives that There is no solution for ω.
Consider next the equation (5) is locally asymptotically stable if all roots of equation (30) are negative real numbers or complex with negative real parts.
In Case 1, system (5) is linear, so the asymptotical stability is global. Equation (30) is very complicated in general, so two special cases are examined.

Case of Symmetric Firms. Assume
Without delay, λ = A + ðn − 1ÞB < 0. Stability switching might occur with λ = iωðω > 0Þ, and substituting it into (31) yields which implies that By adding the squares of these equations, we have If ðn − 1ÞB ≥ A, then the right-hand side is nonpositive with no solution and without stability switch. Otherwise, From (33), it is clear that cos ωτ < 0 and sin ωτ > 0; furthermore, the critical values of the delays are 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 where equation (31) is used again. At the critical value λ = iω, Abstract and Applied Analysis the real part of which has the same sign as Proposition 2. In the symmetric case, the equilibrium is locally asymptotically stable if ðn − 1ÞB ≥ 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, Furthermore with K k = K, which is nonnegative if One numerical example is given to confirm Proposition 2. Only for analytical simplicity, γ k = γ = 1 (that is, no product differentiation) and K k = K = 1 are assumed. From (36), the critical value of the delay at which stability is lost is given by The down-sloping τ * 0 ðKÞ curve 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 ðKÞ curve shifts downward, leading to a smaller stability (shaded) region.

General Duopolies.
In the case of n = 2, equation (22) shows that which is a single-delay equation with τ = τ 1 + τ 2 : Without delay, τ = 0, and (46) becomes Since jA k j > jB k j 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 implying that By adding the squares of these equations and arranging the terms, we have Since jA k j > jB k j 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:

Two-Delay Stability
Assume next that firm k has a delay τ k 1 in its output and delay τ k 2 in the outputs of its competitors. Then, the dynamic equation (9) modifies as follows: Similarly to equations (22) and (25), the characteristic equation can be derived as ð52Þ 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.

Case of Symmetric Firms.
Assume now that A k = A, B k = B,τ k 1 = τ 1 and τ k 2 = τ 2 . From (52), we have to consider two cases. In the first case, we examine the equation 5 Abstract and Applied Analysis which can be rewritten as We will apply the method introduced by Gu et al. (2005) and discussed in detail in Matsumoto and Szidarovszky [4]. Notice that The range of ω is determined by conditions, which simplify in our case as So ω runs through interval ½B − A,−ðA + BÞ. Moreover, by the law of cosine, The stability switching curves are given by pairs ðτ ±k 1 , τ ∓m 2 Þ with In the second case, we consider the other equation which can be rewritten as with   Abstract and Applied Analysis Then, The range of ω is determined again based on conditions (57), which are the following in this case: so a range of ω is in the interval, By the law of cosines, The stability switching curves are given by pairs ðτ ±k 1 , τ ∓m 2 Þ withτ We numerically confirm Proposition 5. 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 ðτ ±0 1 ðωÞ, τ ∓0 2 ðωÞÞ 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 ðτ ±0 1 ðωÞ,τ ∓0 2 ðωÞÞ, 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.
The directions of stability switching in the first case can be assessed by computing the following expressions: with real and imaginary parts, and finally, 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 with real and imaginary parts, sin ωτ 2 , 7 Abstract and Applied Analysis and therefore, 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 S 1 ðor S 2 Þ is positive, then at least one pair of eigenvalues changes the sign of the real part from negative to positive. If S 1 ðor S 2 Þ is negative, then the sign change is in the opposite direction. This equation is analytically intractable since it has four delays, τ 1 1 , τ 2 1 ,τ 1 1 + τ 2 1 , τ 1 2 + τ 2 2 : Therefore, we make the following simplifying assumption: when (76) becomes By multiplying both sides by e λτ , we have Without delay, a quadratic equation is obtained by λ 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), By separating the real and imaginary parts, we have We have to consider now two possibilities from the first equation of (82). then, and from the second equation of (82),