Finite Difference Methods for Option Pricing under Lévy Processes: Wiener-Hopf Factorization Approach

Abstract

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 𝐻^𝑠-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–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 𝑂(𝑚²) operations (𝑚 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 𝑂(𝑚ln𝑚) 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 𝐹(𝑑𝑦). A Lévy process can be completely specified by its characteristic exponent, 𝜓, definable from the equality 𝐸[𝑒^{𝑖𝜉𝑋(𝑡)}]=𝑒^{−𝑡𝜓(𝜉)} (we confine ourselves to the one-dimensional case). The characteristic exponent is given by the Lévy-Khintchine formula: 𝜓(𝜉)=𝜎²2𝜉²−𝑖𝜇𝜉+^+∞_−∞1−𝑒^𝑖𝜉𝑦+𝑖𝜉𝑦𝟏_|𝑦|≤1𝐹𝑑𝑦, (1) where 𝜎² and 𝜇 are the variance and drift coefficient of the Gaussian component and 𝐹(𝑑𝑦) satisfies _𝐑⧵{0}min1,𝑦²𝐹𝑑𝑦<+∞. (2)

Assume that the riskless rate 𝑟 is constant, and, under a risk-neutral measure chosen by the market, the underlying asset evolves as 𝑆_𝑡=𝑆₀𝑒^𝑋𝑡, 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]. Further, the following condition (the EMM requirement) must hold: 𝐸[𝑒^𝑋𝑡]=𝑒^𝑟𝑡. Equivalently, 𝑟+𝜓(−𝑖)=0, (3) where 𝑟 is instantaneous interest rate. The latter condition determines the drift via the other parameters of the Lévy process: 𝜎𝜇=𝑟−²2+^+∞_−∞1−𝑒^𝑦+𝑦𝟏_|𝑦|≤1𝐹𝑑𝑦. (4) Hence, the characteristic exponent may be rewritten as follows: 𝜓(𝜉)=𝜎²2𝜉²𝜎−𝑖𝑟−²2𝜉+^+∞_−∞1−𝑒^𝑖𝜉𝑦−𝑖𝜉1−𝑒^𝑦𝐹𝑑𝑦. (5) Then the infinitesimal generator of 𝑋 (denote it by) 𝐿 is an integrodifferential operator which acts as follows: (𝑥)=𝜎𝐿𝑢²2𝑢^′′(𝑥)+𝜎𝑟−²2𝑢^′(𝑥)+^+∞_−∞𝑢𝑥+𝑦−𝑢(𝑥)−𝑒^𝑦−1𝑢^′(𝑥)𝐹𝑑𝑦. (6)

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 𝜓(𝜉)=−𝑖𝜇𝜉+𝑐Γ(−𝜈)𝜆^𝜈₊−𝜆₊+𝑖𝜉^𝜈+−𝜆₋^𝜈−−𝜆₋−𝑖𝜉^𝜈, (7) where 𝑐>0, 𝜇∈𝐑, 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, 𝐹_∓(𝑑𝑦), given by 𝐹_∓𝑑𝑦=𝑐𝑒^𝜆±𝑦𝑦^−𝜈−1𝑑𝑦. (8)

Example 2. In DEJD model, 𝐹_∓(𝑑𝑦) are given by exponential functions on negative and positive axis, respectively: 𝐹_∓𝑑𝑦=𝑐_±±𝜆_±𝑒^𝜆±𝑦, (9) where 𝜎>0, 𝜇∈𝐑, 𝑐_±>0, and 𝜆₋<−1<0<𝜆₊. Then the characteristic exponent is of the form 𝜓(𝜉)=𝜎²2𝜉²−𝑖𝜇𝜉+𝑖𝑐₊𝜉𝜆₊++𝑖𝜉𝑖𝑐₋𝜉𝜆₋.+𝑖𝜉 (10)

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: 𝑞⁻¹𝑔(𝑥)(𝑥)𝑔(𝑥)𝑞−𝐿=𝐺,𝑥>0,=0,𝑥≤0, (11) where 𝑞>0.

Choose a space step Δ𝑥, and set 𝑥_𝑙=𝑙Δ𝑥, 𝑙∈𝐙. 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]): 𝐿𝑔𝑥_𝑘=_𝑙≠0𝛼_𝑙𝑔𝑥_𝑘+𝑙−_𝑙≠0𝛼_𝑙𝑔𝑥_𝑘, (12) where _𝑙≠0𝛼_𝑙<∞,𝛼_𝑙>0,𝑙≠0. (13) Then we can approximate 𝑞⁻¹(𝑞−𝐿) as follows: 𝑞⁻¹𝑞−𝐿𝑔𝑥_𝑘=_𝑙∈𝐙𝑎_𝑙𝑔𝑥_𝑘−𝑙,𝑎_𝑙=−𝑞⁻¹𝛼_−𝑙,𝑙≠0,𝑎₀=1+𝑞⁻¹_𝑙≠0𝛼_𝑙. (14)(15)

