Finite difference methods for option pricing under Lévy processes: Wiener-Hopf factorization approach.

In the paper, we consider the problem of pricing options in wide classes of Lévy processes. We propose a general approach to the numerical methods based on a finite difference approximation for the generalized Black-Scholes equation. The goal of the paper is to incorporate the Wiener-Hopf factorization into finite difference methods for pricing options in Lévy models with jumps. The method is applicable for pricing barrier and American options. The pricing problem is reduced to the sequence of linear algebraic systems with a dense Toeplitz matrix; then the Wiener-Hopf factorization method is applied. We give an important probabilistic interpretation based on the infinitely divisible distributions theory to the Laurent operators in the correspondent factorization identity. Notice that our algorithm has the same complexity as the ones which use the explicit-implicit scheme, with a tridiagonal matrix. However, our method is more accurate. We support the advantage of the new method in terms of accuracy and convergence by using numerical experiments.

1. Introduction

In recent years more and more attention has been given to stochastic models of financial markets which depart from the traditional Black-Scholes model. We concentrate on one-factor non-Gaussian exponential Lévy models. These models provide a better fit to empirical asset price distributions that typically have fatter tails than Gaussian ones and can reproduce volatility smile phenomena in option prices. For an introduction to applications of these models applied to finance, we refer to [1, 2].

Option valuation under Lévy processes has been dealt with by a host of researchers; therefore, an exhaustive list is virtually impossible. However, the pricing of barrier options in exponential Lévy models still remains a mathematical and computational challenge (see, e.g., [3-6] for recent surveys of the state of the art of exotic option pricing in Lévy models).

The most general method to price barrier or American options deals with solving the corresponding partial integrodifferential equation (the generalized Black-Scholes equation) with appropriate boundary conditions. Note that in the case of American options free boundary problem arises. There are four main numerical methods for solving the partial integrodifferential equation (PIDE): multinomial trees, finite difference schemes, Galerkin methods, and numerical Wiener-Hopf factorization methods.

In [7], A family of Markov chain approximations of jump-diffusion models is constructed. Multinomial trees can be considered as special cases of explicit finite difference schemes. The main advantage of the method is simplicity of implementation; the drawbacks are inaccurate representation of the jumps and slow convergence.

Galerkin methods are based on the variational formulation of PIDE. While implementation of finite difference methods requires only a moderate programming knowledge, Galerkin methods use specialized toolboxes. Finite difference schemes use less memory than Galerkin methods, since there is no overhead for managing grids, but a refinement of the grid is more difficult. A wavelet Galerkin method for pricing American options under exponential Lévy processes is constructed in [8]. A general drawback of variational methods is that, for processes of finite variation, the convergence can be proved in the H -norm only, where s < 1/2; hence, the convergence in C-norm is not guaranteed.

In a finite difference scheme, derivatives are replaced by finite differences. In the presence of jumps, one needs to discretize the integral term as well. Finite difference schemes were applied to pricing of continuous barrier options in [9] and to pricing of American options in [10-12].

A construction of any finite difference scheme involves discretization in space and time, truncation of large jumps, and approximation of small jumps. Truncation of large jumps is necessary because an infinite sum cannot be calculated; approximation of small jumps is needed when Lévy measure diverges at zero. The result is a linear system that needs to be solved at each time step, starting from payoff function. In the general case, solution of the system on each time step by a linear solver requires O(m 2) operations (m is a number of space points), which is too time consuming. In [9-11], the integral part is computed using the solution from the previous time step, while the differential term is treated implicitly. This leads to the explicit-implicit scheme, with tridiagonal system which can be solved in O(mln⁡⁡m) operations. The paper [12] uses the implicit scheme and the iteration method at each time step. The methods in [10-12] are applicable to processes of infinite activity and finite variation; the part of the infinitesimal generator corresponding to small jumps is approximated by a differential operator of first order (additional drift component). The paper [9] uses an approximation by a differential operator of second order (additional diffusion component).

It follows from the analysis of the above methods for option pricing that in general case finite difference schemes seem to be the best choice. However, the essential disadvantage of the existing methods is speed and/or accuracy.

In [5], the fast and accurate numerical method for pricing barrier option in a wide class of Lévy processes was developed. The fast Wiener-Hopf factorization method (FWHF method) constructed in the paper is based on an efficient approximation of the Wiener-Hopf factors in the exact formula for the solution and the fast Fourier transform algorithm. In contrast to finite difference methods where the application entails an analysis of the underlying Lévy model, the FWHF method deals with the characteristic exponent of the process.

