A Fitted Multi-Point Flux Approximation Method for Pricing two options

TL;DR

This paper introduces novel numerical schemes based on Multi-Point Flux Approximation (MPFA) for accurately pricing two-asset options by effectively handling PDE degeneracy and improving upon existing methods.

Contribution

The paper develops a fitted MPFA method combined with upwinding and $ heta$-Euler time discretization, offering enhanced accuracy over current fitted finite volume approaches.

Findings

Fitted MPFA scheme outperforms existing methods in accuracy.

Numerical simulations validate the effectiveness of the new schemes.

The methods handle PDE degeneracy efficiently on non-uniform grids.

Abstract

In this paper, we develop novel numerical methods based on the Multi-Point Flux Approximation (MPFA) method to solve the degenerated partial differential equation (PDE) arising from pricing two-assets options. The standard MPFA is used as our first method and is coupled with a fitted finite volume in our second method to handle the degeneracy of the PDE and the corresponding scheme is called fitted MPFA method. The convection part is discretized using the upwinding methods (first and second order) that we have derived on non uniform grids. The time discretization is performed with $\theta$- Euler methods. Numerical simulations show that our new schemes can be more accurate than the current fitted finite volume method proposed in the literature.

Figures

Tables

Table 1. Table 1: Table of errors. The computational domain of the problem is Ω = [ 0 ; 300 ] × [ 0 ; 300 ] × [ 0 ; T ] Ω 0 300 0 300 0 𝑇 \Omega=[0;300]\times[0;300]\times[0;T] with T = 1 / 6 𝑇 1 6 T=1/6 , K = 100 𝐾 100 K=100 , the volatilities σ 1 = σ 2 = 0.3 subscript 𝜎 1 subscript 𝜎 2 0.3 \sigma_{1}=\sigma_{2}=0.3 . The correlation coefficient is ρ = 0.5 𝜌 0.5 \rho=0.5 , the risk free interest r = 0.1 𝑟 0.1 r=0.1 and Δ ​ τ = 1 / 100 Δ 𝜏 1 100 \Delta\tau=1/100 .

	Fitted fin vol	MPFA-$1^{st}$ upw	MPFA-$2^{nd}$ upw	fitted MPFA-$1^{st}$ upw	fitted MPFA -$2^{nd}$ upw
$50\times 50$	0.0134	0.0060	0.0059	0.0060	0.0060
$70\times 70$	0.0133	0.0044	0.0044	0.0044	0.0044
$85\times 85$	0.0132	0.0037	0.0037	0.0037	0.0037
$100\times 100$	0.0132	0.0032	0.0032	0.0032	0.0032
$150\times 150$	0.0131	0.0024	0.0023	0.0023	0.0023

Table 2. Table 2: Table of errors. The computational domain of the problem is Ω = [ 0 ; 300 ] × [ 0 ; 300 ] × [ 0 ; T ] Ω 0 300 0 300 0 𝑇 \Omega=[0;300]\times[0;300]\times[0;T] with T = 1 / 6 𝑇 1 6 T=1/6 , K = 100 𝐾 100 K=100 , the volatilities σ 1 = σ 2 = 0.3 subscript 𝜎 1 subscript 𝜎 2 0.3 \sigma_{1}=\sigma_{2}=0.3 . The correlation coefficient is ρ = 0.5 𝜌 0.5 \rho=0.5 , the risk free interest r = 0.08 𝑟 0.08 r=0.08 and Δ ​ τ = 1 / 100 Δ 𝜏 1 100 \Delta\tau=1/100 .

	Fitted fin vol	MPFA-$1^{st}$ upw	MPFA-$2^{nd}$ upw	fitted MPFA-$1^{st}$ upw	fitted MPFA -$2^{nd}$ upw
$50\times 50$	0.0134	0.0060	0.0059	0.0060	0.0060
$100\times 100$	0.0104	0.0064	0.0063	0.0064	0.0063
$150\times 150$	0.0131	0.0056	0.0055	0.0056	0.0055

Table 3. Table 3: Table of errors. The computational domain of the problem is Ω = [ 0 ; 300 ] × [ 0 ; 300 ] × [ 0 ; T ] Ω 0 300 0 300 0 𝑇 \Omega=[0;300]\times[0;300]\times[0;T] with T = 1 / 6 𝑇 1 6 T=1/6 , K = 100 𝐾 100 K=100 , the volatilities σ 1 = σ 2 = 0.3 subscript 𝜎 1 subscript 𝜎 2 0.3 \sigma_{1}=\sigma_{2}=0.3 . The correlation coefficient is ρ = 0.5 𝜌 0.5 \rho=0.5 , the risk free interest r = 0 𝑟 0 r=0 and Δ ​ τ = 1 / 100 Δ 𝜏 1 100 \Delta\tau=1/100 .

	Fitted fin vol	MPFA-$1^{st}$ upw	MPFA-$2^{nd}$ upw	fitted MPFA-$1^{st}$ upw	fitted MPFA -$2^{nd}$ upw
$100\times 100$	0.0152	0.0239	0.0235	0.0240	0.0229
$150\times 150$	0.0151	0.0231	0.0228	0.0232	0.0229

Table 4. Table 4: Table of errors. The computational domain of the problem is Ω = [ 0 ; 4 ] × [ 0 ; 4 ] × [ 0 ; T ] Ω 0 4 0 4 0 𝑇 \Omega=[0;4]\times[0;4]\times[0;T] with T = 2 𝑇 2 T=2 , K = 1 𝐾 1 K=1 , the volatilities σ 1 = σ 2 = 1 subscript 𝜎 1 subscript 𝜎 2 1 \sigma_{1}=\sigma_{2}=1 . The correlation coefficient is ρ = 0.3 𝜌 0.3 \rho=0.3 , the risk free interest r = 0.5 𝑟 0.5 r=0.5 and Δ ​ τ = 1 / 100 Δ 𝜏 1 100 \Delta\tau=1/100 .

	Fitted fin vol	MPFA-$1^{st}$ upw	MPFA-$2^{nd}$ upw	fitted MPFA-$1^{st}$ upw	fitted MPFA -$2^{nd}$ upw
$50\times 50$	0.1208	0.0631	0.0669	0.0623	0.0659
$100\times 100$	0.1203	0.0572	0.0648	0.0559	0.0629

