A fitted L-Multi-point Flux Approximation method for pricing options

In this paper, we introduce a special kind of finite volume method called Multi-Point Flux Approximation method (MPFA) to price European and American options in two dimensional domain. We focus on the L-MPFA method for space discretization of the diffusion term of Black-Scholes operator. The degeneracy of the Black Scholes operator is tackled using the fitted finite volume method. This combination of fitted finite volume method and L-MPFA method coupled to upwind methods gives us a novel scheme called the fitted L-MPFA method. Numerical experiments show the accuracy of the novel fitted L-MPFA method comparing to well known schemes for pricing options.

where L is the following Black-Scholes operator with r is the risk free interest, V is the option value at time τ , V is the payoff, τ = T − t with t and T respectively the instantaneous and maturity time. For i, j = 1, . . . , n, S i represents the asset i price, σ i represents the volatility of asset i, ρ ij represents the correlation between the assets i and j. Furthermore,  proposed a power penalty method to solve the linear complementary problem for pricing American options. The power penalty problem is formulated as follows: where β is penalty parameter and k is the power of the method. Let us notice that, when we take the penalty parameter β = 0 in (3), we get the Black-scholes Equation for pricing European options, with the operator L defined in (2). However, the power penalty problem (3) can not be solved analytically, therefore numerical methods are required for its resolution. Nevertheless, the Black-Scholes operator (2) is degenerated when the stock price approaches zero. This degeneracy can affect the accuracy of the numerical method used for the resolution. To tackle this problem, several methods have been proposed. The fitted finite volume method, proposed by S. Wang in Wang [2004] whereby a rigorous proof of convergence is provided, appears to be more attractive. Moreover, the fitted finite volume method has been used for the resolution of the two dimensional second order Black Scholes PDE followed by the convergence proof in Huang et al. [2006]. In spite of the fact that the fitted finite volume methods perform well for the resolution of the Black Scholes PDE, they are only of order 1 with respect to asset price variable. Besides, the fitted O-Multi-Point Flux Approximation (O-MPFA) method has been proposed in Koffi and Tambue [2019] to overcome the degeneracy problem of the Black-Scholes PDE. It has been shown that the O-MPFA is more accurate than the classical fitted finite volume method by Wang [2004]. However, the O-MPFA is heavy (9 points stencil method) and for more general grids, the convergence rate of the O-MPFA method may reduce (see Aavatsmark [2007]).
In this paper, we focus on the L-MPFA method which is based on the approximation of a linear function gradient defined over a given triangle and the continuity of flux through the edges of this triangle. Indeed, the L-MPFA method is a 7 points stencil method while the O-MPFA is a 9 points stencil method. This shows that the O-MPFA method can be computationally more expensive than the L-MPFA method. Moreover, for more general grids, the order reduction in convergence rate is larger for the O-MPFA than the L-MPFA (see Aavatsmark [2002]). Thereby, to approximate the solution of the second order Black Scholes operator, we couple the L-MPFA method with the upwind methods (first and second order). Besides, the degeneracy of the Black Scholes operator (2) is handled by the fitted finite volume , Wang [2004], when the stock price is approaching zero. The L-MPFA method coupled with the upwind methods (1 st and 2 nd order) is used to approximate the solution of (3) when the Black operator is not degenerated. We call fitted L-MPFA method that combination of the L-MPFA method and the fitted finite volume method. Numerical simulations show that the new fitted L-MPFA method is more accurate than the fitted O-MPFA method developed in Koffi and Tambue [2019] and the standard fitted finite volume method developed in Huang et al. [2006].
The paper is structured as follows. In section 2, we present the power penalty problem with the corresponding initial and boundary conditions. The spatial discretization of the linear operator is developed in the section 3. Details on the L-MPFA method of the diffusion term discretization are provided. The convection term is discretized using the upwind methods (1 st and 2 nd method). At the end of the section 3, the novel fitted MPFA method is provided. The θ− Euler method is used for the time discretization method in the section 4. Numerical experiments are presented for the different numerical methods are presented in the section 5. The conclusions of our study are drawn in the last section.
2 Formulation of the problem 2.1 Option with two underlying assets Pricing an American option with 2 underlying assets leads to solve the following power penalty problem: where the Black-Scholes operator is defined as: with the following initial and boundary conditions with K the strike price, V * the payoff for basket options and α i , i = 1, 2, are weights such that α 1 + α 2 = 1.
When the penalty parameter β = 0 in (4), we get the Black-Scholes Partial Differential Equation for pricing European options with the corresponding initial and boundary conditions In order to apply the finite volume method, it is convenient to re-write the Black-Scholes operator (5) in the following divergence form

