Efficient solution of a partial integro-differential equation in finance

Applied Numerical Mathematics 58 (2008) 1687–1703 www.elsevier.com/locate/apnum

Efficient solution of a partial integro-differential equation in finance E.W. Sachs a,b,∗ , A.K. Strauss c a Universität Trier, Fachbereich IV, Abteilung Mathematik, 54286 Trier, Germany b Interdisciplinary Center for Applied Mathematics, Virginia Tech, Blacksburg, VA 24061, USA c Department of Management Science, Lancaster University, Lancaster, LA1 4YX, UK

Available online 23 November 2007

Abstract Jump-diffusion models for the pricing of derivatives lead under certain assumptions to partial integro-differential equations (PIDE). Such a PIDE typically involves a convection term and a non-local integral. We transform the PIDE to eliminate the convection term, discretize it implicitly, and use finite differences on a uniform grid. The resulting dense linear system exhibits so much structure that it can be solved very efficiently by a circulant preconditioned conjugate gradient method. Therefore, this fully implicit scheme requires only on the order of O(n log n) operations. Second order accuracy is obtained numerically on the whole computational domain for Merton’s model. Published by Elsevier B.V. on behalf of IMACS. Keywords: Lévy process; Partial integro-differential equations; Conjugate Gradient method; Toeplitz matrices

1. Introduction Black and Scholes [4] developed a formula to compute the price of a European call which is used in many practical implementations since 1975 in the finance industry. Extensions of this model include among others stochastic volatility models [23,20,3], Lévy models [12] and jump-diffusion models [29,2,25]. In jump-diffusion models the log-price of the underlying commodity follows a diffusion process, similar to the Black–Scholes model, but jumps can occur at random time points which are modeled by a compound Poisson process. Under certain assumptions the jump-diffusion models for the pricing of European call options leads to a partial integro-differential equation (PIDE) involving a non-local integral term. In order to solve the PIDE problem numerically, Andersen and Andreasen [2] use an unconditionally stable ADI finite difference method and accelerate it using the fast Fourier transform (FFT). In [26,28] and [27] wavelet methods are applied to infinite activity jump-diffusion models. Several researchers used implicit-explicit finite difference methods, see for example the recent papers [5,6] and [14]. They avoid the solve of dense linear systems that arise by an implicit integral discretization, and use implicit schemes only for the second order derivative term in the PIDE. Briani [5] also approximates the convection term explicitly, whereas Cont and Voltchkova [14] treat the whole differential part implicitly. However, Cont and Voltchkova [14] use backward or forward difference quotients of only first order depending on the sign of the coefficient of the convection term in order to avoid oscillations. Implicit finite difference approaches in combination with * Corresponding author.

E-mail addresses: [email protected], [email protected] (E.W. Sachs), [email protected] (A.K. Strauss). 0168-9274/$30.00 Published by Elsevier B.V. on behalf of IMACS. doi:10.1016/j.apnum.2007.11.002

1688

E.W. Sachs, A.K. Strauss / Applied Numerical Mathematics 58 (2008) 1687–1703

FFT have been investigated in [16] extended by a penalty method to price American options and in [17] where the correlation integral is evaluated using an FFT method. Almendral and Oosterlee [1] present both an implicit discretization of the PIDE on a uniform grid using finite differences as well as a finite element approach, where a splitting technique combined with FFT is used to accelerate the dense matrix-vector product. Ikonen and Toivanen [24] use a finite difference approach on a non-uniform grid refined around the strike price using the Rannacher time-stepping scheme and a splitting technique like in [1]. They obtain fast convergence but cannot use FFT due to the non-uniform grid. In cases where the volatility is driven by a stochastic process, small values of the volatility can occur for certain time instances, giving more weight for the convection term and thus causing instabilities in the numerics. Therefore, we transform the PIDE analytically in such a way that the convection term disappears completely. The resulting discretized linear system still has a dense coefficient matrix due to the discretized integral term, but since the convection terms are missing, it exhibits a lot of structure: It is a symmetric and positive definite coefficient matrix, which is, in addition, also Toeplitz if we use a uniform grid. This allows to utilize a fully implicit scheme, because the application of a circulant preconditioned conjugate gradient method for the solution of the systems requires an effort of n log n. The numerical results show that the resulting method is fast and of second order. In particular, for Merton’s model second order accuracy is obtained on the whole computational domain improving the results of the numerical experiments in [1]. The paper is organized as follows: First, we introduce the PIDE in a general form and then specifically for Merton’s model. We present the elimination of the convection term in this specific PIDE by an appropriate transformation and discretize the resulting PIDE in the following third section. Then we briefly review some facts about circulant matrices that play an important role in speeding up the numerical solution process. After an analysis of the eigenvalues of the coefficient matrix in Section 5 we improve the spectrum by means of circulant preconditioning and employ a preconditioned conjugate gradient method in Section 7 to solve the arising linear systems. In Section 8 we present the numerical results and a conclusion. 2. Partial integro-differential equation This section contains a short presentation of the stochastic process and the corresponding partial integro-differential equation (PIDE) which describes the option price for general Lévy processes. In the framework of this paper we confine ourselves to the model by Merton [29]. In the first subsection we restate a rather general form of a PIDE formulation, in the second we derive from it a more specific version for the considered model and, finally, in the third subsection we derive an equivalent formulation without convection term. For the purpose of a clear distinction between the two equivalent PIDE versions in Sections 2.2 and 2.3 we call them PIDE W (With convection term) and PIDE N (No convection term), respectively. 2.1. General PIDE for Lévy processes The option price depends on the time to maturity and the price of the underlying St as variables and on some given parameters. We describe the movement of the price of the underlying asset, for example a stock, for which we use a stochastic process (St )t0 on a probability space (Ω, F , Ft , P) with filtration F . We define the price process of the underlying asset by an exponential Lévy model, that is St = S0 exp(rt + Xt ), where r is the riskfree interest rate, S0 the stock price at time zero and Xt a Lévy process given by   Nt  σ2 t + σ Wt + Xt := μ − r − Yi , 2 i=1  where we set μ + λη = r, η := E(ex − 1) = R (ex − 1) dF (x) and F is the distribution of the independent identically distributed random variables (Yi ). Furthermore, (Nt ) is a Poisson process with intensity λ, i.e. the mean number of jumps per unit time, and (Wt ), (Yi ), (Nt ) are mutually independent. As usual, Wt denotes the Brownian motion, σ the volatility and μ the drift. After plugging in μ = r − λη, the price process becomes St = S0 exp(Lt ), where  t Lt := (r − λη − σ 2 /2)t + σ Wt + N 1 Yi .

