Semi-implicit FEM for the valuation of American options under the Heston model

Abstract

In this paper, we present an efficient numerical method for the valuation of American put options under the Heston model. Firstly, by adding a penalty term, the pricing model, which is a linear complementary problem on an unbounded domain, is transformed into a nonlinear parabolic partial differential equation. Then, the perfectly matched layer technique is applied to truncate the solution domain. To deal with the challenge arising from the nonlinearity, a semi-implicit finite element scheme is utilized, which makes the numerical implementation feasible and efficient. Furthermore, we shall prove that the full-discrete matrix is an M-matrix under some moderate assumptions, which implies the nonnegativity of the numerical solutions. Finally, some numerical simulations are carried out to test the performance of the proposed method.

Appendix: The proof of Theorem 5

In this part, we present the proof of Theorem 5 in detail. Without loss of generality, we apply the equidistant partition with $\Delta x=\Delta y=h$, and only take the interior point for example to discuss the properties of the discretized matrix. For the other cases, we can obtain the same conclusions by applying the similar techniques. The partition elements around the interior point k are given in Fig. 5.

Fig. 5

The spatial partition element

Since the weak form (11) can be reorganized as

$$\begin{aligned}&(1+r\Delta \tau )(u^{m+1}_h,\omega _h)+\Delta \tau \Big \{\frac{1}{2}\Big (y(u^{m+1}_h)_x,(\omega _h)_x\Big )+\rho \xi \Big (y(u^{m+1}_h)_x,(\omega _h)_y\Big )\\&\qquad +\frac{\xi ^2}{2}(y(u^{m+1}_h)_y,(\omega _h)_y)+\Big (\Big (\frac{y}{2}-\mu -y\varpi +\rho \xi \Big )(u^{m+1}_h)_x,\omega _h\Big )\\&\qquad +\Big (\big (\kappa (y-\gamma )-\rho \xi y\varpi +\frac{\xi ^2}{2}\big )(u^{m+1}_h)_y,\omega _h\Big )\\&\qquad +\Big (\big (-\frac{y}{2}({\varpi }'+\varpi ^2) +(\frac{y}{2}-\mu )\varpi \big )u^{m+1}_h,\omega _h\Big )\\&\qquad -\rho \xi y_{\max }\int _{-L-\delta }^{L+\delta }(u^{m+1}_h)_x(\tau ,x,y_{\max })\omega _h(x,y_{\max })\mathrm{d}x\Big \}\\&\quad =(u^{m}_h-\Delta \tau f(u^{m}_h),\omega _h),\quad \forall ~ \omega _h\in V_h^0(\Sigma _{pml}), \end{aligned}$$

the k-th equation of the discretized system (11) could be rewritten as

$$\begin{aligned} A(k,k)u_k^{m+1}+\sum _{i=1}^6A(k,k_i)u_{k_i}^{m+1}=F_k^{m+1}. \end{aligned}$$

Now, we will verify that the discretized matrix is an M-matrix. For simplify, we only discuss the case where $\delta =0$, other cases are similar. For the first triangle in Fig. 5, the basis functions can be defined as:

$$\begin{aligned} \begin{aligned}&L_k^{(1)}=\frac{h}{2S}(-x+x_{i+1}),\\&L_{k_2}^{(1)}=\frac{h}{2S}(y-y_j),\\&L_{k_1}^{(1)}=\frac{h}{2S}[(x-x_i)-(y-y_j)]. \end{aligned} \end{aligned}$$

Then we have

$$\begin{aligned} A^{(1)}(k,k) =~&(1+r\Delta \tau )(L_k^{(1)},L_k^{(1)})+\Delta \tau \Big \{\frac{1}{2}\Big (y({L_k^{(1)}}\big )_x,({L_k^{(1)}})_x\Big ) +\rho \xi \Big (y({L_k^{(1)}})_{x},({L_k^{(1)}})_y\Big )&\\&+\frac{\xi ^2}{2}\Big (y({L_k^{(1)}})_{y},({L_k^{(1)}})_y\Big ) +\frac{1}{2}\Big (y({L_k^{(1)}})_{x},{L_k^{(1)}}\Big )+\kappa \Big (y({L_k^{(1)}})_{y},{L_k^{(1)}}\Big )&\\&+(\rho \xi -\mu )\Big (({L_k^{(1)}})_{x},{L_k^{(1)}}\Big ) +\Big (\frac{\xi ^2}{2}-\kappa \gamma \Big )\Big (({L_k^{(1)}})_{y},{L_k^{(1)}}\Big )\Big \},\\ =~&\frac{S}{6}(1+r\Delta \tau )+\Delta \tau \Big (\frac{h^4}{48S^2}(h+3y_j)-\frac{h^5}{192S^2}(h+4y_{j})-\frac{h^5}{24S^2}(\rho \xi -\mu )\Big ),&\\ A^{(1)}(k,k_1)=~&(1+r\Delta \tau )(L_{k_1}^{(1)},L_k^{(1)})+\Delta \tau \Big \{\frac{1}{2}\big (y(L_{k_1}^{(1)}\big )_x,({L_k^{(1)}})_x\big ) +\rho \xi \Big (y(L_{k_1}^{(1)})_{x},({L_k^{(1)}})_y\Big )&\\&+\frac{\xi ^2}{2}\Big (y(L_{k_1}^{(1)})_{y},({L_k^{(1)}})_y\Big ) +\frac{1}{2}\Big (y(L_{k_1}^{(1)})_{x},{L_k^{(1)}}\Big )+\kappa \Big (y(L_{k_1}^{(1)})_{y},{L_k^{(1)}}\Big ),&\\&+(\rho \xi -\mu )\Big (({L_{k_1}})_{x},{L_k^{(1)}}\Big ) +\Big (\frac{\xi ^2}{2}-\kappa \gamma \Big )\Big ((L_{k_1}^{(1)})_{y},{L_k^{(1)}}\Big )\Big \}&\\ =~&\frac{S}{12}(1+r\Delta \tau )+\Delta \tau \Big (-\frac{h^4}{48S^2}(h+3y_j)+\frac{h^5}{192S^2}(1-2\kappa )(h+4y_{j})&\\&+\frac{h^5}{24S^2}(\rho \xi -\mu -\frac{\xi ^2}{2}+\kappa \gamma )\Big ),&\\ A^{(1)}(k,k_2)=~&(1+r\Delta \tau )(L_{k_2}^{(1)},L_k^{(1)})+\Delta \tau \Big \{\frac{1}{2}\big (y(L_{k_2}^{(1)}\big )_x,({L_k^{(1)}})_x\big ) +\rho \xi \Big (y({L_{k_2}^{(1)})_{x}},({L_k^{(1)}})_y\Big )&\\&+\frac{\xi ^2}{2}\Big (y(L_{k_2}^{(1)})_{y},({L_k^{(1)}})_y\Big ) +\frac{1}{2}\Big (y(L_{k_2}^{(1)})_{x},{L_k^{(1)}}\Big )+\kappa \Big (y(L_{k_2}^{(1)})_{y},{L_k^{(1)}}\Big )&\\&+(\rho \xi -\mu )\Big ((L_{k_2}^{(1)})_{x},{L_k^{(1)}}\Big ) +\Big (\frac{\xi ^2}{2}-\kappa \gamma \Big )\Big ((L_{k_2}^{(1)})_{y},{L_k^{(1)}}\Big )\Big \}&\\ =~&\frac{S}{12}(1+r\Delta \tau )+\Delta \tau \Big (\frac{\kappa h^5}{96S^2}(h+4y_{j})+\frac{h^5}{24S^2}(\frac{\xi ^2}{2}-\kappa \gamma )\Big ). \end{aligned}$$

For the second triangle, the basis functions can be defined as:

$$\begin{aligned} \begin{aligned}&L_k^{(2)}=\frac{h}{2S}(-y+y_{j+1}),\\&L_{k_2}^{(2)}=\frac{h}{2S}(x-x_i),\\&L_{k_3}^{(2)}=\frac{h}{2S}[(y-y_j)-(x-x_i)], \end{aligned} \end{aligned}$$

and the corresponding coefficients are

$$\begin{aligned} A^{(2)}(k,k)=~&(1+r\Delta \tau )(L_k^{(2)},L_k^{(2)})+\Delta \tau \Big \{\frac{1}{2}\big (y({L_k^{(2)}}\big )_x,({L_k^{(2)}})_x\big ) +\rho \xi \Big (y({L_k^{(2)}})_{x},({L_k^{(2)}})_y\Big )&\\&+\frac{\xi ^2}{2}(y({L_k^{(2)}})_{y},({L_k^{(2)}})_y) +\frac{1}{2}\Big (y({L_k^{(2)}})_{x},{L_k^{(2)}}\Big )+\kappa \Big (y({L_k^{(2)}})_{y},{L_k^{(2)}}\Big )&\\&+(\rho \xi -\mu )\Big (({L_k^{(2)}})_{x},{L_k^{(2)}}\Big ) +\Big (\frac{\xi ^2}{2}-\kappa \gamma \Big )\Big (({L_k^{(2)}})_{y},{L_k^{(2)}}\Big )\Big \}&\\ =~&\frac{S}{6}(1+r\Delta \tau )+\Delta \tau \Big (\frac{\xi ^2h^4}{48S^2}(2h+3y_j)-\frac{\kappa h^5}{48S^2}(h+2y_{j})-\frac{h^5}{24S^2}(\frac{\xi ^2}{2}-\kappa \gamma )\Big ),&\\ A^{(2)}(k,k_2)=~&(1+r\Delta \tau )(L_{k_2}^{(2)},L_k^{(2)})+\Delta \tau \Big \{\frac{1}{2}\big (y(L_{k_2}^{(2)}\big )_x,({L_k^{(2)}})_x\big )+\rho \xi (y(L_{k_2}^{(2)})_{x},({L_k^{(2)}})_y)&\\&+\frac{\xi ^2}{2}\Big (y(L_{k_2}^{(2)})_{y},({L_k^{(2)}})_y\Big ) +\frac{1}{2}\Big (y(L_{k_2}^{(2)})_{x},{L_k^{(2)}}\Big )+\kappa \Big (y(L_{k_2}^{(2)})_{y},{L_k^{(2)}}\Big )\\&+(\rho \xi -\mu )\Big (({L_{k_2}^{(2)}})_{x},{L_k^{(2)}}\Big )+\Big (\frac{\xi ^2}{2}-\kappa \gamma \Big )\Big ((L_{k_2}^{(2)})_{y},{L_k^{(2)}}\Big )\Big \}\\ =~&\frac{S}{12}(1+r\Delta \tau )\\&+\Delta \tau \Big (-\frac{\rho \xi h^4}{24S^2}(2h+3y_j)+\frac{h^5}{96S^2}(h+2y_{j})+\frac{h^5}{24S^2}(\rho \xi -\mu )\Big ),\\ A^{(2)}(k,k_3)=~&(1+r\Delta \tau )(L_{k_3}^{(2)},L_k^{(2)})+\Delta \tau \Big \{\frac{1}{2}\big (y(L_{k_3}^{(2)}\big )_x,({L_k^{(2)}})_x\big ) +\rho \xi (y({L_{k_3}^{(2)})_{x}},({L_k^{(2)}})_y)\\&+\frac{\xi ^2}{2}\Big (y(L_{k_3}^{(2)})_{y},({L_k^{(2)}})_y\Big ) +\frac{1}{2}\Big (y(L_{k_3}^{(2)})_{x},{L_k^{(2)}}\Big )+\kappa \Big (y(L_{k_3}^{(2)})_{y},{L_k^{(2)}}\Big )&\\&+(\rho \xi -\mu )\Big ((L_{k_3}^{(2)})_{x},{L_k^{(2)}}\Big ) +\Big (\frac{\xi ^2}{2}-\kappa \gamma \Big )\Big ((L_{k_3}^{(2)})_{y},{L_k^{(2)}}\Big )\Big \}\\ =~&\frac{S}{12}(1+r\Delta \tau )+\Delta \tau \Big (\frac{\xi h^4}{48S^2}(2h+3y_{j})(2\rho -\xi )-\frac{h^5}{96S^2}(h+2y_j)(1-2\kappa )\\&+\frac{h^5}{24S^2}(\mu -\rho \xi +\frac{\xi ^2}{2}-\kappa \gamma )\Big ). \end{aligned}$$

For the third triangle, the basis functions can be defined as:

$$\begin{aligned} \begin{aligned}&L_k^{(3)}=\frac{h}{2S}[(x-x_i)-(y-y_{j+1})],\\&L_{k_3}^{(3)}=\frac{h}{2S}(y-y_j),\\&L_{k_4}^{(3)}=\frac{h}{2S}(-x+x_i), \end{aligned} \end{aligned}$$

and the corresponding coefficients are

$$\begin{aligned} \begin{aligned} A^{(3)}(k,k)=~&(1+r\Delta \tau )(L_k^{(3)},L_k^{(3)})+\Delta \tau \Big \{\frac{1}{2}\big (y({L_k^{(3)}}\big )_x,({L_k^{(3)}})_x\big ) +\rho \xi (y({L_k^{(3)}})_{x},({L_k^{(3)}})_y)\\&+\frac{\xi ^2}{2}\Big (y({L_k^{(3)}})_{y},({L_k^{(3)}})_y\Big ) +\frac{1}{2}\Big (y({L_k^{(3)}})_{x},{L_k^{(3)}}\Big )+\kappa \Big (y({L_k^{(3)}})_{y},{L_k^{(3)}}\Big )\\&+(\rho \xi -\mu )\Big (({L_k^{(3)}})_{x},{L_k^{(3)}}\Big ) +\Big (\frac{\xi ^2}{2}-\kappa \gamma \Big )\Big (({L_k^{(3)}})_{y},{L_k^{(3)}}\Big )\Big \}\\ =~&\frac{S}{6}(1+r\Delta \tau )+\Delta \tau \Big (\frac{h^4}{48S^2}(h+3y_j)(1-2\rho \xi +\xi ^2)\\&+\frac{h^5}{192S^2}(1-2\kappa )(h+4y_{j}) +\frac{h^5}{24S^2}(\rho \xi -\mu -\frac{\xi ^2}{2}+\kappa \gamma )\Big ),\\ A^{(3)}(k,k_3)=~&(1+r\Delta \tau )(L_{k_3}^{(3)},L_k^{(3)})+\Delta \tau \Big \{\frac{1}{2}\big (y(L_{k_3}^{(3)}\big )_x,({L_k^{(3)}})_x\big )\\&+\rho \xi (y(L_{k_3}^{(3)})_{x},({L_k^{(3)}})_y)\\&+\frac{\xi ^2}{2}(y(L_{k_3}^{(3)})_{y},({L_k^{(3)}})_y)\\&+\frac{1}{2}(y(L_{k_3}^{(3)})_{x},{L_k^{(3)}})+\kappa (y(L_{k_3}^{(3)})_{y},{L_k^{(3)}})\\&+(\rho \xi -\mu )(({L_{k_3}})_{x},{L_k^{(3)}})\\&+(\frac{\xi ^2}{2}-\kappa \gamma )((L_{k_3}^{(3)})_{y},{L_k^{(2)}})\Big \}\\ =~&\frac{S}{12}(1+r\Delta \tau )+\Delta \tau \Big (-\frac{\xi ^2 h^4}{48S^2}(h+3y_j)\\&+\frac{\kappa h^5}{96S^2}(h+4y_{j})+\frac{h^5}{24S^2}(\frac{\xi ^2}{2}-\kappa \gamma )\Big ),\\ A^{(3)}(k,k_4)=~&(1+r\Delta \tau )(L_{k_4}^{(3)},L_k^{(3)})+\Delta \tau \Big \{\frac{1}{2}\big (y(L_{k_4}^{(3)}\big )_x,({L_k^{(3)}})_x\big ) +\rho \xi \Big (y({L_{k_4}^{(3)})_{x}},({L_k^{(3)}})_y\Big )&\\&+\frac{\xi ^2}{2}\Big (y(L_{k_4}^{(3)})_{y},({L_k^{(3)}})_y\Big ) +\frac{1}{2}\Big (y(L_{k_4}^{(3)})_{x},{L_k^{(3)}}\Big )+\kappa \Big (y(L_{k_4}^{(3)})_{y},{L_k^{(3)}}\Big )&\\&+(\rho \xi -\mu )\Big ((L_{k_4}^{(3)})_{x},{L_k^{(3)}}\Big ) +\Big (\frac{\xi ^2}{2}-\kappa \gamma \Big )\Big ((L_{k_4}^{(3)})_{y},{L_k^{(3)}}\Big )\Big \}&\\ =~&\frac{S}{12}(1+r\Delta \tau )+\Delta \tau \Big (-\frac{ h^4}{48S^2}(h+3y_{j})(1-2\rho \xi )&\\&-\frac{h^5}{192S^2}(h+4y_j)+\frac{h^5}{24S^2}(\mu -\rho \xi )\Big ). \end{aligned} \end{aligned}$$

For the fourth triangle, the basis functions can be defined as:

$$\begin{aligned} \begin{aligned}&L_k^{(4)}=\frac{h}{2S}(x-x_{i-1}),\\&L_{k_5}^{(4)}=\frac{h}{2S}(-y+y_j),\\&L_{k_4}^{(4)}=\frac{h}{2S}[(y-y_j)-(x-x_i)], \end{aligned} \end{aligned}$$

and the corresponding coefficients are

$$\begin{aligned} \begin{aligned} A^{(4)}(k,k)=~&(1+r\Delta \tau )(L_k^{(4)},L_k^{(4)})+\Delta \tau \Big \{\frac{1}{2}\big (y({L_k^{(4)}}\big )_x,({L_k^{(4)}})_x\big ) +\rho \xi \Big (y({L_k^{(4)}})_{x},({L_k^{(4)}})_y\Big )&\\&+\frac{\xi ^2}{2}\Big (y({L_k^{(4)}})_{y},({L_k^{(4)}})_y\Big ) +\frac{1}{2}\Big (y({L_k^{(4)}})_{x},{L_k^{(4)}}\Big )+\kappa \Big (y({L_k^{(4)}})_{y},{L_k^{(4)}}\Big )&\\&+(\rho \xi -\mu )\Big (({L_k^{(4)}})_{x},{L_k^{(4)}}\Big ) +\Big (\frac{\xi ^2}{2}-\kappa \gamma \Big )\Big (({L_k^{(4)}})_{y},{L_k^{(4)}}\Big )\Big \}&\\ =~&\frac{S}{6}(1+r\Delta \tau )+\Delta \tau \Big (\frac{h^4}{48S^2}(3y_j-h)+\frac{h^5}{192S^2}(4y_{j}-h)+\frac{h^5}{24S^2}(\rho \xi -\mu )\Big ),&\\ A^{(4)}(k,k_4)=~&(1+r\Delta \tau )(L_{k_4}^{(4)},L_k^{(4)})+\Delta \tau \Big \{\frac{1}{2}\big (y(L_{k_4}^{(4)}\big )_x,({L_k^{(4)}})_x\big ) +\rho \xi (y(L_{k_4}^{(4)})_{x},({L_k^{(4)}})_y)&\\&+\frac{\xi ^2}{2}\Big (y(L_{k_4}^{(4)})_{y},({L_k^{(4)}})_y\Big ) +\frac{1}{2}\Big (y(L_{k_4}^{(4)})_{x},{L_k^{(4)}}\Big )+\kappa \Big (y(L_{k_4}^{(4)})_{y},{L_k^{(4)}}\Big )&\\&+(\rho \xi -\mu )\Big (({L_{k_4}})_{x},{L_k^{(4)}}\Big ) +(\frac{\xi ^2}{2}-\kappa \gamma )((L_{k_4}^{(4)})_{y},{L_k^{(4)}})\Big \}&\\ =~&\frac{S}{12}(1+r\Delta \tau )+\Delta \tau \Big (\frac{h^4}{48S^2}(h-3y_j)+\frac{h^5}{192S^2}(4y_{j}-h)(2\kappa -1)&\\&+\frac{h^5}{24S^2}(\frac{\xi ^2}{2}-\kappa \gamma -\rho \xi +\mu )\Big ),&\\ A^{(4)}(k,k_5)=~&(1+r\Delta \tau )(L_{k_5}^{(4)},L_k^{(4)})+\Delta \tau \Big \{\frac{1}{2}\big (y(L_{k_5}^{(4)}\big )_x,({L_k^{(4)}})_x\big ) +\rho \xi (y({L_{k_5}^{(4)})_{x}},({L_k^{(4)}})_y)&\\&+\frac{\xi ^2}{2}\Big (y(L_{k_5}^{(4)})_{y},({L_k^{(4)}})_y\Big ) +\frac{1}{2}\Big (y(L_{k_5}^{(4)})_{x},{L_k^{(4)}}\Big )+\kappa \Big (y(L_{k_5}^{(4)})_{y},{L_k^{(4)}}\Big )&\\&+(\rho \xi -\mu )\Big ((L_{k_5}^{(4)})_{x},{L_k^{(4)}}\Big ) +\Big (\frac{\xi ^2}{2}-\kappa \gamma \Big )\Big ((L_{k_5}^{(4)})_{y},{L_k^{(4)}}\Big )\Big \}&\\ =~&\frac{S}{12}(1+r\Delta \tau )+\Delta \tau \Big (\frac{ \kappa h^5}{96S^2}(h-4y_{j})-\frac{h^5}{24S^2}(\frac{\xi ^2}{2}-\kappa \gamma )\Big ). \end{aligned} \end{aligned}$$

For the fifth triangle, the basis functions can be defined as:

$$\begin{aligned} \begin{aligned}&L_k^{(5)}=\frac{h}{2S}(y-y_{j-1}),\\&L_{k_5}^{(5)}=\frac{h}{2S}(-x+x_i),\\&L_{k_6}^{(5)}=\frac{h}{2S}[(x-x_i)-(y-y_j)], \end{aligned} \end{aligned}$$

and the corresponding coefficients are

$$\begin{aligned} A^{(5)}(k,k)=~&(1+r\Delta \tau )(L_k^{(5)},L_k^{(5)})+\Delta \tau \big \{\frac{1}{2}\big (y({L_k^{(5)}}\big )_x,({L_k^{(5)}})_x\big ) +\rho \xi (y({L_k^{(5)}})_{x},({L_k^{(5)}})_y)&\\&+\frac{\xi ^2}{2}(y({L_k^{(5)}})_{y},({L_k^{(5)}})_y) +\frac{1}{2}(y({L_k^{(5)}})_{x},{L_k^{(5)}})+\kappa (y({L_k^{(5)}})_{y},{L_k^{(5)}})&\\&+(\rho \xi -\mu )(({L_k^{(5)}})_{x},{L_k^{(5)}}) +(\frac{\xi ^2}{2}-\kappa \gamma )(({L_k^{(5)}})_{y},{L_k^{(5)}})\big \}&\\ =~&\frac{S}{6}(1+r\Delta \tau )+\Delta \tau \Big (\frac{\xi ^2h^4}{48S^2}(3y_j-2h)+\frac{\kappa h^5}{48S^2}(2y_{j}-h)+\frac{h^5}{24S^2}(\frac{\xi ^2}{2}-\kappa \gamma )\Big ),&\\ A^{(5)}(k,k_5)=~&(1+r\Delta \tau )(L_{k_5}^{(5)},L_k^{(5)})+\Delta \tau \Big \{\frac{1}{2}\big (y(L_{k_5}^{(5)}\big )_x,({L_k^{(5)}})_x\big ) +\rho \xi (y(L_{k_5}^{(5)})_{x},({L_k^{(5)}})_y)&\\&+\frac{\xi ^2}{2}(y(L_{k_5}^{(5)})_{y},({L_k^{(5)}})_y) +\frac{1}{2}(y(L_{k_5}^{(5)})_{x},{L_k^{(5)}})+\kappa (y(L_{k_5}^{(5)})_{y},{L_k^{(5)}})&\\&+(\rho \xi -\mu )(({L_{k_5}})_{x},{L_k^{(5)}}) +(\frac{\xi ^2}{2}-\kappa \gamma )((L_{k_5}^{(5)})_{y},{L_k^{(5)}})\Big \}&\\ =~&\frac{S}{12}(1+r\Delta \tau )\\&+\Delta \tau \Big (-\frac{\rho \xi h^4}{24S^2}(3y_j-2h)-\frac{h^5}{96S^2}(2y_{j}-h)-\frac{h^5}{24S^2}(\rho \xi -\mu )\Big ),&\\ A^{(5)}(k,k_6)=~&(1+r\Delta \tau )(L_{k_6}^{(5)},L_k^{(5)})+\Delta \tau \Big \{\frac{1}{2}\big (y(L_{k_6}^{(5)}\big )_x,({L_k^{(5)}})_x\big ) +\rho \xi (y({L_{k_6}^{(5)})_{x}},({L_k^{(5)}})_y)&\\&+\frac{\xi ^2}{2}(y(L_{k_6}^{(5)})_{y},({L_k^{(5)}})_y) +\frac{1}{2}(y(L_{k_6}^{(5)})_{x},{L_k^{(5)}})+\kappa (y(L_{k_6}^{(5)})_{y},{L_k^{(5)}})&\\&+(\rho \xi -\mu )((L_{k_6}^{(5)})_{x},{L_k^{(5)}}) +(\frac{\xi ^2}{2}-\kappa \gamma )((L_{k_6}^{(5)})_{y},{L_k^{(5)}})\Big \}&\\ =~&\frac{S}{12}(1+r\Delta \tau )+\Delta \tau \Big (\frac{h^4}{48S^2}(2\rho \xi -\xi ^2)(3y_{j}-2h)&\\&+\frac{h^5}{96S^2}(1-2\kappa )(2y_j-h) +\frac{h^5}{24S^2}(\rho \xi -\mu -\frac{\xi ^2}{2}+\kappa \gamma )\Big ). \end{aligned}$$

For the sixth triangle, the basis functions can be defined as:

$$\begin{aligned} \begin{aligned}&L_k^{(6)}=\frac{h}{2S}[(y-y_{j-1})-(x-x_i)],\\&L_{k_6}^{(6)}=\frac{h}{2S}(-y+y_j),\\&L_{k_1}^{(6)}=\frac{h}{2S}(x-x_i), \end{aligned} \end{aligned}$$

and the corresponding coefficients are

$$\begin{aligned} A^{(6)}(k,k)=~&(1+r\Delta \tau )(L_k^{(6)},L_k^{(6)})+\Delta \tau \big \{\frac{1}{2}\big (y({L_k^{(6)}}\big )_x,({L_k^{(6)}})_x\big ) +\rho \xi (y({L_k^{(6)}})_{x},({L_k^{(6)}})_y)&\\&+\frac{\xi ^2}{2}(y({L_k^{(6)}})_{y},({L_k^{(6)}})_y) +\frac{1}{2}(y({L_k^{(6)}})_{x},{L_k^{(6)}})+\kappa (y({L_k^{(6)}})_{y},{L_k^{(6)}})&\\&+(\rho \xi -\mu )(({L_k^{(6)}})_{x},{L_k^{(6)}}) +(\frac{\xi ^2}{2}-\kappa \gamma )(({L_k^{(6)}})_{y},{L_k^{(6)}})\big \}&\\ =~&\frac{S}{6}(1+r\Delta \tau )+\Delta \tau \Big (\frac{h^4}{48S^2}(1-2\rho \xi +\xi ^2)(3y_j-h)&\\&+\frac{h^5}{192S^2}(4y_j-h)(2\kappa -1) +\frac{h^5}{24S^2}(\frac{\xi ^2}{2}-\kappa \gamma -\rho \xi +\mu )\Big ),&\\ A^{(6)}(k,k_6)=~&(1+r\Delta \tau )(L_{k_6}^{(6)},L_k^{(6)})+\Delta \tau \Big \{\frac{1}{2}\big (y(L_{k_6}^{(6)}\big )_x,({L_k^{(6)}})_x\big ) +\rho \xi (y(L_{k_6}^{(6)})_{x},({L_k^{(6)}})_y)&\\&+\frac{\xi ^2}{2}(y(L_{k_6}^{(6)})_{y},({L_k^{(6)}})_y) +\frac{1}{2}(y(L_{k_6}^{(6)})_{x},{L_k^{(6)}})+\kappa (y(L_{k_6}^{(6)})_{y},{L_k^{(6)}})&\\&+(\rho \xi -\mu )(({L_{k_6}^{(6)}})_{x},{L_k^{(6)}}) +(\frac{\xi ^2}{2}-\kappa \gamma )((L_{k_6}^{(6)})_{y},{L_k^{(6)}})\Big \}&\\ =~&\frac{S}{12}(1+r\Delta \tau )+\Delta \tau \Big (\frac{\xi ^2 h^4}{48S^2}(h-3y_j)&\\&-\frac{\kappa h^5}{96S^2}(4y_{j}-h)-\frac{h^5}{24S^2}(\frac{\xi ^2}{2}-\kappa \gamma )\Big ),&\\ A^{(6)}(k,k_1)=~&(1+r\Delta \tau )(L_{k_1}^{(6)},L_k^{(6)})+\Delta \tau \Big \{\frac{1}{2}\big (y(L_{k_1}^{(6)}\big )_x,({L_k^{(6)}})_x\big ) +\rho \xi (y({L_{k_1}^{(6)})_{x}},({L_k^{(6)}})_y)&\\&+\frac{\xi ^2}{2}(y(L_{k_1}^{(6)})_{y},({L_k^{(6)}})_y) +\frac{1}{2}(y(L_{k_1}^{(6)})_{x},{L_k^{(6)}})+\kappa (y(L_{k_1}^{(6)})_{y},{L_k^{(6)}})&\\&+(\rho \xi -\mu )((L_{k_1}^{(6)})_{x},{L_k^{(6)}}) +(\frac{\xi ^2}{2}-\kappa \gamma )((L_{k_1}^{(6)})_{y},{L_k^{(6)}})\Big \}&\\ =~&\frac{S}{12}(1+r\Delta \tau )+\Delta \tau \Big (\frac{h^4}{48S^2}(2\rho \xi -1)(3y_{j}-h)&\\&+\frac{h^5}{192S^2}(4y_j-h) +\frac{h^5}{24S^2}(\rho \xi -\mu )\Big ).&\end{aligned}$$

From the above, we can obtain that

$$\begin{aligned} A(k,k)=~&S(1+r\Delta \tau )+\Delta \tau \Big (\frac{h^4 y_j}{4S^2}(1-\rho \xi +\xi ^2)-\frac{\kappa h^6}{16S^2}\Big ),\\ A(k,k_1)=~&\frac{S}{6}(1+r\Delta \tau )+\Delta \tau \Big (\frac{h^4}{24S^2}(3\rho \xi y_j-\rho \xi h-3y_j)+\frac{h^5}{96S^2}(4y_j-4\kappa y_j- \kappa h)\\&+\frac{h^5}{24S^2}(2\rho \xi -2\mu -\frac{\xi ^2}{2}+\kappa \gamma )\Big ),\\ A(k,k_2)=~&\frac{S}{6}(1+r\Delta \tau )+\Delta \tau \Big (\frac{h^5}{96S^2}(\kappa h+4\kappa y_j+h+2 y_j)&\\&+\frac{h^5}{24S^2}(\frac{\xi ^2}{2}-\kappa \gamma +\rho \xi -\mu )-\frac{\rho \xi h^4}{24S^2}(2h+3y_j)\Big ),\\ A(k,k_3)=~&\frac{S}{6}(1+r\Delta \tau )+\Delta \tau \Big (\frac{\xi h^4}{48S^2}(2\rho (2h+3y_j)-3\xi (h+2y_j))+\frac{h^5}{96S^2}((3\kappa -1)h\\&+2y_j(4\kappa -1))+\frac{h^5}{24S^2}(\mu -\rho \xi +\xi ^2-2\kappa \gamma )\Big ),\\ A(k,k_4)=~&\frac{S}{6}(1+r\Delta \tau )+\Delta \tau \Big (\frac{h^4}{24S^2}(\rho \xi h+3\rho \xi y_j-3y_j)+\frac{h^5}{96S^2}(4\kappa y_j-4 y_j-\kappa h)\\&+\frac{h^5}{24S^2}(\frac{\xi ^2}{2}-\kappa \gamma -2\rho \xi +2\mu )\Big ),\\ A(k,k_5)=~&\frac{S}{6}(1+r\Delta \tau )+\Delta \tau \Big (-\frac{\rho \xi h^4}{24 S^2}(3y_j-2h)+\frac{h^5}{96S^2}((\kappa +1)h-2y_j(2\kappa +1))\\&-\frac{h^5}{24S^2}(\frac{\xi ^2}{2}-\kappa \gamma +\rho \xi -\mu )\Big ),\\ A(k,k_6)=~&\frac{S}{6}(1+r\Delta \tau )+\Delta \tau \Big (\frac{h^4}{48S^2}(2\rho \xi (3y_j-2h)-3\xi ^2(2y_j-h))+\frac{h^5}{96S^2}(2y_j-h\\&+\kappa (3h-8y_j))+\frac{h^5}{24S^2}(\rho \xi -\mu -\xi ^2+2\kappa \gamma )\Big ). \end{aligned}$$

Under the assumptions of Theorem 5, it is not difficult to verify

$$\begin{aligned} \begin{aligned}&A(k,k)>0,\\&A(k,k_i)<0,~~i=1,\ldots ,6,\\&A(k,k)+\sum _{i=1}^6A(k,k_i)>0. \end{aligned} \end{aligned}$$

Hence, the coefficient matrix of system (11) is an M-matrix, whose inverse is nonnegative. Since $F_k^m\ge 0$, $u^m_k(k=1,2,\ldots ,(N_1+1)(N_2+1))$ are all nonnegative.