Transboundary pollution control and environmental absorption efficiency management

Abstract

In this paper, we suggest a two-player differential game model of transboundary pollution that accounts for time-dependent environmental absorption efficiency, which allows the biosphere to switch from a carbon sink to a source. We investigate the impact of negative externalities resulting from a transboundary pollution non-cooperative game wherein countries are dynamically involved. Based on a linear-quadratic specification for the instantaneous revenue function, we assess differences related to both transient path and steady state between cooperative solution, open-loop and Markov perfect Nash equilibria (MPNE). Regarding the methodological contribution of the paper, we suggest a particular structure of the conjectured value function to solve MPNE problems with multiplicative interaction between state variables in one state equation, so that third-order terms that arise in the Hamilton–Jacobi–Bellman equation are made negligible. Using a collocation procedure, we confirm the validity of the particular structure of the conjectured value function. The results suggest unexpected contrasts in terms of pollution control and environmental absorption efficiency management: (i) in the long run, an open-loop Nash equilibrium (OLNE) allows equivalent emissions to the social optimum but requires greater restoration efforts; (ii) although an MPNE is likely to end up with lower emissions and greater restoration efforts than an OLNE, it has a much greater chance of falling in the emergency area; (iii) the absence of cooperation and or precommitment becomes more costly as the initial absorption efficiency decreases; (iv) more heavily discounted MPNE strategies are less robust than OLNE to prevent irreversible switching of the biosphere from a carbon sink to a source.

Appendix

A1. To determine the (non-trivial) cooperative steady state, we solve the canonical system (8)–(11) in the state-costate space. Using (6) and (7) leads to the system’s steady state in (12), $i = 1,\;2$. Note that since:

$$ \begin{aligned} \mathop {\lim }\limits_{t \mapsto + \infty } e^{ - rt} \eta_{1} (t)P(t) & = \mathop {\lim }\limits_{t \mapsto + \infty } \left[ { - {{r\gamma \left( {r^{2} \gamma + h} \right)e^{ - rt} } \mathord{\left/ {\vphantom {{r\gamma \left( {r^{2} \gamma + h} \right)e^{ - rt} } {4\left( {4c + \gamma^{2} } \right)}}} \right. \kern-0pt} {4\left( {4c + \gamma^{2} } \right)}}} \right] = 0 \\ \mathop {\lim }\limits_{t \mapsto + \infty } e^{ - rt} \eta_{2} (t)A(t) & = \mathop {\lim }\limits_{t \mapsto + \infty } \left[ {{{\gamma \left( {2a - r\gamma } \right)e^{ - rt} } \mathord{\left/ {\vphantom {{\gamma \left( {2a - r\gamma } \right)e^{ - rt} } 2}} \right. \kern-0pt} 2}} \right] = 0 \\ \end{aligned} $$

the cooperative steady-state solution satisfies the limiting transversality conditions. □

A2. Plugging the expressions of $e_{i\infty }$, $P_{\infty }$ and $w_{i\infty }$ from (12) in the integrand of (4) and simplifying gives (13). □

A3. Writing the variables in the order $P,A,\eta_{1} ,\eta {}_{2}$, the Jacobian matrix associated with the canonical system (8)–(11) is:

$$ J = \left[ {\begin{array}{*{20}c} { - A} & { - P} & 2 & 0 \\ { - \gamma } & 0 & 0 & 2 \\ {2c} & {\eta_{1} } & {r + A} & \gamma \\ {\eta_{1} } & 0 & P & r \\ \end{array} } \right] $$

The determinant of the Jacobian matrix is:

$$ \left| J \right| = \left( {4c + \gamma^{2} } \right)P_{\infty }^{2} + r\gamma A_{\infty } P_{\infty } = {{\left( {r^{2} \gamma + h} \right)^{2} } \mathord{\left/ {\vphantom {{\left( {r^{2} \gamma + h} \right)^{2} } {4\left( {4c + \gamma^{2} } \right)}}} \right. \kern-0pt} {4\left( {4c + \gamma^{2} } \right)}} + r\gamma \left( {2a - r\gamma } \right) $$

Using Dockner’s formula (Dockner 1985), the sum of the principal minors of $J$ of order 2 minus the squared discount rate, denoted by $K$, is:

$$ K = - A_{\infty } \left( {r + A_{\infty } } \right) - 4c - 2\gamma P_{\infty } = - \frac{{\left[ {2\left( {2a - r\gamma } \right)\left( {4c + \gamma^{2} } \right) + r\left( {r^{2} \gamma + h} \right)} \right]^{2} }}{{\left( {r^{2} \gamma + h} \right)^{2} }} - \frac{{\gamma h + 4c\left( {4c + \gamma^{2} - r^{2} } \right)}}{{\left( {4c + \gamma^{2} } \right)}} $$

Because $\left| J \right| > 0$ and $K < 0$, the cooperative steady state is a saddle-point. We compute:

$$ \varOmega = K^{2} - 4\left| J \right| = A_{\infty } \left( {r + A_{\infty } } \right)\left[ {A_{\infty } \left( {r + A_{\infty } } \right) + 8c} \right] + 4\gamma A_{\infty }^{2} P_{\infty } + 16c\left( {c + \gamma P_{\infty } } \right) - 16cP_{\infty }^{2} $$

The sign of $\varOmega$ can be either negative or positive. In the first case, the steady state is a saddle-focus with transient oscillations. Otherwise, it is a saddle-node and the optimal solution monotonically converges to the steady state. □

A4. The linear approximation of the system (8)–(11) around its steady state is:

$$ \begin{aligned} \dot{\eta }_{1} &= 2cP + \eta_{1} \left( {r + A_{\infty } } \right) + \eta_{1\infty } \left( {A - A_{\infty } } \right) + \gamma \eta_{2} \\ \dot{\eta }_{2} & = r\eta_{2} + \eta_{1\infty } P + P_{\infty } \left( {\eta_{1} - \eta_{1\infty } } \right) \\ \dot{P} & = 2\left( {a + \eta_{1} } \right) - P_{\infty } \left( {A - A_{\infty } } \right) - A_{\infty } P \\ \dot{A} & = 2\eta_{2} - \gamma P \\ \end{aligned} $$

The four eigenvalues associated with the Jacobian matrix of the canonical system are:

$$ \begin{aligned} &{}_{1}^{3} \chi_{2}^{4} = \frac{r}{2} \pm \sqrt {\frac{{r^{2} }}{4} - \frac{K}{2} \pm \frac{1}{2}\sqrt {K^{2} - 4\left| J \right|} } \\ &\quad = \frac{r}{2} \pm \sqrt {\frac{{\left( {r + A_{\infty } } \right)^{2} + A_{\infty }^{2} }}{4} + \gamma P_{\infty } + \frac{8ac}{{r^{2} \gamma^{2} }} \pm \sqrt {\left( {A_{\infty }^{2} + \frac{8ac}{{r^{2} \gamma^{2} }}} \right)\left( {\gamma P_{\infty } + \frac{8ac}{{r^{2} \gamma^{2} }}} \right) + \left[ {\frac{{\left( {r + A_{\infty } } \right)^{2} A_{\infty } }}{4} + \frac{8ac}{{r\gamma^{2} }}} \right]A_{\infty } - 4cP_{\infty }^{2} } } \\ \end{aligned} $$

where $A_{\infty }$ and $P_{\infty }$ are given in (12). As expected, two eigenvalues, $\chi_{1}$ and $\chi_{3}$, have negative real part and two have positive real part, $\chi_{2}$ and $\chi_{4}$. Choosing the roots with negative real part for convergence, the time paths of the costate and state variables are written as:

$$ \begin{aligned} \eta_{1} (t) & = \eta_{1\infty } + B_{1} e^{{\chi_{1} t}} + B_{2} e^{{\chi_{3} t}} \\ \eta_{2} (t) & = \eta_{2\infty } + B_{3} e^{{\chi_{1} t}} + B_{4} e^{{\chi_{3} t}} \\ P(t) & = P_{\infty } + B_{5} e^{{\chi_{1} t}} + B_{6} e^{{\chi_{3} t}} \\ A(t) & = A_{\infty } + B_{7} e^{{\chi_{1} t}} + B_{8} e^{{\chi_{3} t}} \\ \end{aligned} $$

These equations involve 10 unknowns ($B_{1} , \ldots ,B_{8}$,$\eta_{1} (0),\;\eta_{2} (0)$) that can be solved with the 10 following equations, which are drawn from the above expressions and the linearized versions of (8)–(11):

$$ \begin{aligned} \eta_{1} \left( 0 \right) & = \eta_{1\infty } + B_{1} + B_{2} \\ \eta_{2} \left( 0 \right) & = \eta_{2\infty } + B_{3} + B_{4} \\ \end{aligned} $$