E.W. Sachs, A.K. Strauss / Applied Numerical Mathematics 58 (2008) 1687–1703

1689

It is shown in the stochastic finance literature that, after the change of variable x := ln S, the definition of a price for an European call option is the expected value   v(t, x) := e−r(T −t) EQ (ex+LT −t − K)+ , where EQ (·) is the expectation operator under the equivalent martingale measure Q (which is assumed to be already chosen) and leads to the following partial differential equation including an integral term. Theorem 2.1. (See [32].) Let v(t, x) ∈ C 1,2 ([0, T ) × R) ∩ C 0 ([0, T ] × R) and assume satisfies the following partial integro-differential equation a vt + vxx + γ vx − rv + 2

∞

  v(t, x + y) − v(t, x) − vx (t, x)y dν(y) = 0,



R |x| dν(x)

< ∞. Then v

∀(t, x) ∈ [0, T ) × R,

−∞

with final condition v(T , x) = H (ex ),

∀x ∈ R.

The only parameters entering are the risk-free interest rate r and the reduced Lévy triplet (a, γ , ν) under the riskneutral measure Q. The function H (·) denotes the payoff at t = T . 2.2. PIDE W with convection term In this short section we reformulate the PIDE using Merton’s model, where the expected value E(·) is taken with respect to the normal density f (x) = √

1

e−(x−μJ )

2 /(2σ 2 ) J

2πσJ with mean μJ and standard deviation σJ (jump volatility). The cumulative normal distribution is given by 1 Φ(x) = √ 2π

x

e−z

2 /2

dz.

(1)

(2)

−∞

For η := E(ex − 1) we obtain in Merton’s model: 1 2 dx 2 2 x x η := E(e − 1) = (e − 1)f (x) dx = √ ex−(x−μJ ) /(2σJ ) − 1 = eμJ +σJ /2 − 1. σJ 2π R

(3)

R

Using the reduced Lévy triplet

(σ 2 , r

−

σ2 2

− λη + λE(y), λf ) in the general PIDE stated in Theorem 2.1 we get:

  ∞   σ2 σ2 vxx + r − − λη + λE(y) vx − rv + λ v(t, x + y) − v(t, x) − vx (t, x)y f (y) dy = 0, vt + 2 2

(4)

−∞

and after some calculating   ∞ σ2 σ2 vxx + r − − λη vx − (r + λ)v + λ vt + v(t, x + y)f (y) dy = 0. 2 2

(5)

−∞

After reversing the time direction τ := T − t and setting u(τ, ·) := v(T − τ, ·) we change variables in the integral term z = x + y, y = z − x, dy = dz and arrive at the final formulation of what we will call PIDE W: uτ −

  ∞ σ2 σ2 uxx − r − − λη ux + (r + λ)u − λ u(τ, z)f (z − x) dz = 0, 2 2

on (0, T ] × R,

−∞

u(0, x) = H (ex ),

∀x ∈ R.

(6)

1690

E.W. Sachs, A.K. Strauss / Applied Numerical Mathematics 58 (2008) 1687–1703

2.3. Transformation to PIDE N It is well known, that the convection term in (6) can cause numerical instability, especially in cases where the volatility is itself modeled by a stochastic process or is a function of z and τ . If this phenomenon occurs, additional measures have to be taken like adding stabilizing terms (e.g. upwinding, artificial diffusion) which themselves often cause a drop in accuracy of the solution. Another disadvantage when using a second order central difference approximation for the convection term is that the resulting linear system has a non-symmetric coefficient matrix. Instead, we apply a variable transformation to the PIDE W, as has been done in the Black–Scholes framework, in order to arrive at a PIDE formulation in such a way that the convection term ux disappears. This is in particular useful for Merton’s model because the discretized PIDE yields a linear equation with a symmetric Toeplitz coefficient matrix which gives second order accuracy. For convenience of notation, let ζ := r − 12 σ 2 − λη. Lemma 2.2. Let w ∈ C 1,2 ((0, T ] × R) ∩ C 0 ([0, T ] × R) satisfy 1 wτ − σ 2 wξ ξ + (r + λ)w − λ 2

∞ w(τ, z)f (z − ξ ) dz = 0,

on (0, T ] × R,

−∞

w(0, ξ ) = H (e ), ξ

∀ξ ∈ R.

(7)

Then u(τ, x) := w(τ, x + ζ τ ) solves the original equation (6). Proof. With the new variable ξ := x + ζ τ we have uτ = wτ + ζ wξ ,

ux = wξ ,

uxx = wξ ξ .

Note also that for fixed τ we can rewrite the integral term as ∞

∞ w(τ, z)f (z − ξ ) dz =

−∞

  w(τ, z + ζ τ )f z − (ξ − ζ τ ) dz =

−∞

∞ u(τ, z)f (z − x) dz

−∞

Restating (7) in terms of u we obtain σ2 uτ − ζ ux − uxx + (r + λ)u − λ 2

∞ u(τ, z)f (z − x) dz = 0,

∀(τ, x) ∈ (0, T ] × R,

−∞

which completes the proof.

2

This is called PIDE N (since it has No convection term) to distinguish it from its equivalent form W (With convection term) stated in the last section. 3. Discretization 3.1. Domain truncation The domain for x of the PIDE W is usually restricted to (−x, ˆ x) ˆ for some sufficiently large x. ˆ The corresponding restricted domain for PIDE N is chosen to be Ω := (ξ− , ξ+ ) with ξ− := −xˆ + ζ T and ξ+ := xˆ + ζ T , since then in the final time T the function u(T , x) on (−x, ˆ x) ˆ is equal to w(T , ξ ) on (ξ− , ξ+ ). As shown in Cont and Voltchkova [14], the truncation error decreases exponentially with the size of the domain if we use the discounted payoff function as a boundary condition for ξ → ∞. For a European call we get w(τ, ξ ) → 0, w(τ, ξ ) ∼ e

ξ −ζ τ

(ξ → −∞), − Ke−rτ ,

(ξ → +∞).

E.W. Sachs, A.K. Strauss / Applied Numerical Mathematics 58 (2008) 1687–1703

1691

Hence the integral term R for the region outside of Ω +∞ R(τ, ξ, ξ+ ) = (ez−ζ τ − Ke−rτ )f (z − ξ ) dz, ξ+

yields, using the cumulative normal distribution Φ, in case of Merton’s model     ξ − ξ+ + σJ2 ξ − ξ+ 2 − Ke−rτ Φ R(τ, ξ, ξ+ ) = e−ζ τ eξ +σJ /2 Φ σJ σJ which is used in the discretized boundary condition in (9). 3.2. Finite difference discretization We set similarly to [1] ξi := ξ− + (i − 1)h τm := (m − 1)k

