**3. Best Reply Dynamics**

Dynamic interpretation of the oligopoly model depends on how to define a learning process on how each firm observes its competitors' choices. Theocharis (1960) constructs the best reply dynamics with naive expectations in discrete time scales,

$$x\_i(t+1) = a - \beta \sum\_{j \neq i}^n x\_j(t)$$

where the adjustment to the optimal output in each period is perfect. His provocative result shows that the stability of the Cournot equilibrium is determined only by the number of the firms in an industry as mentioned in the Introduction. McManus and Quandt (1961) makes two reasonable modifications of Theocharis' assumptions: the discrete-time scales are replaced with continuous-time scales and the imperfect adjustment assumption is adopted in which the direction of output change is proportional to the discrepancy between the optimal and actual values,

$$\dot{x}\_i(t) = k\_i \left[ a - \beta \sum\_{j \neq i}^n x\_j(t) - x\_i(t) \right] \text{ with } k\_i > 0.$$

It is demonstrated that the Cournot equilibrium is always stable when the adjustment speeds are the same (i.e., *ki* = *k*). Their result is in sharp contrast to Theocharis' result. We also note that this result remains true if all adjustment speeds are positive.

#### *3.1. Stability Switching*

In this study, we move one step forward from the McManus and Quandt model and introduce implementation delays (i.e., *τ*<sup>1</sup> > 0) on the firm's own production and information delays (*τ*<sup>2</sup> > 0) on the competitors' productions,

$$\dot{\mathbf{x}}\_i(t) = k \left[ \mathbf{a} - \beta \sum\_{j \neq i}^n \mathbf{x}\_j(t - \tau\_2) - \mathbf{x}\_i(t - \tau\_1) \right] \text{ for } i = 1, 2, \dots, n. \tag{1}$$

*Mathematics* **2020**, *8*, 1615

Notice that dynamic system (1) has the Cournot equilibrium as the steady-state and its homogeneous part is

$$\dot{\mathbf{x}}\_{i}(t) = k \left[ -\mathbf{x}\_{i}(t-\tau\_{1}) - \beta \sum\_{j \neq i}^{n} \mathbf{x}\_{j}(t-\tau\_{2}) \right] \text{ for } i = 1, 2, \dots, n. \tag{2}$$

The characteristic equation is

$$g(\boldsymbol{\lambda}) = \begin{pmatrix} \lambda + k e^{-\lambda \tau\_1} & k \beta e^{-\lambda \tau\_2} & \cdots & k \beta e^{-\lambda \tau\_2} \\ k \beta e^{-\lambda \tau\_2} & \lambda + k e^{-\lambda \tau\_1} & \cdots & k \beta e^{-\lambda \tau\_2} \\ \cdot & \cdot & \cdots & \cdots \\ k \beta e^{-\lambda \tau\_2} & k \beta e^{-\lambda \tau\_2} & \cdots & \lambda + k e^{-\lambda \tau\_1} \end{pmatrix} = 0.$$

With new notation,

$$\mathbf{D} = \operatorname{diag} \left( \lambda + k e^{-\lambda \tau\_1} - k \beta e^{-\lambda \tau\_2}, \dots, \lambda + k e^{-\lambda \tau\_1} - k \beta e^{-\lambda \tau\_2} \right)\_{(n,n)}$$

