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.
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 payoff , depending on two assets and and strike price . 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: will denote a supremum norm for a given matrix .
2. Spread Option Pricing Problem and Boundary Conditions
The price 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 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. (1)For the payoff suggests Dirichlet’s condition 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 we take the value of the closed form solution of the basic Black-Scholes equation of a call for with given strike , where where is the standard normal cumulative distribution function. (3)As is large versus , then the behaviour of the solution looks like the asymptotic value of a one-dimensional vanilla call option with the strike for (see [22], p. 157), (4)For large values of we assume that the values are almost constants when changes; therefore we can use Neuman’s boundary condition:
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 by means of the substitutionobtaining
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 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 equationwhere and are the conjugate roots of . Solving (12) one getswhere one can relate the integration constant 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 , are small positive values; about and about . 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 stencilwhere
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 with small value of , we have transformed Black-Scholes solution, sowhere
Side corresponds to large values of , so behaviour of the option there leads to
Boundary corresponds with the boundary for large values of . Condition (7) means that component is constant and the null first derivative is approximated by the forward difference ,
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 .
From (25) and (27) it is easy to check that coefficients , , , and of the scheme (24) are positive under the condition
Coefficient is positive, if
Under conditions (35) and (36), values in interior points of the numerical domain , , calculated by the scheme (24), preserve nonnegativity from the previous time level at all interior and boundary points , .
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 according to (38) for the transformed problem for any .
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. Let be the matrix involving the numerical solution of PDE problem (19) computed by scheme (24), (29), (30), (31), (33), and (34),The numerical scheme (24), (29), (30), (31), (33), and (34) is said to be uniformly -stable in the domain , if for every partition with , as defined above, conditionholds true, where is independent of , . When the scheme is uniformly stable under a condition imposed to the step size discretization, then the scheme is said to be conditionally uniformly stable.
Let us find the maximum value of the with respect to , for each fixed . For 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 :
For the next time levels let us consider the numerical scheme (24). Under conditions (35) and (36) the coefficients , . Moreover, it is easy to check that
Therefore the values in interior points at the th level, , 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
The values on the boundary 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 . 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
This value is transformed call option price for the asset price , 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.
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 .
Now consistency of the numerical scheme (24) with PDE (19) will be studied [26].
Let be approximating difference equation (24). The scheme (24) is consistent with (19) if the local truncation error satisfies where denotes the theoretical solution of (19) evaluated at the point and is the operator
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 thatwhere
For the spatial discretization one gets where
Analogously, where thus
The discretization of the second spatial derivatives is as follows: where
Analogously, where
Then the local truncation error takes the following form:
From (54), (57), (60), (63), and (66) one gets
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
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 the problem is reduced to the European put option on asset and strike . Therefore, the option value on this boundary is where and are defined by (5). (2)If , we consider as a put option price on asset with the strike : (3)For large values of we suppose that the slope of the curve tends to 1; that is, (4)For we consider as a put option price with large values of , 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)–(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 is reached at a corner (in this case , 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
Therefore, the transformed function also can be bounded as
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):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 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 levelwhere 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 Example 1 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 parametersand taking and large enough, in particular, ; .
For parameters (84) the stability conditions (35) and (36) on space and time steps are as follows:
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 . 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].
Example 2. Consider the European spread option with the data of Example 1 (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 readswhere 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 and 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 small and 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 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].
Example 4. Let us consider the American spread call option problem with parameters:Table 2 shows the option price at fixed asset prices and 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]:
In [9], three methods are considered to evaluate the spread option price for data 7.4 and fixed values and : 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 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 . 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.
Competing Interests
The authors declare that there is no conflict of interests regarding the publication of this article.
Acknowledgments
This work has been partially supported by the European Union in the FP7- PEOPLE-2012-ITN Program under Grant Agreement no. 304617 (FP7 Marie Curie Action, Project Multi-ITN STRIKE-Novel Methods in Computational Finance) and the Ministerio de Economía y Competitividad Spanish Grant MTM2013-41765-P.