The sequence {𝑎_𝑙}^+∞_𝑙=−∞ generates doubly infinite Laurent matrix 𝐿(𝑎) which is constant along the diagonals: 𝐿(𝑎)=⋯⋯⋯⋯⋯⋯⋯𝑎₀𝑎₋₁𝑎₋₂𝑎₋₃⋯⋯𝑎₁𝑎₀𝑎₋₁𝑎₋₂⋯⋯𝑎₂𝑎₁𝑎₀𝑎₋₁⋯⋯𝑎₃𝑎₂𝑎₁𝑎₀⋯⋯⋯⋯⋯⋯⋯. (16)

After the discretization, the function 𝑔∈𝐿₂(𝐑) turns into a piecewise constant function. Thus, we may consider 𝑔=(…,𝑔(𝑥₋₂),𝑔(𝑥₋₁),𝑔(𝑥₀),𝑔(𝑥₁),𝑔(𝑥₂),…) as an element of 𝑙₂(𝐙). Then we may rewrite (14) as follows: 𝑞⁻¹𝑞−𝐿𝑔𝑥_𝑘=𝐿(𝑎)𝑔_𝑘,𝑘∈𝐙. (17)

Let 𝐓 stand for the complex unit circle. Since {𝑎_𝑙} belongs to 𝑙₁(𝐙), we may introduce the function 𝑎(𝑡)=∑_𝑘𝑎_𝑘𝑡^𝑘, 𝑡∈𝐓, which is known as the symbol of the Laurent matrix or of the Laurent operator 𝐿(𝑎). Recall that the family of all functions with absolutely converging Fourier series is the Wiener algebra 𝑊∶=𝑊(𝐓) (see details in [29]), which is a Banach algebra with respect to pointwise multiplication of functions and the norm ‖𝑎‖_𝑊=∑_𝑘|𝑎_𝑘|.

Further, the sequence {𝑎_𝑘}^+∞_𝑘=−∞ is the sequence of the Fourier coefficients of 𝑎(𝑡): 𝑎_𝑘=12𝜋₀^2𝜋𝑎𝑒^𝑖𝜑𝑒^{−𝑖𝑘𝜑}𝑑𝜑,𝑘∈𝐙. (18) Denote by 𝐹∶𝐿₂(𝐓)→𝑙₂(𝐙) the operator which maps a function 𝑎(𝑡) to the sequence of its Fourier coefficients {𝑎_𝑘}^+∞_𝑘=−∞ (see (18)).

It is well known that the Laurent matrix 𝐿(𝑎) is the matrix representation of the multiplication by 𝑎(𝑡) operator on 𝐿₂(𝐓) with respect to the orthonormal basis {(1/√2𝜋)𝑒^𝑖𝑘𝜑}. Hence, we have 𝐿(𝑎)=𝐹𝑎𝐹⁻¹. (19) It follows from (19) that 𝐿𝑎₁𝐿𝑎₂=𝐿𝑎₁𝑎₂,∀𝑎₁,𝑎₂∈𝐿_∞(𝐓). (20)

According to the Wiener theorem, (see, e.g., [29]) if 𝑎∈𝑊 and 𝑎(𝑡)≠0, ∀𝑡∈𝐓, then 𝑎⁻¹=1/𝑎∈𝑊. Let us denote by 𝐺𝑊 the set of all invertible elements of the algebra 𝑊. It follows from (20) that if 𝑎∈𝐺𝑊 then the matrix 𝐿(𝑎) is invertible with the inverse 𝐿(𝑎⁻¹).

The sequence {𝑎_𝑙}^+∞_𝑙=−∞ also generates the infinite Toeplitz matrix 𝑇(𝑎): 𝑇(𝑎)=𝑎₀𝑎₋₁𝑎₋₂𝑎₋₃⋯𝑎₁𝑎₀𝑎₋₁𝑎₋₂⋯𝑎₂𝑎₁𝑎₀𝑎₋₁⋯𝑎₃𝑎₂𝑎₁𝑎₀⋯⋯⋯⋯⋯⋯. (21)