Finite volume method
Let us consider the new domain Ω of study by truncating D such that and I y = [0, y max ]. In the sequel of this work, the Black-Scholes operator (5) is considered over the truncated domain Ω.
At x = x max and y = y max , the linear boundary condition will be applied Huang et al. [2006]. The intervals I x and I y will be subdivided into N in the following way Huang et al. [2006] Huang et al. [2009 without loss the generality.
Let us set the mid-points x i− 1 2 and y j− 1 2 as follows Note that for i, j = 1, . . . , N the control volume C ij is the area surrounding the grid point (x i , y j ).
Our goal is to approximate the option function V at (x i , y j ) 1 by a function denoted V.
The matrix M in (8) will be replaced by its average value in each control volume Mdxdy, i, j = 1, ..., N.
1 center of the control volume C ij where meas(C ij ) is the measure of C ij . Thereby, we have Now let us consider the divergence form given in (8). According to the finite volume method, we integrate the partial differential equation (8) over each control volume C ij . For i, j = 1, ..., N , we have: The next section will be dedicated to spatial discretization of equation (12). For the term in the left hand side of the equality sign and for the last one in the right hand side of (12), we use the mid-point quadrature rule for their approximation as follows: the convection term Cij ∇(f V)dC (16) of (12) will be approximated using the upwind methods (first or second order).

The diffusion term
Cij ∇ · (M∇V)dC (17) of (12) will be approximated using the Multi-Point Flux Approximation (MPFA) L-method or the fitted multi-point flux approximation Lmethod. More details about these methods will be given in the next section.

Space discretization
The spatial discretization of (8) consists in approximating all terms in (12) over the control volumes of the study domain.

Discretization of the diffusion term
Let us start by applying the divergence theorem to the diffusion term (17) as follows, for i, j = 1, ..., N : where n is the outward vector from the control volume. Now, we can apply the so-called L-Multi-Point Flux Aprroximation(MPFA) method to approximate the integral defined in (18).

L -Multi-Point Flux Approximation (L-MPFA) method
The L-MPFA method takes its name from the fact that the curve connecting the three control volume centres considered for the application of the method, constitutes a stylized "L". Here, we follow the description of the L-method given in Aavatsmark [2002].
Let us consider the triangle x 1 x 2 x 3 (see Figure 2), and a linear function g defined over this triangle. we define Thereby, the gradient of the linear function g may be expressed as follows: where ν 2 , ν 3 are respectively the normal vector to x 2 − x 1 and x 3 − x 1 defined by Rν 3 Let's notice that the matrix R is a rotation of angle − π 2 . Thereby the vector ν 2 and ν 3 have same length with the vectors x 2 − x 1 and x 3 − x 1 .
x 1 Let us called interaction volume R ij a cell grid defined as follows We denote respectively by .We denote also byx 1 ,x 2 ,x 3 andx 4 the midpoints of the segment x 1 x 2 , x 3 x 4 , x 1 x 3 and x 2 x 4 . We may notice that an interaction volume R ij , for i, j = 1, . . . , N + 1, is covering an area in the intersection of the control volumes C ij , C i+1,j , C i,j+1 and C i+1,j+1 . An interaction volume can be divided into 2 triangles such that the half edges 1, 2 are in the triangle T 1 = x 1 x 2 x 3 and the half edges 3, 4 are in the triangle T 2 = x 1 x 3 x 4 (see Figure 3). Here, We follow the procedure in Aavatsmark [2007].

Fig. 3: Interaction volume
In an interaction volume, we aim to compute the flux through the half edges 1, 2, 3 and 4 (See 3).
Thereby, using (18), the flux f ij p through the half edge p seen from the centre of the control volume C ij is expressed as follows: where n p is the vector normal to the half edge p with the same length.
Let us consider the triangle T 1 = x 1 x 2 x 3 from the interaction volume R ij In the triangle T 12 =x 1 x 2x2 , using the flux expression (23) and following the gradient expression in (21) we get Moreover, using the property of the matrix R, we have Using also the equation (19) and the expression of gradient (21) in the triangle T 12 lead to where In the triangle T 11 = x 1x1x5 , we have ReplacingV 5 by its expression (26) in (27),we have where ReplacingV 5 by its expression (26) in (29),we have where Since the flux is continuous through edges, then using (24), (28) and (30) we have by setting The system of equations (31) can be written as Using the expressions at both sides of the second equalities of system equations (31) we have Thereby, by solving (34) with respect to V and replacing in (33) we get where Now considering the triangle T 2 (see figure 3) and applying the above procedure used in the triangle T 1 , we are able to the fluxes through the half edges 3 and 4 as follows: where For simplicity, in an interaction volume R ij , the flux through the half edges 1, 2, 3 and 4 are given by and T ij is 4 × 4 matrix coming from R ij and S ij defined in (35),(36). T ij is called the transmissibility matrix of the interaction volume R ij .
Let us notice that the flux through a full edge will be the addition of the fluxes through its 2 half edges.
Let us recall that, from (18), our aim is to compute the flux through the edges of the control volume C ij . In order to cover the boundary of a control volume, we need to four interaction volumes Let us denote, for the volume control C ij , by E f ij l the flux through lower half eastern edge, by E f ij u the flux through the upper half eastern edge. The flux E f ij through the east edge of the control volume C ij is calculated as follows: The lower half eastern edge is in position 3 in the triangle T 2 of the interaction volume R i+1,j (See figure 4). So by using (37) we have: Similarly, the upper half eastern edge is in position 1 in the triangle T 1 of the interaction volume R i+1,j+1 and it is in position 1 in the interaction volume. Thereby using ((37)) we have: Finally the flux through the eastern edge of the control volume C ij will be the addition of E f ij l and E f ij u . Thereby we have: The same method is applied to calculate the flux through the northern, western and southern edges of the control volume C ij . The flux through the edges of the control volume C ij is obtained by summing up the flux through the 4 edges. This gives : where This leads to a system of equations which can be written as follows: where are tridiagonal matrix and F mp is a N 2 vector coming from the boundary conditions. The diffusion matrix A mp is under the folowing form: Fig. 6: A structure of the diffusion matrix using L-MPFA method.
As we can see on Fig.6, the L-MPFA method is a 7 points stencil method, unlikely to the O-MPFA method (see Koffi and Tambue [2019] ) which is a 9 points stencil method.

Discretisation of the convection term
The integral of convection term will be approximated using the upwind methods (1 st and 2 nd order ). We start by applying the divergence theorem, and we have for i, j = 1, ..., N : with n an outward unit normal vector.

First order upwind method
The first order upwind method discussed in [LeVeque, 2004, chapter 4.8] will be applied to evaluate the second term of (12).
I ij is calculated by summing up the flux through the edges of the control volume C ij . The flux through an edge using the first order upwind will depend on the sign of f · n on this edge. If the sign of f · n is positive, V ij will be used to approximate (f · nV) otherwise we will use the value of V in other side of the edge.
By doing so, we get ∀i, j = 1, ..., N where Equation (43) will lead to a system of equations which be written as follows: where A up is a M × N matrix with 0 N is N × N zeros 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. The advection matrix A up using the first order upwind method is under following form Therefore, combining the L-MPFA method (41) and the first order upwind (44), we get with where A L is a diagonal matrix of size N × N coming from the discretisation of (14). The diagonal elements of A L are A ii = h i l i λ for i = 1, ..., N with λ given in (8). The matrix L is also a diagonal matrix of size N × N whose diagonal elements are L ii = h i l i for i = 1, . . . , N

Second order upwind method
A second order approximation is used to calculated the flux defined in (42). For instance, the flux E J ij through the eastern edge of the control volume C ij is approximated as follows: We will start by giving an approximation of the gradient in the integral expression (18). where is the outward unit normal vector to the eastern edge E ij of the control volume C ij . We set f x = f · n E and we have Then we get We use the same argument to calculate the flux N J ij , W J ij , S J ij through the northern, western and southern edges of the control volume C ij and after sum them up. We get then For the control volumes near the boundary of the study domain, the first order upwind method is used for the approximation of the flux through edges directly connected to the boundary.
Equation (49) leads to a system of equations which can be written as: where for i = 1, . . . , N K i is penta-diagonal matrice and R i , G i , S i , H i are diagonal matrices. F 2up is a vector coming from the boundary conditions. Therefore, combining the L-MPFA method (41) and the second order upwind method (50), we have where A L is a diagonal matrix of size N × N coming from the discretisation of (14). The diagonal elements of A L are A ii = h i l i λ for i = 1, ..., N with λ given in (8). The matrix L is also a diagonal matrix of size N × N whose diagonal elements are L ii = h i l i for i = 1, . . . , N Besides, the ellipticity condition for the PDE (2) is not satisfied when the stocks price (x → 0 and/or y → 0) is near to zero. This may cause some oscillations of the numerical solution when the PDE is degenerate. Nevertheless, Wang [2004] suggested a fitted finite volume method to deal with the degeneracy of the PDE. Thereby, the fitted finite volume method will be applied in the degeneracy region (x → 0 and/or y → 0) in the next section.

Fitted finite volume
The fitted finite volume method is used to approximated the flux through edges which are (fully) in the degeneracy region i.e the western edge of the control volume C 1,j j = 1, .., N and the southern edge of the control volume C i,1 i = 1, . . . , N .
We want to approximate g(V) = ax ∂V ∂x + bV by a constant over I x1 = (0, x 1 ) satisfying the following two-point boundary value problem By solving this problem we get Thereby, by using (52),(53), (54), (55) and the forward difference to approximate the first partial derivative ∂V ∂y we get Similarly, for the flux through the southern edge of the control volume C i,1 i = 1, ..., N , we have The fitted L-Multi-Point Flux Approximation method ( with the 1 st order upwind method ) 1. Fitted L-MPFA method (with 1 st order upwind method) Here the fitted finite volume method is combined with the first order upwind method . Thereby we have: For the control volume C 11 , the western and southern edges of are (fully) in the degeneracy region. The integrals over the western and southern edges of the control volume C 11 are then approximated using the fitted finite volume (56) and (57). The integrals over the eastern and northern edges of the control C 11 , which are not in the degeneracy region, are approximated using the L-MPFA method coupled to the upwind method (1 st and 2 nd order).

C11
∇k(V)dC 11 ≈ aa 11 V 11 + bb 11 V 21 + cc 11 V 22 + dd 11 V 12 + ee 11 V 01 + ββ 11 V 10 where aa 11 = T 22 11 + T 21 34 + T 22 Similarly, for the control volume C 1,j j = 1, . . . , N , only the southern edge is (fully) in the degeneracy region then the integral over it, is approximated using the fitted finite volume method (57). The integrals over the eastern, northern and western edges are approximated using the L-MPFA method coupled to the upwind methods(1 st and 2 nd order) where aa 1,j = T 2,j+1 11 + T 2,j 34 + T 2,j+1 bb 1,j = T 2,j 33 + T 2,j+1 12 − T 2,j 43 + l j min(f 2 x , 0) cc 1,j = T 2,j+1 13 + T 2,j+1 43 dd 1,j = T 1,j+1 Using the same argument as above, for the control volume C i,1 i = 2, .., N , the integral over the southern edge is approximated using the fitted finite volume (57). The integrals over the eastern, northern and western edges are approximated using the L-MPFA method combined with the upwind methods (1 st and 2 nd order) where Besides, for the control volume C ij , i, j = 2, ..., N , the L-MPFA method is used to approximate the diffusion term and the upwind to approximate the advection term. This leads to the following semi-discrete equation with F the vector of boundary conditions, A L a diagonal matrix of size N ×N coming from the discretisation of (14). The diagonal elements of A L are A ii = h i l i λ for i = 1, ..., N with λ given in (8). The matrix L is also a diagonal matrix of size N × N whose diagonal elements are L ii = h i l i for i = 1, . . . , N and The elements of matrix Z are matrices. 0 N is a zeros matrix of size N . The matrices D i , K i , L i are tridiagonal matrices defined as follows: where the elements aa ij , bb ij , cc ij , dd ij , ee ij , ββ ij are defined in (58),(58),(59), and the elements a ij , b ij , c ij , d ij , e ij , Ω ij ∆ ij , β ij , φ ij , α ij , µ ij , η ij are defined in (40) and (43).

Fitted Multi-Point Flux Approximation (2 nd order upwind)
Similarly, the fitted MPFA-L method deriving from the combination of the L-MPFA method and the 2 nd order upwind method leads to the following equation : with F the vector of boundary conditions. A L is a diagonal matrix of size N × N coming from the discretisation of (14). The diagonal elements of A L are A ii = h i l i λ for i = 1, ..., N with λ given in (8). The matrix L is also a diagonal matrix of size N × N whose diagonal elements are L ii = h i l i for i = 1, . . . , N and The elements of matrix Y are matrices. 0 N is a zeros matrix of size N . The matrices H i i = 1, .., N are a penta-diagonal matrices and the matrices P i , R i , Q i , W i are diagonal matrices.

Time discretization
Let us consider the ODE stemming from the spatial discretization and given by (45), (51), (60) and (61) dV dτ = AV + G(V) + F By using the θ-Euler method for the time discretization, we have At every time iteration, the nonlinear system where V m+1 is the solution is solving using the Newton method. Note that where the time step is ∆τ = T M , T being the maturity time.

Numerical experiments
In this section, we perform some numerical simulations for the L-MPFA method combined to the upwind methods (first and second order) and for the fitted L-MPFA method combined to the upwind methods (first and second order).

Errors for European call options
The computational domain of the problem is Ω = [0; 300] × [0; 300] × [0; T ] with T=1/12. The numerical experiments are performed with the strike price E = 100, the volatilities σ 1 = σ 2 = 0.3, the correlation coefficient ρ = 0.3 and the risk free interest r = 0.08 Here, by taking β = 0 in (4), the L-MPFA method illustrated in the previous sections will be compared to the fitted finite volume method, Wang [2004], and the fitted O-MPFA methods for pricing multi-asset options for pricing options introduced in Koffi The solution using the L-MPFA coupled to the 2 nd order upwind method is illustrated as below The L 2 -norm is used to compute the error is where V is the numerical solution, V ana the analytical solution and meas(C ij ) is the measure of the control volume C ij . This gives the following tables: As we can observe in Table 1 and    Since there is no analytical solution for the power penalty problem (4) for pricing American put options, and the numerical solution given by the fitted L-MPFA coupled to 2 nd order upwind method is more accurate when pricing European options (see Table 1 and Table 2), we have chosen for reference solution or "exact solution" the numerical solution given by the fitted L-MPFA coupled to 2 nd order upwind method with dt = T /256. The relative error of all the numerical methods used in this this study will be performed with respect to this reference solution.   Again we can observe in Table 3 and Table 4, the novel fitted L-MPFA coupled to the 2 nd order upwind method is more accurate than the standard fitted finite volume by Huang et al. [2006].

Conclusion
In this paper, the L-MPFA methods have been introduced to approximate the diffusion term of the Black-Scholes PDE. The upwind methods (1 st and 2 nd order) are used for space discretization of the convection term appearing in the two dimensional Black-Scholes PDE. We have provided a novel scheme called the fitted L-MPFA method to handle the degeneracy of the Black Scholes PDE by combining the fitted finite volume and the L-MPFA method coupled to the upwind methods. Numerical experiments are performed and comparison between the L-MPFA methods, the O-MPFA methods by Koffi and Tambue [2019] and the fitted finite methods by Huang et al. [2006] are performed. The results have shown that the fitted L-MPFA method coupled to the 2 nd order upwind method is more accurate than the fitted finite volume by Huang et al. [2006] and the O-MPFA method by Koffi and Tambue [2019] for pricing Europeans and American options.