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

Authors: R. Company; V. N. Egorova; L. Jódar
Journal: Abstract and Applied Analysis
 (2016)

---

## Abstract

This paper deals with numerical analysis and computing of spread option pricing problem described by a two-spatial variables partial differential equation. Both European and American cases are treated. Taking advantage of a cross derivative removing technique, an explicit difference scheme is developed retaining the benefits of the one-dimensional finite difference method, preserving positivity, accuracy, and computational time efficiency. Numerical results illustrate the interest of the approach.


---

## Body

## 1. 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–15].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 payoffg(S1,S2)=S1-S2-E+, depending on two assets S1 and S2 and strike price E≥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.Throughout the paper, we use the following notation:A∞=max⁡aij:1≤i≤N,1≤j≤M will denote a supremum norm for a given matrix A=(aij)∈RN×M.


## 2. Spread Option Pricing Problem and Boundary Conditions

The priceU(S1,S2,τ) of a spread option is the solution of the PDE (1)∂U∂τ=12σ12S12∂2U∂S12+12σ22S22∂2U∂S22+ρσ1σ2S1S2∂2U∂S1∂S2+r-q1S1∂U∂S1+r-q2S2∂U∂S2-rU,0<S1<∞,0<S2<∞,0<τ≤T,where τ=T-t is the time to maturity T, r is the interest rate, σi is the volatility of the asset Si,  qi is the dividend yield of asset Si, 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:(2)US1,S2,0=gS1,S2=max⁡S1-S2-E,0=S1-S2-E+.The caseE=0 corresponds to exchange options that can be reduced to 1D problem by the change of variables x1=S1/S2,  c(x1,τ)=U(S1,S2,τ)/S2.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.(1)
ForS1=0 the payoff S1-S2-E+ suggests Dirichlet’s condition(3)U0,S2,τ=0,0<S2<∞,0<τ≤T,as one-dimensional call option.(2)
Taking the ideas developed by some authors, for instance, Duffy for the case of basket option (see [21], p. 270), for S2=0 we take the value of the closed form solution of the basic Black-Scholes equation of a call for S1 with given strike E,(4)US1,0,τ=e-rτFNϕ1-ENϕ2,where(5)F=S1er-q1τ,ϕ1=1σ1τlog⁡FE+12σ12τ,ϕ2=ϕ1-σ1τ,whereN(x) is the standard normal cumulative distribution function.(3)
AsS1 is large versus S2+E, then the behaviour of the solution looks like the asymptotic value of a one-dimensional vanilla call option with the strike (S2+E) for S1≫S2+E (see [22], p. 157),(6)US1,S2,τ≈e-q1τS1-e-rτS2+E,S1≫S2+E.(4)
For large values ofS2 we assume that the values are almost constants when S2 changes; therefore we can use Neuman’s boundary condition:(7)∂U∂S2=0.These boundary conditions will be validated with the numerical examples.