$$\mathbf{a} = \left( k \beta e^{-\lambda \tau\_2} \right)\_{(n,1)} \text{ and } \mathbf{b} = (1)\_{(n,1)'}$$

the characteristic equation can be written as

$$\begin{aligned} \varphi(\lambda) &= \det\left(D + ab^{\top}\right), \\ &= \det D \det\left(I + D^{-1}ab^{\top}\right), \\ &= \det D \left[1 + b^{\top}D^{-1}a\right]. \end{aligned}$$

Hence

$$\begin{split} \varrho(\lambda) &= \left[\lambda + k e^{-\lambda \tau\_1} - k \beta e^{-\lambda \tau\_2}\right]^n \left[1 + \frac{nk \beta e^{-\lambda \tau\_2}}{\lambda + k e^{-\lambda \tau\_1} - k \beta e^{-\lambda \tau\_2}}\right] \\ &= \left(\lambda + k e^{-\lambda \tau\_1} - k \beta e^{-\lambda \tau\_2}\right)^{n-1} \left(\lambda + k e^{-\lambda \tau\_1} + k \beta (n-1) e^{-\lambda \tau\_2}\right) .\end{split}$$

It follows that we have two possibilities to solve *ϕ*(*λ*) = 0,

$$\begin{aligned} \text{(i) } \lambda + ke^{-\lambda \tau\_1} - k\beta e^{-\lambda \tau\_2} &= 0, \\ \text{(ii) } \lambda + ke^{-\lambda \tau\_1} + k\beta (n - 1)e^{-\lambda \tau\_2} &= 0. \end{aligned}$$

Without delays *τ*<sup>1</sup> = *τ*<sup>2</sup> = 0, the eigenvalues are negative,

$$
\lambda\_1 = -k(1 - \beta) < 0 \text{ and } \lambda\_2 = -k\left[1 + \beta(n - 1)\right] < 0,
$$

implying that the equilibrium is asymptotically stable.

For positive delays, we follow the method discussed in Matsumoto and Szidarovszky [13] based on Gu et al. [14]. Consider equation (i) first. As *λ* = 0 does not solve equation (i), it can be rewritten as

$$1 + a\_1(\lambda)e^{-\lambda\tau\_1} + a\_2(\lambda)e^{-\lambda\tau\_2} = 0\tag{3}$$

where

$$a\_1(\lambda) = \frac{k}{\lambda} \text{ and } a\_2(\lambda) = -\frac{k\beta}{\lambda}.$$

Equation (3) must have a pair of pure imaginary solutions when a stability switch occurs. Hence let *λ* = *iω* with *ω* > 0 (It is possible to take its conjugate with *ω* < 0. Even so, we can arrive at the same result.) and we may consider the three terms in (3) as three vectors in the complex plane with the magnitudes 1, |*a*1(*iω*)| and |*a*2(*iω*)|, respectively. Equation (3) means that if we put these vectors head to tail, they form a triangle with the internal angles *θ*<sup>1</sup> and *θ*<sup>2</sup> as illustrated in Figure 1.

**Figure 1.** Triangle conditions.

These vectors form a triangle if and only if the sum of the lengths of any two adjacent line segments is not shorter than the length of the remaining line segment:

$$|a\_1(i\omega)| + |a\_2(i\omega)| \ge 1$$

and

$$-1 \le |a\_1(i\omega)| - |a\_2(i\omega)| \le 1.$$

For *λ* = *iω*,

$$a\_1(i\omega) = -i\frac{k}{\omega} \text{ and } a\_2(\lambda) = i\frac{k\beta}{\omega}$$

where the absolute values are

$$|a\_1(i\omega)| = \frac{k}{\omega} \text{ and } |a\_2(i\omega)| = \frac{k\beta}{\omega}$$

and the arguments are

$$\arg\left[a\_1(i\omega)\right] = \frac{3\pi}{2} \text{ and } \arg\left[a\_2(i\omega)\right] = \frac{\pi}{2}.$$

From the triangle conditions, we have the interval of *ω* for which *λ* = *iω* can be a solution of equation (i) for some *τ*<sup>1</sup> and *τ*2,

$$
\omega \in I = \left[\frac{1}{2}k, \frac{3}{2}k\right].
$$

The internal angles of *θ*<sup>1</sup> and *θ*<sup>2</sup> are calculated by the law of cosine as

$$\theta\_1(\omega) = \cos^{-1}\left[\frac{4\omega^2 + 3k^2}{8k\omega}\right]$$

and

$$\theta\_2(\omega) = \cos^{-1}\left[\frac{4\omega^2 - 3k^2}{4k\omega}\right].$$

For any *ω* ∈ *I*, we may find all pairs of (*τ*1, *τ*2) satisfying (3) as follows:

$$
\pi\_1^{\pm}(\omega, \ell\_1) = \frac{1}{\omega} \left[ \frac{3}{2} \pi + (2\ell\_1 - 1)\,\pi \pm \theta\_1(\omega) \right] \tag{4}
$$

and

$$
\pi\_2^{\frac{\pi}{\pi}}(\omega, \ell\_2) = \frac{1}{\omega} \left[ \frac{1}{2}\pi + (2\ell\_2 - 1)\,\pi \mp \theta\_2(\omega) \right].\tag{5}
$$

Since a symmetric triangle can be formed below the horizontal axis in Figure 1, four inner angles are defined, ±*θ*1(*ω*) and ∓*θ*2(*ω*) (double-sign correspondence). By the definitions of the interior angles, we have the followings:

$$\arg \left( a\_1(\omega) e^{-i\omega \tau\_1} \right) + 2\ell\_1 \pi \pm \theta\_1(\omega) = \pi i$$

and

$$\arg \left( a\_2(\omega) e^{-i\omega \tau\_2} \right) + 2\ell\_2 \pi \mp \theta\_2(\omega) = \pi$$

for -<sup>1</sup> = 0, 1, 2, ... and -<sup>2</sup> = 0, 1, 2, ... Solving these equations for *τ*<sup>1</sup> and *τ*<sup>2</sup> yields (4) and (5). So we have two sets of line segments,

$$\mathcal{C}\_1^+(\ell\_1, \ell\_2) = \left\{ \left( \pi\_1^+(\omega, \ell\_1), \pi\_2^-(\omega, \ell\_2) \right) \mid \omega \in I, \left( \ell\_1, \ell\_2 \right) \in \mathbf{Z} \right\}$$

and

$$\mathcal{L}\_1^-(\ell\_1, \ell\_2) = \left\{ \left( \mathfrak{r}\_1^-(\omega, \ell\_1), \mathfrak{r}\_2^+(\omega, \ell\_2) \right) \mid \omega \in I, \left( \ell\_1, \ell\_2 \right) \in \mathbf{Z} \right\}.$$

As -<sup>1</sup> is the horizontal shift parameter and -<sup>2</sup> is the vertical shift parameter, changing these values shifts these segments accordingly. Connecting these segments creates the stability switching curve (SSC, henceforth) under equation (i).

We now turn attention to equation (ii) that can be written as

$$1 + b\_1(\lambda)e^{-\lambda\tau\_1} + b\_2(\lambda)e^{-\lambda\tau\_2} = 0\tag{6}$$

where

$$b\_1(\lambda) = \frac{k}{\lambda} \text{ and } b\_2(\lambda) = \frac{k\beta(n-1)}{\lambda}.$$

With *λ* = *iω*,

$$b\_1(i\omega) = -i\frac{k}{\omega} \text{ and } b\_2(i\omega) = -i\frac{k\beta(n-1)}{\omega}.$$

their absolute values are

$$|b\_1(i\omega)| = \frac{k}{\omega} \text{ and } |b\_2(i\omega)| = \frac{k\beta(n-1)}{\omega}\omega$$

and their arguments are

$$\arg\left(b\_1(i\omega)\right) = \frac{3}{2}\pi \text{ and } \arg\left(b\_2(i\omega)\right) = \frac{3}{2}\pi.$$

By the triangle conditions, the domains of *ω* are defined, respectively, by

$$I\_2 = \left[\frac{1}{2}k, \frac{3}{2}k\right] \text{ if } n=2.$$

and

$$I\_{\mathfrak{N}} = \left[\frac{n-3}{2}k, \frac{n+1}{2}k\right] \text{ if } n \ge 3.$$

As in the same way, the internal angles denoted as ¯ *θ*<sup>1</sup> and ¯ *θ*<sup>2</sup> generated under equation (ii) are obtained as

$$\begin{aligned} \bar{\theta}\_1(\omega) &= \cos^{-1} \left[ \frac{4\omega^2 - k^2 \left(n - 3\right) \left(n + 1\right)}{8k\omega} \right] \\\\ \bar{\theta}\_2(\omega) &= \cos^{-1} \left[ \frac{4\omega^2 + k^2 \left(n - 3\right) \left(n + 1\right)}{4k(n - 1)\omega} \right] \end{aligned}$$

and

$$\eta\_{6}$$

4*k*(*n* − 1)*ω*

.

For any *ω* ∈ *I*<sup>2</sup> or *In*, we may find all pairs of (*τ*1, *τ*2) satisfying (6) as follows:

$$\pi\_1^{\pm}(\omega, m\_1) = \frac{1}{\omega} \left[ \frac{3}{2}\pi + (2m\_1 - 1)\,\pi \pm \bar{\theta}\_1(\omega) \right] \tag{7}$$

and

$$
\pi\_2^{\mp}(\omega, m\_2) = \frac{1}{\omega} \left[ \frac{3}{2}\pi + (2m\_2 - 1)\,\pi \mp \bar{\theta}\_2(\omega) \right]. \tag{8}
$$

As before, we have again two sets of line segments,

$$C\_2^+(m\_1, m\_2) = \left\{ \left( \mathfrak{r}\_1^+(\omega, m\_1), \mathfrak{r}\_2^-(\omega, m\_2) \right) \mid \omega \in I\_2 \text{ or } I\_{n\_1} \ (m\_1, m\_2) \in \mathbb{Z} \right\}.$$

and

$$\mathbb{C}\_2^-(m\_1, m\_2) = \left\{ \left( \mathbb{H}\_1^-(\omega, m\_1), \mathbb{H}\_2^+(\omega, m\_2) \right) \mid \omega \in I \text{ or } I\_{n\prime} \ (m\_1, m\_2) \in \mathbb{Z} \right\}$$

which are shifted horizontally and vertically by changing the values of *m*<sup>1</sup> and *m*2. Connecting these segments creates again the stability switching curves under equation (ii).

#### *3.2. Equal Delays*

Having found the delays' critical values, we may draw attention to the equal delay case before proceeding further with the different delay case. When the delays are equal, conditions (i) and (ii) are changed to

$$\text{(i)}'\,\lambda + k\varepsilon^{-\lambda\tau} - k\beta\varepsilon^{-\lambda\tau} = 0,$$

$$\text{(ii)}'\,\lambda + k\varepsilon^{-\lambda\tau} + k\beta(n-1)\varepsilon^{-\lambda\tau} = 0.$$

For *λ* = *iω* with *ω* > 0, equation (i)' is

$$i\omega + k(1 - \beta) \left(\cos\omega\tau - i\sin\omega\tau\right) = 0.$$

Separating the real and imaginary parts gives the equations,

$$k(1 - \beta)\cos\omega\tau = 0$$

$$k(1 - \beta)\sin\omega\tau = \omega$$

from which

$$
\cos\omega\,\tau = 0,\\
\sin\omega\,\tau = 1 \text{ and } \omega = \frac{k}{2}.
$$

Hence the critical values of *τ* for equation (i)' are determined as

$$
\pi\_\ell^\* = \frac{2}{k} \left( \frac{\pi}{2} + 2\ell\pi \right) \text{ for } \ell = 0, 1, 2, \dots \tag{9}
$$

Similarly, for equation (ii)', we have

$$\cos\omega\tau = 0, \sin\omega\tau = 1 \text{ and } \omega = \frac{k\left(n+1\right)}{2}.$$

Hence the critical values of *τ* are determined as

$$\pi\_m^\*(n) = \frac{2}{k(n+1)} \left(\frac{\pi}{2} + 2m\pi\right) \text{ for } m = 0, 1, 2, \dots \tag{10}$$

It is confirmed that

$$
\tau\_0^\*(n) < \tau\_0^\* \text{ for any } n \ge 2.
$$

Therefore stability switching occurs when *τ* = *τ*∗ <sup>0</sup> (*n*). To check the direction of stability switches, we select *τ* as the bifurcation parameter and consider the eigenvalues as functions of *τ*, *λ* = *λ*(*τ*). Then we differentiate equation (ii)' with respect to *τ*,

$$
\lambda' + k(1 + (n-1)\beta)e^{\lambda \tau} \left(-\lambda'\tau - \lambda\right) = 0
$$