with i = 1, . . . , n, h = (ξ+ − ξ− )/(n − 1) = 2x/(n ˆ − 1),

with m = 1, . . . , q, k = T /q.

The solution wim of the finite difference scheme approximates the true solution at the points w(τm , ξi ). The integral is approximated by the composite trapezoidal rule on Ω and the estimate R on Ω c using fi−j := f (ξj − ξi ) = f ((j − i)h) for i = 2, . . . , n − 1: 

n−1  h w(τm , z)f (z − ξi ) dz ≈ wjm fi−j + wnm fi−n + R(τm , ξi , ξ+ ). w1m fi−1 + 2 2 R

j =2

Due to second order accuracy and stability in time, we use a backward difference formula of second order (BDF2) for m  2 and the second order central difference scheme for the spatial derivatives: ( 32 wim − 2wim−1 + 12 wim−2 )/k, for m  2, wτ (τm , ξi ) ≈ for m = 1, (wim − wim−1 )/k, m m wξ ξ (τm , ξi ) ≈ (wi+1 − 2wim + wi−1 )/ h2 .

We use a fully implicit scheme and arrive at the following linear system with symmetric coefficient matrix:

1, m = 1, m m Aw = b , A := γ0 I + C + D, γ0 = 3 2 , m  2, I is the identity matrix and C, D are matrices defined by ⎧ −kσ 2 /2h2 , i = j − 1, 2  i  n − 1, ⎪ ⎨ 2 /h2 , i = j, 2  i  n − 1, (r + λ)k + kσ cij = 2 2 ⎪ i = j + 1, 2  i  n − 1, ⎩ −kσ /2h , 0, otherwise, ⎧ 1 ⎨ − 2 khλfi−j , 2  i  n − 1, j ∈ {1, n}, dij = −khλfi−j , 2  i  n − 1, 2  j  n − 1, ⎩ 0, otherwise, m bi = kλR(τm , ξi , ξ+ ) + γ1 wim−1 + γ2 wim−2 , i = 2, . . . , n − 1,

0, m = 1, 1, m = 1, γ1 = γ2 = −1/2, m  2. 2, m  2,

(8)

(9)

The first and last component of the right hand side b are given by the boundary conditions b1 = 0 and bn = γ0 (eξ+ −ζ τm − Ke−rτm ). The initial solution vector w 0 = [w10 , . . . , wn0 ]T = [H (eξ1 ), . . . , H (eξn )]T represents the payoff function H evaluated at n asset prices. Reformulating the system by using the known values for w0m and wnm from the boundary conditions one obtains a (n − 2) × (n − 2) system with the coefficient matrix T .

1692

E.W. Sachs, A.K. Strauss / Applied Numerical Mathematics 58 (2008) 1687–1703

Definition 3.1. T ∈ Cn×n is called Toeplitz, if T is determined by the 2n − 1 scalars t−(n−1) , . . . , tn−1 with Tij = tj −i for all i and j . Lemma 3.2. The coefficient matrix T = (Aij )2i,j n−1

(10)

of the linear system Aw m = bm is dense and belongs to the class of symmetric Toeplitz matrices. The entries of the matrix T are given by ⎛ γ0 + α + d 0 β + d−1 ··· . ⎜ ⎜ β + d1 γ0 + α + d 0 . . ⎜ .. .. .. T =⎜ . . ⎜ . ⎜ .. ⎝ . dn−4 dn−3 dn−4 ···

d4−n

d3−n

..

d4−n .. .

.

..

. β + d1

β + d−1 γ0 + α + d 0

⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠

(11)

with α = σ2

k + (r + λ)k, h2

β =−

σ2 k , 2 h2

di = −khλfi ,

1 2 2 fi = f (−ih) = √ e−(ih) /(2σJ ) . 2π σJ

(12)

Note that this property requires the choice of a uniform grid, but the Toeplitz property allows us to significantly reduce the effort needed to solve the linear system in each time step as we shall see in the next section. 4. Circulant matrices In this section we review some definitions and key facts for matrices with the Toeplitz property and the circulant property. A special case of Toeplitz matrices are circulant matrices: Definition 4.1. C ∈ Cn×n is called circulant if it is a Toeplitz matrix where each column is a circular shift of its preceding column or c−i = cn−i for i = 1, . . . , n − 1. Fourier matrices are used in the spectral decomposition of circulant matrices. Definition 4.2. The unitary and symmetric matrix Fn ∈ Cn×n is called Fourier matrix if it is of the form ⎞ ⎛ 1 1 ··· 1 1 ⎜1 ω ωn−2 ωn−1 ⎟ ⎟ 1 ⎜ .. .. .. ⎟, ⎜ Fn = √ ⎜ . . . ⎟ n⎝ ω(n−2)(n−1) ⎠ 1 ωn−2 2 1 ωn−1 · · · ω(n−1)(n−2) ω(n−1)

(13)

with ω = e−2πi/n . A matrix-vector product with a Fourier matrix like x=

√ ∗ nFn (a0 , . . . , an−1 )T

or xj =

n−1 

ak e2πij k/n ,

j = 0, 1, . . . , n − 1,

k=0

can be computed efficiently via an inverse Fast Fourier Transform (FFT), i.e. in O(n log n) instead of O(n2 ). With regard to the close relationship between Fourier matrix and FFT, the following result is instrumental for all considered methods in this work:

E.W. Sachs, A.K. Strauss / Applied Numerical Mathematics 58 (2008) 1687–1703

1693

Theorem 4.3. (See [15].) Let Cn ∈ Rn×n be circulant. Then it has the decomposition Cn = Fn∗ ΛFn ,

(14)

where Λ = diag(λ0 , . . . , λn−1 ), λj being the j th eigenvalue of Cn , j = 0, . . . , n − 1. The theorem has in particular two very important consequences: First, the eigenvalues λj of Cn can be easily found using this decomposition. λj =

n−1 

ck e−2πij k/n ,

j = 0, . . . , n − 1.

k=0

Second, the matrix-vector product x = Cn a = Fn∗ ΛFn a can be computed efficiently in four steps, see for example Golub and Van Loan [18]. This way one can compute the circulant matrix-vector product Cn a by three FFTs and one vector multiply, i.e. in O(n log n) operations. Circulant matrices belong to a larger class of matrices called {ω}-circulant matrices: Definition 4.4. Let ω = eiθ0 with θ0 ∈ [−π, π]. An n × n matrix W is said to be an {ω}-circulant matrix if it has the spectral decomposition Wn = Ωn∗ Fn∗ Λn Fn Ωn .