## 3. 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 rU by means of the substitution(8)V=erτU,obtaining (9)∂V∂t=12σ12S12∂2V∂S12+12σ22S22∂2V∂S22+ρσ1σ2S1S2∂2V∂S1∂S2+r-q1S1∂V∂S1+r-q2S2∂V∂S2.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 (10)a122-4a11a22=σ12σ22S12S22ρ2-1<0,where(11)a11=12σ12S12;a12=ρσ1σ2S1S2;a22=12σ22S22.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(12)dS2dS1+σ2S2σ1S1-ρ±iρ~=0,where (13)ρ~=1-ρ2,and σ2S2/σ1S1-ρ±iρ~ are the conjugate roots of a11x2+a12x+a22=0. Solving (12) one gets(14)1σ2log⁡S2+-ρ±iρ~σ1log⁡S1=C0,where one can relate the integration constant C0 to the new variables by(15)x=ImC0,y=-ReC0.From (14) and (15) it follows the expression of the new variables(16)x=ρ~σ1log⁡S1;y=ρσ1log⁡S1-1σ2log⁡S2.Note that(17)y-mx=-1σ2log⁡S2,m=ρρ~.By denoting(18)Wx,y,τ=VS1,S2,τ, (9) takes the following form without cross derivative term(19)∂W∂τ=ρ~22∂2W∂x2+∂2W∂y2+ρ~σ1r-q1-σ122∂W∂x+ρσ1r-q1-σ122-1σ2r-q2-σ222∂W∂y,x,y∈R2,0<τ≤T.In accordance with [23, 24] and (6) we choose the rectangular domain (S1,S2)∈a1,b1×a2,b2=Ω, where ai>0,i=1,2, are small positive values; b2 about 3E and b1 about 3(b2+E). For the transformed problem (19) due to (16), the rectangular domain Ω becomes a rhomboid with vertices ABCD (see Figure 1). From (16) let us denote(20)c1=ρ~σ1log⁡a1,c2=log⁡a2σ2,d1=ρ~σ1log⁡b1,d2=log⁡b2σ2.Figure 1
Transformed numerical domain.By (17) the rhomboid domain has the sides:(21)AD¯=x,y:x=c1,mc1-d2≤y≤mc1-c2,AB¯=x,y:c1≤x≤d1,y=mx-d2,BC¯=x,y:x=d1,md1-d2≤y≤md1-c2,DC¯=x,y:c1≤x≤d1,y=mx-c2.The mesh points in the rhomboid domain are(xi,yj) such that (22)xi=c1+ih,xi∈c1,d1,0≤i≤Nx,Δx=h=d1-c1Nx,Δy=mh,Ny=d2-c2mh,yj=mxi-d2+j-imh,i≤j≤Ny+i.The approximations of the derivatives appearing in (19) at the point (xi,yj,τn) are as follows: (23)∂W∂x~wi+1,jn-wi-1,jn2h,∂W∂y~wi,j+1n-wi,j-1n2mh,∂2W∂x2~wi+1,jn-2wi,jn+wi-1,jnh2,∂2W∂y2~wi,j+1n-2wi,jn+wi,j-1nm2h2,∂W∂τ~wi,jn+1-wi,jnk,where wi,jn~W(xi,yj,τn), τn=nk.Substituting these finite difference approximations of the derivatives into (19) is approximated by the explicit difference scheme with five-point stencil(24)wi,jn+1=α1wi-1,jn+α2wi+1,jn+α3wi,jn+α4wi,j-1n+α5wi,j+1n,where(25)α1,2=ρ~22kh2∓ρ~σ1r-q1-σ122k2h,(26)α3=1-ρ~2kh21+1m2,(27)α4,5=ρ~2k2m2h2∓ρ~r-q1-σ12/2σ1-r-q2-σ22/2σ2k2mh.The discretization of the spatial boundary is given by(28)PAD¯=x0,yj:0≤j≤Ny,PAB¯=xi,yi:0≤i≤Nx,PBC¯=xNx,yj:Nx≤j≤Nx+Ny,PDC¯=xi,yNy+i:0≤i≤Nx.Both initial and boundary conditions (2), (3), (4), (6), and (7) are transformed throughout (8), (16), (18). For the initial condition we have(29)wi,j0=eσ1xi/ρ~-e-σ2yj-mxi-E+,0≤i≤Nx,i≤j≤i+Ny.Boundary condition for sideAD¯ is as follows:(30)w0,jn=0,0≤j≤Ny,1≤n≤Nτ.ForDC¯ with small value of S2, we have transformed Black-Scholes solution, so(31)wi,Ny+in=FinNϕ1,in-ENϕ2,in,0<i<Nx,1≤n≤Nτ,where (32)Fin=exp⁡σ1ρ~xi+r-q1τn,ϕ1,in=1σ1τnlog⁡FiE+12σ12τn,ϕ2,in=ϕ1,in-σ1τn.SideBC¯ corresponds to large values of S1, so behaviour of the option there leads to(33)wNx,jn=exp⁡σ1ρ~xNx+r-q1τn-exp⁡σ2mxNx-yj+E,Nx+1≤j≤Nx+Ny,1≤n≤Nτ.BoundaryAB¯ corresponds with the boundary for large values of S2. Condition (7) means that component S2 is constant and the null first derivative is approximated by the forward difference wi,i+1n-wi,in/mh=0,(34)wi,in=wi,i+1n,1≤i≤Nx,1≤n≤Nτ.