$$ \begin{aligned} P_{0} & = P_{\infty } + B_{5} + B_{6} \\ A_{0} & = A_{\infty } + B_{7} + B_{8} \\ \end{aligned} $$

$$ \begin{aligned} \dot{\eta }_{1} \left( 0 \right) & = 2cP_{0} + \eta_{1} \left( 0 \right)\left( {r + A_{\infty } } \right) + \eta_{1\infty } \left( {A_{0} - A_{\infty } } \right) + \gamma \eta_{2} \left( 0 \right) = B_{1} \chi_{1} + B_{2} \chi_{3} \\ \dot{\eta }_{2} \left( 0 \right) & = r\eta_{2} \left( 0 \right) + \eta_{1\infty } \left( {P_{0} - P_{\infty } } \right) + \eta_{1} \left( 0 \right)P_{\infty } = B_{3} \chi_{1} + B_{4} \chi_{3} \\ \dot{P}\left( 0 \right) & = 2\left( {a + \eta_{1} \left( 0 \right)} \right) - P_{\infty } \left( {A_{0} - A_{\infty } } \right) - A_{\infty } P_{0} = B_{5} \chi_{1} + B_{6} \chi_{3} \\ \dot{A}\left( 0 \right) & = 2\eta_{2} \left( 0 \right) - \gamma P_{0} = B_{7} \chi_{1} + B_{8} \chi_{3} \\ & 2\left( {B_{3} \chi_{1} + B_{4} \chi_{3} } \right) - \gamma \left( {B_{5} \chi_{1} + B_{6} \chi_{3} } \right) = B_{7} \chi_{1}^{2} + B_{8} \chi_{3}^{2} \end{aligned} $$

$$ 2\left( {B_{1} \chi_{1} + B_{2} \chi_{3} } \right) - A_{\infty } \left( {B_{5} \chi_{1} + B_{6} \chi_{3} } \right) - P_{\infty } \left( {B_{7} \chi_{1} + B_{8} \chi_{3} } \right) = B_{5} \chi_{1}^{2} + B_{6} \chi_{3}^{2} $$

This system of equations can be solved numerically to determine real-valued solutions. Finally, using (6)–(7) yields (14)–(17). □

A5. To compute the (non-trivial) OLNE steady state, we solve the system in the state-costate space (21)–(24). Using (19) and (20) leads us to derive the steady state in (25), $i = 1,\;2$. Note that since:

$$ \begin{aligned} \mathop {\lim }\limits_{t \mapsto + \infty } e^{ - rt} \lambda_{1} (t)P(t) & = \mathop {\lim }\limits_{t \mapsto + \infty } \left[ { - {{r\gamma \left( {r^{2} \gamma + k} \right)e^{ - rt} } \mathord{\left/ {\vphantom {{r\gamma \left( {r^{2} \gamma + k} \right)e^{ - rt} } {4\left( {2c + \gamma^{2} } \right)}}} \right. \kern-0pt} {4\left( {2c + \gamma^{2} } \right)}}} \right] = 0 \\ \mathop {\lim }\limits_{t \mapsto + \infty } e^{ - rt} \lambda_{2} (t)A(t) & = \mathop {\lim }\limits_{t \mapsto + \infty } \left[ {{{\gamma \left( {2a - r\gamma } \right)e^{ - rt} } \mathord{\left/ {\vphantom {{\gamma \left( {2a - r\gamma } \right)e^{ - rt} } 2}} \right. \kern-0pt} 2}} \right] = 0 \\ \end{aligned} $$

the OLNE steady-state solution satisfies the limiting transversality conditions. □

A6. Plugging the expressions of $e_{i\infty }^{on}$, $P_{\infty }^{on}$ and $w_{i\infty }^{on}$ from (25) in the integrand of (3) and simplifying gives (26). □

A7. Writing the variables in the order $P,A,\lambda_{1} ,\lambda {}_{2}$, the Jacobian matrix associated with the canonical system (21)–(24) is:

$$ J^{on} = \left[ {\begin{array}{*{20}c} { - A} & { - P} & 2 & 0 \\ { - \gamma } & 0 & 0 & 2 \\ c & {\lambda_{1} } & {r + A} & \gamma \\ {\lambda_{1} } & 0 & P & r \\ \end{array} } \right] $$

The determinant of $J^{on}$ is:

$$ \left| {J^{on} } \right| = \left( {2c + \gamma^{2} } \right)P_{\infty }^{on2} + r\gamma A_{\infty }^{on} P_{\infty }^{on} = \frac{{\left( {r^{2} \gamma + k} \right)^{2} }}{{4\left( {2c + \gamma^{2} } \right)}} + r\gamma \left( {2a - r\gamma } \right) $$

Further, the sum of the principal minors of $J^{on}$ of order 2 minus the squared discount rate is:

$$ K^{on} = - A_{\infty }^{on} \left( {r + A_{\infty }^{on} } \right) - 2c - 2\gamma P_{\infty }^{on} = - \frac{{2\left( {2a - r\gamma } \right)\left( {2c + \gamma^{2} } \right)}}{{r^{2} \gamma + k}}\left[ {r + \frac{{2\left( {2a - r\gamma } \right)\left( {2c + \gamma^{2} } \right)}}{{r^{2} \gamma + k}}} \right] - \frac{{\gamma \left( {r^{2} \gamma + k} \right)}}{{\left( {2c + \gamma^{2} } \right)}} - 2c $$

which is strictly negative. Hence, the OLNE steady state is a saddle-point. We now compute:$\varOmega^{on} = \left( {K^{on} } \right)^{2} - 4\left| {J^{on} } \right| = A_{\infty }^{on2} \left[ {\left( {r + A_{\infty }^{on} } \right)^{2} + 4\gamma P_{\infty }^{on} } \right] + \frac{16ac}{{r^{2} \gamma^{2} }}\left[ {A_{\infty }^{on} \left( {r + A_{\infty }^{on} } \right) + \frac{4ac}{{r^{2} \gamma^{2} }}} \right] + 8c\left( {\frac{4a}{{r^{2} \gamma }} - P_{\infty }^{on} } \right)P_{\infty }^{on}$ where $A_{\infty }^{on}$ and $P_{\infty }^{on}$ are given in (25). The sign of this expression, which is unclear, determines whether the saddle path converges monotonically or oscillates towards the steady state. □

A8. The linear approximation of system (21)–(24) around its steady state is:

$$ \begin{aligned} \dot{\lambda }_{1} & = cP + \lambda_{1} \left( {r + A_{\infty } } \right) + \lambda_{1\infty } \left( {A - A_{\infty } } \right) + \gamma \lambda_{2} \\ \dot{\lambda }_{2} & = r\lambda_{2} + \lambda_{1\infty } P + P_{\infty } \left( {\lambda_{1} - \lambda_{1\infty } } \right) \\ \dot{P} & = 2\left( {a + \lambda_{1} } \right) - P_{\infty } \left( {A - A_{\infty } } \right) - A_{\infty } P \\ \dot{A} & = 2\lambda_{2} - \gamma P \\ \end{aligned} $$

The eigenvalues associated with the Jacobian of the non-cooperative canonical system are:

$$\begin{aligned} &{}_{1}^{3} \xi_{2}^{4} = {r \mathord{\left/ {\vphantom {r 2}} \right. \kern-0pt} 2} \pm \\ &\quad \sqrt {{{\left[ {\left( {r + A_{\infty }^{on} } \right)^{2} + \left( {A_{\infty }^{on} } \right)^{2} } \right]} \mathord{\left/ {\vphantom {{\left[ {\left( {r + A_{\infty }^{on} } \right)^{2} + \left( {A_{\infty }^{on} } \right)^{2} } \right]} 4}} \right. \kern-0pt} 4} + \gamma P_{\infty }^{on} + {{4ac} \mathord{\left/ {\vphantom {{4ac} {r^{2} \gamma^{2} }}} \right. \kern-0pt} {r^{2} \gamma^{2} }} \pm \sqrt {\left[ {{{A_{\infty }^{on} \left( {r + A_{\infty }^{on} } \right)} \mathord{\left/ {\vphantom {{A_{\infty }^{on} \left( {r + A_{\infty }^{on} } \right)} 2}} \right. \kern-0pt} 2} + {{2ac} \mathord{\left/ {\vphantom {{2ac} {r^{2} \gamma^{2} }}} \right. \kern-0pt} {r^{2} \gamma^{2} }} + \gamma P_{\infty }^{on} } \right]^{2} - \left( {2c + \gamma^{2} } \right)\left( {P_{\infty }^{on} } \right)^{2} - 4a} } \end{aligned}$$