When a finite difference approximation is applied (see (14)), (11) may be rewritten as follows: 𝐿(𝑎)𝑔_𝑘=𝐺_𝑘𝑔,𝑘∈𝐍,_𝑘=0,𝑘∈𝐙,𝑘≤0, (22) where 𝐺_𝑘=𝐺(𝑥_𝑘). Let 𝑃 denote the orthogonal projection of 𝑙₂(𝐙) onto 𝑙₂(𝐍): 𝑃𝑢_𝑘=𝑢_𝑘,𝑘>0,0,𝑘≤0. (23) Then, taking into account that 𝑇(𝑎)=𝑃𝐿(𝑎)𝑃, we rewrite (22) in terms of Toeplitz matrices: 𝑇(𝑎)𝑔=𝑃𝐺, (24) where 𝐺=(…,𝐺(𝑥₋₂),𝐺(𝑥₋₁),𝐺(𝑥₀),𝐺(𝑥₁),𝐺(𝑥₂),…) is considered as an element of 𝑙₂(𝐙).

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 𝑎=exp(𝑏),𝑏∈𝑊. (25) Then the operator 𝑇(𝑎) is invertible and there exist 𝑎₊,𝑎₋∈𝐺𝑊 such that 𝑎₊(𝑡)=_𝑘≥0𝑎⁺_𝑘𝑡^𝑘𝑎,𝑡∈𝐓,₋(𝑡)=_𝑘≤0𝑎⁻_𝑘𝑡^𝑘,𝑡∈𝐓,𝑎=𝑎₊𝑎₋,𝑇(𝑎)=𝑇𝑎₋𝑇𝑎₊,𝑇(𝑎)⁻¹=𝑇𝑎₊⁻¹𝑇𝑎₋⁻¹. (26)(27)(28)(29)(30)

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 {𝑎_𝑙}^+∞_𝑙=−∞ defined by (13) and (15) be the sequence of the Fourier coefficients of a 𝑎(𝑡). Then 𝑎(𝑡) satisfies conditions of Theorem 3.

Proof. Set 𝑎(𝑡)𝑎(𝑡)=1−𝑎₀,𝑡∈𝐓, (31) then we have 𝑎(𝑡)=_𝑙≠0𝑎_𝑙𝑡^𝑙,𝑎_𝑙𝑎=−^𝑙𝑎₀. (32) According to (15) and (32), there exists a positive number 𝑟₀<1 such that ‖𝑎(𝑡)‖_𝑊<𝑟₀,∀𝑡∈𝐓. (33) Hence, the Taylor series for ln(1−𝑎(𝑡)) (ln(⋅) is the principal branch of the logarithm), −_𝑛>0𝑎(𝑡)^𝑛𝑛 (34) converges at every point 𝑡∈𝐓.
Set 𝑏(𝑡)=ln(1−𝑎(𝑡))+ln𝑎₀,𝜙_𝑘(𝑡)=−^𝑘_𝑛=1𝑎(𝑡)^𝑛𝑛+ln𝑎₀,𝑘∈𝐍. (35)(36)
Because 𝑎(𝑡)∈𝑊 and the condition (33) is satisfied, we conclude that 𝜙_𝑘(𝑡) is a fundamental sequence which is contained in the Wiener algebra 𝑊. In fact, for any 𝜖>0, there exists a natural number 𝑘 such that 𝜙_𝑘+𝑚−𝜙_𝑘𝑊≤^𝑘+𝑚_𝑛=𝑘+1𝑟^𝑛₀𝑛≤𝑟₀^𝑘+1(𝑘+1)1−𝑟₀<𝜖,∀𝑚∈𝐍. (37)
Thus, we have proved that the function 𝑏(𝑡) belongs to the Wiener algebra 𝑊. Obviously, 𝑎(𝑡)=exp(𝑏(𝑡)), ∀𝑡∈𝐓.

Let a function 𝑎(𝑡) satisfy the conditions of Proposition 4. Set 𝑝(𝑡)=(𝑎(𝑡))⁻¹,𝑝⁺(𝑡)=𝑎₊(𝑡)⁻¹,𝑝⁻(𝑡)=𝑎₋(𝑡)⁻¹. (38) From Theorem 3 and Proposition 4, we deduce 𝑝(𝑡)=_𝑙∈𝐙𝑝_𝑙𝑡^𝑙𝑝,𝑝∈𝑊,₋(𝑡)=^𝑙=0_𝑙=−∞𝑝⁻_𝑙𝑡^𝑙,𝑝⁻𝑝∈𝑊,₊(𝑡)=^𝑙=+∞_𝑙=0𝑝⁺_𝑙𝑡^𝑙,𝑝⁺∈𝑊. (39)(40)(41)