(15)

Here Ωn = diag[1, ω−1/n , . . . , ω−(n−1)/n ], Fn denotes the Fourier matrix and Λn is a diagonal matrix containing the eigenvalues of Wn . In this paper we confine ourselves to the special cases ω = 1 and ω = −1, that means the circulant matrices and the so-called skew-circulant matrices respectively. Also in the skew-circulant case one needs only O(n log n) operations to find Λn and the matrix-vector product needs the same amount of work as in the circulant case. It is possible to embed a Toeplitz matrix in a circulant one. For an n × n Toeplitz matrix Tn , the Toeplitz matrixvector multiplication Tn y can be computed with three FFTs by first embedding Tn into a 2n × 2n circulant matrix by



Tn Bn

Bn Tn

    y Tn y = , 0 ∗

⎛ 0 ⎜ t1−n ⎜ . ⎜ Bn = ⎜ .. ⎜ ⎝ t−2 t−1

tn−1 0 t1−n t−2

···

t2

t1 t2 .. .

tn−1 0 .. . ···

..

.

..

. t1−n

tn−1 0

⎞ ⎟ ⎟ ⎟ ⎟. ⎟ ⎠

(16)

The importance of this embedding as noted by Strang [33] lies in the fact that the conjugate gradient method (CG) applied to a Toeplitz system needs only O(n log n) operations per iteration. There exists another interesting fact that we use in the preconditioned CG framework to reduce the number of FFTs needed in the computation of the matrix-vector product: According to Pustyl’nikov [31] one can easily partition any Toeplitz matrix Tn into a sum of a circulant matrix Un and a skew-circulant matrix Vn with ⎛

t0

⎜ ⎜ t1 + t−n+1 1⎜ Un = ⎜ 2⎜ ⎜ ⎝ tn−1 + t−1

t−1 + tn−1 t0 .. .

t−(n−1) + t1 ..

.

..

. ..

. t1 + t−n+1

t−1 + tn−1 t0

⎞ ⎟ ⎟ ⎟ ⎟, ⎟ ⎟ ⎠

(17)

1694

E.W. Sachs, A.K. Strauss / Applied Numerical Mathematics 58 (2008) 1687–1703

and

⎛

t0

⎜ ⎜ t1 − t−n+1 1⎜ Vn = ⎜ 2⎜ ⎜ ⎝

t−1 − tn−1

t−(n−1) − t1

t0 .. .

..

.

..

. ..

. t1 − t−n+1

tn−1 − t−1

t−1 − tn−1 t0

⎞ ⎟ ⎟ ⎟ ⎟. ⎟ ⎟ ⎠

(18)

5. Eigenvalues of the coefficient matrix Before we implement the conjugate gradient method, we need to ensure that the coefficient matrix T is both symmetric and positive definite. For Merton’s model the symmetry is clear by construction of the discretization scheme. In the convergence behavior of the CG method, the location of the eigenvalues plays an important role. Since T is Toeplitz we can use results from the theory of Toeplitz matrices. The definition of a generating function g which generates a unique infinite Toeplitz matrix T∞ is defined by the Fourier series ∞  tj e−ij x , x ∈ [−π, π], (19) g(x) = j =−∞

where tj is the entry of T∞ on the j th subdiagonal or, in other words, the (i, j )th entry of T∞ is given by ti−j . This function generates Toeplitz matrices in the sense that the diagonals and subdiagonals are given by 1 tj = 2π

π g(x)eij x dx,

j = 0, ±1, ±2, . . . ,

−π

and we denote the so-defined matrices by Tl [g] ∈ Rl×l , l ∈ N. If g is real-valued, we have t−j

1 = 2π

π

g(x)e−ij x dx = t¯j ,

∀j,

−π

and hence Tl [g] must be Hermitian for all l. If g is also an even function, that is g(−x) = g(x), then Tl [g] is real and symmetric for all l. For generating functions there exists a close relationship between g and the spectra of Tl . Not only forms the minimum (maximum) of g a lower (upper) bound for the eigenvalues of Tl , also one can deduce information about the eigenvalue distribution from g, see Theorem 5.7. Theorem 5.1. (See [19].) Let g be a real-valued function in L1 [−π, π]. Then the spectrum σ (Tl ) of Tl satisfies σ (Tl ) ⊆ [gmin , gmax ],

∀l  1,

where gmin and gmax are the essential infimum and the essential supremum of g respectively. Moreover, if gmax > gmin , then gmin < λmin (Tl )  λmax (Tl ) < gmax . In particular, if gmin > 0, then Tl is positive definite for all l ∈ N. The generating function g(x) can be derived explicitly from the coefficient matrix T arising from the discretization scheme for Merton’s model. Using (11) and (12) we obtain for m  2 khλ k 3 + (r + λ)k + − √ , 2 h2 2π σJ σ2 k khλ −h2 /(2σ 2 ) J , t1 = t−1 = − −√ e 2 2 h 2π σJ

t0 = σ 2

E.W. Sachs, A.K. Strauss / Applied Numerical Mathematics 58 (2008) 1687–1703

khλ −(ih)2 /(2σ 2 ) J , ti = t−i = − √ e 2π σJ

1695

i  2.

(20)

Hence by using the definition of a generating function (19) we obtain 

∞  khλ −(j h)2 /(2σJ2 ) g(x) = − √ e cos(j x) + 2t1 cos x + t0 . 2 2π σJ j =2

(21)

Next we show that g is well defined. Lemma 5.2. The real-valued and even function g defined in (21) is in Wiener class, i.e. its Fourier coefficients tj are absolutely summable. Proof. With tj defined in (20) for j ∈ {0, 1, 2, . . .} we have ∞ khλ  −(j h)2 /(2σ 2 ) J <∞ |tj |  √ e 2π σ J j =2 j =2

∞ 

2

by the ratio test. Hence it follows that g is in Wiener class.

In order to show the positive definiteness we establish a lower bound for g. Lemma 5.3. Let g be the function defined in (21) and λ, r > 0. Then a lower bound for g on [−π, π] is    3 λ h g(x)  + k r + λ − √ 1+ . 2 σJ 2π

(22)

Proof.



∞  khλ −(j h)2 /(2σJ2 ) g(x) = − √ e cos(j x) + 2t1 cos(x) + t0 2 2π σJ j =2 

∞  k k khλ 3 −(j h)2 /(2σJ2 ) = −√ e cos(j x) + 1 − σ 2 2 cos(x) + σ 2 2 + (r + λ)k + . 2 2 h h 2π σJ j =1

For convenience define c := h2 /(2σJ2 ). The sum can be estimated by ∞  j =1