The method in [5] uses the interpretation of the factors as the expected present value operators (EPV operators)—integral operators suggested in [13] which calculate the (discounted) expected present values of streams of payoffs under supremum and infimum processes. This interpretation allows one to guess the optimal exercise boundary quite naturally and give a simple proof of optimality; see details in [14].

The goal of the paper is to incorporate the Wiener-Hopf factorization into finite difference methods for pricing options in Lévy models with jumps in terms of Laurent and Toeplitz matrices. The theory of Laurent and Toeplitz operators allows for solving linear algebraic systems related to the finite difference schemes sufficiently fast and accurate. Moreover, the correspondent matrix operators also admit probabilistic interpretation as expectation operators and they have similar properties to the ones of EPV-operators in [14]. It allows for developing effective methods for solving many standard problems of option pricing.

The method presented in the paper combines speed, simplicity, and accuracy. As our numerical examples show that it is rather faster than existing finite difference schemes. We generalize accurate finite difference scheme developed in [12] on processes of order more than 1 and describe the outline of the solution to the standard problems of option pricing.

The rest of the paper is organized as follows. In Section 2, we give necessary definitions of the theory of Lévy processes. In Section 3 we consider model problems related to the option pricing which can be reduced to solving Toeplitz systems. We provide the formulas for Wiener-Hopf factorization in terms of Laurent matrices and give the probabilistic interpretation to the factors. Section 4 incorporates the Wiener-Hopf factorization of Toeplitz matrices into finite difference methods for pricing barrier and American options. In Section 5, we produce numerical examples and compare several methods for pricing barrier and American options. Section 6 concludes the paper.

2. Lévy Processes: General Definitions

A Lévy process is a stochastically continuous process with stationary independent increments (for general definitions, see, e.g., [15]). A Lévy process may have a Gaussian component and/or pure jump component. The latter is characterized by the density of jumps, which is called the Lévy density. We denote it by F(dy). A Lévy process can be completely specified by its characteristic exponent, ψ, definable from the equality E[e ] = e − (we confine ourselves to the one-dimensional case). The characteristic exponent is given by the Lévy-Khintchine formula: where σ 2 and μ are the variance and drift coefficient of the Gaussian component and F(dy) satisfies

Assume that the riskless rate r is constant, and, under a risk-neutral measure chosen by the market, the underlying asset evolves as S = S 0 e , where X is a Lévy process. Then we must have E[e ]<+∞, and, therefore, ψ must admit the analytic continuation into the strip Im⁡ ξ ∈ (−1,0) and continuous continuation into the closed strip Im⁡ ξ ∈ [−1,0]. Further, the following condition (the EMM requirement) must hold: E[e ] = e . Equivalently, where r is instantaneous interest rate. The latter condition determines the drift via the other parameters of the Lévy process: Hence, the characteristic exponent may be rewritten as follows: Then the infinitesimal generator of X (denote it by) L is an integrodifferential operator which acts as follows:

In empirical studies of financial markets, the following classes of Lévy processes are popular: the Merton model [16], double-exponential jump-diffusion model (DEJD) introduced to finance by Lipton [17] and Kou [18], generalization of DEJD model constructed by Levendorskiǐ [19] and labeled later hyperexponential jump-diffusion model (HEJD), variance gamma processes (VGP) introduced to finance by Madan with coauthors (see, e.g., [20]), hyperbolic processes constructed in [21, 22], normal inverse Gaussian processes constructed by Barndorff-Nielsen [23] and generalized in [24], and extended Koponen's family introduced in [25, 26] and labeled KoBoL model in [1]. Koponen [27] introduced a symmetric version; Boyarchenko and Levendorskiǐ [25, 26] gave a nonsymmetric generalization; later, in [28], a subclass of this model appeared under the name CGMY model.

Example 1

The characteristic exponent of a pure jump KoBoL process of order ν ∈ (0,2), ν ≠ 1, is given by where c > 0, μ ∈ R, and λ − < −1 < 0 < λ +. Formula (7) is derived in [1, 25] from the Lévy-Khintchine formula with the Lévy densities of negative and positive jumps, F ∓(dy), given by

Example 2

In DEJD model, F ∓(dy) are given by exponential functions on negative and positive axis, respectively: where σ > 0, μ ∈ R, c ± > 0, and λ − < −1 < 0 < λ +. Then the characteristic exponent is of the form

3. Wiener-Hopf Factorization for Finite Difference Schemes

3.1. Wiener-Hopf Factorization for Finite Difference Schemes: Problems with a Barrier

Notice that many option pricing problems with a barrier can be reduced to the family of the following problems: where q > 0.

Choose a space step Δx, and set x = lΔx, l ∈ Z. Fix q > 0 and apply any finite difference scheme to (11) (see, e.g., [9, 11, 12]), which may be reduced to approximate the infinitesimal generator L as follows (the example of the reduction can be found in [12]): where Then we can approximate q −1(q − L) as follows:

The sequence {a } + generates doubly infinite Laurent matrix L(a) which is constant along the diagonals:

After the discretization, the function g ∈ L 2(R) turns into a piecewise constant function. Thus, we may consider g = (…, g(x −2), g(x −1), g(x 0), g(x 1), g(x 2),…) as an element of l 2(Z). Then we may rewrite (14) as follows:

Let T stand for the complex unit circle. Since {a } belongs to l 1(Z), we may introduce the function a(t) = ∑ a t , t ∈ T, which is known as the symbol of the Laurent matrix or of the Laurent operator L(a). Recall that the family of all functions with absolutely converging Fourier series is the Wiener algebra W : = W(T) (see details in [29]), which is a Banach algebra with respect to pointwise multiplication of functions and the norm ||a|| = ∑ | a |.

Further, the sequence {a } + is the sequence of the Fourier coefficients of a(t): Denote by F : L 2(T) → l 2(Z) the operator which maps a function a(t) to the sequence of its Fourier coefficients {a } + (see (18)).

It is well known that the Laurent matrix L(a) is the matrix representation of the multiplication by a(t) operator on L 2(T) with respect to the orthonormal basis . Hence, we have It follows from (19) that

According to the Wiener theorem, (see, e.g., [29]) if a ∈ W and a(t) ≠ 0, ∀t ∈ T, then a −1 = 1/a ∈ W. Let us denote by GW the set of all invertible elements of the algebra W. It follows from (20) that if a ∈ GW then the matrix L(a) is invertible with the inverse L(a −1).

The sequence {a } + also generates the infinite Toeplitz matrix T(a):

When a finite difference approximation is applied (see (14)), (11) may be rewritten as follows: where G = G(x ). Let P denote the orthogonal projection of l 2(Z) onto l 2(N): Then, taking into account that T(a) = PL(a)P, we rewrite (22) in terms of Toeplitz matrices: where G = (…, G(x −2), G(x −1), G(x 0), G(x 1), G(x 2),…) is considered as an element of l 2(Z).

The standard theory of Toeplitz matrices (see details in Section 1.5, [29]) leads us to the following theorem.

Theorem 3

Let a function a ∈ W be able to be represented in the form Then the operator T(a) is invertible and there exist a +, a − ∈ GW such that

Notice that the identities (29) and (28) are called a Wiener-Hopf factorization for Toeplitz matrices and Wiener functions, respectively.

In the context of the finite difference schemes under consideration one can prove the following proposition.

Proposition 4

Let {a } + defined by (13) and (15) be the sequence of the Fourier coefficients of a a(t). Then a(t) satisfies conditions of Theorem 3.

Proof

Set then we have According to (15) and (32), there exists a positive number r 0 < 1 such that Hence, the Taylor series for (ln⁡⁡(·) is the principal branch of the logarithm), converges at every point t ∈ T.

Set

Because and the condition (33) is satisfied, we conclude that ϕ (t) is a fundamental sequence which is contained in the Wiener algebra W. In fact, for any ϵ > 0, there exists a natural number k such that

Thus, we have proved that the function b(t) belongs to the Wiener algebra W. Obviously, a(t) = exp⁡(b(t)), ∀t ∈ T.

Let a function a(t) satisfy the conditions of Proposition 4. Set From Theorem 3 and Proposition 4, we deduce

Further, we will describe an algorithm for finding coefficients {p },  {p ±}, based on the theory from [29]. By Theorem 3, It follows that {p } is the sequence of the Fourier coefficients of a −1. We have due to (18).

Wiener-Hopf factorization formula (28) gives The factors p ± can be found as follows. From Proposition 4, there exists a function such that where the sequence of the Fourier coefficients of this function is defined as

Notice that p(t) = e −,  t ∈ T. Next we define b ± as Further, we set Finally, we have Obviously, p + = 0 as k < 0, and p − = 0 as k > 0.

Remark 5

Notice that the Wiener-Hopf factorization (44) is also satisfied if we substitute b −(t) + C and b +(t) − C into (49) (with any constant C) instead b −(t) and b +(t), respectively.

Hence, to solve the problem (24), one needs to construct the inverse Toeplitz operator T(a)−1 by using the above algorithm. Thus, an approximate solution to the problem (11) can be written as

From a practical point of view, it is more convenient to rewrite (51) in terms of Laurent operators L(p ±): An efficient numerical realization of (52) is available by means of fast Fourier transform (FFT) due to (19). The complexity of the method is O(Mln⁡M), where M is the number of space discretization points.