where $A_{\infty }^{on}$ and $P_{\infty }^{on}$ are given in (25). Choosing roots with negative real part, $\xi_{1}$ and $\xi_{3}$, and using similar methods as in the cooperative case leads us to find $D_{1} , \ldots ,D_{8}$, and $\lambda_{1} (0),\lambda_{2} (0)$. □

A9. Assuming a second order Taylor-polynomial approximation of the value function centered on $\left( {A_{\infty }^{mp} ,P_{\infty }^{mp} } \right)$, i.e.,

$$ V = \alpha_{1} + \alpha_{2} \left( {A - A_{\infty }^{mp} } \right) + \alpha_{3} \left( {P - P_{\infty }^{mp} } \right) + \alpha_{4} \left( {A - A_{\infty }^{mp} } \right)^{2} + \alpha_{5} \left( {P - P_{\infty }^{mp} } \right)^{2} + \alpha_{6} \left( {A - A_{\infty }^{mp} } \right)\left( {P - P_{\infty }^{mp} } \right) $$

we obtain constant second order derivatives for $V$. We assume that its unknown coefficients $\alpha_{1} ,..,\alpha_{6}$ and the problem parameters are such that the steady state is positive. Plugging (37) and (38) in (1)–(2), respectively, gives the MPNE controlled system, that is:

$$ \begin{aligned} \dot{P} & = 2\left[ {a + \alpha_{3} + 2\alpha_{5} \left( {P - P_{\infty }^{mp} } \right) + \alpha_{6} \left( {A - A_{\infty }^{mp} } \right)} \right] - AP \\ \dot{A} & = 2\left[ {\alpha_{2} + 2\alpha_{4} \left( {A - A_{\infty }^{mp} } \right) + \alpha_{6} \left( {P - P_{\infty }^{mp} } \right)} \right] - \gamma P \\ \end{aligned} $$

The resulting approximation of the HJB Eq. (35) is:

$$ \begin{aligned} & r\left[ {\alpha_{1} + \alpha_{2} \left( {A - A_{\infty }^{mp} } \right) + \alpha_{3} \left( {P - P_{\infty }^{mp} } \right) + \alpha_{4} \left( {A - A_{\infty }^{mp} } \right)^{2} + \alpha_{5} \left( {P - P_{\infty }^{mp} } \right)^{2} + \alpha_{6} \left( {A - A_{\infty }^{mp} } \right)\left( {P - P_{\infty }^{mp} } \right)} \right] \\ & = \left[ {\alpha_{3} + 2\alpha_{5} \left( {P - P_{\infty }^{mp} } \right) + \alpha_{6} \left( {A - A_{\infty }^{mp} } \right)} \right]\left\{ {\frac{3}{2}\left[ {\alpha_{3} + 2\alpha_{5} \left( {P - P_{\infty }^{mp} } \right) + \alpha_{6} \left( {A - A_{\infty }^{mp} } \right)} \right] - AP + 2a} \right\} \\ & \quad + \left[ {\alpha_{2} + 2\alpha_{4} \left( {A - A_{\infty }^{mp} } \right) + \alpha_{6} \left( {P - P_{\infty }^{mp} } \right)} \right]\left\{ {3{{\left[ {\alpha_{2} + 2\alpha_{4} \left( {A - A_{\infty }^{mp} } \right) + \alpha_{6} \left( {P - P_{\infty }^{mp} } \right)} \right]} \mathord{\left/ {\vphantom {{\left[ {\alpha_{2} + 2\alpha_{4} \left( {A - A_{\infty }^{mp} } \right) + \alpha_{6} \left( {P - P_{\infty }^{mp} } \right)} \right]} 2}} \right. \kern-0pt} 2} - \gamma P} \right\}\\ &\quad- {{cP^{2} } \mathord{\left/ {\vphantom {{cP^{2} } 2}} \right. \kern-0pt} 2} + {{a^{2} } \mathord{\left/ {\vphantom {{a^{2} } 2}} \right. \kern-0pt} 2} \\ \end{aligned} $$

Rearranging its terms, this can be expressed as a third-order polynomial in $\left( {A - A_{\infty }^{mp} } \right)$ and $\left( {P - P_{\infty }^{mp} } \right)$. Then, the unknown parameters of the approximation of the value function are found by setting to 0 the coefficients, up to the second order, of such a third-order polynomial (its third-order terms are locally negligible, by construction). To do so, one solves the system:

$$ \begin{aligned} r\alpha_{1} - 3\alpha_{2}^{2} /2 - \alpha_{3} \left( {3\alpha_{3} /2 + 2a} \right) - a^{2} /2 + \gamma \alpha_{2} P_{\infty }^{mp} + \alpha_{3} A_{\infty }^{mp} P_{\infty }^{mp} + c\left( {P_{\infty }^{mp} } \right)^{2} /2 = 0 \hfill \\ \alpha_{2} \left( {\gamma - 3\alpha_{6} } \right) + r\alpha_{3} - 2\alpha_{5} \left( {3\alpha_{3} + 2a} \right) + \alpha_{3} A_{\infty }^{mp} + \left( {c + \gamma \alpha_{6} } \right)P_{\infty }^{mp} + 2\alpha_{5} A_{\infty }^{mp} P_{\infty }^{mp} = 0 \hfill \\ {c \mathord{\left/ {\vphantom {c 2}} \right. \kern-0pt} 2} + \alpha_{5} \left( {r - 6\alpha_{5} + 2A_{\infty }^{mp} } \right) + \alpha_{6} \left( {\gamma - {{3\alpha_{6} } \mathord{\left/ {\vphantom {{3\alpha_{6} } 2}} \right. \kern-0pt} 2}} \right) = 0 \hfill \\ \alpha_{2} \left( {r - 6\alpha_{4} } \right) - \alpha_{6} \left( {3\alpha_{3} + 2a} \right) + \left( {\alpha_{3} + 2\gamma \alpha_{4} } \right)P_{\infty }^{mp} + \alpha_{6} A_{\infty }^{mp} P_{\infty }^{mp} = 0 \hfill \\ \alpha_{3} + 2\gamma \alpha_{4} - 6\alpha_{6} \left( {\alpha_{4} + \alpha_{5} } \right) + \alpha_{6} \left( {r + A_{\infty }^{mp} } \right) + 2\alpha_{5} P_{\infty }^{mp} = 0 \hfill \\ \alpha_{4} \left( {r - 6\alpha_{4} } \right) - \alpha_{6} \left( {3\alpha_{6} /2 - P_{\infty }^{mp} } \right) = 0 \hfill \\ \end{aligned} $$

We have 6 equations in 8 unknowns $\alpha_{1} , \ldots ,\alpha_{6} ,A_{\infty }^{mp} ,P_{\infty }^{mp}$. To get two other equations, we use (37)–(38), imposing a steady state solution for the MPNE:

$$ 2\left( {a + \alpha_{3} } \right) - A_{\infty }^{mp} P_{\infty }^{mp} = 0 $$

$$ 2\alpha_{2} - \gamma P_{\infty }^{mp} = 0 $$

Finally, we obtain a non-linear system of 8 equations in 8 unknowns, which in general might be solved numerically, but not analytically. When it can be solved, its solution provides:

$$ \left[ {P_{\infty }^{mp} ,A_{\infty }^{mp} ,e_{\infty }^{mp} ,w_{\infty }^{mp} } \right]^{T} = \left[ {2\alpha_{2} /\gamma ,\gamma \left( {a + \alpha_{3} } \right)/\alpha_{2} ,a + \alpha_{3} ,\alpha_{2} } \right]^{T} $$

Clearly, the necessary conditions for the positivity of the steady state values and the positivity of each player’s instantaneous revenue are $\alpha_{2} > 0$ and $\alpha_{3} > - a$. In general, the non-linear system above may admit multiple solutions (or even no solutions at all). To solve the system, we use the Matlab function fsolve, using default options, and starting with a random initialization of the 8 unknowns, initially chosen as realizations of independent and uniformly distributed random variables on the interval $[0,10]$. To check whether a numerically found solution of the non-linear system above is likely also a steady state solution for the MPNE, we perform the following additional checks:

• (1) Positivity of the state and control variables, and positivity of each player’s instantaneous revenue: $\alpha_{2} > 0$ and $- a < \alpha_{3} < 0$;

• (2) Global asymptotical stability for the linearized dynamical system: the eigenvalues of the Jacobian matrix (40) should have negative real parts;