Further, we will describe an algorithm for finding coefficients {𝑝_𝑙},  {𝑝^±_𝑙}, based on the theory from [29]. By Theorem 3, 𝐿𝑝=𝐿(𝑎)⁻¹,𝑇𝑝⁺=𝑇𝑎₊⁻¹,𝑇𝑝⁻=𝑇𝑎₋⁻¹. (42) It follows that {𝑝_𝑙} is the sequence of the Fourier coefficients of 𝑎⁻¹. We have 𝑝_𝑘=12𝜋^𝜋_−𝜋𝑎𝑒^𝑖𝜑−1𝑒^{−𝑖𝑘𝜑}𝑑𝜑,𝑘∈𝐙, (43) due to (18).

Wiener-Hopf factorization formula (28) gives 𝑝(𝑡)=𝑝⁺(𝑡)𝑝⁻(𝑡),∀𝑡∈𝐓. (44) The factors 𝑝^± can be found as follows. From Proposition 4, there exists a function 𝑏(𝑡)=^+∞_𝑘=−∞𝑏_𝑘𝑡^𝑘,^+∞_𝑘=−∞𝑏_𝑘<∞, (45) such that 𝑏(𝑡)=ln𝑎(𝑡),𝑡∈𝐓, (46) where the sequence of the Fourier coefficients of this function is defined as 𝑏_𝑘=12𝜋^𝜋_−𝜋ln𝑎𝑒^𝑖𝜑𝑒^{−𝑖𝑘𝜑}𝑑𝜑,𝑘∈𝐙. (47)

Notice that 𝑝(𝑡)=𝑒^{−𝑏(𝑡)},  𝑡∈𝐓. Next we define 𝑏_± as 𝑏₋(𝑡)=^𝑘=−1_𝑘=−∞𝑏_𝑘𝑡^𝑘,𝑏−1₊(𝑡)=^𝑘=+∞_𝑘=1𝑏_𝑘𝑡^𝑘−1,𝑡∈𝐓, (48) Further, we set 𝑝₊(𝑡)=𝑒^{−𝑏+(𝑡)},𝑝₋(𝑡)=𝑒^{−𝑏−(𝑡)},𝑡∈𝐓. (49) Finally, we have 𝑝^±_𝑘=12𝜋^𝜋_−𝜋𝑝_±𝑒^𝑖𝜑𝑒^{−𝑖𝑘𝜑}𝑑𝜑,𝑘∈𝐙. (50) Obviously, 𝑝⁺_𝑘=0 as 𝑘<0, and 𝑝⁻_𝑘=0 as 𝑘>0.

Remark 5. Notice that the Wiener-Hopf factorization (44) is also satisfied if we substitute 𝑏₋(𝑡)+𝐶 and 𝑏₊(𝑡)−𝐶 into (49) (with any constant 𝐶) instead 𝑏₋(𝑡) and 𝑏₊(𝑡), respectively.

Hence, to solve the problem (24), one needs to construct the inverse Toeplitz operator 𝑇(𝑎)⁻¹ by using the above algorithm. Thus, an approximate solution to the problem (11) can be written as 𝑔=𝑇𝑎₊⁻¹𝑇𝑎₋⁻¹𝐺. (51)

From a practical point of view, it is more convenient to rewrite (51) in terms of Laurent operators 𝐿(𝑝^±): 𝑔=𝐿𝑝₊𝑃𝐿𝑝₋𝐺. (52) 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.

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 𝐺(𝑥) be a monotonically increasing function, and it changes sign on the real line. Consider the following problem:𝑞⁻¹𝑔(𝑥)(𝑥)𝑔(𝑥)𝑞−𝐿=𝐺,𝑥>ℎ,=0,𝑥≤ℎ, (53) where the continuous function 𝑔(𝑥) is maximized over barriers ℎ.

Applying a finite difference scheme to (53) which approximates 𝑞⁻¹(𝑞−𝐿) by formulas (14) and (15), we obtain the following discrete equation on the half line: 𝐿(𝑎)𝑔_𝑘=𝐺^𝑘,𝑘>𝑘₀,𝑔_𝑘=0,𝑘≤𝑘₀, (54) where 𝑔_𝑘=𝑔(𝑥_𝑘), 𝐺^𝑘=(𝑞Δ𝑡)⁻¹𝐺(𝑥_𝑘), 𝑘₀ maximizes 𝑔, and 𝑎 is the symbol of a Laurent 𝐿(𝑎) (see (12)–(16)).