and solving this for *λ* gives

$$
\lambda' = \frac{-\lambda^2}{1 + \lambda\tau}.
$$

The sign of the real part for *λ* = *iω* is positive,

$$\operatorname{Re}\left[\left(\lambda'\right)\_{\lambda=i\omega}\right] = \frac{\omega^2}{1 + \left(\omega\tau\right)^2} > 0.$$

As equations (i)' is obtained from (ii)' with *n* = 0, this derivation also applies to equation (i)'. Hence we have the following result when the delays are equal:

**Theorem 1.** *The Cournot equilibrium is locally asymptotically stable for τ* < *τ*∗ <sup>0</sup> (*n*), *loses its stability at τ* = *τ*∗ <sup>0</sup> (*n*) *and stability cannot be regained for τ* > *τ*<sup>∗</sup> <sup>0</sup> (*n*) *where*

$$
\pi\_0^\*(n) = \frac{\pi}{k(n+1)}.
$$

Theorem 1 is numerically confirmed when

$$\kappa = 10, \ k = 0.5.$$

We perform simulations with three different values of *n*, *n* = 2, *n* = 3, and *n* = 4. The simulations are done with Mathematica, ver. 12.1. The corresponding critical values of *τ* are

$$
\tau\_0^\*(2) = \frac{2\pi}{3}, \ \tau\_0^\*(3) = \frac{\pi}{2} \text{ and } \tau\_0^\*(4) = \frac{2\pi}{5}
$$

which imply that the stability region becomes smaller as *n* increases. This is also clear from the form of *τ*∗ <sup>0</sup> (*n*) in Theorem 1. In each simulation below, we take *τ* = *τ*<sup>∗</sup> <sup>0</sup> (*n*) − 0.2 for the red convergent curve and *τ* = *τ*∗ <sup>0</sup> (*n*) + 0.1 for the divergent green curve and assume constant functions for *t* ≤ 0. In duopoly, the initial functions are defined as

$$
\varphi\_1(t) = \mathfrak{x}\_1^\varepsilon - \mathfrak{z} \text{ and } \varphi\_2(t) = \mathfrak{x}\_2^\varepsilon - 2 \text{ for } t \le 0.
$$

In tiropoly and quartopoly, the appropriate functions are similarly defined and the initial values are selected from the neighborhood of the equilibrium point. Although it is clear that the simulation results strongly depend on the model's specification, we can see the followings from those simulations illustrated in Figure 2A–C:


Results (3) and (4) are inevitable because the best reply functions are linear and the resultant dynamical system does not have enough nonlinearities to prevent the trajectories from becoming negative. We also have essentially the same results in the case of different delays.

**Figure 2.** (**A**) Oscillatory dynamics in duopoly. (**B**) Oscillatory dynamics in triopoly. (**C**) Oscillatory dynamics in quartopoly.