## 4. 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αi.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 (35)h<min⁡ρ~σ1r-q1-σ12/2;ρ~2mρ~r-q1-σ12/2/σ1-r-q2-σ22/2/σ2.Coefficientα3 is positive, if(36)k<m2h2ρ~2m2+1.Under conditions (35) and (36), values in interior points of the numerical domain wi,jn+1, 1≤i≤Nx-1;i+1≤j≤Ny+i-1, calculated by the scheme (24), preserve nonnegativity from the previous time level n at all interior and boundary points wi,jn≥0, 0≤i≤Nx;i≤j≤Ny+i.Let us pay attention now to the positivity of the numerical solution at the boundary of the numerical domain. On the boundaryAD¯ from (30) we have zero values. On the boundary DC¯ formula (31) is used preserving the positivity since it is the transformed Black-Scholes formula throughout expression (16). Along the line AB¯ we use Neumann boundary condition (34). Therefore, the positivity at the neighbour interior point guarantees the positivity at the boundary. Finally, on the line BC¯ 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(37)S1max>S2max+Eexp⁡q1-rτn,1≤n≤Nτ.Following the ideas of Kangro and Nicolaides [24], the computational domain has to be large enough to translate the boundary conditions; therefore(38)S1max≥max⁡S2max+Eexp⁡q1-rT,3S2max+E.Therefore, choosingS1max according to (38) for the transformed problem wNx,jn≥0 for any Nx+1≤j≤Nx+Ny,1≤n≤Nτ.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.As authors manage several stability criteria in literature, we recall the concept of stability used in this paper.Definition 1.
Letwn∈R(Nx+1)×(Ny+1) be the matrix involving the numerical solution {wi,jn} of PDE problem (19) computed by scheme (24), (29), (30), (31), (33), and (34),(39)wn=wi,jn0≤i≤Nx,i≤j≤Ny+i,0≤n≤Nτ.The numerical scheme (24), (29), (30), (31), (33), and (34) is said to be uniformly ·∞-stable in the domain c1,d1×c2,d2×0,T, if for every partition with k, h as defined above, condition(40)wn∞≤A,0≤n≤Nτ,holds true, where A is independent of h, k. When the scheme is uniformly stable under a condition imposed to the step size discretization, then the scheme is said to beconditionally uniformly stable.Let us find the maximum value of thewi,jn with respect to i, j for each fixed n. For n=0 the values wi,jn are defined by the initial condition (29). Let us consider the following function:(41)fx,y=exp⁡σ1ρ~x-exp⁡-σ2y-mx-Eand find the maximum value on the domain [c1,d1]×[c2,d2]. The first derivatives are(42)∂f∂x=σ1ρ~exp⁡σ1ρ~x+σ2mexp⁡-σ2y-mx>0,∂f∂y=σ2exp-σ2y-mx>0.Analysing (42) one concludes that the function f(x,y) is increasing with respect to x and y. Therefore, the maximum value is reached at the point (d1,d2):(43)maxx,y∈c1,d1×c2,d2⁡fx,y=fd1,d2=exp⁡σ1ρ~d1-exp⁡-σ2d2-md1-E=C=const.For the next time levelsn≥1 let us consider the numerical scheme (24). Under conditions (35) and (36) the coefficients αi>0,  i=1,…,5. Moreover, it is easy to check that (44)∑i=15αi=1.Therefore the values in interior points at the(n+1)th level, n≥0, can be bounded by the following:(45)wi,jn+1≤max⁡wi-1,jn,wi+1,jn,wi,jn,wi,j-1n,wi,j+1n≤maxi,j⁡wi,jn.This maximum can be reached on the boundary. Therefore, let us study the boundary conditions. From (30) one gets that (46)maxAD¯⁡wi,jn+1=0.The values on the boundaryAB¯ are equal to the interior points; therefore the inequality (45) can be applied.The values on the boundaryBC¯ are increasing with respect to the index j and reach the maximum at the point (Nx,Nx+Ny-1). From the other side, this value can be bounded by the value at the point (Nx,Nx+Ny) 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 S and, moreover, it is always greater than its asymptotic function (see [25], p. 268, formula (13.1)). The proposed transformation preserves monotonicity; therefore, (47)wNx,Nx+Nyn+1=FNxn+1Nα1,Nxn+1-ENα2,Nxn+1≥eσ1/ρ~xNx+r-q1τn+1-eσ2mxNx-yNx+Ny-1+E=wNx,Nx+Ny-1n+1.In summary we can conclude that(48)maxi,j⁡wi,jn+1=wNx,Nx+Nyn+1.This value is transformed call option price for the asset priceS1max, 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,(49)US1,0,τ=e-q1τS1Nd1-e-rτENd2≤e-qτS1Nd1≤e-qτS1,wNx,Nx+Nyn+1≤erτn+1e-q1τn+1S1max≤er-q1TS1max,0≤n≤Nτ-1.Then the following result can be stated.Theorem 2.
With previous notations, under conditions (35) and (36), numerical scheme (24), (29), (30), (31), (33), and (34) is conditionally uniformly ·∞-stable with upper bound A=er-q1TS1max.Now consistency of the numerical scheme (24) with PDE (19) will be studied [26].LetF(wi,jn)=0 be approximating difference equation (24). The scheme (24) is consistent with (19) if the local truncation error Ti,jn satisfies (50)Ti,jn=FWi,jn-LWi,jn⟶0ash⟶0,k⟶0,where Wi,jn denotes the theoretical solution of (19) evaluated at the point (xi,yj,τn) and L is the operator (51)LW=∂W∂τ-ρ~22∂2W∂x2+∂2W∂y2-ρ~σ1r-q1-σ122∂W∂x-ρσ1r-q1-σ122-1σ2r-q2-σ222∂W∂y.Assuming that theoretical solutionW of the PDE problem is four times continuously differentiable with respect to x and y and twice with respect to τ and using Taylor expansion about xi,yj,τn it is easy to check that(52)Wi,jn+1-Wi-1,jnk=∂W∂τxi,yj,τn+kEi,jn1,where (53)Ei,jn1=12∂2W∂τ2xi,yj,δ,nk<δ<n+1k,(54)Ei,jn1≤12Di,jn1max=12max∂2W∂τ2xi,yj,δ;0<δ<T.For the spatial discretization one gets(55)Wi+1,jn-Wi-1,jn2h=∂W∂xxi,yj,τn+h2Ei,jn2,where (56)Ei,jn2=16∂3W∂x3ξ,yj,τn,xi-1<ξ<xi+1,(57)Ei,jn2≤16Di,jn2max=16max⁡∂3W∂x3ξ,yj,τn;c1<ξ<d1.Analogously,(58)Wi,j+1n-Wi,j-1n2mh=∂W∂yxi,yj,τn+h2Ei,jn3, where (59)Ei,jn2=m26∂3W∂y3xi,η,τn,yj-1<η<yj+1;thus(60)Ei,jn3≤m26max⁡∂3W∂y3xi,η,τn;c2<η<d2.The discretization of the second spatial derivatives is as follows:(61)Wi+1,jn-2Wi,jn+Wi-1,jnh2=∂2W∂x2xi,yj,τn+h4Ei,jn4,where (62)Ei,jn4=112∂4W∂x4ξ,yj,τn,xi-1<ξ<xi+1,(63)Ei,jn4≤112max⁡∂4W∂x4ξ,yj,τn;c1<ξ<d1.Analogously,(64)wi,j+1n-2wi,jn+wi,j-1nm2h2=∂2W∂y2xi,yj,τn+h4Ei,jn5,where (65)Ei,jn5=m412∂4W∂y4xi,η,τn,yj-1<η<yj+1,(66)Ei,jn5≤m412max⁡∂3W∂y3xi,η,τn;c2<η<d2.Then the local truncation error takes the following form:(67)Ti,jn=FWi,jn-LWi,jn=kEi,jn1-ρ~22h4Ei,jn4+Ei,jn5-ρ~σ1r-q1-σ122h2Ei,jn2-ρσ1r-q1-σ122-1σ2r-q2-σ222h2Ei,jn3.From (54), (57), (60), (63), and (66) one gets(68)Ti,jn=Ok+Oh2.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).


