An Efficient Method for Solving Spread Option Pricing Problem : Numerical Analysis and Computing

and Applied Analysis 3 3. Removing the Correlation Term and Problem Discretization Webegin this section by transforming the PDE problem (1) in order to remove the cross derivative term. Firstlywe eliminate the reaction term rU by means of the substitution V = erτU, (8) obtaining ∂V ∂t = 1 2σ 1S2 1 ∂ 2V ∂S2 1 + 1 2σ 2S2 2 ∂ 2V ∂S2 2 + ρσ1σ2S1S2 ∂2V ∂S1∂S2 + (r − q1) S1 ∂V ∂S1 + (r − q2) S2 ∂V ∂S2 . (9) In order to remove the cross derivative term, note that the right-hand side of (9) is a linear differential operator of two variables and classical techniques to obtain the canonical form of second-order linear PDEs can be applied; see, for instance, chapter 3 of [20]. Under the assumption of correlated variables with −1 < ρ < 1 the sign of the discriminant is a2 12 − 4a11a22 = σ2 1σ2 2S2 1S2 2 (ρ2 − 1) < 0, (10) where a11 = 1 2σ 1S2 1; a12 = ρσ1σ2S1S2; a22 = 1 2σ 2S2 2. (11) Therefore, right-hand side of (9) becomes of elliptic type and a convenient substitution to obtain the canonical form is given by solving the ordinary differential equation dS2 dS1 + σ2S2 σ1S1 (−ρ ± i?̃?) = 0, (12)


Introduction
Multiasset option pricing problem has two main challenges arising from their high dimensions and the correlations between assets.This is a type of option that derives its value from the difference between the prices of two or more assets.Spread options can be written on all types of financial products including equities, bonds, and currencies.Spread options are frequently traded in the energy market [1].Two examples are as follows.
Crack Spreads.Crack spreads are the options on the spread between refined petroleum products and crude oil.The spread represents the refinement margin made by the oil refinery by "cracking" the crude oil into a refined petroleum product.
Spark Spreads.Spark spreads are the options on the spread between electricity and some type of fuel.The spread represents the margin of the power plant, which takes fuel to run its generator to produce electricity.
Since the closed form solution is not available for such problems, many numerical methods are employed.Probably the most popular methods of solving multiasset option pricing problems are those around Monte Carlo method and its modifications [2,3].However, their simulation are usually time consuming and slow convergent [4].
Analytical-numerical methods are also existing in the literature as it occurs in the one asset case.So, Kirk and Aron [5] provide an approximation method to price a bivariate spread option problem, but it is inaccurate when the strike prices are high [1,6,7].Alexander and Venkatramanan improve Kirk's approach in [6] based on the hypothesis that the spread option price is the sum of the prices of two compound exchange options which avoids any strike convention.
Fourier Transform approach and numerical integration techniques are used by Chiarella and Ziveyi in [8,9] for numerical evaluation of American options written on two underlying assets.
Among the numerical methods we mention the so-called lattice methods [10,11].Recently authors in [4] improve the results of [10,11] and give fast and accurate results although it requires some minor correlation restrictions.
Dealing with pure finite difference methods the existence of the cross derivative term in the partial differential equation (PDE) problem due to asset correlation introduces additional challenges in the numerical solution due to the existence of negative coefficient terms in the numerical scheme as well as involving expensive stencil schemes.Both facts may produce numerical drawbacks coming from the dominance of the convection versus the diffusion as well as more expensive computational cost [12][13][14][15].

Abstract and Applied Analysis
In order to overcome these difficulties the authors have recently proposed several strategies dealing with finite difference approach: high order compact schemes are proposed in [16] using central difference approximations for stochastic volatility Heston model.The authors in [17] use ADI methods and clever discretizations to obtain stable numerical schemes under minor coefficients restrictions in the PDE.The authors of [14,18] use special finite difference approximations of the cross derivative term based on a seven-point stencil instead of the nine-point stencil scheme resulting from the central finite difference approach.
In this paper we continue removing cross derivative approach based on the transformation of the PDE problem initiated in [12] for the case of one asset stochastic volatility Heston model and applied to the Bates stochastic volatility jump diffusion model in [19].Apart from avoiding negative coefficients in the proposed scheme, it has the additional computational advantage that uses only a five-point stencil.
Here we consider spread option pricing problems for both European and American cases.The paper is organized as follows.The first part of the paper deals with European spread call option with payoff ( 1 ,  2 ) = ( 1 −  2 − ) + , depending on two assets  1 and  2 and strike price  ≥ 0. In Section 2 the PDE with the boundary conditions is established for the spread option pricing problem.In Section 3 we remove the cross derivative term of the PDE using the canonical form transformation method [20].In that section we also develop the discretization of the continuous European spread call option problem achieving an explicit finite difference scheme.Section 4 is devoted to the study of the properties of the method like consistence and conditional stability and positivity.In Section 5 we extend the study of previous sections to the European spread put option summarizing the main results including the treatment of the boundary conditions.Section 6 deals with American spread option pricing problem using the LCP technique based on the proposed difference scheme.Finally, in Section 7 we illustrate the numerical results computing, simulating, and comparing accuracy and CPU time with other well acknowledge methods.