e−cj cos(j x)  2

∞  j =1

e−cj  2

∞ 0

1 2 e−cy dy = √ 2c

∞

e−y˜

2 /2

dy˜ =

σJ . 2h

0

The estimate above gives     k khλ σJ 3 g(x)  √ − − 1 + 1 − cos(x) σ 2 2 + (r + λ)k + h 2 h 2π σJ kλ 3 khλ 2  −√ + (r + λ)k + . −√ 2 2π σJ 2π Using a rough estimate, this results in the following lemma: Lemma 5.4. The matrices Tl [g], l ∈ N, are symmetric and their smallest eigenvalue is bounded from below by if the quantity h used in the definition of g satisfies √ 0 < h < ( 2π − 1)σJ .

3 2

+ rk, (23)

1696

E.W. Sachs, A.K. Strauss / Applied Numerical Mathematics 58 (2008) 1687–1703

This estimate is obvious from the previous lemma due to   3 kλ √ h g(x)  + rk + √ . 2π − 1 − 2 σJ 2π Clearly, one could refine this estimate using (22) directly. We apply this result to the matrix T as defined in (11). If we choose n ∈ N such that 2xˆ n>1+ √ , (24) ( 2π − 1)σJ then h = 2x/(n ˆ − 1) as defined in the finite difference discretization satisfies (23) in Lemma 5.4. Hence, for the function g as defined in (21) using h as discussed before, we obtain from Lemma 5.4 that all Tl [g], l ∈ N, are symmetric and positive definite. Since we have for T as defined in (11) that T = Tn−2 [g], we obtain that T is positive definite. This results in the following Remark 5.5. The matrix T defined in (11) is symmetric and positive definite if n satisfies (24). We also obtain information on the distribution of the eigenvalues of Tl . (l)

Definition 5.6. Let g be a real-valued function in L1 [−π, π]. A sequence (αk ) is said to be equally distributed as g(x) if π l   1 1 (l) F (αk ) = F g(x) dx lim l→∞ l 2π k=1

−π

for any continuous function F with bounded support. Theorem 5.7. (See [19].) Let g ∈ L2 [−π, π]. Then the singular values of the matrices Tl generated by g are equally distributed as |g(x)|. In particular, for real-valued g, the eigenvalues of Tl are equally distributed as g(x). To illustrate this theorem, let us consider the generating function g with the parameters h = 0.125, k = 0.1, r = 0, σJ = 0.5, σ = 0.2 and λ = 0.1. Then we consider two cases: T63 [g] and T127 [g], obtained by adjusting the computational domain (−x, ˆ x) ˆ such that h remains constant.

Fig. 1. Generating function and eigenvalues of generated Toeplitz matrices for Merton’s model and m  2, l = 63 and l = 127. Parameters: h = 0.125, k = 0.1, r = 0, σJ = 0.5, σ = 0.2 and λ = 0.1.

E.W. Sachs, A.K. Strauss / Applied Numerical Mathematics 58 (2008) 1687–1703

1697

Both plots in Fig. 1 contain the graph of the same generating function g(x) as well as the separately computed eigenvalues of the matrices T63 and T127 that are generated by g(x). These eigenvalues are plotted against a uniform grid on the interval [0, π] and apparently follow the curve g independently of the dimension l; they are equally distributed as g. The larger l, the more densely the eigenvalues are located along the graph of g, especially near the endpoints of the interval. 6. Preconditioning Improving the eigenvalue distribution of the Toeplitz system Tn x = b can be done efficiently by preconditioning. In order to compete with the tridiagonal splitting approach of Almendral and Oosterlee [1] we need a fast matrix-vector product evaluation. We search for circulant preconditioners in order to exploit the use of FFT. The first circulant preconditioner we consider is Strang’s preconditioner Sn [33]. Several other circulant preconditioners have been pron [35] posed, among them T. Chan’s optimal preconditioner c(Tn ) [13], Tyrtyshnikov’s super-optimal preconditioner C and R. Chan’s preconditioner Rn [7]. We find that all these preconditioners perform very well and give a more detailed account of the properties of Sn . 6.1. Strang’s preconditioner Gilbert Strang first proposed a circulant preconditioner in [33]. His idea was to keep the central diagonals of the Toeplitz matrix Tn and to wrap them around to obtain a circulant matrix. The diagonals of the Strang preconditioner Sn = (sk−l )0k,l 0, there exists a constant c() > 0 such that the error vector ek of the preconditioned conjugate gradient method at the kth iteration satisfies ek  c() k , e0 −1/2

where v 2 = v ∗ Sn

−1/2

Tn S n

v.

1698

E.W. Sachs, A.K. Strauss / Applied Numerical Mathematics 58 (2008) 1687–1703

It follows in particular that the number of iterations required for convergence is independent of the size of the matrix Tn when n is large. To sum up, Sn seems to be indeed an adequate choice. Theorem 6.1 shows that Sn is an optimal approximation of Tn , it copies the symmetry for real Tn , it will be positive definite under the assumptions of Theorem 6.2 and its construction is simple. 6.2. Other circulant-type preconditioners Tony Chan proposed in [13] a circulant preconditioner c(Tn ) that is optimal with respect to the Frobenius norm, that means it is defined as the minimizer of Cn − Tn F over the space of all n × n circulant matrices Cn . It is called the optimal circulant preconditioner and has the nice feature that c(Tn ) is s.p.d. whenever Tn is. Theorem 6.5. (See [9].) Let Tn ∈ Cn×n and c(Tn ) be the minimizer of Cn −Tn F over all circulant n×n matrices Cn . Then c(Tn ) is uniquely determined by Tn . Moreover, (i) c(Tn ) is given by c(Tn ) =

n−1  1 j =0

n



 tpq Π j ,

p−q≡j mod n

where Π is the n × n circulant “push” matrix defined by ⎞ ⎛ 0 1 ⎟ ⎜1 0 ⎟. Π =⎜ .. .. ⎠ ⎝ . . 0 1 (ii) c(Tn ) is also given by

0

c(Tn ) = Fn∗ δ(Fn Tn Fn∗ )Fn ,

(26)

where Fn is the Fourier matrix defined in (13) and δ(·) denotes the diagonal matrix with the same diagonal as the matrix in the argument. Note that c(Tn ) can be defined for any square matrix, unlike the Strang preconditioner. The eigenvalues are the entries of the diagonal matrix δ(Fn Tn Fn∗ ). Although c(Tn ) inherits symmetry and positive definiteness from Tn [9], it is not necessarily regular. The non-singularity of Tn cannot guarantee δ(Fn Tn Fn∗ ) to be non-singular. For a Toeplitz matrix Tn the diagonals of c(Tn ) are given by

