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

Authors: Oleg Kudryavtsev
Journal: The Scientific World Journal
 (2013)

---

## 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.


---

## Body

## 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 Hs-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 requiresO(m2) 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(mlnm) 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 Section2, 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[eiξX(t)]=e-tψ(ξ) (we confine ourselves to the one-dimensional case). The characteristic exponent is given by the Lévy-Khintchine formula:
(1)ψ(ξ)=σ22ξ2-iμξ+∫-∞+∞‍(1-eiξy+iξy1|y|≤1)F(dy),
where σ2 and μ are the variance and drift coefficient of the Gaussian component and F(dy) satisfies
(2)∫R∖{0}‍min{1,y2}F(dy)<+∞.Assume that the riskless rater is constant, and, under a risk-neutral measure chosen by the market, the underlying asset evolves as St=S0eXt, where Xt is a Lévy process. Then we must have E[eXt]<+∞, 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[eXt]=ert. Equivalently,
(3)r+ψ(-i)=0,
where r is instantaneous interest rate. The latter condition determines the drift via the other parameters of the Lévy process:
(4)μ=r-σ22+∫-∞+∞‍(1-ey+y1|y|≤1)F(dy).
Hence, the characteristic exponent may be rewritten as follows:
(5)ψ(ξ)=σ22ξ2-i(r-σ22)ξ+∫-∞+∞‍(1-eiξy-iξ(1-ey))F(dy).
Then the infinitesimal generator of X (denote it by) L is an integrodifferential operator which acts as follows:
(6)Lu(x)=σ22u′′(x)+(r-σ22)u′(x)+∫-∞+∞‍(u(x+y)-u(x)-(ey-1)u′(x))F(dy).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)ψ(ξ)=-iμξ+cΓ(-ν)[λ+ν-(λ++iξ)νhhhh+(-λ-)ν-(-λ--iξ)ν],
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
(8)F∓(dy)=ceλ±y|y|-ν-1dy.Example 2.
In DEJD model,F∓(dy) are given by exponential functions on negative and positive axis, respectively:
(9)F∓(dy)=c±(±λ±)eλ±y,
where σ>0, μ∈R, c±>0, and λ-<-1<0<λ+. Then the characteristic exponent is of the form
(10)ψ(ξ)=σ22ξ2-iμξ+ic+ξλ++iξ+ic-ξλ-+iξ.