Spread Option Pricing Problem and Boundary Conditions
The price ( 1 ,  2 , ) of a spread option is the solution of the PDE where  =  −  is the time to maturity ,  is the interest rate,   is the volatility of the asset   ,   is the dividend yield of asset   , and  is the correlation parameter; see [1].Since we consider backwards time, it changes the signs of the coefficients of (1) and instead of a terminal condition the following initial condition is considered: The case  = 0 corresponds to exchange options that can be reduced to 1D problem by the change of variables Since the numerical solution of European multiasset options requires the selection of a bounded numerical domain and the translation of the boundary conditions to the boundary of the domain, it is important to pay attention to such conditions.For the initial problem (1)-( 2) suitable boundary conditions areas follows.
(2) Taking the ideas developed by some authors, for instance, Duffy for the case of basket option (see [21], p. 270), for  2 = 0 we take the value of the closed form solution of the basic Black-Scholes equation of a call for  1 with given strike ,  ( 1 , 0, ) =  − ( ( 1 ) −  ( 2 )) , where where () is the standard normal cumulative distribution function.
(3) As  1 is large versus  2 + , then the behaviour of the solution looks like the asymptotic value of a onedimensional vanilla call option with the strike ( 2 +) for  1 ≫  2 +  (see [22], p. 157), (4) For large values of  2 we assume that the values are almost constants when  2 changes; therefore we can use Neuman's boundary condition: These boundary conditions will be validated with the numerical examples.

Removing the Correlation Term and Problem Discretization
We begin this section by transforming the PDE problem (1) in order to remove the cross derivative term.Firstly we eliminate the reaction term  by means of the substitution obtaining In order to remove the cross derivative term, note that the right-hand side of ( 9) is a linear differential operator of two variables and classical techniques to obtain the canonical form of second-order linear PDEs can be applied; see, for instance, chapter 3 of [20].Under the assumption of correlated variables with −1 <  < 1 the sign of the discriminant is where Therefore, right-hand side of (9) becomes of elliptic type and a convenient substitution to obtain the canonical form is given by solving the ordinary differential equation where and where one can relate the integration constant  0 to the new variables by From ( 14) and ( 15) it follows the expression of the new variables Note that By denoting (9) takes the following form without cross derivative term In accordance with [23,24] and (6) we choose the rectangular domain where   > 0,  = 1, 2, are small positive values;  2 about 3 and  1 about 3( 2 + ).For the transformed problem (19) due to (16), the rectangular domain Ω becomes a rhomboid with vertices  (see Figure 1).From ( 16) let us denote By (17) the rhomboid domain has the sides: The mesh points in the rhomboid domain are (  ,   ) such that The approximations of the derivatives appearing in (19) at the point (  ,   ,   ) are as follows: where   , ∼ (  ,   ,   ),   = .Substituting these finite difference approximations of the derivatives into (19) is approximated by the explicit difference scheme with five-point stencil where The discretization of the spatial boundary is given by Both initial and boundary conditions ( 2), ( 3), ( 4), (6), and ( 7) are transformed throughout ( 8), ( 16), (18).For the initial condition we have Boundary condition for side  is as follows: For C with small value of  2 , we have transformed Black-Scholes solution, so where Side  corresponds to large values of  1 , so behaviour of the option there leads to Boundary  corresponds with the boundary for large values of  2 .Condition (7)