((n − j )tj + j tj −n )/n, 0  j < n, cj = (27) −n < j < 0. cn+j , Therefore, in our case the preconditioner can be constructed in only O(n) operations. It can also be shown [8] that the spectrum of the preconditioned matrix will be clustered around 1 for Toeplitz matrices being generated by a positive function in the Wiener class. Another possible preconditioner would be the minimizer of In − Cn−1 Tn F over all non-singular circulant matrices Cn . Tismenetsky [34] and Tyrtyshnikov [35] propose this preconditioner independently and Tyrtyshnikov [35] called it the super-optimal circulant preconditioner. n be the super-optimal circulant preconditioner Theorem 6.6. (See [9,10].) Let Tn ∈ Cn×n be positive definite. Let C for Tn , defined by n Tn F = min I − Cn−1 Tn F , I − C where the minimum is taken over all n × n non-singular circulant matrices Cn . Then   n = c(Tn∗ ) c(Tn Tn∗ ) −1 , C where c(A) denotes T. Chan’s preconditioner constructed from A.

(28)

E.W. Sachs, A.K. Strauss / Applied Numerical Mathematics 58 (2008) 1687–1703

1699

Tyrtyshnikov [35] also showed that it inherits both symmetry and positive definiteness from Tn , just like the one from T. Chan. In this regard these two preconditioners are better than Strang’s preconditioner. n is uniquely defined. Theorem 6.7. (See [35].) Let Tn ∈ Rn×n non-singular. Then the super-optimal preconditioner C n is non-singular and, moreover, also positive definite. If Tn = TnT , then C n = C nT . If Tn is positive definite, then C n will cause the eigenvalues of the preconditioned matrix C n Tn to be Like the previous two preconditioners also C clustered around 1 [10]. R. Chan [7] proposed a preconditioner Rn whose diagonals are defined by

+ t , 0  k < n, t rk = k−n k 0 < −k < n, r−k , where t−n is set equal to zero and Tn assumed to be real. Note that Rn is symmetric by construction. He also showed that the eigenvalue distribution of Rn is asymptotically the same as for Strang’s preconditioner Sn for large n and limn→∞ Sn − Rn 2 = 0. An optimal preconditioner was also proposed by T. Huckle [21]. This preconditioner is optimal in the sense that for the solution the norm −1/2

I − Cn

−1/2

T n Cn

F

(29)

is minimal. It is shown that this preconditioner and T. Chan’s preconditioner are asymptotically equivalent with respect to the spectral norm. 6.3. Comparison of the preconditioners To compare the effect of various preconditioners, we graph the eigenvalue distribution of the preconditioned coefficient matrix in Fig. 2. To test the performance of the preconditioners for the problem under investigation we applied them to the coefficient matrix and measured how many iterations the conjugate gradient method needs for the preconditioned system. The outcome is presented in Table 1. Note that without preconditioning the number of iterations increases steadily, whereas for the preconditioned systems it stays almost the same. All preconditioners are obviously effective and it may be surprising that Rn is equally successful as Sn ; this is because for our specific system Rn turns out to be nearly equal to Sn . For a comparison of various preconditioners based on analytical methods one can use different convolution kernels. It has been shown in [11] that both the R. Chan preconditioner and Strang preconditioner are obtained by convoluting with a Dirichlet kernel, however of different length. Recalling that for Sn the construction comes for free we conclude that Sn will be the best choice for our system. For this reason we will confine our numerical tests to Strang preconditioned systems.

Table 1 Number of CG iterations for different preconditioners. Parameters: k = T /40, T = 1, truncation point xˆ = 4, r = 0, σ = 0.2, σJ = 0.5, λ = 0.1 and strike K = 1 n

Unprecond.

Strang Sn

T. Chan c(Tn )

Tyrtyshnikov Yn

R. Chan Rn

17 33 65 129 257 513 1025 2049

3 3 4 6 9 18 36 73

3 3 3 4 4 4 4 4

3 3 3 4 4 5 5 5

3 3 3 4 4 5 6 8

3 3 3 4 4 4 4 4

1700

E.W. Sachs, A.K. Strauss / Applied Numerical Mathematics 58 (2008) 1687–1703

Fig. 2. Eigenvalue distribution of the preconditioner coefficient matrix for several preconditioners. Parameters r = 0, T = 1, xˆ = 4, σ = 0.2, σJ = 0.5, λ = 0.1, n = 65, m = 40 and K = 1.

7. CG solver For the matrix-vector product two FFTs of length 2n are needed by using the circulant embedding (16). Note that two FFTs are sufficient as long as the coefficient matrix does not change over the time. We simply compute the eigenvalues of its circulant embedding once at the beginning, thus saving one of the usual three necessary FFTs. Additionally, we would need two FFTs of length n to solve the circulant system Mzj +1 = rj +1 . In fact, this overhead of the PCG versus the unpreconditioned CG where we only need two FFTs of length 2n can be reduced to zero by using an idea of Huckle [21]. If we let PCG run in the Fourier domain, then we get the preconditioning virtually for free. The key ingredient to this implementation is the decomposition Tn = Un + Vn defined in (17) and (18). For leftpreconditioned systems an implementation can be found in Ng [30]. We will implement the same trick equivalently for −1/2 −1/2 = a split preconditioned framework with Strang’s preconditioner Sn to guarantee symmetry. Since Sn−1 = Sn Sn −1/2 −1/2 Sn (Sn )T this preconditioned system is −1/2

Sn

−1/2

Tn S n

−1/2

y = Sn

b,

1/2

y = Sn x.

(30)

˜ n . Then Sn−1/2 = Fn∗ Λ1 Fn with Λ1 = Λ˜ 1/2 . The matrix Tn can be decomposed in a sum of a circuLet Sn−1 = Fn∗ ΛF lant Un and a skew-circulant matrix Vn by (17) and (18) with Un = Fn∗ Λ2 Fn and Vn = Ωn∗ Fn∗ Λ3 Fn Ωn , where Λ1 , Λ2 , Λ3 and Ωn are diagonal matrices. Substituting this into Eq. (30) yields −1/2

−1/2 1/2 −1/2 −1/2 1/2 −1/2 Sn x = Sn (Un + Vn )Sn Sn x = Sn b, Fn∗ Λ1 Fn (Fn∗ Λ2 Fn + Ωn∗ Fn∗ Λ3 Fn Ωn )(Fn∗ Λ1 Fn )(Fn∗ Λ−1 1 Fn )x

Sn

Tn S n

Λ1 (Λ2 + Fn Ωn∗ Fn∗ Λ3 Fn Ωn Fn∗ )Λ1 x˜ = b˜