## 5. 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(69)US1,S2,0=S2+E-S1+,E>0.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)
IfS2=0 the problem is reduced to the European put option on asset S1 and strike E. Therefore, the option value on this boundary is(70)US1,0,τ=Ee-rτN-ϕ2-S1eq1τN-ϕ1,whereϕ1 and ϕ2 are defined by (5).(2)
IfS1=0, we consider U(0,S2,τ) as a put option price on asset S1=0 with the strike S2+E:(71)U0,S2,τ=S2+Ee-rτ.(3)
For large values ofS2 we suppose that the slope of the curve tends to 1; that is,(72)∂U∂S2S1,S2⟶∞,τ=e-q2τ.(4)
ForS1≫S2+E  we consider U(S1,S2,τ) as a put option price with large values of S1, and therefore(73)US1,S2,τ=0.Under transformations (8) and (16) the boundary conditions (70)–(73) are translated to the rhomboid numerical domain ABCD as follows: (74)DC¯:wi,Ny+in=EN-α2,i+eσ1/ρ~xi-q1τnN-α1,i;(75)AD¯:w0,jn=E+e-σ2yj-mx0;(76)AB¯:wi,in=wi,i+1n+mhσ2e-σ2yi+1-mxi-rτ;(77)BC¯:wNx+1,Nx+jn=0.Positivity. As the numerical scheme (24) for the interior points is the same and the new boundary conditions (74)–(77) are nonnegative, then under conditions (35) and (36) the numerical solution for the European spread put option is nonnegative.Stability. Stability of the scheme is preserved with new boundary conditions (74)–(77), because analogously to the call option case the maximum value of the transformed function W(x,y,τ) is reached at a corner (in this case (c1,d2), where we use the boundary condition (75)). At the corresponding point of the original domain we use boundary condition (71). In the original coordinates it can be bounded as (78)U0,S2,τ=S2+Ee-rτ≤S2max+E=4E=const.Therefore, the transformed function also can be bounded as(79)maxi,j⁡wi,jn=4EerT.