3.2. Wiener-Hopf Method for Finite Difference Schemes: Optimization Problems

Optimal stopping problems play a very important role in the mathematical finance and they are connected with pricing American, Bermudan, and other types of options. Pricing such options can be typically reduced to the sequence of the following problems.

Let G(x) be a monotonically increasing function, and it changes sign on the real line. Consider the following problem:   where the continuous function g(x) is maximized over barriers h.

Applying a finite difference scheme to (53) which approximates q −1(q − L) by formulas (14) and (15), we obtain the following discrete equation on the half line: where g = g(x ), G = (qΔt)−1 G(x ), k 0 maximizes g, and a is the symbol of a Laurent L(a) (see (12)–(16)).

Introduce the orthogonal projection P as follows: Further, we factorize the corresponding Toeplitz operator T(a) (see Theorem 3, Proposition 4, and (38)–(50)). The factorization formulas (48) are chosen in such a way that functions p, p +, and p − are characteristic functions of infinitely divisible distributions.

Theorem 6

Let sequences {p }, {p ±} be as defined by (43) and (47)–(50), and assume that the conditions of Proposition 4 are satisfied. Set Then P, P +, and P − are characteristic functions of infinitely divisible lattice distributions supported on {x = kΔx | k ∈ Z}, {x = kΔx | k ∈ Z, k ≥ 0}, and {x = kΔx | k ∈ Z, k ≤ 0}, respectively.

Recall that the characteristic function of an infinitely divisible lattice distribution with the maximal step Δx has the following form (see, e.g., [30]): where c ≥ 0, ∀k ∈ Z,  ∑ c < +∞.

From (35), we have Since b ∈ W and b(1) = 0, the function −b can be written as follows: where b are defined by (47).

Notice that by the definition (see (32)) , k ≠ 0. It follows that the Fourier coefficients of are also positive for every n due to Cauchy's series product theorem. Since , . Hence, all the coefficients −b in the formula (59) are positive.

Clearly, the functions f(t) ∈ W (t = e ∈ T) are continuous on T and, when regarded as functions f(e −), ξ ∈ R, they are (2π/Δx)-periodic continuous functions. Hence, the function p(t) (see (39) and (43)) can be rewritten in the form (57): Analogously, we rewrite P ±(ξ) in the form (57): It follows that P, P +, and P − are characteristic functions of infinitely divisible lattice distributions with the maximal step Δx supported on {x = kΔx | k ∈ Z}, {x = kΔx | k ∈ Z, k ≥ 0}, and {x = kΔx | k ∈ Z, k ≤ 0}, respectively.

From Theorem 6 we have that there exist discrete random variables X, X +, and X − taking on values of the form x = kΔx, k ∈ Z, such that where {p } and {p ±} are defined by (43) and (47)–(50).

It follows that the corresponding Laurent operators L(p), L(p +), and L(p −) can be interpreted as expectation operators conditioned on current values of X, X −, and X +, respectively:

The following simple properties are immediate from the interpretation of L(p ±) as expectation operators.

Proposition 7

Laurent operators L(p ±) enjoy the following properties.

If g = 0,  ∀k ≥ k 0, then, ∀k ≥ k 0, (L(p −)g) = 0.

If g = 0,  ∀k ≤ k 0, then, ∀k ≤ k 0, (L(p +)g) = 0.

If g ≥ 0,  ∀k, then, (L(p −)g) ≥ 0, ∀k. If, in addition, there exists k 0 such that g > 0  ∀k > k 0, then (L(p −)g) > 0,  ∀k.

If g ≥ 0,  ∀k, then (L(p +)g) ≥ 0, ∀x. If, in addition, there exists k 0 such that g > 0  ∀k < k 0, then (L(p +)g) > 0,  ∀k.

If g = {g } is monotone, then {(L(p −)g)} and {(L(p +)g)} are also monotone.

Proposition 7 is a direct analog of the properties of the expected present value operators introduced in [14]; see Proposition 6.2.1.

Taking into account that G = {G } in (54) is a monotonically increasing sequence and it changes the sign, then from Proposition 7 an approximate solution g = (…, g(x −2), g(x −1), g(x 0), g(x 1), g(x 2),…) to the problem (53) can be written in terms of Laurent operators L(p ±): where the only number k 0 can be found from the following conditions:   We remark that (64) includes the requirement that the series are convergent. An efficient numerical realization of (64) is based on (19) and fast Fourier transform.

3.3. Wiener-Hopf Factorization for Finite Difference Schemes: Algorithm

In the subsection, we give an algorithm of the construction of an approximate Wiener-Hopf factorization for finite difference schemes.

Wiener-Hopf Factorization