## 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:(11)q-1(q-L)g(x)=G(x),x>0,g(x)=0,x≤0,
where q>0.Choose a space stepΔx, and set xl=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]):
(12)Lg(xk)=∑l≠0‍αlg(xk+l)-∑l≠0‍αlg(xk),
where
(13)∑l≠0‍αl<∞,αl>0,l≠0.
Then we can approximate q-1(q-L) as follows:
(14)q-1(q-L)g(xk)=∑l∈Z‍alg(xk-l),(15)al=-q-1α-l,l≠0,a0=1+q-1∑l≠0‍αl.The sequence{al}l=-∞+∞ generates doubly infinite Laurent matrix L(a) which is constant along the diagonals:
(16)L(a)=(⋯⋯⋯⋯⋯⋯⋯a0a-1a-2a-3⋯⋯a1a0a-1a-2⋯⋯a2a1a0a-1⋯⋯a3a2a1a0⋯⋯⋯⋯⋯⋯⋯).After the discretization, the functiong∈L2(R) turns into a piecewise constant function. Thus, we may consider g=(…,g(x-2),g(x-1),g(x0),g(x1),g(x2),…) as an element of l2(Z). Then we may rewrite (14) as follows:
(17)q-1(q-L)g(xk)=(L(a)g)k,k∈Z.LetT stand for the complex unit circle. Since {al} belongs to l1(Z), we may introduce the function a(t)=∑k‍aktk, 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∥W=∑k‍|ak|.Further, the sequence{ak}k=-∞+∞ is the sequence of the Fourier coefficients of a(t):
(18)ak=12π∫02π‍a(eiφ)e-ikφdφ,k∈Z.
Denote by F:L2(T)→l2(Z) the operator which maps a function a(t) to the sequence of its Fourier coefficients {ak}k=-∞+∞ (see (18)).It is well known that the Laurent matrixL(a) is the matrix representation of the multiplication by a(t) operator on L2(T) with respect to the orthonormal basis {(1/2π)eikφ}. Hence, we have
(19)L(a)=FaF-1.
It follows from (19) that
(20)L(a1)L(a2)=L(a1a2),∀a1,a2∈L∞(T).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{al}l=-∞+∞ also generates the infinite Toeplitz matrix T(a):
(21)T(a)=(a0a-1a-2a-3⋯a1a0a-1a-2⋯a2a1a0a-1⋯a3a2a1a0⋯⋯⋯⋯⋯⋯).When a finite difference approximation is applied (see (14)), (11) may be rewritten as follows:
(22)(L(a)g)k=Gk,k∈N,gk=0,k∈Z,k≤0,
where Gk=G(xk). Let P denote the orthogonal projection of l2(Z) onto l2(N):
(23)Puk={uk,k>0,0,k≤0.
Then, taking into account that T(a)=PL(a)P, we rewrite (22) in terms of Toeplitz matrices:
(24)T(a)g=PG,
where G=(…,G(x-2),G(x-1),G(x0),G(x1),G(x2),…) is considered as an element of l2(Z).The standard theory of Toeplitz matrices (see details in Section 1.5, [29]) leads us to the following theorem.Theorem 3.
Let a functiona∈W be able to be represented in the form
(25)a=exp(b),b∈W.
Then the operator T(a) is invertible and there exist a+,a-∈GW such that
(26)a+(t)=∑k≥0‍ak+tk,t∈T,(27)a-(t)=∑k≤0‍ak-tk,t∈T,(28)a=a+a-,(29)T(a)=T(a-)T(a+),(30)T(a)-1=T(a+-1)T(a--1).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{al}l=-∞+∞ 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(31)a~(t)=1-a(t)a0,t∈T,
then we have
(32)a~(t)=∑l≠0‍a~ltl,a~l=-ala0.
According to (15) and (32), there exists a positive number r0<1 such that
(33)∥a~(t)∥W<r0,∀t∈T.
Hence, the Taylor series for ln(1-a~(t)) (ln(·) is the principal branch of the logarithm),
(34)-∑n>0‍a~(t)nn
converges at every point t∈T.
Set(35)b(t)=ln(1-a~(t))+ln(a0),(36)ϕk(t)=-∑n=1k‍a~(t)nn+ln(a0),k∈N.
Becausea~(t)∈W and the condition (33) is satisfied, we conclude that ϕk(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
(37)∥ϕk+m-ϕk∥W≤∑n=k+1k+m‍r0nn≤r0k+1(k+1)(1-r0)<ϵ,∀m∈N.
Thus, we have proved that the functionb(t) belongs to the Wiener algebra W. Obviously, a(t)=exp(b(t)), ∀t∈T.Let a functiona(t) satisfy the conditions of Proposition 4. Set
(38)p(t)=(a(t))-1,p+(t)=(a+(t))-1,p-(t)=(a-(t))-1.
From Theorem 3 and Proposition 4, we deduce
(39)p(t)=∑l∈Z‍pltl,p∈W,(40)p-(t)=∑l=-∞l=0‍pl-tl,p-∈W,(41)p+(t)=∑l=0l=+∞‍pl+tl,p+∈W.Further, we will describe an algorithm for finding coefficients{pl},  {pl±}, based on the theory from [29]. By Theorem 3,
(42)L(p)=L(a)-1,T(p+)=T(a+)-1,T(p-)=T(a-)-1.
It follows that {pl} is the sequence of the Fourier coefficients of a-1. We have
(43)pk=12π∫-ππ‍a(eiφ)-1e-ikφdφ,k∈Z,
due to (18).Wiener-Hopf factorization formula (28) gives
(44)p(t)=p+(t)p-(t),∀t∈T.
The factors p± can be found as follows. From Proposition 4, there exists a function
(45)b(t)=∑k=-∞+∞‍bktk,∑k=-∞+∞‍|bk|<∞,
such that
(46)b(t)=lna(t),t∈T,
where the sequence of the Fourier coefficients of this function is defined as
(47)bk=12π∫-ππ‍lna(eiφ)e-ikφdφ,k∈Z.Notice thatp(t)=e-b(t),  t∈T. Next we define b± as
(48)b-(t)=∑k=-∞k=-1‍bk(tk-1),b+(t)=∑k=1k=+∞‍bk(tk-1),t∈T,
Further, we set
(49)p+(t)=e-b+(t),p-(t)=e-b-(t),t∈T.
Finally, we have
(50)pk±=12π∫-ππ‍p±(eiφ)e-ikφdφ,k∈Z.
Obviously, pk+=0 as k<0, and pk-=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
(51)g=T(a+-1)T(a--1)G.From a practical point of view, it is more convenient to rewrite (51) in terms of Laurent operators L(p±):
(52)g=L(p+)PL(p-)G.
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(MlnM), 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.LetG(x) be a monotonically increasing function, and it changes sign on the real line. Consider the following problem:(53)q-1(q-L)g(x)=G(x),x>h,g(x)=0,x≤h,
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:
(54)(L(a)g)k=Gk,k>k0,gk=0,k≤k0,
where gk=g(xk), Gk=(qΔt)-1G(xk), k0 maximizes g, and a is the symbol of a Laurent L(a) (see (12)–(16)).Introduce the orthogonal projectionPl as follows:
(55)Pluk={uk,k>l,0,k≤l.
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{pl}, {pl±} be as defined by (43) and (47)–(50), and assume that the conditions of Proposition 4 are satisfied. Set
(56)P(ξ)=∑l∈Z‍plexp(-ilξΔx),ξ∈R,P+(ξ)=∑l=-∞l=0‍pl-exp(-ilξΔx),ξ∈R,P-(ξ)=∑l=0l=+∞‍pl+exp(-ilξΔx),ξ∈R.
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.Proof.
Recall that the characteristic function of an infinitely divisible lattice distribution with the maximal stepΔx has the following form (see, e.g., [30]):
(57)c(ξ)=exp(∑k∈Z‍ck(exp(ikξΔx)-1)),ξ∈R,
where ck≥0,∀k∈Z,  ∑k∈Z‍ck<+∞.
From (35), we have
(58)-b(t)=∑n=1+∞‍a~(t)nn-ln(a0).
Since b∈W and b(1)=0, the function -b can be written as follows:
(59)-b(t)=∑k=-∞+∞(-bk)(tk-1),
where bk are defined by (47).
Notice that by the definition (see (32)) a~k>0, k≠0. It follows that the Fourier coefficients of a~(t)n/n are also positive for every n due to Cauchy’s series product theorem. Since ∥a~∥W<r0<1, ∥a~n/n∥W<r0n/n. Hence, all the coefficients -bk in the formula (59) are positive.
Clearly, the functionsf(t)∈W (t=eiϕ∈T) are continuous on T and, when regarded as functions f(e-iΔxξ), ξ∈R, they are (2π/Δx)-periodic continuous functions. Hence, the function p(t) (see (39) and (43)) can be rewritten in the form (57):
(60)P(ξ)=exp(∑l≠0‍(-b-l)(exp(ilξΔx)-1)).
Analogously, we rewrite P±(ξ) in the form (57):
(61)P-(ξ)=exp(∑l<0‍(-b-l)(exp(ilξΔx)-1)),P+(ξ)=exp(∑l>0‍(-b-l)(exp(ilξΔx)-1)).
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 Theorem6 we have that there exist discrete random variables X, X+, and X- taking on values of the form xk=kΔx, k∈Z, such that
(62)P(X=xk)=p-k,k∈Z,P(X-=xk)=p-k+,k∈Z,P(X+=xk)=p-k-,k∈Z,
where {pk} and {pk±} are defined by (43) and (47)–(50).It follows that the corresponding Laurent operatorsL(p), L(p+), and L(p-) can be interpreted as expectation operators conditioned on current values of X, X-, and X+, respectively:(63)L(p)g(xk)=∑l∈Z‍plg(xk-l)=E[g(xk+X)],L(p+)g(xk)=∑l=0l=+∞‍pl+g(xk-l)=E[g(xk+X-)],L(p-)g(xk)=∑l=-∞l=0‍pl-g(xk-l)=E[g(xk+X+)].The following simple properties are immediate from the interpretation ofL(p±) as expectation operators.Proposition 7.
Laurent operatorsL(p±) enjoy the following properties.(a)
Ifgk=0,  ∀k≥k0, then, ∀k≥k0, (L(p-)g)k=0.(b)
Ifgk=0,  ∀k≤k0, then, ∀k≤k0, (L(p+)g)k=0.(c)
Ifgk≥0,  ∀k, then, (L(p-)g)k≥0,∀k. If, in addition, there exists k0 such that gk>0  ∀k>k0, then (L(p-)g)k>0,  ∀k.(d)
Ifgk≥0,  ∀k, then (L(p+)g)k≥0,∀x. If, in addition, there exists k0 such that gk>0  ∀k<k0, then (L(p+)g)k>0,  ∀k.(e)
Ifg={gk} is monotone, then {(L(p-)g)k} and {(L(p+)g)k} are also monotone.Proposition7 is a direct analog of the properties of the expected present value operators introduced in [14]; see Proposition 6.2.1.Taking into account thatG={Gk} 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(x0),g(x1),g(x2),…) to the problem (53) can be written in terms of Laurent operators L(p±):
(64)g=L(p+)Pk0L(p-)G,
where the only number k0 can be found from the following conditions:(65)(L(p-)G)k>0,k>k0,(L(p-)G)k≤0,k≤k0.
We remark that (64) includes the requirement that the series
(66)wk=∑l=-∞l=0‍pl-Gk-l,k∈Z
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 FactorizationStep 1.
Input the interest rater 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 generatorL.Step 4.
Choose desired truncation errorϵ for coefficients {αl} in (12) (as a rule, the choice ϵ=10-6 is optimal). Due to the FDS, calculate αl, l=-1,-2,…,l-, where l-=max{l<0∣|αl|<ϵ/2}; calculate αl, l=1,2,…,l+, where l+=min{l>0∣|αl|<ϵ/2}. Set αl=0, as l<l- or l>l+.Step 5.
Input the terminal dateT 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.
Inputxmin and xmax, the lower and upper bounds for the space variable x. As a rule, the choice xmin=ln(0.4) and xmax=ln(2.5) is optimal.Step 7.
Define the number of space pointsm as follows. Set l0=max{-l-;l+;(xmax-xmin)/2Δx}. We find integer number k0 such that 2k0-1<l0≤2k0 and set m=2k0. 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 coefficientsal,  l=-m+1,…,m, by the formula (15).Step 9.
Denote byτk=exp(iπk/m),  k=-m+1,…,m. Find a(τk)=∑l=-m+1l=m‍alτkl,  k=-m+1,…,m, using FFT.Step 10.
We find the symbol ofL(p)=(L(a-1)): p(τk)=a-1(τk),  k=-m+1,…,m  (see (39)).Step 11.
We findb(τk):=ln(a(τk)), k=-m+1,…,m (see (46)). Using inverse FFT, we obtain the sequence of coefficients bk,  k=-m+1,…,m, for decomposition of b(τ) to the series
(67)b(τ)=∑l=-m+1l=m‍blτlStep 12.
Setb0-=-∑l=-m+1l=-1‍bl and b-(τ)=∑l=-m+1l=-1‍blτl+b0-; set b0+=b0-b0- and b+(τ)=∑l=1l=m‍blτl+b0+ (see (48)). Using FFT we obtain b±(τk),  k=-m+1,…,m.Step 13.
We find the symbols ofL(p±): p±(τk)=exp(-b±(τk)),  k=-m+1,…,m (see (49)).


## 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:(11)q-1(q-L)g(x)=G(x),x>0,g(x)=0,x≤0,
where q>0.Choose a space stepΔx, and set xl=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]):
(12)Lg(xk)=∑l≠0‍αlg(xk+l)-∑l≠0‍αlg(xk),
where
(13)∑l≠0‍αl<∞,αl>0,l≠0.
Then we can approximate q-1(q-L) as follows:
(14)q-1(q-L)g(xk)=∑l∈Z‍alg(xk-l),(15)al=-q-1α-l,l≠0,a0=1+q-1∑l≠0‍αl.The sequence{al}l=-∞+∞ generates doubly infinite Laurent matrix L(a) which is constant along the diagonals:
(16)L(a)=(⋯⋯⋯⋯⋯⋯⋯a0a-1a-2a-3⋯⋯a1a0a-1a-2⋯⋯a2a1a0a-1⋯⋯a3a2a1a0⋯⋯⋯⋯⋯⋯⋯).After the discretization, the functiong∈L2(R) turns into a piecewise constant function. Thus, we may consider g=(…,g(x-2),g(x-1),g(x0),g(x1),g(x2),…) as an element of l2(Z). Then we may rewrite (14) as follows:
(17)q-1(q-L)g(xk)=(L(a)g)k,k∈Z.LetT stand for the complex unit circle. Since {al} belongs to l1(Z), we may introduce the function a(t)=∑k‍aktk, 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∥W=∑k‍|ak|.Further, the sequence{ak}k=-∞+∞ is the sequence of the Fourier coefficients of a(t):
(18)ak=12π∫02π‍a(eiφ)e-ikφdφ,k∈Z.
Denote by F:L2(T)→l2(Z) the operator which maps a function a(t) to the sequence of its Fourier coefficients {ak}k=-∞+∞ (see (18)).It is well known that the Laurent matrixL(a) is the matrix representation of the multiplication by a(t) operator on L2(T) with respect to the orthonormal basis {(1/2π)eikφ}. Hence, we have
(19)L(a)=FaF-1.
It follows from (19) that
(20)L(a1)L(a2)=L(a1a2),∀a1,a2∈L∞(T).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{al}l=-∞+∞ also generates the infinite Toeplitz matrix T(a):
(21)T(a)=(a0a-1a-2a-3⋯a1a0a-1a-2⋯a2a1a0a-1⋯a3a2a1a0⋯⋯⋯⋯⋯⋯).When a finite difference approximation is applied (see (14)), (11) may be rewritten as follows:
(22)(L(a)g)k=Gk,k∈N,gk=0,k∈Z,k≤0,
where Gk=G(xk). Let P denote the orthogonal projection of l2(Z) onto l2(N):
(23)Puk={uk,k>0,0,k≤0.
Then, taking into account that T(a)=PL(a)P, we rewrite (22) in terms of Toeplitz matrices:
(24)T(a)g=PG,
where G=(…,G(x-2),G(x-1),G(x0),G(x1),G(x2),…) is considered as an element of l2(Z).The standard theory of Toeplitz matrices (see details in Section 1.5, [29]) leads us to the following theorem.Theorem 3.
Let a functiona∈W be able to be represented in the form
(25)a=exp(b),b∈W.
Then the operator T(a) is invertible and there exist a+,a-∈GW such that
(26)a+(t)=∑k≥0‍ak+tk,t∈T,(27)a-(t)=∑k≤0‍ak-tk,t∈T,(28)a=a+a-,(29)T(a)=T(a-)T(a+),(30)T(a)-1=T(a+-1)T(a--1).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{al}l=-∞+∞ 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(31)a~(t)=1-a(t)a0,t∈T,
then we have
(32)a~(t)=∑l≠0‍a~ltl,a~l=-ala0.
According to (15) and (32), there exists a positive number r0<1 such that
(33)∥a~(t)∥W<r0,∀t∈T.
Hence, the Taylor series for ln(1-a~(t)) (ln(·) is the principal branch of the logarithm),
(34)-∑n>0‍a~(t)nn
converges at every point t∈T.
Set(35)b(t)=ln(1-a~(t))+ln(a0),(36)ϕk(t)=-∑n=1k‍a~(t)nn+ln(a0),k∈N.
Becausea~(t)∈W and the condition (33) is satisfied, we conclude that ϕk(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
(37)∥ϕk+m-ϕk∥W≤∑n=k+1k+m‍r0nn≤r0k+1(k+1)(1-r0)<ϵ,∀m∈N.
Thus, we have proved that the functionb(t) belongs to the Wiener algebra W. Obviously, a(t)=exp(b(t)), ∀t∈T.Let a functiona(t) satisfy the conditions of Proposition 4. Set
(38)p(t)=(a(t))-1,p+(t)=(a+(t))-1,p-(t)=(a-(t))-1.
From Theorem 3 and Proposition 4, we deduce
(39)p(t)=∑l∈Z‍pltl,p∈W,(40)p-(t)=∑l=-∞l=0‍pl-tl,p-∈W,(41)p+(t)=∑l=0l=+∞‍pl+tl,p+∈W.Further, we will describe an algorithm for finding coefficients{pl},  {pl±}, based on the theory from [29]. By Theorem 3,
(42)L(p)=L(a)-1,T(p+)=T(a+)-1,T(p-)=T(a-)-1.
It follows that {pl} is the sequence of the Fourier coefficients of a-1. We have
(43)pk=12π∫-ππ‍a(eiφ)-1e-ikφdφ,k∈Z,
due to (18).Wiener-Hopf factorization formula (28) gives
(44)p(t)=p+(t)p-(t),∀t∈T.
The factors p± can be found as follows. From Proposition 4, there exists a function
(45)b(t)=∑k=-∞+∞‍bktk,∑k=-∞+∞‍|bk|<∞,
such that
(46)b(t)=lna(t),t∈T,
where the sequence of the Fourier coefficients of this function is defined as
(47)bk=12π∫-ππ‍lna(eiφ)e-ikφdφ,k∈Z.Notice thatp(t)=e-b(t),  t∈T. Next we define b± as
(48)b-(t)=∑k=-∞k=-1‍bk(tk-1),b+(t)=∑k=1k=+∞‍bk(tk-1),t∈T,
Further, we set
(49)p+(t)=e-b+(t),p-(t)=e-b-(t),t∈T.
Finally, we have
(50)pk±=12π∫-ππ‍p±(eiφ)e-ikφdφ,k∈Z.
Obviously, pk+=0 as k<0, and pk-=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
(51)g=T(a+-1)T(a--1)G.From a practical point of view, it is more convenient to rewrite (51) in terms of Laurent operators L(p±):
(52)g=L(p+)PL(p-)G.
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(MlnM), 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.LetG(x) be a monotonically increasing function, and it changes sign on the real line. Consider the following problem:(53)q-1(q-L)g(x)=G(x),x>h,g(x)=0,x≤h,
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:
(54)(L(a)g)k=Gk,k>k0,gk=0,k≤k0,
where gk=g(xk), Gk=(qΔt)-1G(xk), k0 maximizes g, and a is the symbol of a Laurent L(a) (see (12)–(16)).Introduce the orthogonal projectionPl as follows:
(55)Pluk={uk,k>l,0,k≤l.
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{pl}, {pl±} be as defined by (43) and (47)–(50), and assume that the conditions of Proposition 4 are satisfied. Set
(56)P(ξ)=∑l∈Z‍plexp(-ilξΔx),ξ∈R,P+(ξ)=∑l=-∞l=0‍pl-exp(-ilξΔx),ξ∈R,P-(ξ)=∑l=0l=+∞‍pl+exp(-ilξΔx),ξ∈R.
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.Proof.
Recall that the characteristic function of an infinitely divisible lattice distribution with the maximal stepΔx has the following form (see, e.g., [30]):
(57)c(ξ)=exp(∑k∈Z‍ck(exp(ikξΔx)-1)),ξ∈R,
where ck≥0,∀k∈Z,  ∑k∈Z‍ck<+∞.
From (35), we have
(58)-b(t)=∑n=1+∞‍a~(t)nn-ln(a0).
Since b∈W and b(1)=0, the function -b can be written as follows:
(59)-b(t)=∑k=-∞+∞(-bk)(tk-1),
where bk are defined by (47).
Notice that by the definition (see (32)) a~k>0, k≠0. It follows that the Fourier coefficients of a~(t)n/n are also positive for every n due to Cauchy’s series product theorem. Since ∥a~∥W<r0<1, ∥a~n/n∥W<r0n/n. Hence, all the coefficients -bk in the formula (59) are positive.
Clearly, the functionsf(t)∈W (t=eiϕ∈T) are continuous on T and, when regarded as functions f(e-iΔxξ), ξ∈R, they are (2π/Δx)-periodic continuous functions. Hence, the function p(t) (see (39) and (43)) can be rewritten in the form (57):
(60)P(ξ)=exp(∑l≠0‍(-b-l)(exp(ilξΔx)-1)).
Analogously, we rewrite P±(ξ) in the form (57):
(61)P-(ξ)=exp(∑l<0‍(-b-l)(exp(ilξΔx)-1)),P+(ξ)=exp(∑l>0‍(-b-l)(exp(ilξΔx)-1)).
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 Theorem6 we have that there exist discrete random variables X, X+, and X- taking on values of the form xk=kΔx, k∈Z, such that
(62)P(X=xk)=p-k,k∈Z,P(X-=xk)=p-k+,k∈Z,P(X+=xk)=p-k-,k∈Z,
where {pk} and {pk±} are defined by (43) and (47)–(50).It follows that the corresponding Laurent operatorsL(p), L(p+), and L(p-) can be interpreted as expectation operators conditioned on current values of X, X-, and X+, respectively:(63)L(p)g(xk)=∑l∈Z‍plg(xk-l)=E[g(xk+X)],L(p+)g(xk)=∑l=0l=+∞‍pl+g(xk-l)=E[g(xk+X-)],L(p-)g(xk)=∑l=-∞l=0‍pl-g(xk-l)=E[g(xk+X+)].The following simple properties are immediate from the interpretation ofL(p±) as expectation operators.Proposition 7.
Laurent operatorsL(p±) enjoy the following properties.(a)
Ifgk=0,  ∀k≥k0, then, ∀k≥k0, (L(p-)g)k=0.(b)
Ifgk=0,  ∀k≤k0, then, ∀k≤k0, (L(p+)g)k=0.(c)
Ifgk≥0,  ∀k, then, (L(p-)g)k≥0,∀k. If, in addition, there exists k0 such that gk>0  ∀k>k0, then (L(p-)g)k>0,  ∀k.(d)
Ifgk≥0,  ∀k, then (L(p+)g)k≥0,∀x. If, in addition, there exists k0 such that gk>0  ∀k<k0, then (L(p+)g)k>0,  ∀k.(e)
Ifg={gk} is monotone, then {(L(p-)g)k} and {(L(p+)g)k} are also monotone.Proposition7 is a direct analog of the properties of the expected present value operators introduced in [14]; see Proposition 6.2.1.Taking into account thatG={Gk} 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(x0),g(x1),g(x2),…) to the problem (53) can be written in terms of Laurent operators L(p±):
(64)g=L(p+)Pk0L(p-)G,
where the only number k0 can be found from the following conditions:(65)(L(p-)G)k>0,k>k0,(L(p-)G)k≤0,k≤k0.
We remark that (64) includes the requirement that the series
(66)wk=∑l=-∞l=0‍pl-Gk-l,k∈Z
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 FactorizationStep 1.
Input the interest rater 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 generatorL.Step 4.
Choose desired truncation errorϵ for coefficients {αl} in (12) (as a rule, the choice ϵ=10-6 is optimal). Due to the FDS, calculate αl, l=-1,-2,…,l-, where l-=max{l<0∣|αl|<ϵ/2}; calculate αl, l=1,2,…,l+, where l+=min{l>0∣|αl|<ϵ/2}. Set αl=0, as l<l- or l>l+.Step 5.
Input the terminal dateT 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.
Inputxmin and xmax, the lower and upper bounds for the space variable x. As a rule, the choice xmin=ln(0.4) and xmax=ln(2.5) is optimal.Step 7.
Define the number of space pointsm as follows. Set l0=max{-l-;l+;(xmax-xmin)/2Δx}. We find integer number k0 such that 2k0-1<l0≤2k0 and set m=2k0. 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 coefficientsal,  l=-m+1,…,m, by the formula (15).Step 9.
Denote byτk=exp(iπk/m),  k=-m+1,…,m. Find a(τk)=∑l=-m+1l=m‍alτkl,  k=-m+1,…,m, using FFT.Step 10.
We find the symbol ofL(p)=(L(a-1)): p(τk)=a-1(τk),  k=-m+1,…,m  (see (39)).Step 11.
We findb(τk):=ln(a(τk)), k=-m+1,…,m (see (46)). Using inverse FFT, we obtain the sequence of coefficients bk,  k=-m+1,…,m, for decomposition of b(τ) to the series
(67)b(τ)=∑l=-m+1l=m‍blτlStep 12.
Setb0-=-∑l=-m+1l=-1‍bl and b-(τ)=∑l=-m+1l=-1‍blτl+b0-; set b0+=b0-b0- and b+(τ)=∑l=1l=m‍blτl+b0+ (see (48)). Using FFT we obtain b±(τk),  k=-m+1,…,m.Step 13.
We find the symbols ofL(p±): p±(τk)=exp(-b±(τk)),  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 rater>0 is constant and, under a risk-neutral measure chosen by the market, the log price of the stock Xt=logSt 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 amountG(ST) 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 somerebate. 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,St) 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(Hex), and v(t,x)=V(t,Hex). Then,
(68)(∂t+L-r)v(t,x)=0,0≤t≤T,x>0,(69)v(T,x)=g(x),x>0,(70)v(t,x)=0,0≤t≤T,x≤0.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 tj=jΔt,j=0,1,…,n, where Δt=T/n, and denote by vj(x) the approximation to v(x,tj). Then vn(x)=g(x), and by discretizing the derivative ∂t in (68), we obtain, for j=n-1,n-2,…,0,
(71)vj+1(x)-vj(x)Δt-(r-L)vj(x)=0,x>0.
Equation (70) assumes the form
(72)vj(x)=0,x≤0.
Set q=Δt-1+r; then (71) can be rewritten as follows:
(73)q-1(q-L)vj(x)=(qΔt)-1vj+1(x),x>0.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).LetV(t,St) 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-ex), and v(t,x)=V(t,Kex). Assume that the optimal stopping time is of the form τB′∧T, where τB′ is the hitting time of a closed set B⊂R×(-∞,T] by the two-dimensional process X^t=(Xt,t). Set 𝒞=R×[0,T)∖B (this is the continuation region, where the option remains alive), and consider the following boundary value problem:
(74)(∂t+L-r)v(t,x)=0,(t,x)∈𝒞,(75)v(t,x)=g(x),(t,x)∈Bort=T,(76)v(t,x)≥(x)+,t≤T,x∈R,(77)(∂t+L-r)v(t,x)≤0,t<T,(t,x)∉𝒞¯,
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 Section6.2.2 of [1] for the American put. Normalize the strike price to 1; divide [0,T] into n subperiods by points tj=jΔt,j=0,1,…,n, where Δt=T/n; and denote by vj(x) the approximation to v(x,tj); hj denotes the approximation to the early exercise boundary at time tj. Then vn(x)=K(1-ex)+, and by discretizing the derivative ∂t in (74), we obtain, for j=n-1,n-2,…,0,
(78)vj+1(x)-vj(x)Δt-(r-L)vj(x)=0,x>hj.
Equation (75) assumes the form
(79)vj(x)=g(x),x≤hj.
The approximation hj to the early exercise boundary is found so that the vj is maximal.Introducev~j(x)=vj(x)-g(x) and substitute vj(x)=v~j(x)+g(x) into (78) and (79):
(80)v~j+1(x)-v~j(x)Δt-(r-L)v~j(x)=r,x>hj,(81)v~j(x)=0,x≤hj.Setq=Δt-1+r and Gj=(qΔt)-1v~j+1-q-1(r-L)g=(qΔt)-1v~j+1-q-1Kr; then (80) can be rewritten as follows:
(82)q-1(q-L)v~j(x)=Gj(x),x>hj.Notice that the sequence of problems (81) and (82) has the form (53), where Gj(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.


## 4.1. Barrier Options

Consider a contract which pays the specified amountG(ST) 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 somerebate. 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,St) 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(Hex), and v(t,x)=V(t,Hex). Then,
(68)(∂t+L-r)v(t,x)=0,0≤t≤T,x>0,(69)v(T,x)=g(x),x>0,(70)v(t,x)=0,0≤t≤T,x≤0.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 tj=jΔt,j=0,1,…,n, where Δt=T/n, and denote by vj(x) the approximation to v(x,tj). Then vn(x)=g(x), and by discretizing the derivative ∂t in (68), we obtain, for j=n-1,n-2,…,0,
(71)vj+1(x)-vj(x)Δt-(r-L)vj(x)=0,x>0.
Equation (70) assumes the form
(72)vj(x)=0,x≤0.
Set q=Δt-1+r; then (71) can be rewritten as follows:
(73)q-1(q-L)vj(x)=(qΔt)-1vj+1(x),x>0.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).LetV(t,St) 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-ex), and v(t,x)=V(t,Kex). Assume that the optimal stopping time is of the form τB′∧T, where τB′ is the hitting time of a closed set B⊂R×(-∞,T] by the two-dimensional process X^t=(Xt,t). Set 𝒞=R×[0,T)∖B (this is the continuation region, where the option remains alive), and consider the following boundary value problem:
(74)(∂t+L-r)v(t,x)=0,(t,x)∈𝒞,(75)v(t,x)=g(x),(t,x)∈Bort=T,(76)v(t,x)≥(x)+,t≤T,x∈R,(77)(∂t+L-r)v(t,x)≤0,t<T,(t,x)∉𝒞¯,
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 Section6.2.2 of [1] for the American put. Normalize the strike price to 1; divide [0,T] into n subperiods by points tj=jΔt,j=0,1,…,n, where Δt=T/n; and denote by vj(x) the approximation to v(x,tj); hj denotes the approximation to the early exercise boundary at time tj. Then vn(x)=K(1-ex)+, and by discretizing the derivative ∂t in (74), we obtain, for j=n-1,n-2,…,0,
(78)vj+1(x)-vj(x)Δt-(r-L)vj(x)=0,x>hj.
Equation (75) assumes the form
(79)vj(x)=g(x),x≤hj.
The approximation hj to the early exercise boundary is found so that the vj is maximal.Introducev~j(x)=vj(x)-g(x) and substitute vj(x)=v~j(x)+g(x) into (78) and (79):
(80)v~j+1(x)-v~j(x)Δt-(r-L)v~j(x)=r,x>hj,(81)v~j(x)=0,x≤hj.Setq=Δt-1+r and Gj=(qΔt)-1v~j+1-q-1(r-L)g=(qΔt)-1v~j+1-q-1Kr; then (80) can be rewritten as follows:
(82)q-1(q-L)v~j(x)=Gj(x),x>hj.Notice that the sequence of problems (81) and (82) has the form (53), where Gj(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 asCV 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 Table1, 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.Table 1
Convergence of the down-and-out put prices in KoBoL model,ν<1: FDS & WH versus CV.
(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.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 Table2, 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.
(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.


### 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 stepsNτ
ϵ
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.
ϵ
p and ϵb are the maximums of the relative differences between correspondent prices and boundaries, respectively, in the region S≤1.3K.


## 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 asCV 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 Table1, 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.Table 1
Convergence of the down-and-out put prices in KoBoL model,ν<1: FDS & WH versus CV.
(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.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 Table2, 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.
(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.


## 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 stepsNτ
ϵ
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.
ϵ
p and ϵb 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(MlnM) 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).



---

*Source: 963625-2013-12-30.xml*