## 6. 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):(80)∂U∂τ-LUU-g=0,∂U∂τ-LU≥0,U-gS1,S2≥0,where L represents the spatial differential operator of (1).If we rewrite (19) in the form (81)∂W∂τ-LW=0,then under transformation (16) LCP (80) has the following form:(82)∂W∂τ-LWW-fx,y=0,∂W∂τ-LW≥0;W-fx,y≥0,where L is the right-hand side of (19) and f(x,y) 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(83)w~i,jn=max⁡wi,jn,fxi,yj≥0,where wi,jn is defined by the finite difference equation (24).


## 7. Numerical Examples

In this section we illustrate numerical results for both European and American spread options.In Example1 we show that the stability conditions (35) and (36) for the European spread option can not be disregarded.Example 1.
We tested the algorithm for the European spread call option pricing problem with the parameters(84)T=1;E=100;σ1=0.2;σ2=0.1;r=0.1;q1=0.01;q2=0.05;ρ=0.1;and taking S1 and S2 large enough, in particular, S2max=400; S1max=3(S2max+E)=1500.For parameters (84) the stability conditions (35) and (36) on space and time steps are as follows:(85)h<min2.8428;23.7358=2.8428,k<0.0112.These conditions are crucial for the qualitative properties such as positivity and stability of the scheme. In Figure2 the numerical solution with h and k satisfying (85) is presented. On the other hand, in Figure 3 there is a numerical simulation with time step k=0.012>0.0112. The solution is unstable and it has negative values.Figure 2
Parametersh and k satisfy (85).Figure 3
Stability condition is broken.Example2 illustrates the computational time as well as the comparison with analytical approximation, provided by [6].Example 2.
Consider the European spread option with the data of Example1 (84) and choosing step sizes satisfying the stability condition.
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(86)US1,S2,T=S1e-q1TNd1-Ee-rT+S2e-r-r~-q~2TNd2,where (87)d1=1σTln⁡S1S2+Ee-rT+r-r~+q~2-q1+σ22T,d2=d1-σT,σ=σ12-2ρσ1σ2S2S2+Ee-rT+S2S2+Ee-rT2σ22,r~=S2S2+Ee-rTr,q~2=S2S2+Ee-rTq2.
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 S1 and S2 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.018 sec. It is the main part of computational time. In the case of FDM cdf takes less than half of computational time - 48.867 sec.
Furthermore, unsuitable negative values appear in analytical approximation for smallS1 and S2 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 Figure4.Table 1
Comparison proposed method (FDM) with analytical approximation (86).
S 2 S 1 FDM Approximation 100 400 210.5978 209.7261 200 24.2895 23.4543 100 0.1150 0.0123 50 3.1179 e - 05 - 4.7741 e - 10  90 400 220.0678 220.1714 200 29.8745 29.8724 100 0.1691 0.0372 50 5.2150 e - 05 - 1.9574 e - 09  CPU time (sec) 106.188 929.467Figure 4
European spread put option with parameters (84).Example 3.
With the parameters in (84) of Example 1, but with payoff (S2+E-S1)+ in Figure 4, we show the simulation of the numerical solution obtained by the proposed numerical scheme.Finally, in Example4 we compare the numerical solution obtained with our scheme (83) using LCP approach with other tested efficient methods presented in [28].Example 4.
Let us consider the American spread call option problem with parameters:(88)T=1;σ1=0.1;σ2=0.2;r=0.1;q1=0.05;q2=0.05;ρ=0.5.
Table2 shows the option price at fixed asset prices S1=96 and S2=100 for different values of strike price E 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)T=0.5;E=100;σ1=0.25;σ2=0.3;r=0.03;q1=0.06;q2=0.02;ρ=0.5.Table 2
Comparing several methods for American option pricing.
E Analytic FFT MC FDM 0 8.513201 8.513079 8.516613 8.508198 2 7.542296 7.542242 7.545496 7.534064 4 6.653060 6.652976 6.656364 6.648792In [9], three methods are considered to evaluate the spread option price for data 7.4 and fixed values S1=100 and S2=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 h=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.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).
Nτ NIM MOL Monte Carlo FDM 50 296.598 59.20 47 970 26.547 0.003 0.001 0.1735 0.0036  32 127.069 55.12 19 811 17.911 0.0035 0.006 0.1148 0.0122  16 33.493 51.718 5 149 9.432 0.0012 0.0236 0.132 0.0371  8 10.360 49.01 1 350 4.787 0.0105 0.074 0.1316 0.0907  4 4.117 45.47 355.75 2.815 0.0364 0.2512 0.1351 0.2254



---

*Source: 1549492-2016-12-07.xml*