Step 1

Input the interest rate r and the parameters of the Lévy exponents (5).

Step 2

Input the space step Δx.

Step 3

Choose a finite difference scheme (FDS) for an approximation of the infinitesimal generator L.

Step 4

Choose desired truncation error ϵ for coefficients {α } in (12) (as a rule, the choice ϵ = 10−6 is optimal). Due to the FDS, calculate α , l = −1, −2,…, l −, where l − = max⁡⁡{l < 0 | |α | < ϵ/2}; calculate α , l = 1,2,…, l +, where l + = min⁡⁡{l > 0 | |α | < ϵ/2}. Set α = 0, as l < l − or l > l +.

Step 5

Input the terminal date T and define the number of time steps n (the choice of n typically depends on the finite difference scheme). Set space step Δt = T/n and q = (Δt)−1 + r.

Step 6

Input x min⁡ and x max⁡, the lower and upper bounds for the space variable x. As a rule, the choice x min⁡ = ln⁡⁡(0.4) and x max⁡ = ln⁡⁡(2.5) is optimal.

Step 7

Define the number of space points m as follows. Set l 0 = max⁡⁡{−l −; l +; (x max⁡ − x min⁡)/2Δx}. We find integer number k 0 such that 2 < l 0 ≤ 2 and set m = 2. We will use fast Fourier transform for real-valued functions (FFT); see details in [5, 31]. That is why we choose the number of space points as a power of 2.

Step 8

Find coefficients a ,  l = −m + 1,…, m, by the formula (15).

Step 9

Denote by τ = exp⁡⁡(iπk/m),  k = −m + 1,…, m. Find a(τ ) = ∑ a τ ,  k = −m + 1,…, m, using FFT.

Step 10

We find the symbol of L(p) = (L(a −1)): p(τ ) = a −1(τ ),  k = −m + 1,…, m  (see (39)).

Step 11

We find b(τ ): = ln⁡⁡(a(τ )), k = −m + 1,…, m (see (46)). Using inverse FFT, we obtain the sequence of coefficients b ,  k = −m + 1,…, m, for decomposition of b(τ) to the series

Step 12

Set b 0 − = −∑ b and b −(τ) = ∑ b τ + b 0 −; set b 0 + = b 0 − b 0 − and b +(τ) = ∑ b τ + b 0 + (see (48)). Using FFT we obtain b ±(τ ),  k = −m + 1,…, m.

Step 13

We find the symbols of L(p ±): p ±(τ ) = exp⁡(−b ±(τ )),  k = −m + 1,…, m (see (49)).

4. Implementation of Wiener-Hopf Method for Solution to Standard Problems on Option Pricing

We assume that the riskless rate r > 0 is constant and, under a risk-neutral measure chosen by the market, the log price of the stock X = log⁡S follows a Lévy process with the infinitesimal generator L (see (6)) and characteristic exponent ψ (see (5)).

4.1. Barrier Options

Consider a contract which pays the specified amount G(S ) at the terminal date T, provided during the life-time of the contract, the price of the stock does not cross a specified constant barrier H from above (down-and-out barrier options) or from below (up-and-out barrier options). When the barrier is crossed, the option expires worthless or the option owner is entitled to some rebate. We restrict ourselves to the case of down-and-out barrier options without rebate; the generalization to the cases of a up-and-out barrier options and barrier options with rebate is straightforward. The price V(t, S ) of such barrier option can be found as the solution to the following integrodifferential equation with initial and boundary conditions (see [1]). Set x = ln⁡⁡(S/H), g(x) = G(He ), and v(t, x) = V(t, He ). Then,

The most numerical methods start with a time discretization (the method of lines); see, for example, [5]. Divide [0, T] into n subperiods by points t = jΔt, j = 0,1,…, n, where Δt = T/n, and denote by v (x) the approximation to v(x, t ). Then v (x) = g(x), and by discretizing the derivative ∂ in (68), we obtain, for j = n − 1, n − 2,…, 0, Equation (70) assumes the form Set q = Δt −1 + r; then (71) can be rewritten as follows:

Notice that the sequence of problems (72) and (73) has the form (11). Applying a finite difference scheme to (72) and (73) which approximates q −1(q − L) by formulas (14) and (15), we obtain the discrete problem of the form (22) which can be easily solved by using (52). See details in Section 3.1. One can speed up the calculations by using real-valued FFT and similar tricks as in [5].

4.2. American Options

We consider the American put on a stock which pays no dividends; the generalization to the case of a dividend-paying stock and the American call is straightforward. (Moreover, as it is well-known, changing the direction on the line, the unknown function, the riskless rate, and the process, one can reduce the pricing problem for the American call to the pricing problem for the American put).