Positivity, Stability, and Consistency
In this section we show that the proposed numerical scheme ( 24)-( 27) presents suitable qualitative and computational properties, such as consistency and conditional positivity and stability.
In order to guarantee positivity of the numerical solution let us check the positivity of the coefficients   .
From ( 25) and ( 27) it is easy to check that coefficients  1 ,  2 ,  4 , and  5 of the scheme (24) are positive under the condition Under conditions (35) and (36), values in interior points of the numerical domain  +1 , , 1 ≤  ≤   − 1;  + 1 ≤  ≤   +  − 1, calculated by the scheme (24), preserve nonnegativity from the previous time level  at all interior and boundary points   , ≥ 0, 0 ≤  ≤   ;  ≤  ≤   + .Let us pay attention now to the positivity of the numerical solution at the boundary of the numerical domain.On the boundary  from (30) we have zero values.On the boundary  formula (31) is used preserving the positivity since it is the transformed Black-Scholes formula throughout expression (16).Along the line  we use Neumann boundary condition (34).Therefore, the positivity at the neighbour interior point guarantees the positivity at the boundary.Finally, on the line  boundary conditions are determined by formula (33).This is transformed expression for the asymptotic behaviour of the call option (6) that is positive under the condition Following the ideas of Kangro and Nicolaides [24], the computational domain has to be large enough to translate the boundary conditions; therefore Therefore, choosing  1 max according to (38) for the transformed problem     , ≥ 0 for any   +1 ≤  ≤   +   , 1 ≤  ≤   .
In summary, the nonnegativity of the numerical solution follows from the nonnegativity of the initial values and positivity preservation under the scheme and boundary conditions.
Let us find the maximum value of the   , with respect to ,  for each fixed .For  = 0 the values   , are defined by the initial condition (29).Let us consider the following function: and find the maximum value on the domain The first derivatives are Analysing (42) one concludes that the function (, ) is increasing with respect to  and .Therefore, the maximum value is reached at the point ( 1 ,  2 ): For the next time levels  ≥ 1 let us consider the numerical scheme (24).Under conditions (35) and (36) the coefficients   > 0,  = 1, . . ., 5.Moreover, it is easy to check that Abstract and Applied Analysis Therefore the values in interior points at the (+1)th level,  ≥ 0, can be bounded by the following: This maximum can be reached on the boundary.Therefore, let us study the boundary conditions.From (30) one gets that max The values on the boundary A are equal to the interior points; therefore the inequality (45) can be applied.
The values on the boundary  are increasing with respect to the index  and reach the maximum at the point (  ,   +   − 1).From the other side, this value can be bounded by the value at the point (  ,   +   ) that is calculated by formula (31).From the theory of option pricing, the price of European call is increasing with respect to the asset price  and, moreover, it is always greater than its asymptotic function (see [25], p. 268, formula (13.1)).The proposed transformation preserves monotonicity; therefore, In summary we can conclude that max This value is transformed call option price for asset price  1 max , that is, the solution of the 1D Black-Scholes equation at the fixed point.European call option gives a right to the option holder to purchase one share stock for a certain price.The option itself cannot be worth more than the stock.In other words, Then the following result can be stated.35) and (36), numerical scheme ( 24), ( 29), ( 30), ( 31), (33), and (34) is conditionally uniformly ‖ ⋅ ‖ ∞ -stable with upper bound  =  |− 1 |  1 max .Now consistency of the numerical scheme (24) with PDE (19) will be studied [26].

Theorem 2. With previous notations, under conditions (
Let (  , ) = 0 be approximating difference equation (24).The scheme ( 24) is consistent with (19)  where   , denotes the theoretical solution of ( 19) evaluated at the point (  ,   ,   ) and  is the operator ))   . (51) Assuming that theoretical solution  of the PDE problem is four times continuously differentiable with respect to  and  and twice with respect to  and using Taylor expansion about (  ,   ,   ) it is easy to check that where For the spatial discretization one gets where (60) The discretization of the second spatial derivatives is as follows: where Then the local truncation error takes the following form: Theorem 3. Assuming that the solution of the PDE (19) admits two times continuous partial derivative with respect to time and up to order four with respect to each of space directions solution, the numerical solution computed by scheme ( 24) is consistent with (19).