Introduce the orthogonal projection 𝑃_𝑙 as follows: 𝑃_𝑙𝑢_𝑘=𝑢_𝑘,𝑘>𝑙,0,𝑘≤𝑙. (55) Further, we factorize the corresponding Toeplitz operator 𝑇(𝑎) (see Theorem 3, Proposition 4, and (38)–(50)). The factorization formulas (48) are chosen in such a way that functions 𝑝, 𝑝⁺, and 𝑝⁻ are characteristic functions of infinitely divisible distributions.

Theorem 6. Let sequences {𝑝_𝑙}, {𝑝^±_𝑙} be as defined by (43) and (47)–(50), and assume that the conditions of Proposition 4 are satisfied. Set 𝑃(𝜉)=_𝑙∈𝐙𝑝_𝑙()𝑃exp−𝑖𝑙𝜉Δ𝑥,𝜉∈𝐑,₊(𝜉)=^𝑙=0_𝑙=−∞𝑝⁻_𝑙()𝑃exp−𝑖𝑙𝜉Δ𝑥,𝜉∈𝐑,₋(𝜉)=^𝑙=+∞_𝑙=0𝑝⁺_𝑙exp(−𝑖𝑙𝜉Δ𝑥),𝜉∈𝐑. (56) Then 𝑃, 𝑃⁺, and 𝑃⁻ are characteristic functions of infinitely divisible lattice distributions supported on {𝑥=𝑘Δ𝑥∣𝑘∈𝐙}, {𝑥=𝑘Δ𝑥∣𝑘∈𝐙,𝑘≥0}, and {𝑥=𝑘Δ𝑥∣𝑘∈𝐙,𝑘≤0}, respectively.

Proof. Recall that the characteristic function of an infinitely divisible lattice distribution with the maximal step Δ𝑥 has the following form (see, e.g., [30]): 𝑐(𝜉)=exp_𝑘∈𝐙𝑐_𝑘exp(𝑖𝑘𝜉Δ𝑥)−1,𝜉∈𝐑, (57) where 𝑐_𝑘≥0,∀𝑘∈𝐙,  ∑_𝑘∈𝐙𝑐_𝑘<+∞.
From (35), we have −𝑏(𝑡)=^+∞_𝑛=1𝑎(𝑡)^𝑛𝑛−ln𝑎₀. (58) Since 𝑏∈𝑊 and 𝑏(1)=0, the function −𝑏 can be written as follows: −𝑏(𝑡)=^+∞_𝑘=−∞−𝑏_𝑘𝑡^𝑘−1, (59) where 𝑏_𝑘 are defined by (47).
Notice that by the definition (see (32)) 𝑎_𝑘>0, 𝑘≠0. It follows that the Fourier coefficients of 𝑎(𝑡)^𝑛/𝑛 are also positive for every 𝑛 due to Cauchy’s series product theorem. Since ‖𝑎‖_𝑊<𝑟₀<1, ‖𝑎^𝑛/𝑛‖_𝑊<𝑟^𝑛₀/𝑛. Hence, all the coefficients −𝑏_𝑘 in the formula (59) are positive.
Clearly, the functions 𝑓(𝑡)∈𝑊 (𝑡=𝑒^𝑖𝜙∈𝐓) are continuous on 𝐓 and, when regarded as functions 𝑓(𝑒^{−𝑖Δ𝑥𝜉}), 𝜉∈𝐑, they are (2𝜋/Δ𝑥)-periodic continuous functions. Hence, the function 𝑝(𝑡) (see (39) and (43)) can be rewritten in the form (57): 𝑃(𝜉)=exp_𝑙≠0−𝑏_−𝑙exp(𝑖𝑙𝜉Δ𝑥)−1. (60) Analogously, we rewrite 𝑃_±(𝜉) in the form (57): 𝑃₋(𝜉)=exp_𝑙<0−𝑏_−𝑙(),𝑃exp𝑖𝑙𝜉Δ𝑥−1₊(𝜉)=exp_𝑙>0−𝑏_−𝑙exp(𝑖𝑙𝜉Δ𝑥)−1. (61) It follows that 𝑃, 𝑃⁺, and 𝑃⁻ are characteristic functions of infinitely divisible lattice distributions with the maximal step Δ𝑥 supported on {𝑥=𝑘Δ𝑥∣𝑘∈𝐙}, {𝑥=𝑘Δ𝑥∣𝑘∈𝐙,𝑘≥0}, and {𝑥=𝑘Δ𝑥∣𝑘∈𝐙,𝑘≤0}, respectively.

