Optimal Abatement and Emission Permit Trading Policies in a Dynamic Transboundary Pollution Game

Abstract

We obtain optimal emission levels and abatement expenditures in a finite-horizon transboundary pollution game with emission trading between two regions. We show that emission trading has significant impact on the optimal strategies and profits of the two regions. We find that cooperation between the regions leads to increased abatement and lower emissions, resulting in a lower pollution stock. We also provide a stochastic extension in which the pollution stock and the emission trading price are diffusion processes and solve it numerically.

Appendices

Appendix

A. Proof of Proposition 1

Proof

For $i=1,2$, we guess the form of the value functions $V_{\mathrm{Ni}}(P,t)$ to be linear in P, that is,

$$\begin{aligned} V_{\mathrm{Ni}}(P,t)=l_i(t)P+k_{i}(t), \end{aligned}$$

(20)

where $l_i(t)$ and $k_i(t)$ are functions of t. Then, we obtain these functions to verify that our guess is correct. By substituting the first-order conditions (6) and (20) into the system of HJB Eq. (5), we have

$$\begin{aligned} \left\{ \begin{aligned} \left( l^{'}_1(t)-(r+\theta )l_1(t)-D\right) P+k^{'}_1(t)-rk_1(t)+f_1(t)=0, \\ \left( l^{'}_2(t)-(r+\theta )l_2(t)-\beta D\right) P+k^{'}_2(t)-rk_2(t)+f_2(t)=0, \end{aligned} \right. \end{aligned}$$

(21)

where

$$\begin{aligned} f_1(t)=SE_{10}+\frac{1}{2}\Big (A-S+l_1(t)\Big )^2+l_1(t)\Big (\alpha A-S+\frac{\gamma C+1}{\gamma C}l_2(t)+\frac{1}{2C}l_1(t)\Big ), \end{aligned}$$

and

$$\begin{aligned} f_2(t)=SE_{20}+\frac{1}{2}\Big (\alpha A-S+l_2(t)\Big )^2+l_2(t)\Big (A-S+\frac{C+1}{\gamma C}l_1(t)+\frac{1}{2\gamma C}l_2(t)\Big ). \end{aligned}$$

Noticing that the system (21) should be satisfied for all $P\ge 0$, we can determine $l_1(t)$, $l_2(t)$, $k_1(t)$, and $k_2(t)$ by solving the following system of ordinary differential equations:

$$\begin{aligned} \left\{ \begin{aligned}&l^{'}_1(t)-(r+\theta )l_1(t)-D=0,\\&l^{'}_2(t)-(r+\theta )l_2(t)-\beta D=0,\\&k^{'}_1(t)-rk_1(t)+f_1(t)=0, \\&k^{'}_2(t)-rk_2(t)+f_2(t)=0. \end{aligned} \right. \end{aligned}$$

(22)

Solving the above ODE system, we can obtain

$$\begin{aligned} l_1(t)&=\Big (\frac{D}{r+\theta }-g_1\Big )e^{-(r+\theta )(T-t)}-\frac{D}{r+\theta },\\ l_2(t)&=\Big (\frac{\beta D}{r+\theta }-g_2\Big )e^{-(r+\theta )(T-t)}-\frac{\beta D}{r+\theta },\\ k_1(t)&=G_{N1}+H_{N1}e^{-(r+\theta )(T-t)}+I_{N1}e^{-2(r+\theta )(T-t)}+\frac{SE_{10}e^{-\rho _1 t}}{r+\rho _1}\\&\quad +\Big (g_1{\bar{P}}_1-I_{N1}-H_{N1}-G_{N1}-\frac{SE_{10}e^{-\rho _1 T}}{r+\rho _1}\Big )e^{-r(T-t)},\\ k_2(t)&=G_{N2}+H_{N2}e^{-(r+\theta )(T-t)}+I_{N2}e^{-2(r+\theta )(T-t)}+\frac{SE_{20}e^{-\rho _2 t}}{r+\rho _2}\\&\quad +\Big (g_2{\bar{P}}_2-I_{N2}-H_{N2}-G_{N2}-\frac{SE_{20}e^{-\rho _2 T}}{r+\rho _2}\Big )e^{-r(T-t)}, \end{aligned}$$

where

$$\begin{aligned} G_{N1}&=\frac{1}{r}\Bigg (\frac{S^2}{2C}+\frac{1}{2}(A-S)^2-\frac{D}{r+\theta }\Big ((1+\alpha )A-\frac{2\gamma C+\gamma +1}{\gamma C}S\Big )\\&\quad +\frac{1}{C}\Big (\frac{C+1}{2}+\beta \frac{\gamma C+1}{\gamma }\Big )\Big (\frac{D}{r+\theta }\Big )^2\Bigg ),\\ H_{N1}&=\frac{1}{\theta }\Bigg (-\Big (\frac{D}{r+\theta }-g_1\Big )\Big ((1+\alpha )A-\frac{2\gamma C+\gamma +1}{\gamma C}S\Big )+\frac{C+1}{C}\frac{D}{r+\theta }\Big (\frac{D}{r+\theta }-g_1\Big )\\&\quad +\frac{\gamma C+1}{\gamma C}\frac{D}{r+\theta }\Big (\frac{2\beta D}{r+\theta }-\beta g_1-g_2\Big )\Bigg ),\\ I_{N1}&=-\frac{1}{r+2\theta }\Bigg (\frac{D}{r+\theta }-g_1\Bigg )\Bigg (\frac{C+1}{2C}\left( \frac{D}{r+\theta }-g_1\right) +\frac{\gamma C+1}{\gamma C}\left( \frac{\beta D}{r+\theta }-g_2\right) \Bigg ),\\ G_{N2}&=\frac{1}{r}\Bigg (\frac{S^2}{2\gamma C}+\frac{1}{2}(\alpha A-S)^2-\frac{\beta D}{r+\theta }\Big ((1+\alpha )A-\frac{2\gamma C+\gamma +1}{\gamma C}S\Big )\\&\quad +\frac{\beta }{C}\Big (C+1+\beta \frac{\gamma C+1}{2\gamma }\Big )\frac{ D^2}{(r+\theta )^2}\Bigg ),\\ H_{N2}&=\frac{1}{\theta }\Bigg (-\Big (\frac{\beta D}{r+\theta }-g_2\Big )\Big ((1+\alpha )A-\frac{2\gamma C+\gamma +1}{\gamma C}S\Big )+\frac{\gamma C+1}{\gamma C}\frac{\beta D}{r+\theta }\Big (\frac{\beta D}{r+\theta }-g_2\Big )\\&\qquad +\,\frac{C+1}{ C}\frac{D}{r+\theta }\Big (\frac{2\beta D}{r+\theta }-\beta g_1-g_2\Big )\Bigg ),\\ I_{N2}&=-\frac{1}{r+2\theta }\Bigg (\frac{\beta D}{r+\theta }-g_2\Bigg )\Bigg (\frac{\gamma C+1}{2\gamma C}\Big (\frac{\beta D}{r+\theta }-g_2\Big )+\frac{C+1}{C}\Big (\frac{ D}{r+\theta }-g_1\Big )\Bigg ). \end{aligned}$$

Then, the proof can be completed by substituting the above results into (6). $\square$

B. Proof of Corollary 1

Proof

• (i) Comparing (7b) and (7e) with the two regions’ instantaneous quotas $E_{10}$ and $E_{20}$, we can obtain the above mentioned conclusions.

• (ii) From Eqs. (7a) and (7d) we know

  $$\begin{aligned} \frac{\partial V_{N1}}{\partial P}= & {} \left( \frac{ D}{r+\theta }-g_1\right) e^{-(r+\theta )(T-t)}-\frac{ D}{r+\theta }<0,\\ \frac{\partial V_{N1}}{\partial E_{10}}= & {} \left( e^{\rho _1(T-t)}-e^{-r(T-t)}\right) \frac{e^{-\rho _1T}S}{r+\rho _1}>0,\\ \frac{\partial V_{N2}}{\partial P}= & {} \left( \frac{\beta D}{r+\theta }-g_2\right) e^{-(r+\theta )(T-t)}-\frac{\beta D}{r+\theta }<0, \\ \frac{\partial V_{N2}}{\partial E_{20}}= & {} \left( e^{\rho _2(T-t)}-e^{-r(T-t)}\right) \frac{e^{-\rho _2T}S}{r+\rho _2}>0, \end{aligned}$$

  which completes the proof of results in (ii).

• (iii) Taking derivative on (7a) and (7d) with respect to S, respectively, we have

  $$\begin{aligned} \frac{\partial V_{N1}}{\partial S}&=\frac{1}{r}\Big (\frac{C+1}{C}S-A+\frac{2\gamma C+\gamma +1}{\gamma C}\frac{D}{r+\theta }\Big )\Big (1-e^{-r(T-t)}\Big )\\&\quad +\frac{E_{10}e^{-\rho _1 T}}{1+\rho _1}\Big (e^{\rho _1(T-t)}-e^{-r(T-t)}\Big )\\&\quad +\frac{1}{\theta }\frac{2\gamma C+\gamma +1}{\gamma C}e^{-r(T-t)}\Big (\frac{D}{r+\theta }-g_1\Big )\Big (e^{-\theta (T-t)}-1\Big ), \\ \frac{\partial V_{N2}}{\partial S}&=\frac{1}{r}\Big (\frac{\gamma C+1}{\gamma C}S-A+\frac{2\gamma C+\gamma +1}{\gamma C}\frac{\beta D}{r+\theta }\Big )\Big (1-e^{-r(T-t)}\Big )\\&\quad +\frac{E_{20}e^{-\rho _2 T}}{1+\rho _2}\Big (e^{\rho _2(T-t)}-e^{-r(T-t)}\Big )\\&\quad +\frac{1}{\theta }\frac{2\gamma C+\gamma +1}{\gamma C}e^{-r(T-t)}\Big (\frac{\beta D}{r+\theta }-g_2\Big )\Big (e^{-\theta (T-t)}-1\Big ), \end{aligned}$$

  which imply that $V_{N1}$ and $V_{N2}$ are two quadratic functions concerning S on the real axis at any time before T, and the minimum points should be $\phi -\tau _1(t)E_{10}$ and $\psi -\tau _2(t)E_{20}$, respectively. Also, combining the condition $0<S<\eta$, which must be satisfied to keep the emissions to be positive, and the condition $S>0$, we can conclude the results in (iii).$\square$

C. Proof of Corollary 2

Proof

Substituting (7b), (7c), (7e) and (7f) into (3), we gain the following ordinary differential equation satisfied by the dynamic process of the pollution stock:

$$\begin{aligned} \dot{P_N}+\theta P_N=(1+\alpha )A-\frac{2\gamma C+\gamma +1}{\gamma C}S+\frac{C+1}{C}l_1(t)+\frac{\gamma C+1}{\gamma C}l_2(t). \end{aligned}$$

By solving the above equation, we can obtain the trajectory of the pollution stock (8).

Then, taking the derivative of (8) with respect to t, we have

$$\begin{aligned} P'_N(t)=e^{-\theta t}\left( Y_{N}(r+\theta )e^{-(r+\theta )(T-t)+\theta t}-\theta \left( P_0-X_{N}-Y_{N}e^{-(r+\theta )T}\right) \right) , \end{aligned}$$

from which the conclusions follow. $\square$

D. Proof of Proposition 2

Proof

We first guess that the joint value function $V_C$ of the two regions is of the form

$$\begin{aligned} V_{C}(P,t)=l_C(t)P+k_C(t), \end{aligned}$$

(24)

where $l_C(t)$ and $k_C(t)$ are the functions of t. Then, by substituting (10) and (24) into the HJB Eq. (16) satisfied by the joint value function, we obtain

$$\begin{aligned} \left( l^{'}_C(t)-(r+\theta )l_C(t)-(1+\beta )DP\right) P+k^{'}_C(t)-rk_C(t)+f_C(t)=0, \end{aligned}$$

(25)

where

$$\begin{aligned} f_C(t)=S(E_{10}+E_{20})+\frac{1}{2}\Big (A-S+l_C(t)\Big )^2+\frac{1}{2}\Big (\alpha A-S+l_C(t)\Big )^2+\frac{1}{2}\frac{\gamma C +1}{\gamma C}l^2_C(t). \end{aligned}$$

Since Eq. (25) holds for every $P\ge 0$, we have the following system of ordinary differential equations

$$\begin{aligned} \left\{ \begin{aligned}&l^{'}_C(t)-(r+\theta )l_C(t)-(1+\beta )D=0,\\&k^{'}_C(t)-rk_C(t)+f_C(t)=0. \end{aligned} \right. \end{aligned}$$

(26)

Solving this ODE system, we obtain

$$\begin{aligned} l(t)= & {} \Big (\frac{(1+\beta )D}{r+\theta }-g_1-g_2\Big )e^{-(r+\theta )(T-t)}-\frac{ (1+\beta )D}{r+\theta }, \\ k(t)= & {} G_{C}+H_{C}e^{-(r+\theta )(T-t)}+I_{C}e^{-2(r+\theta )(T-t)}+\frac{SE_{10}e^{-\rho _1 t}}{r+\rho _1}+\frac{SE_{20}e^{-\rho _2 t}}{r+\rho _2}\\&+\,\Big (g_1{\bar{P}}_1+g_2{\bar{P}}_2-I_{C}-H_{C}-G_{C}-\frac{SE_{10}e^{-\rho _1 T}}{r+\rho _1}-\frac{SE_{20}e^{-\rho _2 T}}{r+\rho _2}\Big )e^{-r(T-t)}, \end{aligned}$$

where

$$\begin{aligned} G_{C}=&\frac{1}{r}\Bigg (\frac{1}{2}(A-S)^2+\frac{1}{2}(\alpha A-S)^2-\frac{(1+\beta )D}{r+\theta }\Big ((1+\alpha )A-\frac{2\gamma C+\gamma +1}{\gamma C}S\Big )\\&+\Big (1+\frac{\gamma +1}{2\gamma C}\Big )\Big (\frac{(1+\beta )D}{r+\theta }\Big )^2+\frac{\gamma +1}{2\gamma C}S^2\Bigg ),\\ H_{C}=&\,\frac{1}{\theta }\Bigg (\frac{(1+\beta )D}{r+\theta }-g_1-g_2\Bigg )\Bigg (\frac{2\gamma C+\gamma +1}{\gamma C}S-(1+\alpha )A+2\Big (1+\frac{\gamma +1}{2\gamma C}\Big )\frac{(1+\beta )D}{r+\theta }\Bigg ),\\ I_{C}=&\,-\frac{1}{r+2\theta }\Bigg (1+\frac{\gamma +1}{2\gamma C}\Bigg )\Bigg (\frac{(1+\beta )D}{r+\theta }-g_1-g_2\Bigg )^2.\\ \end{aligned}$$

$\square$

E. Proof of Corollary 3

Proof

• (i) Comparing (11b) and (11d) with the two regions’ instantaneous quotas $E_{10}$ and $E_{20}$, we can obtain the above conclusions.

• (ii) From Eq. (11a), we have

  $$\begin{aligned} \frac{\partial V_{C}}{\partial P}= & {} \Big (\frac{(1+\beta )D}{r+\theta }-g_1-g_2\Big )e^{-(r+\theta )(T-t)}-\frac{ (1+\beta )D}{r+\theta }<0,\\ \frac{\partial V_{C}}{\partial E_{10}}= & {} \left( e^{\rho _1(T-t)}-e^{-r(T-t)}\right) \frac{Se^{-\rho _1T}}{r+\rho _1}>0,\\ \frac{\partial V_{C}}{\partial E_{20}}= & {} \left( e^{\rho _2(T-t)}-e^{-r(T-t)}\right) \frac{Se^{-\rho _2T}}{r+\rho _2}>0, \end{aligned}$$

  which completes the proof of (ii).

• (iii) Taking the derivative of (11a) with respect to S, we have

  $$\begin{aligned} \frac{\partial V_{C}}{\partial S}&=\frac{1}{r}\Big (\frac{2\gamma C+\gamma +1}{2\gamma C}S-(1+\alpha ) A\\ {}&\quad +\frac{(1+\beta )D}{r+\theta }\frac{2\gamma C+\gamma +1}{2\gamma C}\Big )(1-e^{-r(T-t)})\\&\quad +\frac{1}{\theta }\frac{2\gamma C+\gamma +1}{2\gamma C}e^{-r(T-t)}\Big (\frac{(1+\beta ) D}{r+\theta }-g_1-g_2\Big )\left( e^{-\theta (T-t)}-1\right) \\&\quad +\left( e^{\rho _1(T-t)}-e^{-r(T-t)}\right) \frac{E_{10}e^{-\rho _1T}}{r+\rho _1}+\left( e^{\rho _2(T-t)}-e^{-r(T-t)}\right) \frac{E_{20}e^{-\rho _2T}}{r+\rho _1}, \end{aligned}$$

  which implies that $V_{C}$ is a quadratic function of S at any time before T with the minimum point $\frac{1}{2}({\bar{\phi }}-\tau _{C1}(t)E_{10}-\tau _{C2}(t)E_{20})$. Also, combining the condition $0<S<{\bar{\eta }}$, which must be satisfied to keep the emissions nonnegative, and the condition $S>0$, we can conclude (iii). $\square$

F. Proof of Corollary 4

Proof

Substituting (11b), (11c), (11d) and (11e) into (3), we have the following ordinary differential equation for the dynamics of the pollution stock in the cooperative game:

$$\begin{aligned} \dot{P_C}+\theta P_C=(1+\alpha )A-\frac{2\gamma C+\gamma +1}{\gamma C}S+\frac{2\gamma C+\gamma +1}{\gamma C}l_C(t). \end{aligned}$$

Its solution gives the trajectory of the pollution stock (12). Moreover, taking its derivative with respect to t gives

$$\begin{aligned} P'_C(t)=e^{-\theta t}\left( Y_{C}(r+\theta )e^{-(r+\theta )(T-t)+\theta t}-\theta \left( P_0-X_{C}-Y_{C}e^{-(r+\theta )T}\right) \right) , \end{aligned}$$

from which the results follow directly. $\square$

G. Numerical Method for (16)

Since the structure of HJB equations arising from the noncooperative case is similar to that of the cooperative one, here we only discuss the latter to save space.

From the first-order optimality condition, we know that Eq. (16) can split into the following coupled equations:

$$\begin{aligned}&\frac{\partial V_C}{\partial t}+\left( E_{C1}^*(t)+E_{C2}^*(t)-a_{C1}^*(t)-a_{C2}^*(t)-\theta _PP\right) \frac{\partial V_C}{\partial P}+\frac{1}{2}\sigma ^2_PP^2\frac{\partial ^2 V_C}{\partial P^2}+\mu _SS\frac{\partial V_C}{\partial S} \nonumber \\&\quad +\,\frac{1}{2}\sigma ^2_SS^2\frac{\partial ^2 V_C}{\partial S^2}+\rho \sigma _P\sigma _SPS\frac{\partial ^2 V_C}{\partial P \partial S}-rV_C+F_C\left( P,S,E_{C1}^*(t),E_{C2}^*(t),t\right) =0,\nonumber \\ \end{aligned}$$

(28a)

$$\begin{aligned}&E^*_{C1}=A-S+\frac{\partial V_{C}}{\partial P},\quad a^*_{C1}=\frac{1}{C}\Big (S-\frac{\partial V_{C}}{\partial P}\Big ),\nonumber \\&E^*_{C2}=\alpha A-S+\frac{\partial V_{C}}{\partial P},\quad a^*_{C2}=\frac{1}{\gamma C}\left( S-\frac{\partial V_{C}}{\partial P}\right) . \end{aligned}$$

(28b)

1.1 The Fitted Finite Volume Method for Spatial Discretization

A defined mesh for $(P_{\min },P_{\max }) \times (S_{\min },S_{\max })$ is significant in the process of discretization. So, we first divide the intervals $I_{P}$ and $I_{S}$ into $N_P$ and $N_S$ subintervals, respectively:

$$\begin{aligned} I_{P_i}:=(P_i,P_{i+1}),\quad I_{S_j}:=(S_j,S_{j+1}),\quad i=0,1,\ldots ,N_P-1,\quad j=0,1,\ldots ,N_S-1, \end{aligned}$$

in which

$$\begin{aligned} P_{\min }=P_0<P_1<\cdots<P_{N_P}=P_{\max }\quad \text{ and } \quad S_{\min }=S_0<S_1<\cdots <S_{N_S}=S_{\max }. \end{aligned}$$

Thus, a mesh on $I_P \times I_S$, whose all mesh lines are perpendicular to the axes, is defined. Next we define another partition of $I_P \times I_S$ by letting

$$\begin{aligned} P_{i-\frac{1}{2}}=\frac{P_{i-1}+P_{i}}{2},\quad P_{i+\frac{1}{2}}=\frac{P_{i}+P_{i+1}}{2},\quad S_{j-\frac{1}{2}}=\frac{S_{j-1}+S_{j}}{2},\quad S_{j+\frac{1}{2}}=\frac{S_{j}+S_{j+1}}{2} \end{aligned}$$

for any $i=1,2,\ldots ,N_P-1$ and $j=1,2,\ldots ,N_S-1$. To keep completeness, we also define $P_{-\frac{1}{2}}=P_{\min }$, $P_{N_P+\frac{1}{2}}=P_{\max }$, $S_{-\frac{1}{2}}=S_{\min },$ and $S_{N_S+\frac{1}{2}}=S_{\max }$. The step sizes are defined by $h_{P_i}=P_{i+\frac{1}{2}}-P_{i-\frac{1}{2}}$ and $h_{S_j}=S_{j+\frac{1}{2}}-S_{j-\frac{1}{2}}$ for each $i=0,1,\ldots ,N_P$ and $j=0,1,\ldots ,N_S$.

Then, for the purpose of formulating finite volume scheme, we write Eq. (16) in the following divergence form:

$$\begin{aligned} -\frac{\partial V_C}{\partial t}-\nabla \cdot (A\nabla V_C+{\underline{b}}V_C)+cV_C=F_C, \end{aligned}$$

(29)

where

$$\begin{aligned} \begin{aligned} A&=\left( \begin{array}{c@{\quad }c@{\quad }c} a_{11} &{} a_{12} \\ a_{21} &{} a_{22} \\ \end{array} \right) =\left( \begin{array}{c@{\quad }c@{\quad }c} \frac{1}{2}\sigma ^2_{P}P^2 &{} \frac{1}{2}\rho \sigma _P\sigma _S PS \\ \frac{1}{2}\rho \sigma _P\sigma _S PS&{} \frac{1}{2}\sigma ^2_{S}S^2 \\ \end{array} \right) ,\\ {\underline{b}}&=\left( \begin{array}{c@{\quad }c@{\quad }c} b_1\\ b_2\\ \end{array} \right) =\left( \begin{array}{c@{\quad }c@{\quad }c} E_{C1}^*+E_{C2}^*-a_{C1}^*-a_{C2}^*-\theta _PP-\sigma ^2_{P}P-\frac{1}{2}\rho \sigma _P\sigma _S P\\ \mu _S S-\sigma ^2_{S}S-\frac{1}{2}\rho \sigma _P\sigma _S S\\ \end{array} \right) ,\\ c&=r+\mu _S+2\frac{\partial ^2 V_C}{\partial P^2}-\theta _P-\sigma ^2_{P}-\sigma ^2_{S}-\rho \sigma _P\sigma _S. \end{aligned} \end{aligned}$$

(30)

It follows from integrating Eq. (29) over ${\mathcal {R}}_{i,j}=[S_{i-\frac{1}{2}},S_{i+\frac{1}{2}}]\times [\delta _{j-\frac{1}{2}},\delta _{j+\frac{1}{2}}]$ and applying the midpoint quadrature rule to the resulting equation that

$$\begin{aligned} -\frac{\partial V_{C_{i,j}}}{\partial t}R_{i,j}-\int _{{\mathcal {R}}_{i,j}} \nabla \cdot (A \nabla V_C+{\underline{b}}V_C) \hbox {d}P \hbox {d}S + c_{i,j}V_{C_{i,j}}R_{i,j}=F_{C_{i,j}}R_{i,j} \end{aligned}$$

(31)

for $i=1,2,\ldots ,N_P-1$, $j=1,2,\ldots ,N_S-1$, where $R_{i,j}=(P_{i+\frac{1}{2}}-P_{i-\frac{1}{2}})\times (S_{j+\frac{1}{2}}-S_{j-\frac{1}{2}})$, $c_{i,j}=c(P_i,S_j,t)$, $V_{C_{i,j}}=V_C(P_i,S_j,t),$ and $F_{C_{i,j}}=F_C(P_i,S_j,E^*_{C1},E_{C2}^*,a^*_{C1},a_{C2}^*,t)$.

The approximation of the second term in Eq. (31) is the key and difficult point of this numerical method. According to the definition of flux $A \nabla V_C+{\underline{b}}V_C$ and integrating by parts, we have

$$\begin{aligned} \int _{{\mathcal {R}}_{i,j}} \nabla \cdot \left( A \nabla V_C+{\underline{b}}V_C\right) dSd\delta&=\int _{\partial {\mathcal {R}}_{i,j}} \left( A \nabla V_C+{\underline{b}}V_C\right) \cdot l \hbox {d}s \nonumber \\&=\int ^{\left( P_{i+\frac{1}{2}},S_{j+\frac{1}{2}}\right) }_{\left( P_{i+\frac{1}{2}},S_{j-\frac{1}{2}}\right) } \left( a_{11}\frac{\partial V_C}{\partial P}+a_{12}\frac{\partial V_C}{\partial S}+b_1V_C\right) \hbox {d}S \nonumber \\&\quad -\int ^{\left( P_{i-\frac{1}{2}},S_{j+\frac{1}{2}}\right) }_{\left( P_{i-\frac{1}{2}},S_{j-\frac{1}{2}}\right) } \left( a_{11}\frac{\partial V_C}{\partial P}+a_{12}\frac{\partial V_C}{\partial S}+b_1V_C\right) \hbox {d}S \nonumber \\&\quad +\int ^{\left( P_{i+\frac{1}{2}},S_{j+\frac{1}{2}}\right) }_{\left( P_{i-\frac{1}{2}},S_{j+\frac{1}{2}}\right) } \left( a_{21}\frac{\partial V_C}{\partial P}+a_{22}\frac{\partial V_C}{\partial S}+b_2V_C\right) \hbox {d}P \nonumber \\&\quad -\int ^{\left( P_{i-\frac{1}{2}},S_{j-\frac{1}{2}}\right) }_{\left( P_{i+\frac{1}{2}},S_{j-\frac{1}{2}}\right) } \left( a_{21}\frac{\partial V_C}{\partial P}+a_{22}\frac{\partial V_C}{\partial S}+b_2V_C\right) \hbox {d}P, \end{aligned}$$

(32)

where l denotes the unit vector outward-normal to $\partial {\mathcal {R}}_{i,j}$. We approximate the first integral of Eq. (32) by a constant:

$$\begin{aligned}&\int ^{\left( P_{i+\frac{1}{2}},S_{j+\frac{1}{2}}\right) }_{ \left( P_{i+\frac{1}{2}},S_{j-\frac{1}{2}}\right) } \left( a_{11}\frac{\partial V_C}{\partial P}+a_{12}\frac{\partial V_C}{\partial S}+b_1V_C\right) \hbox {d}S\\&\quad \approx \left( a_{11}\frac{\partial V_C}{\partial P}+a_{12}\frac{\partial V_C}{\partial S}+b_1V_C\right) \left| _{\left( P_{i+\frac{1}{2}},S_j\right) } \cdot h_{S_j}.\right. \end{aligned}$$

Now, we are in the position to derive the approximations to $(a_{11}\frac{\partial V_C}{\partial P}+a_{12}\frac{\partial V_C}{\partial S}+b_1V_C)$ at the midpoint, $(P_{i+\frac{1}{2}},S_j)$, of the interval $I_{P_i}$ for any $i=0,1,\ldots ,N_P-1$. To begin with, the term $a_{11}\frac{\partial V_C}{\partial P}+b_1V_C$ can be approximated by a constant, which means that its derivative equals zero, that is,

$$\begin{aligned}&\left( \frac{1}{2}\sigma _P^2P^2\frac{\partial V_C}{\partial P}+(E_{C1}^*+E_{C2}^*-a_{C1}^*-a_{C2}^*-\theta _PP-\sigma ^2_PP-\frac{1}{2}\rho \sigma _P\sigma _SP)V_C\right) '\nonumber \\&\quad \equiv \left( a P^2 \frac{\partial V_C}{\partial P}+b_1^{i+\frac{1}{2},j}V_C\right) '=0,\end{aligned}$$

(33a)

$$\begin{aligned}&V_C(P_i,S_j) = V_{C_{i,j}}, \quad V_C(P_{i+1},S_j) = V_{C_{i+1,j}}, \end{aligned}$$

(33b)

where $a=\frac{1}{2}\sigma _P^2$ and $b_1^{i+\frac{1}{2},j}=b_1(P_{i+\frac{1}{2}},S_j)$, $V_{C_{i,j}}$ and $V_{C_{i+1,j}}$ denote the values of $V_C$ at $(P_i,S_j)$ and $(P_{i+1},C_j)$, respectively. The above first-order ordinary differential equation can be solved by integrating both sides of Eq. (33a):

$$\begin{aligned} aP^2\frac{\partial V_C}{\partial P}+b_1^{i+\frac{1}{2},j}V_C=C_{1}, \end{aligned}$$

where $C_1$ is an arbitrary constant and can be determined by the boundary conditions Eq. (33b) as follows:

$$\begin{aligned} C_1=b_1^{i+\frac{1}{2},j}\frac{V_{C_{i+1,j}}e^{-\frac{\alpha _{i,j}}{P_{i+1}}}-V_{C_{i,j}}e^{-\frac{\alpha _{i,j}}{P_{i}}}}{e^{-\frac{\alpha _{i,j}}{P_{i+1}}}-e^{-\frac{\alpha _{i,j}}{P_{i}}}}, \end{aligned}$$

where $\alpha _{i,j}=\frac{b_1^{i+\frac{1}{2},j}}{a}$. Additionally, the derivative $\frac{\partial V_C}{\partial S}$ can be approximated by a forward difference

$$\begin{aligned} \frac{V_{C_{i,j+1}}-V_{C_{i,j}}}{h_{S_j}}. \end{aligned}$$

As a result, we have

$$\begin{aligned}&\left( a_{11}\frac{\partial V_C}{\partial P}+a_{12}\frac{\partial V_C}{\partial S}+b_1 V_C\right) |{(P_{i+\dfrac{1}{2}},S_j)} \cdot h_{S_j}\nonumber \\&\quad \approx \left( b_1^{i+\frac{1}{2},j}\frac{V_{C_{i+1,j}}e^{-\frac{\alpha _{i,j}}{P_{i+1}}}-V_{C_{i,j}}e^{-\frac{\alpha _{i,j}}{P_{i}}}}{e^{-\frac{\alpha _{i,j}}{P_{i+1}}}-e^{-\frac{\alpha _{i,j}}{P_{i}}}}+d_{i,j}\frac{V_{C_{i,j+1}}-V_{C_{i,j}}}{h_{S_{j}}}\right) \cdot h_{S_j}, \end{aligned}$$

(34)

where $d=\frac{1}{2}\rho \sigma _P\sigma _SPS$ and $d_{i,j}=d(P_i,S_j)$. Applying the similar method to the other three terms of Eq. (32), we get the following results:

$$\begin{aligned}&\left( a_{11}\frac{\partial V_C}{\partial P}+a_{12}\frac{\partial V_C}{\partial S}+b_1 V_C\right) |_{(P_{i-\frac{1}{2}},S_j)} \cdot h_{S_j}\nonumber \\&\quad \approx \left( b_1^{i-\frac{1}{2},j}\frac{V_{C_{i,j}}e^{-\frac{\alpha _{i-1,j}}{P_{i}}}-V_{C_{i-1,j}}e^{-\frac{\alpha _{i-1,j}}{P_{i-1}}}}{e^{-\frac{\alpha _{i-1,j}}{P_{i}}}-e^{-\frac{\alpha _{i-1,j}}{P_{i-1}}}}+d_{i,j}\frac{V_{C_{i,j+1}}-V_{C_{i,j}}}{h_{S_{j}}}\right) \cdot h_{S_j}, \end{aligned}$$

(35)

$$\begin{aligned}&\left( a_{21}\frac{\partial V_C}{\partial P}+a_{22}\frac{\partial V_C}{\partial S}+b_2 V_C\right) |_{(P_{i},S_{j+\frac{1}{2}})} \cdot h_{P_i}\nonumber \\&\quad \approx S_{j+\frac{1}{2}} \left( {\bar{b}}_{i,{j+\frac{1}{2}}}\frac{S_{j+1}^{{\bar{\alpha }}_{i,j}}V_{C_{i+1,j}}-S_{j}^{{\bar{\alpha }}_{i,j}}V_{C_{i,j}}}{S_{j+1}^{{\bar{\alpha }}_{i,j}}-S_{j}^{{\bar{\alpha }}_{i,j}}}+{\bar{d}}_{i,j}\frac{V_{C_{i,j+1}}-V_{C_{i,j}}}{h_{P_i}} \right) \cdot h_{P_i}, \end{aligned}$$

(36)

and

$$\begin{aligned}&\left( a_{21}\frac{\partial V_C}{\partial P}+a_{22}\frac{\partial V_C}{\partial S}+b_2 V_C\right) |_{(P_{i},S_{j-\frac{1}{2}})} \cdot h_{P_i}\nonumber \\&\quad \approx S_{j-\frac{1}{2}} \left( {\bar{b}}_{i,{j-\frac{1}{2}}}\frac{S_{j}^{{\bar{\alpha }}_{i,j-1}}V_{C_{i,j}}-S_{j-1}^{{\bar{\alpha }}_{i,j-1}}V_{C_{i,j-1}}}{S_{j}^{{\bar{\alpha }}_{i,j-1}}-S_{j-1}^{{\bar{\alpha }}_{i,j-1}}}+{\bar{d}}_{i,j}\frac{V_{C_{i,j+1}}-V_{C_{i,j}}}{h_{P_i}} \right) \cdot h_{P_i}, \end{aligned}$$

(37)

where ${\bar{\alpha }}_{i,j}=\frac{{\bar{b}}_{i,j+\frac{1}{2}}}{{\bar{a}}_j}$, ${\bar{a}}=\frac{1}{2}\sigma _S^2$, ${\bar{b}}=\mu -\sigma _S^2-\frac{1}{2}\rho \sigma _P\sigma _S$, and ${\bar{d}}_{i,j}=\frac{1}{2}\rho \sigma _P\sigma _S P_{i}$. Hence, we obtain the following equations by combining Eqs. (31), (32), and (34)–(37):

$$\begin{aligned}&-\frac{\partial V_{C_{i,j}}}{\partial t}R_{i,j}+e^{i,j}_{i-1,j}V_{C_{i-1,j}}+e^{i,j}_{i,j-1}V_{C_{i,j-1}}+e^{i,j}_{i,j}V_{C_{i,j}}+e^{i,j}_{i,j+1}V_{C_{i,j+1}}\nonumber \\&\quad +e^{i,j}_{i+1,j}V_{C_{i+1,j}}=F_{C_{i,j}}R_{i,j}, \end{aligned}$$

(38)

where

$$\begin{aligned} e^{i,j}_{i-1,j}= & {} -\,b_1^{i-\frac{1}{2},j}\frac{e^{-\frac{\alpha _{i-1,j}}{P_{i-1}}}h_{S_j}}{e^{-\frac{\alpha _{i-1,j}}{P_{i}}}-e^{-\frac{\alpha _{i-1,j}}{P_{i-1}}}}, \nonumber \\ e^{i,j}_{i,j-1}= & {} -\,S_{j-\frac{1}{2}}{\bar{b}}_{i,j-\frac{1}{2}}\frac{S^{{\bar{\alpha }}_{i,j-1}}_{j-1}h_{P_i}}{S^{{\bar{\alpha }}_{i,j-1}}_{j}-S^{{\bar{\alpha }}_{i,j-1}}_{j-1}}, \end{aligned}$$

(39)

$$\begin{aligned} e^{i,j}_{i,j}= & {} h_{S_j}\left( \frac{b_1^{i+\frac{1}{2},j}e^{-\frac{\alpha _{i,j}}{P_{i}}}}{e^{-\frac{\alpha _{i,j}}{P_{i+1}}}-e^{-\frac{\alpha _{i,j}}{P_{i}}}}+ \frac{b_1^{i-\frac{1}{2},j}e^{-\frac{\alpha _{i-1,j}}{P_{i}}}}{e^{-\frac{\alpha _{i-1,j}}{P_{i}}}-e^{-\frac{\alpha _{i-1,j}}{P_{i-1}}}}+{\bar{d}}_{i,j}\right) \nonumber \\&+\,h_{P_i}\left( S_{j+\frac{1}{2}}\frac{{\bar{b}}_{i,j+\frac{1}{2}}S_j^{{\bar{\alpha }}_{i,j}}}{S_{j+1}^{{\bar{\alpha }}_{i,j}}-S_j^{{\bar{\alpha }}_{i,j}}} +S_{j-\frac{1}{2}}\frac{{\bar{b}}_{i,j-\frac{1}{2}}S_j^{{\bar{\alpha }}_{i,j-1}}}{S_{j}^{\alpha _{i,j-1}}-S_{j-1}^{{\bar{\alpha }}_{i,j-1}}}\right) +c_{i,j}R_{i,j}, \end{aligned}$$

(40)

$$\begin{aligned} e^{i,j}_{i,j+1}= & {} -\,S_{j+\frac{1}{2}}{\bar{b}}_{i,j+\frac{1}{2}}\frac{S^{{\bar{\alpha }}_{i,j}}_{j+1}h_{P_i}}{S^{{\bar{\alpha }}_{i,j}}_{j+1}-S^{{\bar{\alpha }}_{i,j}}_{j}} ,\nonumber \\ e^{i,j}_{i+1,j}= & {} -\,b_1^{i+\frac{1}{2},j}\frac{e^{-\frac{\alpha _{i,j}}{P_{i+1}}}h_{S_j}}{e^{-\frac{\alpha _{i,j}}{P_{i+1}}}-e^{-\frac{\alpha _{i,j}}{P_{i}}}}-h_{S_j}{\bar{d}}_{i,j}, \end{aligned}$$

(41)

for $i=1,2,\ldots ,N_P-1$, $j=1,2,\ldots ,N_S-1$. The other elements $e^{i,j}_{m,n}$ equal zeros when $m\ne i-1,i,i+1$ and $n \ne j-1,j,j+1$. We can see that system (38) is an $(N_P-1)^2 \times (N_S-1)^2$ linear system of equations for

$$\begin{aligned} V_C=(V_{C_{1,1}},\ldots ,V_{C_{1,N_S-1}},V_{C_{2,1}},\ldots ,V_{C_{2,N_S-1}},\ldots ,V_{C_{N_P-1,1}},V_{C_{N_P-1,2}},\ldots ,V_{C_{N_P-1,N_S-1}})^{\top }. \end{aligned}$$

Note that $V_{C_{0,j}}(t),V_{C_{i,0}}(t),V_{C_{0,N_S}}(t),$ and $V_{C_{N_P,0}}(t)$ for $i=1,2,\ldots ,N_P$ and $j=1,2,\ldots ,N_S$ equal to the given boundary conditions. Obviously, the coefficient matrix of system (38) is pentadiagonal.

1.2 The Implicit Difference Method for Time Discretization

Next we embark on the time discretization of system (38). To this purpose, we first rewrite Eq. (38) as

$$\begin{aligned} -\frac{\partial V_{C_{i,j}}}{\partial t}R_{i,j}+D_{i,j}V_C=F_{C_{i,j}}R_{i,j}, \end{aligned}$$

(42)

where

$$\begin{aligned} D_{i,j}=(0,\ldots ,0,e^{i,j}_{i-1,j},0,\ldots ,0,e^{i,j}_{i,j-1},e^{i,j}_{i,j},e^{i,j}_{i,j+1},0,\ldots ,0,e^{i,j}_{i+1,j},0,\ldots ,0) \end{aligned}$$

for $i=1,2,\ldots ,N_P-1$ and $j=1,2,\ldots ,N_S-1$. We select $M-1$ points numbered from $t^1$ to $t^{M-1}$ between 0 and T,  and let $T=t^0$, $t^M=0$ to form a partition of time $T=t^0>t^1>\cdots >t^M=0$. Then, the full discrete form of Eq. (42) can be obtained by applying the two-level implicit time-stepping method with a splitting parameter $\theta \in [\frac{1}{2},1]$ to it:

$$\begin{aligned} \begin{aligned}&\left( \theta D\left( P,S,E_{C1}^*\left( t^{m+1}\right) ,E_{C2}^*\left( t^{m+1}\right) ,a_{C1}^*\left( t^{m+1}\right) ,a_{C2}^*\left( t^{m+1}\right) ,t^{m+1}\right) + G^{m}\right) V_C^{m+1}\\&\quad =\theta F_C\left( P,S,E_{C1}^*\left( t^{m+1}\right) ,E_{C2}^*\left( t^{m+1}\right) ,a_{C1}^*\left( t^{m+1}\right) ,a_{C2}^*\left( t^{m+1}\right) ,t^{m+1}\right) \\ {}&\quad \quad +\left( 1-\theta \right) F_C\left( P,S,E_{C1}^*\left( t^{m}\right) ,E_{C2}^*\left( t^{m}\right) ,a_{C1}^*\left( t^{m}\right) ,a_{C2}^*\left( t^{m}\right) ,t^{m}\right) \\ {}&\quad \quad +\left( G^m-\left( 1-\theta \right) D\left( P,S,E_{C1}^*\left( t^{m}\right) ,E_{C2}^*\left( t^{m}\right) ,a_{C1}^*\left( t^{m}\right) ,a_{C2}^*\left( t^{m}\right) ,t^{m}\right) \right) V_C^m, \end{aligned} \end{aligned}$$

(43)

where

$$\begin{aligned} V_C^m= & {} \left( V^m_{C_{1,1}},\ldots ,V^m_{C_{1,N_S-1}},V^m_{C_{2,1}},\ldots ,V^m_{C_{2,N_S-1}},\ldots ,V^m_{C_{N_P-1,1}},\ldots ,V^m_{C_{N_P-1,N_S-1}}\right) ^{\top },\nonumber \\ G^m= & {} \text{ diag }\,\left( -R_{1,1}/\Delta t^m,\ldots ,-R_{N_P-1,N_S-1}/\Delta t^m\right) ^{\top }, \end{aligned}$$

(44)

for $m=0,1,\ldots ,M-1$. Note that $\Delta t^m=t^{m+1}-t^m<0$, and $V_C^m$ denotes the approximation of $V_C$ at $t=t^m$. Particularly, when we set $\theta =\frac{1}{2}$, the scheme (43) becomes the famous Crank–Nicolson scheme and is second-order accuracy; when we set $\theta =1$, the scheme (43) becomes the backward Euler scheme and is first-order accuracy.

1.3 Decoupling of the System

In the above discussion, we have assumed that the control variables $E_{C1}^*$, $a_{C1}^*$, $E_{C2}^*$, and $a_{C2}^*$ are known. However, we can see from (43) that they are coupled with $V_C$ when $\theta \ne 0$. To deal with this dilemma, we replace $E_{C1}^*(t^{m+1})$, $a_{C1}^*(t^{m+1})$, $E_{C2}^*(t^{m+1})$, and $a_{C2}^*(t^{m+1})$ by $E_{C1}^*(t^{m})$, $a_{C1}^*(t^{m})$, $E_{C2}^*(t^{m})$, and $a_{C2}^*(t^{m})$, respectively. This method should be reasonable because the control variables are just replaced by their values in the previous time step. The error is small if $\Delta t^m$ is sufficiently small. The resulting system corresponding to (43) is as follows:

$$\begin{aligned} \begin{aligned}&\left( \theta D\left( P,S,E_{C1}^*\left( t^{m}\right) ,E_{C2}^*\left( t^{m}\right) ,a_{C1}^*\left( t^{m}\right) ,a_{C2}^*\left( t^{m}\right) ,t^{m+1}\right) + G^{m}\right) V_C^{m+1}\\&\quad =\theta F\left( P,S,E_{C1}^*\left( t^{m}\right) ,E_{C2}^*\left( t^{m}\right) ,a_{C1}^*\left( t^{m}\right) ,a_{C2}^*\left( t^{m}\right) ,t^{m+1}\right) \\&\quad \quad +\left( 1-\theta \right) F\left( P,S,E_{C1}^*\left( t^{m}\right) ,E_{C2}^*\left( t^{m}\right) ,a_{C1}^*\left( t^{m}\right) ,a_{C2}^*\left( t^{m}\right) ,t^{m}\right) \\&\quad \quad +\left( G^m-\left( 1-\theta \right) D\left( P,S,E_{C1}^*\left( t^{m}\right) ,E_{C2}^*\left( t^{m}\right) ,a_{C1}^*\left( t^{m}\right) ,a_{C2}^*\left( t^{m}\right) ,t^{m}\right) \right) V_C^m. \end{aligned} \end{aligned}$$

(45)