European Spread Put Option Case
For the sake of brevity and in order to avoid repetition of proofs we summarize in this section the results proved in previous sections for the call option case.
Let us consider a European spread put option with payoff Equation ( 1) and all the transformed versions of it do not change for the put option.Like in the study of the previous call option case, due to different nature of the put option, we need to establish the boundary conditions and translate them to the boundary of the numerical domain.
(1) If  2 = 0 the problem is reduced to the European put option on asset  1 and strike .Therefore, the option value on this boundary is where  1 and  2 are defined by ( 5).(2) If  1 = 0, we consider (0,  2 , ) as a put option price on asset  1 = 0 with the strike  2 + : (3) For large values of  2 we suppose that the slope of the curve tends to 1; that is, (4) For  1 ≫  2 +  we consider ( 1 ,  2 , ) as a put option price with large values of  1 , and therefore Under transformations ( 8) and ( 16) the boundary conditions (70)-(73) are translated to the rhomboid numerical domain  as follows: Positivity.As the numerical scheme (24) for the interior points is the same and the new boundary conditions ( 74

American Spread Options
In the case of American type option the holder has the right to exercise option at any moment before expiring.It leads to a free boundary problem [27].In order to simplify the computational procedure this problem can be considered as a linear complementarity problem (LCP): where  represents the spatial differential operator of (1).If we rewrite (19) in the form then under transformation ( 16) LCP (80) has the following form: where L is the right-hand side of (19) and (, ) is given by (41).Numerical solution for problem (82) can be easily obtained by a small modification of the calculating procedure described above in Sections 3 and 4. On each time level where   , is defined by the finite difference equation (24).

Numerical Examples
In this section we illustrate numerical results for both European and American spread options.
In Example 1 we show that the stability conditions (35) and (36) for the European spread option can not be disregarded.These conditions are crucial for the qualitative properties such as positivity and stability of the scheme.In Figure 2 the numerical solution with ℎ and  satisfying (85) is presented.On the other hand, in Figure 3 there is a numerical simulation with time step  = 0.012 > 0.0112.The solution is unstable and it has negative values.
Example 2 illustrates the computational time as well as the comparison with analytical approximation, provided by [6].The price of European spread option can be expressed as the price of a compound exchange option (see [6]).The known derived analytical approximation for the spread option reads where The comparison of our proposed numerical method and the analytical approximation (86) is presented in Table 1.Analytical approximation is calculated for the same arrays of  1 and  2 using MATLAB-function cdf for calculating cumulative distribution function that requires a lot of computational resources, for each point.In the proposed finite difference method (FDM) cdf is only used for the boundary conditions.Last row in the Table 1 presents the CPU time for each method for the same sizes of array.For approximation method cdf-time is 925.018sec.It is the main  part of computational time.In the case of FDM cdf takes less than half of computational time -48.867sec.
Furthermore, unsuitable negative values appear in analytical approximation for small  1 and  2 while the proposed method preserves nonnegativity of the solution.
In the next example we confirm that for the spread put option case the removing technique is fruitful and the selection of the boundary conditions at the boundary of the numerical domain is appropriated because the behaviour of the solution at the boundaries fits the solution at the interior points as it is shown in Figure 4.
Example 3.With the parameters in (84) of Example 1, but with payoff ( 2 +− 1 ) + in Figure 4, we show the simulation of the numerical solution obtained by the proposed numerical scheme.
Finally, in Example 4 we compare the numerical solution obtained with our scheme (83) using LCP approach with other tested efficient methods presented in [28].(88) Table 2 shows the option price at fixed asset prices  1 = 96 and  2 = 100 for different values of strike price  evaluated by one-dimensional integration analytic method (Analytic), Fast Fourier Transform (FFT), Monte Carlo (MC) method, and our proposed method.The accuracy is competitive with other methods, such as Fourier Transform with high number of discretization (4096) and Monte Carlo with big number of time steps (2000).
Another set of parameters can be considered in order to compare it with data in [9]: (89) In [9], three methods are considered to evaluate the spread option price for data 7.4 and fixed values  1 = 100 and  2 = 200: a numerical integration method (NIM), a method of lines (MOL), and a Monte Carlo approach.NIM used in [9] is based in the recursive integration techniques developed by Huang et al. [29].In NIM the option is treated as a Bermudan that is an option that can be exercised at discrete points of time and Simpsons rule is used.For the experiment, the authors used 1024 spatial nodes; see [8,9] for more details.The MOL is implemented by a spatial discretization of 1438 × 200 mesh points and a three-time level scheme is used for the integration in time [9].50 000 simulations have been used in the Monte Carlo method described in [30].In our method FDM, the step size discretization considered is ℎ = 0.1.To compare FDM with the three previous methods, absolute errors are computed in Table 3 with respect to a reference value obtained by the Fourier space time-stepping approach of Surkov [31] with large number of nodes: 4096 spatial nodes and 300 time steps.From Table 3 we can state that the proposed method has the same order of accuracy as MOL, more efficient, and requires less computational time than Monte Carlo method.It is slightly less efficient than the NIM, but we prove that it possesses a very highly desirable qualitative behaviour in finance.
)-(77) are nonnegative, then under conditions (35) and (36) the numerical solution for the European spread put option is nonnegative.

Table 2 :
Comparing several methods for American option pricing.

Table 3 :
CPU time in sec (first row) and absolute difference (second row) for different methods depending on number of time steps for fixed number of space steps for parameters (89).