From Theorem 6 we have that there exist discrete random variables 𝑋, 𝑋⁺, and 𝑋⁻ taking on values of the form 𝑥_𝑘=𝑘Δ𝑥, 𝑘∈𝐙, such that 𝑃𝑋=𝑥_𝑘=𝑝_−𝑘𝑃𝑋,𝑘∈𝐙,₋=𝑥_𝑘=𝑝⁺_−𝑘𝑃𝑋,𝑘∈𝐙,₊=𝑥_𝑘=𝑝⁻_−𝑘,𝑘∈𝐙, (62) where {𝑝_𝑘} and {𝑝^±_𝑘} are defined by (43) and (47)–(50).

It follows that the corresponding Laurent operators 𝐿(𝑝), 𝐿(𝑝⁺), and 𝐿(𝑝⁻) can be interpreted as expectation operators conditioned on current values of 𝑋, 𝑋⁻, and 𝑋⁺, respectively:𝐿𝑝𝑔𝑥_𝑘=_𝑙∈𝐙𝑝_𝑙𝑔𝑥_𝑘−𝑙=𝐸𝑔𝑥_𝑘,𝐿𝑝+𝑋₊𝑔𝑥_𝑘=^𝑙=+∞_𝑙=0𝑝⁺_𝑙𝑔𝑥_𝑘−𝑙=𝐸𝑔𝑥_𝑘+𝑋⁻,𝐿𝑝⁻𝑔𝑥_𝑘=^𝑙=0_𝑙=−∞𝑝⁻_𝑙𝑔𝑥_𝑘−𝑙=𝐸𝑔𝑥_𝑘+𝑋⁺. (63)

The following simple properties are immediate from the interpretation of 𝐿(𝑝_±) as expectation operators.

Proposition 7. Laurent operators 𝐿(𝑝_±) enjoy the following properties.(a)If 𝑔_𝑘=0,  ∀𝑘≥𝑘₀, then, ∀𝑘≥𝑘₀, (𝐿(𝑝₋)𝑔)_𝑘=0.(b)If 𝑔_𝑘=0,  ∀𝑘≤𝑘₀, then, ∀𝑘≤𝑘₀, (𝐿(𝑝₊)𝑔)_𝑘=0.(c)If 𝑔_𝑘≥0,  ∀𝑘, then, (𝐿(𝑝₋)𝑔)_𝑘≥0,∀𝑘. If, in addition, there exists 𝑘₀ such that 𝑔_𝑘>0  ∀𝑘>𝑘₀, then (𝐿(𝑝₋)𝑔)_𝑘>0,  ∀𝑘.(d)If 𝑔_𝑘≥0,  ∀𝑘, then (𝐿(𝑝₊)𝑔)_𝑘≥0,∀𝑥. If, in addition, there exists 𝑘₀ such that 𝑔_𝑘>0  ∀𝑘<𝑘₀, then (𝐿(𝑝₊)𝑔)_𝑘>0,  ∀𝑘.(e)If 𝑔={𝑔_𝑘} is monotone, then {(𝐿(𝑝₋)𝑔)_𝑘} and {(𝐿(𝑝₊)𝑔)_𝑘} 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 𝐺={𝐺^𝑘} in (54) is a monotonically increasing sequence and it changes the sign, then from Proposition 7 an approximate solution 𝑔=(…,𝑔(𝑥₋₂),𝑔(𝑥₋₁),𝑔(𝑥₀),𝑔(𝑥₁),𝑔(𝑥₂),…) to the problem (53) can be written in terms of Laurent operators 𝐿(𝑝^±): 𝑔=𝐿𝑝₊𝑃_𝑘0𝐿𝑝₋𝐺, (64) where the only number 𝑘₀ can be found from the following conditions:𝐿𝑝₋𝐺_𝑘>0,𝑘>𝑘₀,𝐿𝑝₋𝐺_𝑘≤0,𝑘≤𝑘₀. (65) We remark that (64) includes the requirement that the series 𝑤_𝑘=^𝑙=0_𝑙=−∞𝑝⁻_𝑙𝐺^𝑘−𝑙,𝑘∈𝐙 (66) 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 𝑟 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 4. Choose desired truncation error 𝜖 for coefficients {𝛼_𝑙} in (12) (as a rule, the choice 𝜖=10⁻⁶ is optimal). Due to the FDS, calculate 𝛼_𝑙, 𝑙=−1,−2,…,𝑙⁻, where 𝑙⁻=max{𝑙<0∣|𝛼_𝑙|<𝜖/2}; calculate 𝛼_𝑙, 𝑙=1,2,…,𝑙⁺, where 𝑙⁺=min{𝑙>0∣|𝛼_𝑙|<𝜖/2}. Set 𝛼_𝑙=0, as 𝑙<𝑙⁻ or 𝑙>𝑙⁺.

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 𝑞=(Δ𝑡)⁻¹+𝑟.

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 𝑙₀=max{−𝑙₋;𝑙₊;(𝑥_max−𝑥_min)/2Δ𝑥}. We find integer number 𝑘₀ such that 2^𝑘0−1<𝑙₀≤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.