Let V(t, S ) be the price of American put with the strike price K and the terminal date T. Set x = ln⁡⁡(S/K), g(x) = K(1 − e ), and v(t, x) = V(t, Ke ). Assume that the optimal stopping time is of the form τ ′∧T, where τ ′ is the hitting time of a closed set B ⊂ R × (−∞, T] by the two-dimensional process . Set 𝒞 = R × [0, T)∖B (this is the continuation region, where the option remains alive), and consider the following boundary value problem: where g(x)+ : = max⁡{g(x), 0}.

Under certain regularity conditions (see Theorem 6.1 in [1]), the continuous bounded solution to the free boundary problem (74)–(77) gives the optimal early exercise region, B, and the rational option price, v.

We apply the Lévy analog of Carr's randomization procedure developed in Section 6.2.2 of [1] for the American put. Normalize the strike price to 1; divide [0, T] into n subperiods by points t = jΔt, j = 0,1,…, n, where Δt = T/n; and denote by v (x) the approximation to v(x, t ); h denotes the approximation to the early exercise boundary at time t . Then v (x) = K(1 − e )+, and by discretizing the derivative ∂ in (74), we obtain, for j = n − 1, n − 2,…, 0, Equation (75) assumes the form The approximation h to the early exercise boundary is found so that the v is maximal.

Introduce and substitute into (78) and (79):

Set q = Δt −1 + r and ; then (80) can be rewritten as follows:

Notice that the sequence of problems (81) and (82) has the form (53), where G (x) is monotonically increasing function. Applying a finite difference scheme to (81) and (82) which approximates q −1(q − L) by formulas (14) and (15), we obtain the discrete problem of form (54) which can be easily solved by using (64). See details in Section 3.2.

5. Numerical Examples

5.1. The FDS&WH Method and the Method of Cont and Voltchkova (2005)

In this subsection we apply our finite difference scheme with Wiener-Hopf method (we refer to the method FDS&WH) to KoBoL process and compare barrier option prices with the results obtained by the method in Cont and Voltchkova (2005) [9] (we refer to this method as CV method). We study convergence of two methods for processes of order ν < 1 and ν > 1.

Example 1 (Process of order ν < 1). To compare CV method with FDS&WH for processes of order ν < 1, we take KoBoL model with parameters σ = 0, ν = 0.5, λ + = 4.0, λ − = −6.0, and c = 1.0. We choose instantaneous interest rate r = 0.04879, time to expiry T = 0.5 year, strike price K = 100, and the barrier H = 90. As the base finite difference scheme we choose the one developed in [12].

In Table 1, we compare the down-and-out barrier put option prices calculated by FDS&WH and CV methods for spot prices S = 91,101,111,121 (the values are obtained on a PC with characteristics AMD Turion (tm) 64X2 1.6 GHz, 896 Mb, under Windows XP). We see that FDS&WH method demonstrates very fast convergence: in few seconds the accuracy reaches less than 0.5%. In the same time CV method converges very slowly and gives an error in 2-3% only after several hours of calculations. From Table 1 we clearly see that prices computed by FDS&WH stabilize sufficiently fast, while the ones computed by CV method essentially vary from the previous space step. Notice that near the barrier the prices computed by CV method are especially unstable.

Example 2 (Process of order ν > 1). In the case of processes of order ν > 1, we take KoBoL model with parameters σ = 0, ν = 1.2, λ + = 8.8, λ − = −14.5, and c = 1. We choose riskless rate r = 0.04879, time to expiry T = 0.1 year, strike price K = 100, and the barrier H = 80.

In Table 2, we compare the down-and-out barrier put option prices calculated by FDS&WH and CV methods. The FDS method can be extended for the case of infinite variation (ν > 1); see details in [32]. CV method demonstrates better convergence in comparison with the previous example, but FDS&WH method converges faster, especially in the neighborhood of the barrier.

5.2. The FDS&WH Method and FDS (Levendorskiǐ et al. (2006)), [12]

In this subsection we apply our finite difference scheme with Wiener-Hopf method to KoBoL process and compare American option prices with the results obtained by the finite difference method in [12] (we refer to this method as FDS method). We take KoBoL model with parameters σ = 0, ν = 0.2, λ + = 3.2, λ − = −5.4, and c = 1. We choose riskless rate r = 0.03, time to expiry T = 0.5 year, and strike price K = 100. The differences between prices and early exercise boundaries computed by the both methods are insignificant. Table 3 confirms our observation. As we see form Table 3 the time of computation by FDS&WH method is in several times lesser.

Table 3

American put, time of computation: FDS & WH versus FDS.

