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.


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 onefactor 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][4][5][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 2 The Scientific World Journal 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 -norm only, where < 1/2; hence, the convergence in -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][11][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 ( 2 ) operations ( is a number of space points), which is too time consuming. In [9][10][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 ( ln ) operations. The paper [12] uses the implicit scheme and the iteration method at each time step. The methods in [10][11][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.

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 ( ). A Lévy process can be completely specified by its characteristic exponent, , definable from the equality [ ( ) ] = − ( ) (we confine ourselves to the onedimensional 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 ( ) satisfies Assume that the riskless rate is constant, and, under a risk-neutral measure chosen by the market, the underlying asset evolves as = 0 , where is a Lévy process. Then we must have [ ] < +∞, and, therefore, must admit the analytic continuation into the strip Im ∈ (−1, 0) and continuous continuation into the closed strip Im ∈ [−1, 0].
where 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 (denote it by) 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.

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 > 0. Choose a space step Δ , and set = Δ , ∈ Z. Fix > 0 and apply any finite difference scheme to (11) (see, e.g., [9,11,12]), which may be reduced to approximate the infinitesimal generator as follows (the example of the reduction can be found in [12]): where Then we can approximate −1 ( − ) as follows: The sequence { } +∞ =−∞ generates doubly infinite Laurent matrix ( ) which is constant along the diagonals: After the discretization, the function ∈ 2 (R) turns into a piecewise constant function. Thus, we may consider = (. . . , ( −2 ), ( −1 ), ( 0 ), ( 1 ), ( 2 ), . . .) as an element of 2 (Z). Then we may rewrite (14) as follows: Let T stand for the complex unit circle. Since { } belongs to 1 (Z), we may introduce the function ( ) = ∑ , ∈ T, which is known as the symbol of the Laurent matrix 4 The Scientific World Journal or of the Laurent operator ( ). Recall that the family of all functions with absolutely converging Fourier series is the Wiener algebra := (T) (see details in [29]), which is a Banach algebra with respect to pointwise multiplication of functions and the norm ‖ ‖ = ∑ | |.
The standard theory of Toeplitz matrices (see details in Section 1.5, [29]) leads us to the following theorem.

Theorem 3. Let a function ∈ be able to be represented in the form
Then the operator ( ) is invertible and there exist + , − ∈ 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. (13) and (15) be the sequence of the Fourier coefficients of a ( ). Then ( ) satisfies conditions of Theorem 3.
Proof. Set̃( then we havẽ( According to (15) and (32), there exists a positive number 0 < 1 such that Hence, the Taylor series for ln(1 −̃( )) (ln(⋅) is the principal branch of the logarithm), converges at every point ∈ T. Set Becausẽ( ) ∈ and the condition (33) is satisfied, we conclude that ( ) is a fundamental sequence which is The Scientific World Journal 5 contained in the Wiener algebra . In fact, for any > 0, there exists a natural number such that Thus, we have proved that the function ( ) belongs to the Wiener algebra . Obviously, ( ) = exp( ( )), ∀ ∈ T.
Let a function ( ) satisfy the conditions of Proposition 4. Set ( ) = ( ( )) −1 , From Theorem 3 and Proposition 4, we deduce Further, we will describe an algorithm for finding coefficients { }, { ± }, based on the theory from [29]. By Theorem 3, It follows that { } is the sequence of the Fourier coefficients of −1 . We have due to (18). Wiener-Hopf factorization formula (28) gives The factors ± 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 ( ) = − ( ) , ∈ T. Next we define ± as Further, we set Finally, we have Obviously, + = 0 as < 0, and − = 0 as > 0.
Hence, to solve the problem (24), one needs to construct the inverse Toeplitz operator ( ) −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 ( ± ): An efficient numerical realization of (52) is available by means of fast Fourier transform (FFT) due to (19). The complexity of the method is ( ln ), where is the number of space discretization points.

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 ( ) be a monotonically increasing function, and it changes sign on the real line. Consider the following problem: The Scientific World Journal where the continuous function ( ) is maximized over barriers ℎ. Applying a finite difference scheme to (53) which approximates −1 ( − ) by formulas (14) and (15), we obtain the following discrete equation on the half line: where = ( ), = ( Δ ) −1 ( ), 0 maximizes , and is the symbol of a Laurent ( ) (see (12)-(16)).
where the only number 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.

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 and the parameters of the Lévy exponents (5).
Step 2. Input the space step Δ .
Step 3. Choose a finite difference scheme (FDS) for an approximation of the infinitesimal generator .
Step 5. Input the terminal date and define the number of time steps (the choice of typically depends on the finite difference scheme). Set space step Δ = / and = (Δ ) −1 + .
Step 6. Input min and max , the lower and upper bounds for the space variable . As a rule, the choice min = ln(0.4) and max = ln(2.5) is optimal.
Step 7. Define the number of space points as follows. Set 0 = max{− − ; + ; ( max − min )/2Δ }. We find integer number 0 such that 2 0 −1 < 0 ≤ 2 0 and set = 2 0 . 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.

Implementation of Wiener-Hopf Method for Solution to Standard Problems on Option Pricing
We assume that the riskless rate > 0 is constant and, under a risk-neutral measure chosen by the market, the log price of the stock = log follows a Lévy process with the infinitesimal generator (see (6)) and characteristic exponent (see (5)).

Barrier Options.
Consider a contract which pays the specified amount ( ) at the terminal date , provided during the life-time of the contract, the price of the stock does not cross a specified constant barrier from above (downand-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 8 The Scientific World Journal options without rebate; the generalization to the cases of a up-and-out barrier options and barrier options with rebate is straightforward. The price ( , ) of such barrier option can be found as the solution to the following integrodifferential equation with initial and boundary conditions (see [1]). Set = ln( / ), ( ) = ( ), and V( , ) = ( , ). Then, The most numerical methods start with a time discretization (the method of lines); see, for example, [5]. Divide [0, ] into subperiods by points = Δ , = 0, 1, . . . , , where Δ = / , and denote by V ( ) the approximation to V( , ). Then V ( ) = ( ), and by discretizing the derivative in (68), we obtain, for = − 1, − 2, . . . , 0, Equation (70) assumes the form Set = Δ −1 + ; 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 −1 ( − ) 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].

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 ( , ) be the price of American put with the strike price and the terminal date . Set = ln( / ), ( ) = (1 − ), and V( , ) = ( , ). Assume that the optimal stopping time is of the form ∧ , where is the hitting time of a closed set ⊂ R × (−∞, ] by the two-dimensional procesŝ= ( , ). Set C = R × [0, ) \ (this is the continuation region, where the option remains alive), and consider the following boundary value problem: where ( ) + := max{ ( ), 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, , 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, ] into subperiods by points = Δ , = 0, 1, . . . , , where Δ = / ; and denote by V ( ) the approximation to V( , ); ℎ denotes the approximation to the early exercise boundary at time . Then V ( ) = (1 − ) + , and by discretizing the derivative in (74), we obtain, for = − 1, − 2, . . . , 0, Equation (75) assumes the form The approximation ℎ to the early exercise boundary is found so that the V is maximal. (78) and (79): Notice that the sequence of problems (81) and (82) has the form (53), where ( ) is monotonically increasing function. Applying a finite difference scheme to (81) and (82) which approximates −1 ( − ) 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. 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.

The FDS&WH Method and the Method of
Example 1 (Process of order ] < 1). To compare CV method with FDS&WH for processes of order ] < 1, we take KoBoL The Scientific World Journal 9  We choose instantaneous interest rate = 0.04879, time to expiry = 0.5 year, strike price = 100, and the barrier = 90. As the base finite difference scheme we choose the one developed in [12].
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.
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. (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 = 1. We choose riskless rate = 0.03, time to expiry = 0.5 year, and strike price = 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.

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][10][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 ( ln ) operations, where 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 -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).