Step 8. Find coefficients 𝑎_𝑙,  𝑙=−𝑚+1,…,𝑚, by the formula (15).

Step 9. Denote by 𝜏_𝑘=exp(𝑖𝜋𝑘/𝑚),  𝑘=−𝑚+1,…,𝑚. Find 𝑎(𝜏_𝑘)=∑^𝑙=𝑚_{𝑙=−𝑚+1}𝑎_𝑙𝜏^𝑙_𝑘,  𝑘=−𝑚+1,…,𝑚, using FFT.

Step 10. We find the symbol of 𝐿(𝑝)=(𝐿(𝑎⁻¹)): 𝑝(𝜏_𝑘)=𝑎⁻¹(𝜏_𝑘),  𝑘=−𝑚+1,…,𝑚  (see (39)).

Step 11. We find 𝑏(𝜏_𝑘)∶=ln(𝑎(𝜏_𝑘)), 𝑘=−𝑚+1,…,𝑚 (see (46)). Using inverse FFT, we obtain the sequence of coefficients 𝑏_𝑘,  𝑘=−𝑚+1,…,𝑚, for decomposition of 𝑏(𝜏) to the series 𝑏(𝜏)=^𝑙=𝑚_{𝑙=−𝑚+1}𝑏_𝑙𝜏^𝑙 (67)

Step 12. Set 𝑏⁻₀=−∑^𝑙=−1_{𝑙=−𝑚+1}𝑏_𝑙 and 𝑏₋(𝜏)=∑^𝑙=−1_{𝑙=−𝑚+1}𝑏_𝑙𝜏^𝑙+𝑏⁻₀; set 𝑏⁺₀=𝑏₀−𝑏⁻₀ and 𝑏₊(𝜏)=∑^𝑙=𝑚_𝑙=1𝑏_𝑙𝜏^𝑙+𝑏⁺₀ (see (48)). Using FFT we obtain 𝑏_±(𝜏_𝑘),  𝑘=−𝑚+1,…,𝑚.

Step 13. We find the symbols of 𝐿(𝑝^±): 𝑝^±(𝜏_𝑘)=exp(−𝑏_±(𝜏_𝑘)),  𝑘=−𝑚+1,…,𝑚 (see (49)).

4. 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)).