Table 5. Table 5: Table of errors. The computational domain of the problem is Ω = [ 0 ; 4 ] × [ 0 ; 4 ] × [ 0 ; T ] Ω 0 4 0 4 0 𝑇 \Omega=[0;4]\times[0;4]\times[0;T] with T = 2 𝑇 2 T=2 , K = 1 𝐾 1 K=1 , the volatilities σ 1 = σ 2 = 1 subscript 𝜎 1 subscript 𝜎 2 1 \sigma_{1}=\sigma_{2}=1 . The correlation coefficient is ρ = 0.3 𝜌 0.3 \rho=0.3 , the risk free interest r = 0.5 𝑟 0.5 r=0.5 and Δ ​ τ = 1 / 10 Δ 𝜏 1 10 \Delta\tau=1/10 .

	Fitted fin vol	MPFA-$1^{st}$ upw	MPFA-$2^{nd}$ upw	fitted MPFA-$1^{st}$ upw	fitted MPFA -$2^{nd}$ upw
$50\times 50$	0.1196	0.0562	0.0643	0.0555	0.0624
$100\times 100$	0.1201	0.0626	0.0664	0.0618	0.0654

Equations

$$\frac{\partial U}{\partial \tau }=\frac{1}{2}\sum _{i,j=1}^n{\sigma }_i{\sigma }_j{\rho }_{ij}{S}_i{S}_j\frac{{\partial }^{2}U}{\partial S_i\partial S_j}+r\sum _{i=1}^n{S}_i\frac{\partial U}{\partial S_i}-rU$$