• (3) Comparison between the value of $\alpha_{1}$ assumed in $\left( {A_{\infty }^{mp} ,P_{\infty }^{mp} } \right)$ by the local quadratic approximation of the value function, and the one obtained by replacing the resulting steady state values in (3), the latter being equal to $\left[ {e_{\infty }^{mp} \left( {a - e_{\infty }^{mp} /2} \right) - c\left( {P_{\infty }^{mp} } \right)^{2} /2 - \left( {w_{\infty }^{mp} } \right)^{2} /2} \right]/r$: in the case where $P_{\infty }^{mp} ,A_{\infty }^{mp} ,e_{\infty }^{mp} ,w_{\infty }^{mp}$ is really an approximation of a steady state solution for the MPNE, and the local quadratic approximation is a good approximation of the value function, then the two values should be approximately the same. In practice, we check whether the relative error of the first expression with respect to the second one is smaller than 0.01.

In the case where multiple solutions are obtained that satisfy the requirements (1), (2), and (3) above, one can choose the one associated with the largest value of $\left.\left[ {e_{\infty }^{mp} \left( {a - e_{\infty }^{mp} /2} \right) - c\left( {P_{\infty }^{mp} } \right)^{2} /2 - \left( {w_{\infty }^{mp} } \right)^{2} /2} \right]\right/r$. □

A10. To perform an additional check for the validity of the local quadratic approximation of the value function, we also exploit an alternative method to approximate the HJB equation and then we check whether the two obtained approximations are similar. As the alternative approximation method,⁷ we use the collocation method for PDEs, applied to the HJB Eq. (35), on a small neighborhood of the steady state solution determined by the local quadratic approximation (in our implementation, we chose the box $\left[ {0.75A_{\infty }^{mp} ,1.25A_{\infty }^{mp} } \right] \times \left[ {0.75P_{\infty }^{mp} ,1.25P_{\infty }^{mp} } \right]$ as such a neighborhood). The collocation method approximates locally on that neighborhood the unknown value function as a high-order polynomial in $\left( {A - A_{\infty }^{mp} } \right)$ and $\left( {P - P_{\infty }^{mp} } \right)$, with unknown coefficients. These are determined by imposing that the HJB Eq. (35) holds exactly on a finite number of fixed points (“collocation points”) inside the neighborhood (or equivalently, that its Bellman residual error is 0 on such points). The number of these points is equal to the number of unknown coefficients of the high-order polynomial. In this way, to find the values of the coefficients, one is reduced to solving a non-linear system (due to the non-linearity of the HJB equation), where each equation is associated with one collocation point. For the implementation of the collocation method, we follow the procedure proposed in (Doraszelski 2003), that is:

• We construct a polynomial approximation of order $K = 12$, using a tensor-product basis of univariate Chebyshev polynomials. This choice of the approximation order is guided by the need to find a suitable trade-off between the accuracy of the approximation, and the computational effort needed to find it (as in the cited reference, the number of coefficients in the expansion is $(K + 1)^{2}$);

• Likewise in Doraszelski (2003), we employ a two-dimensional expanded Chebyshev array (as defined in Judd 1998) to define the $(K + 1)^{2}$ collocation points.

To initialize the collocation method, at first we fit the initial local quadratic approximation of the value function by a linear combination of the same basis functions used by the collocation method itself (in practice, we apply the collocation method itself to this function approximation problem, instead of to the PDE). Then, we apply the collocation method to the PDE, starting from the resulting vector of coefficients. Finally, to check for the validity of the resulting high-order polynomial approximation and for the validity of the original local quadratic approximation, we compute, for both approximations, the Bellman residual error on a larger uniform grid of $N^{2}$ points in the same neighborhood (with $N = 30$), which are different from the collocation points, checking whether this error is small. We also evaluate, on each grid point, the relative error of the quadratic approximation with respect to the high-order polynomial approximation obtained using the collocation method. The results of all these comparisons are reported in Sect. 5, for a specific choice of the problem parameters. We conclude by mentioning that, in contrast with Doraszelski (2003), it is not necessary to check a posteriori whether the partial derivatives of the approximation of the value function obtained by the collocation method are good approximations of the corresponding partial derivatives of the value function, which are needed to express the optimal controls [see the expressions (37) and (38)]. Indeed, in the case where the local quadratic approximation and the one found by the collocation method are similar, we can use the partial derivatives of the local quadratic approximation itself. In any case, it is worth mentioning that the numerical results have shown good agreement between the approximations of the partial derivatives obtained by the two methods, as one can infer from Fig. 2 in Sect. 5. □

A11. The results in (41)–(45) follow from the respective comparisons of (12) and (25). □