4.1. 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 (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 𝑉(𝑡,𝑆_𝑡) 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 𝑣(𝑡,𝑥)=𝑉(𝑡,𝐻𝑒^𝑥). Then, 𝜕_𝑡𝑣()𝑣()(𝑥)𝑣()+𝐿−𝑟𝑡,𝑥=0,0≤𝑡≤𝑇,𝑥>0,𝑇,𝑥=𝑔,𝑥>0,𝑡,𝑥=0,0≤𝑡≤𝑇,𝑥≤0. (68)(69)(70)

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 𝑣_𝑗(𝑥) the approximation to 𝑣(𝑥,𝑡_𝑗). Then 𝑣_𝑛(𝑥)=𝑔(𝑥), and by discretizing the derivative 𝜕_𝑡 in (68), we obtain, for 𝑗=𝑛−1,𝑛−2,…,0, 𝑣_𝑗+1(𝑥)−𝑣_𝑗(𝑥)−()𝑣Δ𝑡𝑟−𝐿_𝑗(𝑥)=0,𝑥>0. (71) Equation (70) assumes the form 𝑣_𝑗(𝑥)=0,𝑥≤0. (72) Set 𝑞=Δ𝑡⁻¹+𝑟; then (71) can be rewritten as follows: 𝑞⁻¹𝑞−𝐿𝑣_𝑗(𝑥)=𝑞Δ𝑡⁻¹𝑣_𝑗+1(𝑥),𝑥>0. (73)

Notice that the sequence of problems (72) and (73) has the form (11). Applying a finite difference scheme to (72) and (73) which approximates 𝑞⁻¹(𝑞−𝐿) 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 𝑉(𝑡,𝑆_𝑡) be the price of American put with the strike price 𝐾 and the terminal date 𝑇. Set 𝑥=ln(𝑆/𝐾), 𝑔(𝑥)=𝐾(1−𝑒^𝑥), and 𝑣(𝑡,𝑥)=𝑉(𝑡,𝐾𝑒^𝑥). Assume that the optimal stopping time is of the form 𝜏^′_𝐵∧𝑇, where 𝜏^′_𝐵 is the hitting time of a closed set 𝐵⊂𝐑×(−∞,𝑇] by the two-dimensional process 𝑋_𝑡=(𝑋_𝑡,𝑡). Set 𝒞=𝐑×[0,𝑇)⧵𝐵 (this is the continuation region, where the option remains alive), and consider the following boundary value problem: 𝜕_𝑡𝑣()()𝑣()(𝑥),()𝑣()≥(𝑥)+𝐿−𝑟𝑡,𝑥=0,𝑡,𝑥∈𝒞,𝑡,𝑥=𝑔𝑡,𝑥∈𝐵or𝑡=𝑇,𝑡,𝑥⁺𝜕,𝑡≤𝑇,𝑥∈𝐑,_𝑡+𝐿−𝑟𝑣(𝑡,𝑥)≤0,𝑡<𝑇,(𝑡,𝑥)∉𝒞, (74)(75)(76)(77) 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, 𝑣.

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 𝑣_𝑗(𝑥) the approximation to 𝑣(𝑥,𝑡_𝑗); ℎ_𝑗 denotes the approximation to the early exercise boundary at time 𝑡_𝑗. Then 𝑣_𝑛(𝑥)=𝐾(1−𝑒^𝑥)⁺, and by discretizing the derivative 𝜕_𝑡 in (74), we obtain, for 𝑗=𝑛−1,𝑛−2,…,0, 𝑣_𝑗+1(𝑥)−𝑣_𝑗(𝑥)−()𝑣Δ𝑡𝑟−𝐿_𝑗(𝑥)=0,𝑥>ℎ_𝑗. (78) Equation (75) assumes the form 𝑣_𝑗(𝑥)=𝑔(𝑥),𝑥≤ℎ_𝑗. (79) The approximation ℎ_𝑗 to the early exercise boundary is found so that the 𝑣_𝑗 is maximal.

Introduce 𝑣_𝑗̃(𝑥)=𝑣_𝑗(𝑥)−𝑔(𝑥) and substitute 𝑣_𝑗(𝑥)=𝑣_𝑗̃(𝑥)+𝑔(𝑥) into (78) and (79): 𝑣̃_𝑗+1(𝑥)−𝑣_𝑗̃(𝑥)−()𝑣̃Δ𝑡𝑟−𝐿_𝑗(𝑥)=𝑟,𝑥>ℎ_𝑗,𝑣_𝑗̃(𝑥)=0,𝑥≤ℎ_𝑗. (80)(81)

Set 𝑞=Δ𝑡⁻¹+𝑟 and 𝐺_𝑗=(𝑞Δ𝑡)⁻¹𝑣̃_𝑗+1−𝑞⁻¹(𝑟−𝐿)𝑔=(𝑞Δ𝑡)⁻¹𝑣̃_𝑗+1−𝑞⁻¹𝐾𝑟; then (80) can be rewritten as follows: 𝑞⁻¹𝑞−𝐿𝑣_𝑗̃(𝑥)=𝐺_𝑗(𝑥),𝑥>ℎ_𝑗. (82)

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 𝑞⁻¹(𝑞−𝐿) 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 𝑐=1.0. 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].

In Table 1, we compare the down-and-out barrier put option prices calculated by FDS&WH and CV methods for spot prices 𝑆=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.

Table 1 

Convergence of the down-and-out put prices in KoBoL model, 𝜈<1: FDS & WH versus CV.

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 𝑐=1. We choose riskless rate 𝑟=0.04879, time to expiry 𝑇=0.1 year, strike price 𝐾=100, and the barrier 𝐻=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.

Table 2 

Convergence of the down-and-out put prices in KoBoL model, 𝜈>1: FDS & WH versus CV.

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 𝑐=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.

Table 3 

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

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 𝑂(𝑀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).

Conflict of Interests

The author declares that there is no conflict of interests regarding the publication of this paper.