= Fn∗ Λ1 Fn b, (31)

˜ with x˜ = Λ−1 1 Fn x and b = Λ1 Fn b. The cost for solving the preconditioned s.p.d. system (31) with CG is dominated by the matrix-vector product which can be accomplished in essentially four FFTs of length n, so the cost is approximately the same like for the unpreconditioned system where we needed two FFTs of length 2n. Thus the overhead can indeed be approximately reduced to zero and we see that PCG solves the system at cost comparable to that of the splitting method. The cost per iteration is of order O(n log n). Furthermore, by Theorem 6.4 PCG converges within O(1) iterations, so that in each time step just O(n log n) operations are needed. As for the necessary storage, only some vectors with the diagonals of Λ1 , Λ2 , Λ3 , Ωn and the right-hand side b are needed. A careful comparison of numerical techniques for the solution of Toeplitz systems can be found in [22]. It involves the solution by Levinson-type methods, superfast solvers and optimally preconditioned CG solvers. The results include also a comparison of flop counts for these techniques.

E.W. Sachs, A.K. Strauss / Applied Numerical Mathematics 58 (2008) 1687–1703

1701

8. Numerical results There exists an analytical solution for an European call option in form of an infinite series found by Merton [29]. We compute the numerical error by using a version of this infinite series given in Wilmott [36]. For μJ = 0 the analytical solution is then w(t, s) =

∞ −λ(1+η)τ  e (λ(1 + η)τ )m VBS (τ, s, K, rm , σm ), m!

m=0

2

where τ := T − t is the time to maturity, η = eσJ − 1 denotes the expected percentage change in the stock price originating from a jump, σm2 := σ 2 + mσJ2 /τ the volatility, rm := r − λη + m log(1 + η)/τ and VBS the Black–Scholes price of a call [4]. Problems arise if the interest rate r > 0 and the volatility σ the underlying asset is small, because in this case the convection term gains more influence and oscillations arise particularly in a neighborhood of the strike price, see Fig. 3. These oscillations appear in a neighborhood around the strike price, which is an important area for a lot of options.The splitting method is used with the discretization scheme W1 and PCG is used with scheme N 1. The variable x denotes the log-asset price and the parameters are: r = 0.07, T = 1, x¯ = 4, σ = 0.01, σJ = 0.5, λ = 0.1, n = 129, m = 10 and K = 1. In the following we choose data with not so small volatilities and show that the fully implicit scheme presented here is competitive and slightly faster than the semi-implicit scheme in Almendral and Oosterlee [1]. We first apply the method in [1] to the original equation W1 and then we apply our scheme to N 1. (i) Tridiagonal splitting iterations for scheme W1. (ii) PCG for scheme N 1. Let us start with the results for (i). As stopping criterion we use xk+1 − xk ∞ < 10−8 ,

(32)

where xk denotes the kth iterate. The initial guess x0 is chosen to be the solution from the last time step. Observe in Table 2 the second order accuracy both at the strike price and on the whole interval, the latter being mirrored by the l ∞ error. By l ∞ error we mean the infinity norm of the difference between the numerical solution vector and the analytical solution evaluated at the grid points at the final time τ = T . Furthermore, concerning the required workload the Time ratio in Table 2 is of interest. It is the fraction of the run time in the current row divided by the run time of the previous one. Note that although all computations were done in M ATLAB, the Thomas solver was implemented and pre-compiled in C to obtain better comparability since PCG makes heavy use of the fft(·) function in M ATLAB, which in turn uses already compiled code. The run time

Fig. 3. δ(x): difference between the analytical and the numerical solutions.

1702

E.W. Sachs, A.K. Strauss / Applied Numerical Mathematics 58 (2008) 1687–1703

Table 2 FD results with tridiagonal splitting method (using pre-compiled Thomas algorithm) for Merton’s model using discretization scheme W1. Error at the strike price xK = log(K) and in · ∞ norm over the whole interval for T = 1. Truncation point xˆ = 4, r = 0, σ = 0.2, σJ = 0.5, λ = 0.1 and strike K = 1 n

k

Error at xK

l ∞ error

(sec)

Time ratio

65 129 257 513 1025 2049

T/5 T/10 T/20 T/40 T/80 T/160

0.0044 0.0010 2.5e–4 6.3e–5 1.6e–5 4.0e–6

0.0051 0.0013 3.2e–4 7.9e–5 2.0e–5 5.0e–6

0.06 0.14 0.49 1.88 7.50 30.2

2.3 3.5 3.8 4.0 4.0

Table 3 FD results with split preconditioned CG method (using Strang’s preconditioner) for Merton’s model using discretization scheme N 1. Error at the strike price xK = log(K) and in · ∞ norm over the whole interval for T = 1. Truncation point xˆ = 4, r = 0, σ = 0.2, σJ = 0.5, λ = 0.1 and strike K = 1 n

k

Error at xK

l ∞ error

(sec)

Time ratio

65 129 257 513 1025 2049

T/5 T/10 T/20 T/40 T/80 T/160

9.8e–4 8.1e–6 1.7e–4 3.1e–5 3.2e–6 3.7e–8

0.0024 6.0e–4 1.8e–4 3.8e–5 9.6e–6 2.4e–6

0.05 0.12 0.42 1.62 6.65 27.2

2.4 3.5 3.9 4.1 4.1

increases by a factor of about 4 which shows that the approach only needs O(n log n) operations since we double both space and time grid points from row to row. Next we use the transformed PIDE N , leading us to the approach (ii). Here we use the conjugate gradient method with stopping criterion zk ∞ < 10−8 , where zk denotes the kth residual. As an initial guess for CG we use the solution at the previous time step. We already saw in Section 6.3 that CG will just take a few iterations per time step if we use an appropriate preconditioner, e.g. Strang’s preconditioner, which we use in all computations in Table 3. The accuracy is of second order on the whole interval as we see from Table 3, and especially around the strike price xK the error is mostly about one order of magnitude better than for the splitting approach (i) with convection term. With regard to the needed amount of work again the Time ratio in Table 3 shows that the work is of order O(n log n) since the run time only increases by a factor of 4 despite doubling both the number of space and time grid points. The stability is improved owing to the scheme N 1, see Fig. 3. Without the convection term the oscillations disappear giving a far better approximation around the strike price. Comparing the run times in Tables 2 and 3 yields that the PCG approach (ii) is slightly faster than the splitting approach (i), despite the fact that the tridiagonal solver was implemented in a pre-compiled version. Overall we summarize that the passage from PIDE W to N yields improvements of the error behavior as oscillations of the solution around the strike price are avoided, and that using a preconditioned CG method instead of splitting iterations can shorten the run time. 9. Conclusion For Merton’s model we observe globally a second order convergence for the preconditioned CG (PCG) approach using N 1 improving observations about a loss of accuracy on the entire interval by Almendral and Oosterlee [1]. We proposed a new approach to efficiently solve the PIDE describing Merton’s model by first transforming it in order to