$$\displaystyle\left\{\begin{array}[]{l}dx(t)\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ =\leavevmode\nobreak\ \mu_{1}xdt+\sigma_{1}xdW_{1}\\ \\ dy(t)\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ =\leavevmode\nobreak\ \mu_{2}ydt+\sigma_{2}ydW_{2}\\ \\ dW_{1}(t)dW_{2}(t)=\leavevmode\nobreak\ \rho dt\end{array}\right.$$

$$\frac{\partial U}{\partial \tau }=\frac{1}{2}{\sigma }_1^2{x}^2\frac{{\partial }^{2}U}{\partial x^2}+\rho {\sigma }_1{\sigma }_{2}xy\frac{{\partial }^{2}U}{\partial x\partial y}+\frac{1}{2}{\sigma }_2^2{y}^2\frac{{\partial }^{2}U}{\partial y^2}+rx\frac{\partial U}{\partial x}+ry\frac{\partial U}{\partial y}-rU$$

$$\displaystyle\left\{\begin{array}[]{l}U(x,y,0)=\max\left(\max(x,y)-K,0\right)\\ \\ U(0,y,\tau)=0\\ \\ U(x,0,\tau)=0\\ \end{array}\right.$$

$$\frac{\partial U}{\partial \tau }=\nabla \cdot (\mathbf{M}\nabla U)+\nabla (fU)+\lambda U$$

$$\displaystyle\mathbf{M}=\frac{1}{2}\left(\begin{array}[]{lr}\sigma_{1}^{2}x^{2}&\rho\sigma_{1}\sigma_{2}xy\\ &\\ \rho\sigma_{1}\sigma_{2}xy&\sigma_{2}^{2}y^{2}\end{array}\right),f=\left(\begin{array}[]{c}(r-\sigma_{1}^{2}-\frac{1}{2}\rho\sigma_{1}\sigma_{2})x\\ \\ (r-\sigma_{2}^{2}-\frac{1}{2}\rho\sigma_{1}\sigma_{2})y\end{array}\right)$$

$$\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\text{leavevmode}\lambda =-3r+{\sigma }_1^2+{\sigma }_2^2+\rho {\sigma }_1{\sigma }_2$$

$$I_{x_i}=[x_{i-1};x_i],I_{y_j}=[y_{j-1};y_j]i,j=1,...,N+1.$$

$$x_{i-\frac{1}{2}}=\frac{x_{i-1}+x_i}{2},y_{j-\frac{1}{2}}=\frac{y_{j-1}+y_j}{2}i,j=1,...,N,$$

$$x_{-\frac{1}{2}}=x_0=0,x_{N+\frac{3}{2}}=x_{N+1}=x_{\max },y_{-\frac{1}{2}}=y_0=0,y_{N+\frac{3}{2}}=y_{N+1}=y_{\max }.$$

$${\mathbf{M}}^{ij}=\frac{1}{\mathrm{meas}({\mathcal{C}}_{i,j})}{\int }_{{\mathcal{C}}_{i,j}}\mathbf{M}dxdy,i,j=1,...,N.$$

$$M^{i,j}=\left[\begin{array}[]{lcr}\frac{\sigma_{1}^{2}}{6}\frac{x_{i+\frac{1}{2}}^{3}-x_{i-\frac{1}{2}}^{3}}{x_{i+\frac{1}{2}}-x_{i-\frac{1}{2}}}&&\frac{\rho\sigma_{1}\sigma_{2}}{8}(x_{i+\frac{1}{2}}+x_{i-\frac{1}{2}})(y_{j+\frac{1}{2}}+y_{j-\frac{1}{2}})\\ &&\\ \frac{\rho\sigma_{1}\sigma_{2}}{8}(x_{i+\frac{1}{2}}+x_{i-\frac{1}{2}})(y_{j+\frac{1}{2}}+y_{j-\frac{1}{2}})&&\frac{\sigma_{2}^{2}}{6}\frac{y_{j+\frac{1}{2}}^{3}-y_{j-\frac{1}{2}}^{3}}{y_{j+\frac{1}{2}}-y_{j-\frac{1}{2}}}\end{array}\right].$$

$${\int }_{{\mathcal{C}}_{ij}}\frac{\partial U}{\partial \tau }d\mathcal{C}={\int }_{{\mathcal{C}}_{ij}}\nabla \cdot (\mathbf{M}\nabla U)d\mathcal{C}+{\int }_{{\mathcal{C}}_{ij}}\nabla (fU)d\mathcal{C}+{\int }_{{\mathcal{C}}_{ij}}\lambda Ud\mathcal{C}.$$

$${\int }_{{\mathcal{C}}_{ij}}\frac{\partial U}{\partial \tau }d\mathcal{C}\approx \mathrm{meas}({\mathcal{C}}_{ij})\frac{d\mathcal{U}}{d\tau }(x_i,y_j,\tau )$$

$${\int }_{{\mathcal{C}}_{ij}}\lambda Ud\mathcal{C}\approx meas({\mathcal{C}}_{ij})\lambda \mathcal{U}(x_i,y_j,\tau ).$$

$${\int }_{{\mathcal{C}}_{ij}}\nabla \cdot (\mathbf{M}\nabla \mathcal{U})d\mathcal{C}$$

$${\int }_{{\mathcal{C}}_{ij}}\nabla (f\mathcal{U})d\mathcal{C}$$

$${\mathcal{F}}^{ij}={\int }_{{\mathcal{C}}_{ij}}\nabla \cdot ({\mathbf{M}}^{ij}\nabla \mathcal{U})={\int }_{\partial {\mathcal{C}}_{ij}}({\mathbf{M}}^{ij}\nabla \mathcal{U})\cdot \vec{n}d\partial \mathcal{C}$$

$$\nabla f=-\frac{1}{2\mathcal{A}}\left[\Big{(}f(x_{2})-f(x_{1})\Big{)}\nu_{2}+\Big{(}f(x_{3})-f(x_{1})\Big{)}\nu_{3}\right]$$

$$\nabla \mathcal{U}=\frac{1}{2\mathcal{A}}\left[({\bar{\mathcal{U}}}_1-{\mathcal{U}}_{ij}){\omega }_1+({\bar{\mathcal{U}}}_2-{\mathcal{U}}_{ij}){\omega }_2\right]$$

$$\text{leavevmode}\text{leavevmode}\text{leavevmode}{\mathcal{R}}_{ij}=[x_{i-1};x_i]\times [y_{j-1};y_j],i,j=1,\dots ,N+1.$$

$$f_p^{i-1,j-1}={\Gamma }_p{\vec{n}}_p^T{\mathbf{M}}^{i-1,j-1}\nabla \mathcal{U}$$

$$\left[\begin{array}[]{c}f_{1}^{i-1,j-1}\\ \\ f_{3}^{i-1,j-1}\end{array}\right]=G^{i-1,j-1}\left[\begin{array}[]{c}\bar{\mathcal{U}}_{1}-\mathcal{U}_{i-1,j-1}\\ \\ \bar{\mathcal{U}}_{3}-\mathcal{U}_{i-1,j-1}\end{array}\right]$$

$$G^{i-1,j-1}=\left[\begin{array}{ccccc}{\Gamma }_1{n}_1^T{M}^{i-1,j-1}{\omega }_1& & & & {\Gamma }_1{n}_1^T{M}^{i-1,j-1}{\omega }_2\\ & & & & \\ {\Gamma }_2{n}_1^T{M}^{i-1,j-1}{\omega }_2& & & & {\Gamma }_2{n}_2^T{M}^{i-1,j-1}{\omega }_2\end{array}\right]$$

$$\displaystyle\left[\begin{array}[]{c}f_{1}^{i,j-1}\\ \\ f_{4}^{i,j-1}\end{array}\right]=G^{i,j-1}\left[\begin{array}[]{c}\mathcal{U}_{i,j-1}-\bar{\mathcal{U}}_{1}\\ \\ \bar{\mathcal{U}}_{4}-\mathcal{U}_{i,j-1}\end{array}\right]\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \left[\begin{array}[]{c}f_{2}^{i-1,j}\\ \\ f_{3}^{i-1,j}\end{array}\right]=G^{i-1,j}\left[\begin{array}[]{c}\bar{\mathcal{U}}_{2}-\mathcal{U}_{i-1,j}\\ \\ \mathcal{U}_{i-1,j}-\bar{\mathcal{U}}_{3}\end{array}\right]$$

$$\displaystyle\left[\begin{array}[]{c}f_{2}^{ij}\\ \\ f_{4}^{ij}\end{array}\right]=G^{ij}\left[\begin{array}[]{c}\mathcal{U}_{ij}-\bar{\mathcal{U}}_{2}\\ \\ \mathcal{U}_{ij}-\bar{\mathcal{U}}_{4}\end{array}\right]\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\$$

$$\begin{array}[]{ccccc}f_{1}&=&f_{1}^{i-1,j-1}&=&f_{1}^{i-1,j}\\ &&&&\\ f_{2}&=&f_{2}^{ij}&=&f_{2}^{i-1,j}\\ &&&&\\ f_{3}&=&f_{3}^{i-1,j}&=&f_{3}^{i-1,j-1}\\ &&&&\\ f_{4}&=&f_{4}^{i,j-1}&=&f_{4}^{ij}.\end{array}$$

$$f_1$$

$$f_2$$

$$f_3$$

$$f_4$$

$$f=\left[\begin{array}[]{c}f_{1}\\ f_{2}\\ f_{3}\\ f_{4}\end{array}\right],\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \mathcal{U}=\left[\begin{array}[]{c}\mathcal{U}_{i-1,j-1}\\ \mathcal{U}_{i,j-1}\\ \mathcal{U}_{i-1,j}\\ \mathcal{U}_{ij}\end{array}\right],\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \mathcal{V}=\left[\begin{array}[]{c}\bar{\mathcal{U}}_{1}\\ \bar{\mathcal{U}}_{2}\\ \bar{\mathcal{U}}_{3}\\ \bar{\mathcal{U}}_{4}\end{array}\right]$$

$$f=C^{ij}\mathcal{V}+F^{ij}\mathcal{U}$$

Full text

1 Introduction

Pricing multi-assets options is of great interest in the financial industry (see Persson and Sydow [2007]). Multi-asset options are options based on more than one underlying. There are several kinds of multi-assets options, few of them are exchange options, rainbow options, baskets options, best or worst options, quotient options, foreign exchange options, quanto options, spread options, dual-strike options and out-performance options. Pricing these options lead to the resolution of the following second order degenerated Black-Scholes Partial Differential Equations (PDE)(see Persson and Sydow [2007])

[TABLE]

where $r$ is the risk free interest, $U$ is the option value at time $\tau$, $\tau=T-t$ with $t$ and $T$ respectively the instantaneous and maturity time, \,S_{i}\leavevmode\nobreak\ represents the asset $i$ price, $\sigma_{i}$ represents the volatility of asset $i$, $\rho_{ij}$ represents the correlation between the assets $i$ and $j$, where $i,j=1,...,n$. The main difference between multi-assets options is their payoff functions which represent the initial condition of the corresponding backward PDE. The spatial domain of the PDE is infinite, but for its numerical resolution, a truncation is required (see Duffy [2013],Chapter 3). It has been observed that when the stock price $S$ approaches the region near to zero, the Black Scholes PDE is degenerated (see Duffy [2013], chapter 30.3). Moreover, the initial condition of the PDE has a discontinuity in its first derivative when the stock price is equal to the strike $K$. This discontinuity has an adverse impact on the accuracy when the finite difference method is used (see Wilmott [2005], chapter 26). Therefore, for the spatial discretization of the PDE, it is suitable to use non-uniform grids with more points in the region around $S=0$ and $S=K$ in order to handle the degeneracy and the discontinuity. To overcome the above challenges, many methods have been proposed in the literature. Thereby, Wang [2004] proposed a fitted finite volume method for one dimensional Black Scholes PDE and the rigorous convergence proof is provided by Angermann and Wang [2007]. Besides, Huang et al. [2006] adapted the fitted finite volume discretization method for the two-dimensional Black-Scholes PDE and its rigorous convergence proof is analysed by Huang et al. [2009]. Although these two fitted finite volume methods are stable, they are only order 1 with respect to asset price variables.

In this paper, we present two novel discretization methods for the two-dimensional Black Scholes PDE based on a special kind of finite volume method, the so-called Multi-Point Flux Approximation (MPFA) method. This method was introduced by Aavatsmark [2002] and has been used in fluid dynamics for flow and transport equations (see Sandve et al. [2012] and references therein). Actually, the MPFA was designed to give a correct discretization of the flow equation for general grids including fractures (see Aavatsmark [2002], Sandve et al. [2012]). The MPFA method is essentially based on the approximation of a linear function gradient over a triangle, the calculation and the continuity of flux through edges of this triangle. The convergence of MPFA method is usually second order in space domain on rough grids (see Aavatsmark [2007], Stephansen [2012]). Our first numerical method here is the standard MPFA , which is fully used to approximate the second order operator. To the best of our knowledge, this method was not yet used to solve degenerated Black Scholes PDE in finance. To build our new fitted MPFA method, we couple the standard MPFA with the upwind methods (first and second order) to approximate two dimensional options pricing. Besides, the fitted finite volume proposed by Wang [2004] is used to handle the degeneracy of the PDE in the region where the stocks price approach zero (degeneracy region). In the region, where the PDE in not degenerated, we apply the MPFA method. The novel numerical technique from this combination is called fitted MPFA method and will obviously improve the accuracy of the current fitted finite volume in the literature, since more approximations involving are second order in space. Naturally, these two methods are applicable to other types of multi-asset options and also to financial models such as Heston [1993] model and Bates [1996] model on non-uniform grids. Another advantage of our novel fitted MPFA is that it can easily be adapted to more structured commercial or open-source softwares as the standard MPFA (see Lie et al. [2012]).

The rest of the paper is organized as follows. In section 2, we start by introducing the Black Scholes model for option with 2 stocks and the corresponding partial differential equation. Afterwards, we set the frame of the numerical domain of study suitable for the finite volume method application. Section 3 is devoted to the spatial discretization of the PDE. We describe the Multi-Point Flux Approximation method for the discretization of the diffusion term of the PDE. The upwind methods (first and second order) are used for the the convection term discretization. We end the section 3 with the fitted MPFA which is a combination of a fitted finite volume method and the MPFA method. The time discretization is performed using the $\theta-$Euler methods in the section 4. In section 5, we perform numerical experiments. Those numerical simulations show that the two proposed schemes (the standard MPFA method and fitted MPFA method ) can be more accurate than the current fitted finite volume method proposed in the literature. General conclusion is given in section 6.

2 Formulation of the problem

2.1 Black-Scholes model with 2 underlying assets

An option with two underlying assets modeled by the Black Scholes equation is formulated as follows

[TABLE]

where $\mu_{i},\sigma_{i},W_{i}$ are respectively the drift, the volatility and the Wiener process governing the stocks $x,y$ and $\rho$ is the correlation coefficient between the two Wiener processes. By applying the Ito’s formula and using the standard arbitrage argument, it is well known ( see Hull [2003], Kwok [2008], Wilmott et al. [1993] ) that the value of the option $U$ follows the following two-dimensional Black-Scholes Partial differential equation on the domain $D=[0,+\infty)\times[0,+\infty)\times[0,T]$

[TABLE]

where $\tau=T-t$, $T$ is the maturity time, $t$ the current time and $r$ is the risk-free interest. For European rainbow option price on maximum of two risky assets, the following initial and boundary conditions are used

[TABLE]

with $K$ the strike price. But to compare our numerical solution with the existing fitted finite volume method, the exact solution will be used at the boundary. In order to apply the finite volume method, it is convenient to re-write the Partial Differential Equation (8) in the following divergence form

[TABLE]

where

[TABLE]

Note that $\mathbf{M}$ does not satisfying the standard ellipticity condition (see [Tambue, 2016, (3)]), so the PDE (15) is degenerated.

We will assume Dirichlet boundary condition in the entire domain.

2.2 Finite volume method

Let us consider the new domain $\Omega$ of study by truncating $D$ such that $\Omega=I_{x}\times I_{y}\times[0,T]$ where $I_{x}=[0,x_{\max}]$ and $I_{y}=[0,y_{\max}]$. In the sequel of this work, the Black-Scholes partial differential equation (8) is considered over the truncated domain $\Omega$. At $x=x_{\max}$ and $y=y_{\max}$, the linear boundary condition will be applied (see Huang et al. [2006]). The intervals $I_{x}$ and $I_{y}$ will be subdivided into $N$ part in the following way (see Huang et al. [2006, 2009]) without loss the generality as irregular grids such as triangular grids can be used.

[TABLE]

Let us set the mid-points $x_{i-\frac{1}{2}}$ and $y_{j-\frac{1}{2}}$ as follows

[TABLE]

with $h_{i}=x_{i+\frac{1}{2}}-x_{i-\frac{1}{2}},\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ l_{j}=y_{j+\frac{1}{2}}-y_{j-\frac{1}{2}}$     and

[TABLE]

For $i,j=1,\ldots,N$, we denote by $\mathcal{C}_{ij}=[x_{i-\frac{1}{2}};x_{i+\frac{1}{2}}]\times[y_{j-\frac{1}{2}};y_{j+\frac{1}{2}}]$ a control volume associated to our subdivision.

Note that the control volume $\mathcal{C}_{ij}$ is the area surrounding the grid point $(x_{i},y_{j})$. Our goal is to approximate the option function $U$ at $(x_{i},y_{j})$ 111center of the control volume $\mathcal{C}_{i,j}$ by a function denoted $\mathcal{U}$. The matrix $\mathbf{M}$ in (15) will be replaced by its average value within each control volume

[TABLE]

where $\leavevmode\nobreak\ \mathrm{meas}(\mathcal{C}_{ij})$ is the measure of $\mathcal{C}_{ij}$. Thereby, we have

[TABLE]

Now let us consider the divergence form given in (15). Following the finite volume method’s principle, we integrate the partial differential equation (15) over each control volume $\mathcal{C}_{ij}$ and we have

[TABLE]

The next section will be dedicated to spatial discretization of equation (20). For the term in the left hand side of (20) and for the last term in its right hand side, we use the mid-point quadrature rule for their approximations. More precisely

[TABLE]

[TABLE]

The diffusion term

[TABLE]

of (20) will be approximated using the Multi-point flux approximation (MPFA) method or our novel fitted Multi-point flux approximation. More details will be given in the next section. Besides, the convection term

[TABLE]

of (20) will be approximated using the upwind methods (first or second order). Note that the standard two -point flux approximation in Tambue [2016] can only be consistent in the approximation of (23) if and only if the grid is $\mathbf{M}-$orthogonal.

3 Space discretization

The spatial discretization of (15) consists of approximating all terms in (20) over the control volumes of the study domain.

3.1 Discretization of the diffusion term

Let us start by applying the divergence theorem to the diffusion term (23) as follows, for $\leavevmode\nobreak\ \leavevmode\nobreak\ i,j=1,...,N$

[TABLE]

where $\vec{n}$ is the outward vector from the control volume.

Now, we can apply the so-called Multi-Point Flux Aprroximation(MPFA) to approximate the integral defined in (25).

3.1.1 Multi-Point Flux Approximation (MPFA) method

There exists several types of Multi-Point Flux Approximation methods. The most known of MPFA methods are the O-method and the L-method. In our study, we focus on the O-method because it is the classical MPFA method and it is more intuitive comparing to the L-method which is fairly new and less intuitive (see Aavatsmark [2002]). Here, we follow the description of the O-method developed by Aavatsmark [2002].

We will start by giving an approximation of the gradient in the integral expression (25).

• Let us consider a triangle $x_{1}x_{2}x_{3}$, $\nu_{i}$ the outer normal vector of the edge located opposite of vertex $x_{i},\leavevmode\nobreak\ i=1,2,3$ and $f$ a linear function over this triangle (see Figure 2a). The length of $\nu_{i}$ is equal to the length of the edge to which it is normal.

The gradient expression of the function $f$ in the triangle may be written in the form

[TABLE]

where $\mathcal{A}$ is the area of the triangle.

Thereby, assuming that our solution $\mathcal{U}$ is linear over the control volume $\mathcal{C}_{ij}$ with center $x_{1}(x_{i},y_{j})$, and applying (26) in the triangle $x_{1}\bar{x}_{1}\bar{x}_{2}$ (see Figure 2b), we have

[TABLE]

where $\mathcal{U}_{ij}=\mathcal{U}(x_{1})=\mathcal{U}(x_{i},y_{j}),\bar{\mathcal{U}}_{1}=\mathcal{U}(\bar{x}_{1}),\bar{\mathcal{U}}_{2}=\mathcal{U}(\bar{x}_{2})$ and the vectors $\omega_{1}$ and $\omega_{2}$ are respectively inner normal vector to the edge $x_{1}\bar{x}_{1}$ and $x_{1}\bar{x}_{2}$ with the same length with those vectors, and $\mathcal{A}$ is the area of the triangle $x_{1}\bar{x}_{2}\bar{x}_{2}$.

Let us called interaction volume $\mathcal{R}_{ij}$ a cell grid defined as follows

[TABLE]

We may notice that an interaction volume $\mathcal{R}_{ij}$ is covering an area in the intersection of the control volumes $\mathcal{C}_{i-1,j-1},\mathcal{C}_{i-1,j},\mathcal{C}_{i,j-1}$ and $\mathcal{C}_{ij}$. Here, we follow closely Aavatsmark [2007].

• We denote respectively by $x_{1}(x_{i-1},y_{j-1}),x_{2}(x_{i},y_{j-1}),x_{3}(x_{i-1},y_{j})$  and $x_{4}(x_{i},y_{j})$ the centre of the control volume $\mathcal{C}_{i-1,j-1},\mathcal{C}_{i,j-1},\mathcal{C}_{i-1,j}$ and $\mathcal{C}_{i,j}$. We denote also by $\bar{x}_{1},\bar{x}_{2},\bar{x}_{3}\leavevmode\nobreak\ and\leavevmode\nobreak\ \bar{x}_{4}$ the midpoints of the segment $x_{1}x_{2}$, $x_{3}x_{4},x_{1}x_{3}\leavevmode\nobreak\ \text{and}\leavevmode\nobreak\ x_{2}x_{4}$

Our goal in an interaction volume is to compute the flux through the half edges $1,2,3\leavevmode\nobreak\ and\leavevmode\nobreak\ 4$ inside the interaction volume (see Figure 3). The flux through the half edge p seen from the centre $x_{1}=(x_{i-1},y_{j-1})$ of the control volume $\mathcal{C}_{i-1,j-1}$ is denoted $f_{p}^{i-1,j-1}$. By using the expression (25), we have

[TABLE]

where $\Gamma_{p}$ is the length of half edge p, $\vec{n}_{p}$ is the outward unit normal vector to the half edge p. It is convenient to let $\vec{n}_{p}$ point in the direction of increasing global cell indices. In that case, we have two kinds of inner normal vectors. The vertical ones denoted $\omega_{1}$ and the horizontal ones denoted $\omega_{2}$.

By considering the triangle $x_{1}\bar{x}_{1}\bar{x}_{3}$ (Figure 3) in the control volume $\mathcal{C}_{i-1,j-1}$, using the expression of gradient (27) and the flux expression (29), we have for $i,j=1,...,N$

[TABLE]

with

[TABLE]

By applying (30) in the triangles $x_{2}\bar{x}_{1}\bar{x}_{4},x_{3}\bar{x}_{2}\bar{x}_{3}$ and $x_{4}\bar{x}_{4}\bar{x}_{2}$ (see Figure 3), we have

[TABLE]

Since the flux through an edge is continuous, from (30) and (43) we have

[TABLE]

It follows that

[TABLE]

Let us set

[TABLE]

The equation (3.1.1) allows to have

[TABLE]

where

[TABLE]

[TABLE]

From (3.1.1), we can also have

[TABLE]

where

[TABLE]

[TABLE]

Thereby, $\mathcal{V}$ can be eliminated from (55) by solving (58) with respect to $\mathcal{V}$. This gives the following the expression of the flux through the 4 half edges inside the interaction volume $\mathcal{R}_{ij}$

[TABLE]

where

[TABLE]

$T^{ij}$ is called transmissibility matrix of the interaction volume $\mathcal{R}_{ij}$.

From (61), we are now able to get the flux through the half edges 1,2,3 and 4 inside the interaction volume $\mathcal{R}_{ij}$.

Let us recall that to approximate the integral in (25), we need to compute the flux through the edges on a control volume $\mathcal{C}_{ij}$. We might notice that we need four interaction volume with centres the four vertices of the control volumes in order to cover all the edges of the considered control volume (see Figure 4).

For the volume control $\mathcal{C}_{ij}$, we denote by ${}_{\mathcal{E}}f_{l}^{ij}$ the flux through lower half eastern edge, by ${}_{\mathcal{E}}f_{u}^{ij}$ the flux through the upper half eastern edge. The flux ${}_{\mathcal{E}}f^{ij}$ through the east edge of the control volume $\mathcal{C}_{ij}$ is calculated as follows: The lower half eastern edge is contained in the interaction volume $\mathcal{R}_{i+1,j}$ and it is in position 2 in the interaction of volume (see Figure 4). So by using (61) we have:

[TABLE]

Similarly, the upper half eastern edge is contained in the interaction volume $\mathcal{R}_{i+1,j+1}$ and it is in position 1 in the interaction volume. So by using (61) we have:

[TABLE]

Finally the flux through the east edge of the control volume $\mathcal{C}_{i+1,j+1}$ will be the addition of ${}_{\mathcal{E}}f_{l}^{ij}$ and ${}_{\mathcal{E}}f_{u}^{ij}$. Thereby we have

[TABLE]

Similarly, we compute the flux through the northern, western and southern edge of the control volume $\mathcal{C}_{ij}$. Afterwards, we sum up the flux through the 4 edges of the control to get the outflux $\mathcal{F}^{ij}$ through the edges of the control volume $\mathcal{C}_{ij}$. Therefore we have for $i,j=1,...,N$

[TABLE]

where

[TABLE]

Let us notice that for the control volumes near to the boundary of the our domain, some terms from the boundary conditions will be involved in (63) .

Hence (25) becomes

[TABLE]

where $A_{mp}$ is a $N^{2}\times N^{2}$ matrix and

[TABLE]

with $0_{N}$ is $N\times N$ null matrix , $W_{i},Y_{i},X_{i}$ are tridiagonal matrices, and $F_{mp}$ is a $N^{2}$ vector coming from the boundary conditions. The structure of the diffusion matrix $A_{mp}$ can be viewed in Figure 5

3.2 Discretization of the convection term

In this section, the convection term

[TABLE]

with

[TABLE]

will be approximated by the upwind methods (first and second order).

3.2.1 First order upwind

The first order upwind method discussed by [LeVeque, 2004, chapter 4.8] or Tambue [2016] will be applied to approximate the second term of (20). Using the divergence theorem, we have for $i,j=2,...,N$

[TABLE]

Note that $I^{ij}$ is calculated by summing up the flux through the edges of the control volume $\mathcal{C}_{ij}$. The flux through an edge using the first order upwind will depend on the sign of $f\cdot\vec{n}$ on this edge. If the sign of $f\cdot\vec{n}$ is positive, $\mathcal{U}_{ij}$ will be used to approximate $\mathcal{U}$ in the expression $(f\cdot\vec{n}\mathcal{U})$ otherwise we will use the value of $\mathcal{U}$ in other side of the edge. Note that an edge may be the interface of two control volumes. By doing so, we have for $i,j=1,...,N$

[TABLE]

where

[TABLE]

with

[TABLE]

Let us notice that for the control volumes near to the boundary of the our domain, some terms from the boundary conditions will be involved in (66). Hence, (66) gives

[TABLE]

where $A_{up}$ is a $N^{2}\times N^{2}$ matrix

[TABLE]

with $0_{N}$ is $N\times N$ null matrix, $H_{i}$ is a tridiagonal matrix, $P_{i},Q_{i}$ are diagonal matrices and $F_{up}$ is a vector coming from the boundary conditions. Therefore, combining the MPFA method (64) and the first order upwind (67), we have

[TABLE]

with

[TABLE]

where $A_{L}$ is a diagonal matrix of size $N^{2}\times N^{2}$ coming from the discretisation of (22). The diagonal elements of $A_{L}$ are $A_{ii}=h_{i}l_{i}\lambda$ for $i=1,...,N$ with $\lambda$ given in (15). The matrix $L$ is also a diagonal matrix of size $N^{2}\times N^{2}$ whose diagonal elements are $L_{ii}=h_{i}l_{i}$ for $i=1,\ldots,N$

3.2.2 Upwind second order

We start by applying the mid-quadrature rule as follows.

[TABLE]

where $\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ p_{i}=(r-\sigma_{1}^{2}-\frac{1}{2}\rho\sigma_{1}\sigma_{2})x_{i},\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ q_{j}=(r-\sigma_{2}^{2}-\frac{1}{2}\rho\sigma_{1}\sigma_{2})y_{j}$ and $\omega=2r-\sigma_{1}^{2}-\sigma_{2}^{2}-\rho\sigma_{1}\sigma_{2}$. Let us use the second order upwind to approximate the first derivatives in (69) at the point $(x_{i},y_{j})$.

Approximation of the first derivative using a 3 points stencil

Here, we want to express the first derivative $\frac{\partial\mathcal{U}_{ij}}{\partial x}$ in terms of $\mathcal{U}_{i+2,j},\mathcal{U}_{i+1,j}$ and $\mathcal{U}_{ij}$. Set $h=\underset{1\leq i\leq N}{\max}h_{i}$. Let us find $a,b$ and $c$ such that

[TABLE]

Thereby, using a $2^{nd}$ order Taylor expansion at the point $(x_{i},y_{j})$ on $\mathcal{U}_{i+2,j}$ and $\mathcal{U}_{i+1,j}$, we have

[TABLE]

By matching, we have

[TABLE]

Solving (76), we have

[TABLE]

Therefore we have

[TABLE]

**Application to the $2^{nd}$ order upwind method on non uniform grids

**By analogy with the procedure to get the expression in (78), the term $p_{i}\frac{\partial\mathcal{U}_{ij}}{\partial x}$ is approximated as follows:

• (i)

$p_{i}>0$ then

[TABLE]

• (ii)

$p_{i}<0$ then

[TABLE]

Similarly for the first derivative $\frac{\partial\mathcal{U}_{ij}}{\partial y}$, we have

• (iii)

when $q_{j}>0$ then

[TABLE]

• (iv)

when $q_{j}<0$

[TABLE]

By combining $(i),(ii),(iii),(iv)$ in (69), for $i,j=2,\ldots,N-1$, we have

[TABLE]

where

[TABLE]

For the control volumes near the boundary of the study domain, two ghost points or the first order upwind method can be used. Finally, we have the following matrix form

[TABLE]

where

[TABLE]

and

[TABLE]

where $H_{1},H_{N}$ are tridiagonal matrices, for $i=2,\ldots,N-1,\leavevmode\nobreak\ H_{i}$ are penta-diagonal matrices and $P_{i},R_{i},W_{i},Q_{i}$ are diagonal matrices, and $F_{up}$ is a vector coming from the boundary conditions. A structure of the advection matrix using the second order upwind method can be viewed in Figure 6.

As for the first order upwinding, combining the MPFA method (64) and the second order upwind method (80), we have

[TABLE]

[TABLE]

where $A_{L}$ is a diagonal matrix of size $N^{2}\times N^{2}$ coming from the discretisation of (22). The elements of $A_{L}$ are $h_{i}l_{j}\lambda$ for $i,j=1,...,N$ with $\lambda$ given in (15). The matrix $L$ is also a diagonal matrix of size $N^{2}\times N^{2}$ whose diagonal elements are $L_{ii}=h_{i}l_{i}$ for $i=1,\ldots,N$

Actually, the PDE (8) is degenerated when the stock price is approaching zero $(x\rightarrow 0,y\rightarrow 0)$ which has an adverse impact on the accuracy of the numerical method. However, to overcome the degeneracy, we are going to apply a fitted finite volume method in the degeneracy region $(x\rightarrow 0,y\rightarrow 0)$. More details about this fitted method is given in the next section.

3.3 Fitted Multi-Point Flux Approximation

The fitted Multi-Point Flux Approximation is a combination of the fitted finite volume method ( see Huang et al. [2006, 2009]) and the Multi-Point Flux Approximation method. The fitted finite volume helps to deal with the degeneracy of the PDE (8). We approximate simultaneously the diffusion term and the convection term in the degeneracy region by solving a two-points boundary problem. In the region where the PDE is not degenerated, we apply the standard Multi-point flux approximation to the diffusion term as described in the previous section.

Let us set

[TABLE]

where $\mathbf{M}$ and $f$ are defined in (15). Thereby, we have the following decomposition over a control volume $\mathcal{C}_{ij}$, for $\leavevmode\nobreak\ \leavevmode\nobreak\ i,j=1,...,N$

[TABLE]

with $\vec{n}$ is the outward unit normal vector, $m_{11},m_{12},m_{21},m_{22}$ the coefficients of the matrix $\mathbf{M}$ and $p,q$ coefficients of vector $f$ defined in (15).

In their work, Huang et al. [2006, 2009] showed how the fitted finite method is used to approximate each of the integral in (3.3).

3.3.1 Fitted Finite volume method in the degeneracy region

Following Huang et al. [2006], the fitted finite volume method is used to approximate the flux through the edges which are effectively in the degeneracy region notably the western edge of the control volume $\mathcal{C}_{1,j}$ for $j=1,\ldots,N$ and the southern edge of the control volume $\mathcal{C}_{i,1}$ for $i=1,\ldots,N$ .

Thereby, the flux through the southern edge of the control volume $\mathcal{C}_{i,1}$ for $i=1,\ldots,N$ is calculated as follows.

The fitted finite volume method is applied to approximate the integral along the southern edge of control volume $\mathcal{C}_{i,1}$. The idea is to approximate the integral over $[x_{i-\frac{1}{2}};x_{i+\frac{1}{2}}]$ by a constant. We start by applying the mid-quadrature rule as follows:

[TABLE]

Besides we have

[TABLE]

with $e=\frac{1}{2}\sigma_{2}^{2},\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ k=r-\sigma_{2}^{2}-\frac{1}{2}\rho\sigma_{1}\sigma_{2}$ and $h^{\prime}=\frac{1}{2}\rho\sigma_{1}\sigma_{2}x$.

We want to approximate

[TABLE]

by a linear function over $I_{y_{1}}=(0,y_{1})$ satisfying the following two-points boundary value problem

[TABLE]

By solving this problem we get

[TABLE]

Thereby, by using (84), (89), (90) and the forward difference for approximating the first partial derivative $\frac{\partial\mathcal{U}}{\partial x}$ we get

[TABLE]

where

[TABLE]

Similarly, for the western edge of the control volume $\mathcal{C}_{1,j},\leavevmode\nobreak\ \leavevmode\nobreak\ for\leavevmode\nobreak\ \leavevmode\nobreak\ j=1,...,N$, we have

[TABLE]

with

[TABLE]

3.3.2 Fitted Multi-Point Flux Approximation

The fitted Multi-Point Approximation method consists of calculating the flux through the edges which are totally in the degeneracy region using the fitted finite volume method as described in the previous paragraph. For the edges which are not totally in the degeneracy region, the flux is approximated using simultaneously the Multi-point flux approximation and the upwind methods (first order or second order). In the other hand, the MPFA method and the upwind methods are used to approximate respectively the diffusion term and the convection term over the control volumes which are not in the degeneracy region.

Considering (3.3), in fact, in the control volume $\mathcal{C}_{11}$, the southern and western edges are in the degeneracy region, the northern and the eastern edges are not in the degeneracy region. Thereby, the flux through the southern and western edges are approximated using the fitted finite volume method, while the flux through the eastern and northern edges are approximated using simultaneously of the MPFA method and the upwind method. This gives

[TABLE]

with

[TABLE]

Similarly, for the control volume $\mathcal{C}_{1,j}\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ j=1,\ldots,N$, we have

[TABLE]

[TABLE]

For the control $\mathcal{C}_{i,1}\leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ \leavevmode\nobreak\ i=2,\ldots,N$, we have:

[TABLE]

with

[TABLE]

As we already mentioned, for the control volumes which are not in the degeneracy region, we use the multi-Point flux approximation to approximate the diffusion term and the upwind methods (first and second order) to approximate the convection term. So by combining as before, we obtain the following ODE

[TABLE]

where

[TABLE]

with $F$ the vector of boundary conditions, $A_{L}$ is a diagonal matrix of size $N^{2}\times N^{2}$ coming from the discretisation of (22). The elements of $A_{L}$ are $h_{i}l_{j}\lambda$ for $i,j=1,...,N$ with $\lambda$ given in (15). The matrix $L$ is also a diagonal matrix of size $N^{2}\times N^{2}$ whose diagonal elements are $h_{i}l_{j}$ for $i,j=1,\ldots,N$ and

[TABLE]

The fitted matrix $Z$ uses the first order upwind method. The matrices $D_{i},K_{i},L_{i}$ are tri-diagonal matrices defined as follows. For $i=1,N$

[TABLE]

For $i=2,\ldots,N-1$

[TABLE]

where all the elements $a_{i,j}^{1},b_{i,j}^{1},c_{i,j}^{1},d_{i,j}^{1},e_{i,j}^{1},\gamma_{i,j}^{1},\lambda_{i,j}^{1}$ are defined in (93),(94),(95) and the others elements are defined in (63) and (66).

Similarly, combining the fitted finite volume method, the MPFA and the second order upwind method we have

[TABLE]

where

[TABLE]

with $G$ the vector of boundary conditions, $A_{L}$ is a diagonal matrix of size $N^{2}\times N^{2}$ coming from the discretisation of (22). The elements of $A_{L}$ are $h_{i}l_{j}\lambda$ for $i,j=1,...,N$ with $\lambda$ given in (15). The matrix L is also a diagonal matrix of size $N^{2}\times N^{2}$ whose elements are $h_{i}l_{j}$ for $i,j=1,\ldots,N$ and

[TABLE]

The elements of matrix Y are matrices. Indeed $0_{N}$ is a zeros matrix of size $N\times N$. The matrices $H_{i},P_{i},Q$ are tri-diagonal matrices and $W_{i},R_{i}$ are diagonal matrices defined as follows:

[TABLE]

For $i=2,\ldots,N-1$

[TABLE]

and

[TABLE]

where all the elements $a_{i,j}^{1},b_{i,j}^{1},c_{i,j}^{1},d_{i,j}^{1},e_{i,j}^{1},\gamma_{i,j}^{1},\lambda_{i,j}^{1}$ are defined (93),(94),(95), and the others elements are defined in (63) and (79).

4 Time discretization

Let us consider the ODE stemming from the spatial dicretization and given by (68),(81),(96) and (97)

[TABLE]

Using the $\theta$-method for the time discretization, we have

[TABLE]

Hence

[TABLE]

with

[TABLE]

5 Numerical experiments

In this section, we perform some numerical simulations and compare different numerical schemes developed in this work. More precisely, we compare the novel fitted MPFA method combined to the upwind methods, first method (fitted MPFA-$1^{st}$ upw) and second order (fitted MPFA-$2^{nd}$ upw), with the fitted finite volume method by Huang et al. [2006] (fitted FV) and the standard MPFA method combined to the upwind methods, first (MPFA-$1^{st}$ upw) and second order (MPFA-$2^{nd}$ upw). The analytical solution of the PDE (8) is well known (see Haug [2007] ) and given as

[TABLE]

where

[TABLE]

and

[TABLE]

Note that in all our numerical schemes, the Dirichlet Boundary condition is used with the value equal to the analytical solution.

The graphs of option price with different methods are given in Figure 8,Figure 9 and Figure 10

In this paragraph, we consider the four numerical methods illustrated in the previous sections and the fitted finite volume method Huang et al. [2006]. We evaluate the error of these numerical method with respect to the analytical solution (5). The $L^{2}$-norm is used to compute the error as follows:

[TABLE]

where $\mathcal{U}$ is the numerical solution, $U^{ana}$ the analytical solution and $meas(\mathcal{C}_{i,j})$ is the measure of the control volume $\mathcal{C}_{ij}$. This gives the following table:

As we can observe in Table 1-Table 5, the errors from our fitted MPFA and MPFA methods are smaller compared to those of fitted finite volume in Huang et al. [2006]. We can also note that when $r$ become smaller, the gaps between the errors of the fitted finite volume in Huang et al. [2006] and our fitted MPFA and MPFA methods reduce.

6 Conclusion

In this paper, we have presented the Multi-Point Flux Approximation (MPFA) to approximate the diffusion term of Black-Scholes Partial Differential Equation in its divergence form. The MPFA method coupled with the upwind methods (first and second order) have been used to solve numerically the Black-Scholes PDE. To handle the degeneracy of Black Scholes PDE, we have proposed a novel method based on a combination of the MPFA method and fitted finite volume by Huang et al. [2006]. We have performed some numerical simulations, which show that our fitted MPFA method coupled with first or second order upwinding methods are more accurate than the fitted finite volume method by Huang et al. [2006]. Rigorous convergence proof of the fitted MPFA will be our nearest future work.

Acknowledgement

This work was supported by the Robert Bosch Stiftung through the AIMS ARETE Chair programme (Grant No 11.5.8040.0033.0).