Parameters		Relative difference		Time of computation, sec.
Space step Δx	Number of time steps N τ	ϵ p	ϵ b	FDS	FDS & WH
0.002	65	0.21%	0.2%	14	1
0.001	112	0.21%	0.2%	64	3
0.0005	203	0.2%	0.2%	536	13

KoBoL parameters: σ = 0, ν = 0.2, λ + = 3.2, λ − = −5.4, c = 1.

K = 100, r = 0.03, T = 0.5.

ϵ and ϵ are the maximums of the relative differences between correspondent prices and boundaries, respectively, in the region S ≤ 1.3K.

6. Conclusion

Many option pricing problems can be solved by using finite difference schemes. The method is very popular in practice, because, in a diffusion model, the correspondent system has a tridiagonal matrix which can be easily inverted.

In the presence of jumps, we have the additional integral term which can be replaced by a discrete sum. As the result, one needs to invert a dense Toeplitz matrix. To avoid this problem, many authors (see, e.g., [9-11]) suggest computing the integral part by using the solution from the previous time step, while the differential term is treated implicitly. This leads to the explicit-implicit scheme, with tridiagonal system which can be solved in O(Mln⁡M) operations, where M is a number of space discretization points. However, the advantage in speed turns to the drawback in accuracy, especially in the case of barrier options. In the infinite activity case described in [9], the explicit-implicit scheme demonstrates bad convergence near the barrier and hence also becomes time consuming (see details in [5]).

In [12], an accurate implicit finite difference scheme for pricing American options was developed. The procedure of inversion for the dense matrix of the system is iterative, and it requires 5–10 iterations on each time step. Hence, for a fixed space and time steps, modification of the scheme for barrier options is in several times slower than that of the scheme in [9], but more accurate as examples in [5] show.

In [33], the case of discrete monitoring is considered. The usual backward recursion that arises in discrete barrier option pricing is converted into a set of independent integral equations by using a z-transform approach. In order to solve these equations, the rectangle quadrature rule transforms each integral equation into a Toeplitz linear system which is solved by iterative algorithms as in [12].

In the paper, we suggest a new approach which incorporates the Wiener-Hopf factorization method into a finite difference scheme with a Toeplitz system. Notice that our algorithm has the same complexity as the ones which use the explicit-implicit scheme, with a tridiagonal matrix. However, our method is more accurate, because it inverts the whole Toeplitz matrix, not only its tridiagonal part.

We give an important probabilistic interpretation based on the infinitely divisible distributions theory to the Laurent operators in the Wiener-Hopf factorization identity for finite difference schemes. It implies very useful properties which allow developing effective methods for solving many standard problems on option pricing (e.g., European, barrier, first touch digital, and American options).

(a)  

Parameters		CV			FDS & WH
Δx	N τ	Option price	ϵ	CPU time, sec.	Option price	ϵ	CPU time, sec.
0.001	93	0.0577	−59.1%	1	0.1380	−2.2%	1
0.0005	152	0.0716	−49.2%	3	0.1400	−0.8%	2
0.00025	253	0.0873	−38.1%	19	0.1408	−0.2%	8
0.0001	520	0.1073	−24.0%	78	0.1411	0.0%	33
0.00005	926	0.1197	−15.2%	324	0.1411	0.0%	126
0.000025	1688	0.1281	−9.2%	1348	0.1411		568
0.00001	4000	0.1330	−5.7%	14655

(b)  

Parameters		CV			FDS & WH
Δx	N τ	Option price	ϵ	CPU time, sec.	Option price	ϵ	CPU time, sec.
0.001	93	0.2344	−19.8%	1	0.2899	−0.8%	1
0.0005	152	0.2464	−15.7%	3	0.2915	−0.2%	2
0.00025	253	0.2571	−12.0%	19	0.2920	0.0%	8
0.0001	520	0.2679	−8.3%	78	0.2922	0.0%	33
0.00005	926	0.2740	−6.2%	324	0.2922	0.0%	126
0.000025	1688	0.2787	−4.6%	1348	0.2922		568
0.00001	4000	0.2832	−3.1%	14655

(c)  

Parameters		CV			FDS & WH
Δx	N τ	Option price	ϵ	CPU time, sec.	Option price	ϵ	CPU time, sec.
0.001	93	0.2182	−16.7%	1	0.2583	−1.4%	1
0.0005	152	0.2273	−13.3%	3	0.2604	−0.6%	2
0.00025	253	0.2353	−10.2%	19	0.2614	−0.3%	8
0.0001	520	0.2434	−7.1%	78	0.2618	−0.1%	33
0.00005	926	0.2481	−5.3%	324	0.2620	0.0%	126
0.000025	1688	0.2517	−4.0%	1348	0.2621		568
0.00001	4000	0.2552	−2.6%	14655