E.W. Sachs, A.K. Strauss / Applied Numerical Mathematics 58 (2008) 1687–1703

1703

eliminate the convection term, and then using a preconditioned CG method to solve the symmetric Toeplitz systems arising in every time step due to an implicit discretization. The transformation stabilizes the solution which might otherwise contain oscillations in a neighborhood of the strike price. Furthermore, the PCG method uses FFTs throughout and the solution is therefore obtained in only O(n log n) operations. PCG is equally effective as the tridiagonal splitting technique and turns out to be even slightly faster if we use an appropriate circulant preconditioner for the Toeplitz systems. References [1] A. Almendral, C. Oosterlee, Numerical valuation of options with jumps in the underlying, Appl. Numer. Math. 53 (2005) 1–18. [2] L. Andersen, J. Andreasen, Jump-diffusion processes: Volatility smile fitting and numerical methods for option pricing, Rev. Derivatives Res. 4 (2000) 231–262. [3] D. Bates, Jumps and stochastic volatility: Exchange rate processes implicit in Deutsche Mark options, Rev. Financial Studies 9 (1996) 69–107. [4] F. Black, M. Scholes, The pricing of options and corporate liabilities, J. Political Economy 81 (3) (1973) 637–654. [5] M. Briani, Numerical methods for option pricing in jump-diffusion markets, PhD thesis, Università degli Studi di Roma “La Sapienza”, 2003. [6] M. Briani, R. Natalini, G. Russo, Implicit-explicit numerical schemes for jump-diffusion processes. Technical Report 38 (4/2004), IAC Report, 2004. [7] R. Chan, Circulant preconditioners for Hermitian Toeplitz systems, SIAM J. Matrix Anal. Appl. 10 (1989) 542–550. [8] R. Chan, The spectrum of a family of circulant preconditioned Toeplitz systems, SIAM J. Numer. Anal. 26 (1989) 503–506. [9] R. Chan, X. Jin, M. Yeung, The circulant operator in the Banach algebra of matrices, Linear Algebra Appl. 149 (1991) 41–53. [10] R. Chan, X. Jin, M. Yeung, The spectra of super-optimal circulant preconditioned Toeplitz systems, SIAM J. Numer. Anal. 28 (1991) 871–879. [11] R. Chan, M. Yeung, Circulant preconditioners for Toeplitz matrices with positive continuous generating functions, Math. Comp. 58 (1992) 233–240. [12] T. Chan, Pricing contingent claims on stocks driven by Lévy processes, Ann. Appl. Probab. 9 (2) (1999) 504–528. [13] T. Chan, An optimal circulant preconditioner for Toeplitz systems, SIAM J. Sci. Stat. Comput. 9 (1988) 766–771. [14] R. Cont, E. Voltchkova, A finite difference scheme for option pricing in jump diffusion and exponential Lévy models, SIAM J. Numer. Anal. 43 (4) (2005) 1596–1626. [15] P. Davis, Circulant Matrices, John Wiley & Sons, New York, 1979. [16] Y. d’Halluin, P. Forsyth, G. Labahn, A penalty method for American options with jump diffusion processes, Numer. Math. 97 (2004) 321–352. [17] Y. d’Halluin, P. Forsyth, K. Vetzal, Robust numerical methods for contingent claims under jump diffusion processes, IMA J. Numer. Anal. 25 (2005) 87–112. [18] G. Golub, C. Van Loan, Matrix Computations, third ed., The Johns Hopkins University Press, Baltimore, 1996. [19] U. Grenander, G. Szegö, Toeplitz Forms and their Applications, second ed., Chelsea Publishing, New York, 1984. [20] S. Heston, A closed-form solution for options with stochastic volatility with applications to bond and currency options, Rev. Financial Studies 6 (2) (1993) 327–343. [21] T. Huckle, Some aspects of circulant preconditioners, SIAM J. Sci. Comput. 14 (3) (1993) 531–541. [22] T. Huckle, Iterative methods for Toeplitz-like matrices, Technical Report SCCM 94-5, Stanford Scientific Computing/Computational Mathematics, 1994. [23] J. Hull, A. White, The pricing of options on assets with stochastic volatilities, J. Finance 42 (2) (1987) 281–300. [24] S. Ikonen, J. Toivanen, Numerical valuation of European and American options with Kou’s jump-diffusion model, in: Amamef Conference on Numerical Methods in Finance, INRIA-Rocquencourt, France, February 2006. [25] S. Kou, A jump-diffusion model for option pricing, Management Science 48 (8) (August 2002) 1086–1101. [26] A.-M. Matache, P.-A. Nitsche, C. Schwab, Wavelet Galerkin pricing of American options on Lévy driven assets, Quantitative Finance 5 (4) (2005) 403–424. [27] A.-M. Matache, C. Schwab, T. Wihler, Fast numerical solution of parabolic integrodifferential equations with applications in finance, SIAM J. Sci. Comput. 27 (2005) 369–393. [28] A.-M. Matache, T. von Petersdorff, C. Schwab, Fast deterministic pricing of options on Lévy driven assets, Math. Modelling Numer. Anal. 38 (1) (2004) 37–71. [29] R. Merton, Option pricing when underlying stock returns are discontinuous, J. Financial Economics 3 (1976) 125–144. [30] M. Ng, Iterative Methods for Toeplitz Systems, Oxford University Press, New York, 2004. [31] L. Pustyl’nikov, On the algebraic structure of the spaces of Toeplitz and Hankel matrices, Soviet Math. Dokl. 21 (1980) 141–144. [32] S. Raible, Lévy processes in finance: Theory, numerics, and empirical facts, PhD thesis, Albert-Ludwigs-Universität Freiburg i. Br., 2000. [33] G. Strang, A proposal for Toeplitz matrix calculations, Stud. Appl. Math. 74 (1986) 171–176. [34] M. Tismenetsky, A decomposition of Toeplitz matrices and optimal circulant preconditioner, Linear Algebra Appl. 156 (1991) 105–121. [35] E. Tyrtyshnikov, Optimal and superoptimal circulant preconditioners, SIAM J. Matrix Anal. Appl. 13 (2) (1992) 459–473. [36] P. Wilmott, Derivatives – The Theory and Practice of Financial Engineering, John Wiley & Sons, 1998.