(d)  

Parameters		CV			FDS & WH
Δx	N τ	Option price	ϵ	CPU time, sec.	Option price	ϵ	CPU time, sec.
0.001	93	0.1718	−15.4%	1	0.1995	−1.7%	1
0.0005	152	0.1781	−12.2%	3	0.2012	−0.9%	2
0.00025	253	0.1838	−9.4%	19	0.2022	−0.4%	8
0.0001	520	0.1896	−6.6%	78	0.2027	−0.2%	33
0.00005	926	0.1930	−4.9%	324	0.2029	−0.1%	126
0.000025	1688	0.1956	−3.6%	1348	0.2030		568
0.00001	4000	0.1981	−2.4%	14655

KoBoL parameters: σ = 0, ν = 0.5, λ + = 4, λ − = −6, c = 1.

K = 100, H = 90, r = 0.04879, T = 0.5, ϵ: the relative difference between the current option price and the price computed by FDS & WH method for space step Δx = 0.000025.

(a): S = 91; (b): S = 101; (c): S = 111; (d): S = 121.

(a)  

Parameters		CV			FDS & WH
Δx	N τ	Option price	ϵ	CPU time, sec.	Option price	ϵ	CPU time, sec.
0.001	95	0.372	−28.9%	1	0.551	5.3%	1
0.0005	218	0.402	−23.3%	2	0.543	3.6%	3
0.00025	501	0.431	−17.7%	10	0.536	2.3%	15
0.0001	1501	0.464	−11.4%	188	0.529	1.0%	98
0.00005	3440	0.481	−8.2%	1045	0.526	0.4%	482
0.000025	7882	0.492	−6.1%	5520	0.524		2697

(b)  

Parameters		CV			FDS & WH
Δx	N τ	Option price	ϵ	CPU time, sec.	Option price	ϵ	CPU time, sec.
0.001	95	2.314	−5.7%	1	2.565	4.5%	1
0.0005	218	2.344	−4.5%	2	2.522	2.7%	3
0.00025	501	2.369	−3.5%	10	2.494	1.6%	15
0.0001	1501	2.394	−2.5%	188	2.471	0.7%	98
0.00005	3440	2.407	−1.9%	1045	2.461	0.3%	482
0.000025	7882	2.418	−1.5%	5520	2.455		2697

(c)  

Parameters		CV			FDS & WH
Δx	N τ	Option price	ϵ	CPU time, sec.	Option price	ϵ	CPU time, sec.
0.001	95	2.304	−3.4%	1	2.428	1.8%	1
0.0005	218	2.320	−2.7%	2	2.414	1.2%	3
0.00025	501	2.335	−2.1%	10	2.403	0.7%	15
0.0001	1501	2.349	−1.5%	188	2.393	0.3%	98
0.00005	3440	2.357	−1.2%	1045	2.389	0.1%	482
0.000025	7882	2.364	−0.9%	5520	2.386		2697

(d)  

Parameters		CV			FDS & WH
Δx	N τ	Option price	ϵ	CPU time, sec.	Option price	ϵ	CPU time, sec.
0.001	95	1.547	−2.2%	1	1.567	−0.9%	1
0.0005	218	1.554	−1.8%	2	1.575	−0.4%	3
0.00025	501	1.560	−1.4%	10	1.579	−0.2%	15
0.0001	1501	1.566	−1.0%	188	1.581	0.0%	98
0.00005	3440	1.570	−0.7%	1045	1.581	0.0%	482
0.000025	7882	1.573	−0.5%	5520	1.581		2697

(e)  

Parameters		CV			FDS & WH
Δx	N τ	Option price	ϵ	CPU time, sec.	Option price	ϵ	CPU time, sec.
0.001	95	0.862	−1.5%	1	0.849	−3.0%	1
0.0005	218	0.865	−1.2%	2	0.861	−1.7%	3
0.00025	501	0.867	−0.9%	10	0.868	−0.9%	15
0.0001	1501	0.870	−0.6%	188	0.872	−0.3%	98
0.00005	3440	0.872	−0.4%	1045	0.874	−0.1%	482
0.000025	7882	0.873	−0.3%	5520	0.875		2697

KoBoL parameters: σ = 0, ν = 1.2, λ + = 8.8, λ − = −14.5, c = 1.

K = 100, H = 80, r = 0.04879, T = 0.1; N : number of time steps; ϵ: the relative difference between the current option price and the price computed by FDS & WH method for space step Δx = 0.000025.

(a): S = 81; (b): S = 91; (c): S = 101; (d): S = 111; (e): S = 121.