Journal of Computational Finance

Error analysis in Fourier methods for option pricing

Fabián Crocce, Juho Häppölä, Jonas Kiessling and Raúl Tempone

  • We present an error analysis in using Fourier methods for pricing European options when the underlying asset follows an exponential Levy process.
  • The derived bound is minimised to achieve optimal parameters for the numerical method.
  • We propose a scheme to use the error bound in choosing parameters in a systematic fashion to meet a pre-described error tolerance at minimal cost.
  • Using numerical examples, we present results comparable to or superior to relevant points of comparison

We provide a bound for the error committed when using a Fourier method to price European options, when the underlying follows an exponential Lévy dynamic. The price of the option is described by a partial integro-differential equation (PIDE). Applying a Fourier transformation to the PIDE yields an ordinary differential equation (ODE) that can be solved analytically in terms of the characteristic exponent of the Lévy process. Then, a numerical inverse Fourier transform allows us to obtain the option price. We present a bound for the error and use this bound to set the parameters for the numerical method. We analyze the properties of the bound and demonstrate the minimization of the bound to select parameters for a numerical Fourier transformation method in order to solve the option price efficiently.

Lévy processes form a rich field within mathematical finance. They allow the modeling of asset prices with possibly discontinuous dynamics. An early, and probably the best-known, model involving a Lévy process is the Merton (1976) model, which generalizes the Black and Scholes (1973) model. More recently, we have seen more complex models allowing for more general dynamics of the asset price. Examples of such models include the Kou (2002) model (see also Dotsis et al 2007), the normal inverse Gaussian model (Barndorff-Nielsen 1997; Rydberg 1997), the variance gamma (VG) model (Madan and Seneta 1990; Madan et al 1998) and the Carr–Geman–Madan–Yor (CGMY) model (Carr et al 2002, 2003). For a good exposition on jump processes in finance, we refer the reader to Cont and Tankov (2004) (see also Raible 2000; Eberlein 2001).

The prices of European options with underlying assets driven by the Lévy process are solutions to partial integro-differential equations (PIDEs) (Nualart and Schoutens 2001; Briani et al 2004; Almendral and Oosterlee 2005; Kiessling and Tempone 2011); they generalize the Black–Scholes equation by incorporating a nonlocal integral term to account for the discontinuities in the asset price. This approach has also been extended to cases where the option price features path dependence (see, for example, Boyarchenko and Levendorskii 2002; d’Halluin et al 2004; Lord et al 2008).

The Lévy–Khintchine formula provides an explicit representation of the characteristic function of a Lévy process (see Tankov 2004). As a consequence, one can derive an exact expression for the Fourier transform of the solution of the relevant PIDE. Using the inverse fast Fourier transform (iFFT) method, one may efficiently compute the option price for a range of asset prices simultaneously. Further, in the case of European call options, one may use the duality property presented by Dupire (1997) and iFFT to efficiently compute option prices for a wide range of strike prices.

Despite the popularity of Fourier methods for option pricing, few works can be found on error analysis and related parameter selection for these methods. A bound for the error not only provides an interval for the precise value of the option, but also suggests a method for selecting the parameters of the numerical method. An important work in this direction is Lee (2004), in which several payoff functions are considered for a rather general set of models, whose characteristic function is assumed to be known. Feng and Linetsky (2008) presents the framework and theoretical approach for the error analysis and establishes polynomial convergence rates for approximations of the option prices. For a more contemporary review on the error committed in various FT-related methods, we refer the reader to Boyarchenko and Levendorskii (2011). That paper extends the classical flat Fourier methods by deforming the integration contours on the complex plane; it also looks at the discretely monitored barrier options studied in De Innocentis and Levendorskii (2014).

In this work, we present a methodology for studying and bounding the error committed when using FT methods to compute option prices. We also provide a systematic way of choosing the parameters of the numerical method in a way that minimizes the strict error bound, thus guaranteeing adherence to a pre-described error tolerance. We focus on exponential Lévy processes that may be either diffusive or pure jump. Our contribution is to derive a strict error bound for a Fourier transform method when pricing options under risk-neutral Lévy dynamics. We derive a simplified bound that separates the contributions of the payoff and the process in an easily processed and extensible product form; this is independent of the asymptotic behavior of the option price at extreme prices and strike parameters. We also provide a proof for the existence of optimal parameters of the numerical computation that minimize the presented error bound. When comparing our work with that of Lee (2004), we find that Lee’s work is more general than ours, in that he studies a wider range of processes. However, our results apply to a larger class of payoffs. In test examples of practical relevance, we also find that the bound presented produces comparable or better results than those previously presented in the literature, with an acceptable computational cost.

This paper is organized as follows. Section 2 introduces the PIDE setting in the context of risk-neutral asset pricing; we show the Fourier representation of the relevant PIDE for asset pricing with Lévy processes and use that representation for derivative pricing. In Section 3, we derive a representation for the numerical error and divide it into quadrature and cutoff contributions. We also describe the methodology for choosing numerical parameters to obtain minimal error bounds for the FT method. The derivation is supported by numerical examples using relevant test cases with both diffusive and pure-jump Lévy processes in Section 4. Numerics are followed by conclusions in Section 5.

2 Fourier method for option pricing

Consider an asset whose price at time t is modeled by the stochastic process S=(St), defined by St=S0eXt, where X=(Xt) is assumed to be a Lévy process whose jump measure ν satisfies

  {0}min{y2,1}ν(dy)<.   (2.1)

Assuming the risk-neutral dynamic for St, the price at time t=T-τ of a European option with payoff G and maturity time T is given by


where r is the short rate that we assume to be constant and τ:0τT is the time to maturity. Extensions to nonconstant deterministic short rates are straightforward.

The infinitesimal generator of a Lévy process X is given by (see Applebaum 2004)

  Xf(x) limh0E(f(Xt+h)Xt=x)-f(x)h  
      +{0}(f(x+y)-f(x)-y𝟏|y|1f(x))ν(dy),   (2.2)

where (γ,σ2,ν) is the characteristic triple of the Lévy process. The risk-neutral assumption on (St) implies

  |y|>1eyν(dy)<   (2.3)

and fixes the drift term (see Kiessling and Tempone 2011) γ of the Lévy process to

  γ=r-σ22-{0}(ey-1-y𝟏|y|1)ν(dy).   (2.4)

Thus, the infinitesimal generator of X may be written under the risk-neutral assumption as

  Xf(x) =(r-σ22)f(x)+σ22f′′(x)  
      +{0}(f(x+y)-f(x)-(ey-1)f(x))ν(dy).   (2.5)

Consider g as the reward function in log prices (ie, defined by g(x)=G(S0ex)). Now, take f to be defined as


Then, f solves the following PIDE:


Observe that f and Π are related by

  Π(τ,S0ex)=e-rτf(τ,x).   (2.6)

Consider a damped version of f defined by fα(τ,x)=e-αxf(τ,x); we see that τfα=e-αxXf(τ,x).

There are different conventions for the Fourier transform. Here, we consider the operator such that

  [f](ω)eiωxf(x)dx,   (2.7)

defined for functions f for which the previous integral is convergent. We also use f^(ω) as a shorthand notation of [f](ω). To recover the original function f, we define the inverse Fourier transform as


We have that -1[f^](x)=f(x).

Applying to fα, we get fα^(ω)=f^(ω+iα). Observe also that the Fourier transform applied to Xf(τ,x) gives Ψ(-iω)f^(τ,ω), where Ψ() is the characteristic exponent of the process X, which satisfies E(ezXt)=etΨ(z). The explicit expression for Ψ() is

  Ψ(z)=(r-σ22)z+σ22z2+(ezy-1-(ey-1)z)ν(dy).   (2.8)

From the previous considerations, it can be concluded that

  τf^α=Ψ(α-iω)f^(ω-iα).   (2.9)

Now, f^(ω-iα)=f^α(ω), so f^α satisfies the following ODE:

  τf^α(τ,ω)f^α(τ,ω) =Ψ(α-iω),  
  f^α(0,ω) =g^α(ω).   (2.10)

Solving the previous ODE explicitly, we obtain

  f^α(τ,ω)=eτΨ(α-iω)g^α(ω).   (2.11)

Observe that the first factor on the right-hand side of the above equation is E(e(α-iω)Xτ) (ie, φ1(-iα-ω)), where φτ() denotes the characteristic function of the random variable Xτ:

  φτ(ω)E(τΨ(iω)).   (2.12)

Now, we employ the inverse Fourier transformation to obtain the value function:

  fα(τ,x) =-1[fα^](τ,x)=12πe-iωxfα^(τ,ω)dω,   (2.13)
  fα(τ,x) =1π0+Re[e-iωxf^α(τ,ω)]dω.   (2.14)

As it is typically not possible to compute the inverse Fourier transform analytically, we approximate it by discretizing and truncating the integration domain using trapezoidal quadrature (2.13). Consider the following approximation:

  fα,Δω,n(τ,x) =Δω2πk=-nn-1e-i(k+12)Δωxf^α(τ,(k+12)Δω)   (2.15)
    =Δωπk=0n-1Re[e-i(k+12)Δωxf^α(τ,(k+12)Δω)].   (2.16)

Bounding and consequently minimizing the error in the approximation of f(τ,x) by


is the main focus of this paper and will be addressed in the following section.

Remark 2.1.

Although we are mainly concerned with option pricing when the payoff function can be damped in order to guarantee regularity in the L1 sense, we note here that our main results are naturally extendable to include the Greeks of the option. Indeed, we have by (2.11) that

  f(t,x)=12πe(α-iω)xf^α(τ,ω)dω,   (2.17)

so the Delta and Gamma of the option equal

  Δ(t,x) f(t,x)x=12π(α-iω)e(α-iω)xf^α(τ,ω)dω,   (2.18)
  Γ(t,x) 2f(t,x)x2=12π(α-iω)2e(α-iω)xf^α(τ,ω)dω.   (2.19)

Because the expressions involve partial derivatives with respect to only x, the results in this work are applicable to the computation of Δ and Γ through a modification of the payoff function:

  g^α,Δ(ω)= g^α(ω)(α-iω),   (2.20)
  g^α,Γ(ω)= g^α(ω)(α-iω)2.   (2.21)

When the Fourier space payoff function manifests exponential decay, the introduction of a coefficient that is polynomial in ω does not change the regularity of g^ in a way that would significantly change the following analysis. Last, we note that since we do our analysis for PIDEs on a mesh of x’s, we may also compute the option values in one go and obtain the Greeks with little additional effort, using a finite difference approach for the derivatives.

2.1 Evaluation of the method for multiple values of x simultaneously

The fast Fourier transform (FFT) algorithm provides an efficient way of computing (2.15) for an equidistantly spaced mesh of values for x simultaneously. Examples of works that consider this widely extended tool are Lord et al (2008), Jackson et al (2008), Hurd and Zhou (2010) and Schmelzle (2010).

Similarly, one may define the Fourier frequency ω as the conjugate variable of some external parameter on which the payoff depends. In particular, for the practically relevant case of call options, we can denote the log-strike as k, treat x as a constant and write

  f¯k,α(ω)e(α+iω)kfk(x)dk.   (2.22)

Using this convention, the time dependence is given by

  f~k,α(τ,x)=e(iω+α+1)xφτ(ω-i(α+1))(iω+α)(iω+α+1)   (2.23)

contrasted with the x-space solution

  f^k,α(τ,x)=e(iω-α+1)kφτ(ω+iα)(iω+α)(iω+α+1).   (2.24)

We note that for a call option payoff to be in L1, we demand that α in (2.23) is positive. By omitting the exponential factors that contain the x and k dependence in (2.23) and (2.24), respectively, one can get from (2.23) to (2.24) using the mapping α-α-1. Thanks to this, much of the analysis regarding the x-space transformation generalizes in a straightforward manner to the k-space transform.

3 Error bound

The aim of this section is to compute a bound of the error when approximating the option price f(τ,x) by fα,Δω,n(τ,x), which is defined in (2.15). Considering

  fα,Δω(τ,x)=Δω2πk𝒁e-i(k+12)Δωxf^α(τ,(k+12)Δω),   (3.1)

the total error can be split into a sum of two terms: the quadrature and truncation errors. The former is the error from the approximation of the integral in (2.13) by the infinite sum in (3.1), while the latter is due to the truncation of the infinite sum. Using the triangle inequality, we have

  :=|f(τ,x)-fΔω,n(τ,x)|𝒬+   (3.2)


  𝒬 =eαx|fα(τ,x)-fα,Δω(τ,x)|,  

Observe that each , 𝒬 and depends on three kinds of parameters:

  • those underlying the model and payoff, such as volatility and strike price, which we call physical parameters;

  • those relating to the numerical scheme, such as α and n;

  • auxiliary parameters, which will be introduced in the process of deriving the error bound; these parameters do not enter the computation of the option price, but they need to be chosen appropriately to have as tight a bound as possible.

We start by analyzing the quadrature error.

3.1 Quadrature error

Denote by Aa, with a>0, the strip of width 2a around the real line


The following theorem presents conditions under which the quadrature error goes to zero at a spectral rate as Δω goes to zero. Later in this section, we will discuss simpler conditions to verify the hypotheses and analyze in more detail the case where the process X is a diffusive process, or there are “enough small jumps”.

Theorem 3.1.

Assume that, for a>0,

  • (H1)

    the characteristic function of the random variable X1 has an analytic extension to the set

  • (H2)

    the Fourier transform of gα(x) is analytic in the strip Aa ,

  • (H3)

    there exists a continuous function γL1() such that |fα^(τ,ω+iβ)|<γ(ω) for all ω and for all β[-a,a] .

Then, the quadrature error is bounded by


where Mα,a(τ,x) is given by

  Mα,a(τ,x):=β{-a,a}|e-i(ω+iβ)xf^α(τ,ω+iβ)|dω.   (3.3)

Mα,a(τ,x) equals the Hardy norm (defined in (3.4)) of the function


which is finite.

The proof of Theorem 3.1 is an application of Stenger (1993, Theorem 3.2.1), whose relevant parts we include for ease of reading. Using the notation in Stenger (1993), HAa1 is the family of functions w that are analytic in Aa, such that

  wHAa1:=limε0Aa(ε)|w(z)|d|z|<,   (3.4)


Lemma 3.2 (Stenger 1993, Theorem 3.2.1).

Let wHAa1; then, define

  I(w) =w(x)dx,   (3.5)
  J(w,h) =hj=-NNw(jh),   (3.6)
  ζ(w,h) =I(w)-J(w,h),   (3.7)

and then

  |ζ(w,h)|e-πa/hwHAa12sinh(πa/h).   (3.8)
Proof of Theorem 3.1.

First, observe that H1 and H2 imply that the function w(z)=e-ixz+Δω/2f^α(τ,z+Δω/2) is analytic in Aa. H3 allows us to use dominated convergence theorem to prove that wHAa1 is finite and coincides with Mα,a(τ,x). Applying Lemma 3.2, the proof is completed. ∎

Regarding the hypotheses of Theorem 3.1, the next propositions provide simpler conditions that imply H1 and H2, respectively.

Proposition 3.3.

If α, a and ν are such that

  y>1e(α+a)yν(dy)<𝑎𝑛𝑑y<-1e(α-a)yν(dy)<,   (3.9)

and then H1 in Theorem 3.1 is fulfilled.


Denoting by φ1() the characteristic function of X1, we want to prove that zφ1(z+αi) is analytic in Aa. Considering that φ1(z+αi)=eΨ(iz-α), the only nontrivial part of the proof is to verify that

  zp(z,y)ν(dy)   (3.10)

is analytic in Aa, where p:Aa× is given by


To prove this fact, we demonstrate that we can apply the main result and the only theorem in Mattner (2001), which, given a measure space (Ω,𝒜,μ) and an open subset G, ensures the analyticity of f(,ω)dμ(ω), provided that f:G×Ω satisfies the following: f(z,) is 𝒜-measurable for all zG, f(,ω) is holomorphic for all ωΩ and |f(,ω)|dμ(ω) is locally bounded. In our case, we consider the measure space to be , with the Borel σ-algebra and the Lebesgue measure, G=Aa and f=p. It is clear that p(x,) is Borel measurable and p(,y) is holomorphic. It remains for us to verify that


is locally bounded. To this end, we assume that Re[z]<b (and, since zAa, Im[z]<a) and split the integration domain into |y|>1 and 0<|y|1 to prove that both integrals are uniformly bounded.

Regarding the integral in |y|>1, we observe that

  |p(z,y)|ey(α+Im[z])+1+(ey+1)(α+a+b);   (3.11)

for y<-1, we have ey(α+Im[z])<ey(α-a), while for y>1, we have ey(α+Im[z])<ey(α+a). Using the previous bounds and hypotheses together with (2.1) and (2.3), we obtain the needed bound.

For the integral in 0<|y|1, observe that, denoting f(z,y)=|p(z,y)|, we have f(z,0)=0 for every z, yf(z,0)=0 for every z and |yyf(z,y)|<c for zAa, Re[z]<b, |y|<1. From these observations, we get that the Maclaurin polynomial of degree 1 of yf(z,y) is null for every z. We can bound f(z,y) by the remainder term, which, in our region of interest, is bounded by (c/2)y2; thus, we obtain

  0<|y|1|p(z,y)|ν(dy)c20<|y|1y2ν(dy),   (3.12)

which is finite by the hypothesis on ν. This finishes the proof. ∎

Proposition 3.4.

If, for all b<a, the function xeb|x|gα(x) is in L2(), then H2 in Theorem 3.1 is fulfilled.


The proof is a direct application of Reed and Simon (1975, Theorem IX.13). ∎

We now turn our attention to a more restricted class of Lévy processes; namely, processes such that either σ2>0 or there exists λ(0,2) such that C(λ) defined in (3.13) is strictly positive. For this class of processes, we can state our main result explicitly in terms of the characteristic triplet.

Given λ(0,2), define C(λ) as

  C(λ)=infκ>1{κλ0<|y|<1/κy2ν(dy)}.   (3.13)

Observe that C(λ)0, and, by our assumptions on the jump measure ν, C(λ) is finite. Further, if λ(0,2) is such that

  lim infϵ01ϵλ0<|y|<ϵy2ν(dy)>0,   (3.14)

then C(λ)>0. To see this, note that (3.14) implies the existence of ϵ0, such that


If ϵ0<1, observe that


where, for the first inequality, it was taken into account that 1/ϵλ1 and that the integral is increasing with ϵ. By combining the two previous infima and considering |κ|=1/ϵ, we get that C(λ)>0.

Further, we note that for a Lévy model with finite jump intensity, such as the Black–Scholes and Merton models that satisfy the first of our assumptions, C(λ)=0 for all λ(0,2).

Theorem 3.5.

Assume that α and a are such that (3.9) holds, g^αLAa and either σ2>0 or C(λ)>0 for some λ(0,2). Then, the quadrature error is bounded by



  M~α,a(τ,x) =c{-1,1}ecaxeτΨ(ca)|g^α(ca)|  
      ×exp(-τ(σ22ω2+|ω|2-λ4C(λ)𝟏|ω|>1))dω.   (3.15)

Further, if σ2>0, we have

  M~α,a(τ,x)2πστc{-1,1}ecaxeτΨ(ca)|g^α(ca)|.   (3.16)

Considering hα,a(τ,x,ω), defined by

  hα,a(τ,x,ω)=c{-1,1}|e-i(ω+ica)xf^α(τ,ω+ica)|,   (3.17)

we have that


However, for β(-a,a),

  |e-i(ω+iβ)xf^α(τ,ω+iβ)| =eβx|f^α(τ,ω+iβ)|  
    =eβx|eτΨ(α+β-iω)||g^α(ω+iβ)|.   (3.18)

For the factor involving the characteristic exponent, we have

  |eτΨ(α+β-iω)|=eτRe[Ψ(α+β-iω)].   (3.19)

Now, observe that

  Re[Ψ(α+β-iω)] =(α+β)(r-σ22)+σ22((α+β)2-ω2)  
    +{0}(e(α+β)ycos(-yω)-1-(α+β)(ey-1))ν(dy).   (3.20)

If |ω|1, we bound cos(-yω) by 1, getting

  Re[Ψ(α+β-iω)] (α+β)(r-σ22)+σ22((α+β)2-ω2)  
    =Ψ(α+β)-σ22ω2.   (3.21)

Assume that |ω|>1. Using that for |x|<1 it holds that cos(x)<1-x2/4, we can bound the first term of the integral in the following manner:

  {0}e(α+β)ycos(yω)ν(dy) 0<|y|<1/|ω|e(α+β)y(1-ω2y2/4)ν(dy)+|y|1/|ω|e(α+β)yν(dy)  
    {0}e(α+β)yν(dy)-|ω|2-λ4C(λ).   (3.22)

Inserting (3.1) back into (3.1), we get

  Re[Ψ(α+β-iω)] (α+β)(r-σ22)+σ22((α+β)2-ω2)  

Taking the previous considerations and integrating in with respect to ω, we obtain (3.5).

Finally, observing that C(λ)0 and bounding it by 0, the bound (3.16) is obtained by evaluating the integral. ∎

Remark 3.6.

In the case of call options, hypothesis H2 implies a dependence between the strip-width parameter a and the damping parameter α. We have that the damped payoff of the call option is in L1() if and only if α>1; hence, the appropriate choice of strip-width parameter is given by 0<a<α-1. A similar argument holds for the case of put options, for which the Fourier-transformed damped payoff is identical to the calls, with the distinction that α<0. In such a case, we require a<-α.

In the case of binary options with payoffs that have finite support (G(x)=𝟏[x-,x+](x)), we can set any a (ie, no damping is needed at all, and even if such damping is chosen, it has no effect on the appropriate choice of a).

Remark 3.7.

The bound we provide for the quadrature error is naturally positive and increasing in Δω. It decays to zero at a spectral rate as Δω decreases to zero.

3.2 Frequency truncation error

The frequency truncation error is given by


If a function c:(ω0,)(0,) satisfies

  |Re[e-i(k+12)Δωxf^α(τ,(k+12)Δω)]|c((k+12)Δω)   (3.23)

for every natural number k, then we have that


Further, if c is a nonincreasing concave integrable function, we get

  eαxπnΔωc(ω)dω.   (3.24)

When g^αL[ω0,) and either σ2>0 or C(λ)>0, then the function c in (3.23) can be chosen as

  c(ω)=g^αL[ω0,)eτΨ(α)exp(-τ(σ22ω2+|ω|2-λ4C(λ)𝟏|ω|>1)).   (3.25)

To prove that this function satisfies (3.23), we can use the same bound we found in the proof of Theorem 3.5, with β=0, to obtain


from which the result is straightforward.

3.3 Bound for the full error

In this section, we summarize the bounds obtained for the error under different assumptions and analyze their central properties.

In general, the bounds provided in this paper are of the form

  ¯=eαxπ(M¯e2πa/Δω-1+nΔωc(ω)dω),   (3.26)

where M¯ is an upper bound of Mα,a(τ,x) defined in (3.3), and c is nonincreasing, integrable and satisfies (3.23). Both M¯ and c may depend on the parameters of the model and the artificial parameters, but they are independent of Δω and n. Typically, one can remove the dependence of some of the parameters, simplifying the expressions but obtaining less tight bounds.

When analyzing the behavior of the bound, one can observe that the term correspondent with the quadrature error decreases to zero spectrally when Δω goes to zero. The second term goes to zero if nΔω diverges, but we are unable to determine the rate of convergence without further assumptions.

Once an expression for the error bound is obtained, the problem of how to choose the parameters of the numerical method in order to minimize the bound arises, assuming a constraint on the computational effort one is willing to use. The computational effort of the numerical method depends only on n. For this reason, we aim to find the parameters that minimize the bound for a fixed n. The following result shows that the bound obtained, as a function of Δω, has a unique local minimum, which is the global minimum.

Proposition 3.8.

Fix α, a, n and λ, and consider the bound ¯ as a function of Δω. There exists an optimal Δω*[ω0/n,) such that ¯ is decreasing in (ω0/n,Δω*) and increasing in (Δω*,); thus, a global minimum of ¯ is attained at Δω*.

Further, the optimal Δω is either the only point at which Δωp(nΔω,b)-c(nΔω), with p defined in (3.27), changes sign, or Δω=ω0/n if p(ω0,b)-c(ω0)>0.


Let us simplify the notation by calling y=nΔω, b=2πan and ~=πe-αx¯. We want to prove the existence of y*:y*ω0 such that ~(y) is decreasing for ω0<y<y* and increasing for y>y*. We have


The first term is differentiable with respect to y and goes to 0 if y0+. This allows us to express it as an integral of its derivative. We can then express ~(y) as


The first term on the right-hand side of the previous equation is constant. Now, we move on to proving that the integrand is increasing with y and is positive if y is large enough. This can be denoted by

  p(y,b)=bM¯eb/y(eb/y-1)2y2.   (3.27)

Taking into account that c is integrable, we can compute the limit of the integrand in , obtaining


Let us prove that p(y,b) is increasing with y for all b>0, which renders p(y,b)-c(y) also increasing with y. The derivative of p with respect to y is given by


in which the denominator and the first factor in the numerator are clearly positive. To prove that the remaining factor is also positive, observe that xex-2ex+x+2>0 if x>0. ∎

3.4 Explicit error bounds

In the case where either σ2>0 or C(λ)>0 for some λ(0,2), we can give an explicit version of (3.26). Substituting M with M~ (defined in Theorem 3.5) and c with the function given in (3.25), we obtain

  ¯=¯Q+¯F,   (3.28)


  ¯Q =c{-1,1}eαxecaxeτΨ(ca)|gα^(ca)|π(e(2πa/Δω)-1)  
      ×exp(-τ(σ22ω2+|ω|2-λ4C(λ)𝟏|ω|>1))dω,   (3.29)
  ¯F =eαxπgα^LeτΨ(α)nΔωexp(-τ(σ22ω2+|ω|2-λ4C(λ)𝟏|ω|>1))dω.   (3.30)

This reproduces the essential features of Feng and Linetsky (2008, Theorem 6.6); the bound (3.30) can be further improved by substituting gα^L with gα^L[nΔω,).

Remark 3.9.

Observe that the bound of both the quadrature and the cutoff error is given by a product of one factor that depends exclusively on the payoff and another factor that depends on the asset dynamic. This property makes it easy to evaluate the bound for a specific option under different dynamics of the asset price. In Section 4.4, we analyze the terms that depend on the payoff function for the particular case of call options.

Remark 3.10.

From (3.4), it is evident that the speed of the exponential convergence of the trapezoidal rule for analytic functions is dictated by the width of the strip in which the function being transformed is analytic. Thus, in the limit of small error tolerances, it is desirable to set a as large as possible to obtain optimal rates. However, non-asymptotic error tolerances are often practically relevant, and in these cases the trade-off between optimal rates and the constant term |g^α| becomes nontrivial. As an example, for the particular case of the Merton model, we have that any finite value of a will do. However, this improvement of the rate of spectral convergence is more than compensated for by the divergence in the constant term.

The integrals in (3.4) and (3.30) can, in some cases, be computed analytically, or bounded from above by a closed-form expression. Consider, for instance, dissipative models with finite jump intensity. These models are characterized by σ2>0 and C(λ)=0. Thus, the integrals can be expressed in terms of the cumulative normal distribution Φ:

  e-τ(σ2ω2/2)dω =2πτσ2,   (3.31)
  ςe-τ(σ2ω2/2)dω =2πτσ2(1-Φ(ςτσ2)).   (3.32)

Now we consider the case of pure-jump processes (ie, σ2=0) that satisfy the condition C(λ)>0 for some λ(0,2). In this case, the integrals are expressible in terms of the incomplete gamma function γ. First, let us define the auxiliary integral:


for a,b>0. Using this, the integrals become

  exp(-τ|ω|2-λ4C(λ)𝟏|ω|>1) =2(1+I(τC(λ)4,2-λ)),   (3.33)
  ςexp(-τ|ω|2-λ4C(λ)𝟏|ω|>1) ={I(τC(λ)4,2-λ)+1-ς,ς<1,ςI(τς2-λC(λ)4,2-λ),ς1.   (3.34)

An example of a process for which the previous analysis works is the CGMY model, presented in Carr et al (2002, 2003), for the regime Y>0.

Last, when both C(λ) and σ2 are positive, the integrals in (3.4) and (3.30) can be bounded by a simpler expression. Consider the two following auxiliary bounds for the same integral, where ς1:

  ςexp(-τ(σ22ω2+|ω|2-λ4C(λ)))dω e-τ(σ2/2)ς2ςe-τ(|ω|2-λ/4)C(λ)dω  
    =ςe-τ(σ2/2)ς2I(τς2-λC(λ)4,2-λ),   (3.35)
  ςexp(-τ(σ22ω2+|ω|2-λ4C(λ)))dω e-τ(ς2-λ/4)C(λ)ςe-τ(σ2/2)ω2dω  
    =2πτσ2e-τ(ς2-λ/4)C(λ)(1-Φ(ςτσ2)).   (3.36)

We have that b(ς), defined as the minimum of the right-hand sides of the two previous equations, is

  b(ς) =min{ςe-τ(σ2/2)ς2I(τς2-λC(λ)4,2-λ),  

a bound for the integral. Bearing this in mind, we have

  exp(-τ(σ22ω2+|ω|2-λ4C(λ)𝟏|ω|>1))dω2Φ(τσ2)-1+2b(1)   (3.37)
  ςexp(-τ(σ22ω2+|ω|2-λ4C(λ)𝟏|ω|>1))dωb(ς),   (3.38)

provided that ς1.

4 Computation and minimization of the bound

In this section, we present numerical examples on the bound presented in the previous section, using practical models known from the literature. We gauge the tightness of the bound compared with the true error using both dissipative and pure-jump processes. We also demonstrate the feasibility of using the expression of the bound as a tool for choosing numerical parameters for the Fourier inversion.

4.1 Call option in VG model

The VG model provides a test case to evaluate the bound in the pure-jump setting. We note that of the two numerical examples presented, it is the less regular model, in the sense that σ2=0 and C(λ)=0 for 0<λ<2, indicating that Theorem 3.5 in particular is not applicable.

The Lévy measure of the VG model is given by


and the corresponding characteristic function is given by Madan et al (1998, Equation (7)):

  φτω =(1-iθχω+σ2χ2)-τ/χ,  
  K =χ-1,  
  η- =(θ2χ24+σ2ν2-θχ2)-1,  
  η+ =(θ2χ24+σ2ν2+θχ2)-1.  

By Proposition 3.3, we get that

  a<min{η--α,η++α},   (4.1)

which, combined with the requirement that gαL1() (see Remark 3.6), implies

  a <min{η+-α,η-+α,α-1},  
  a <min{η+-α,η-+α,-α}   (4.2)

for calls and puts, respectively. We note that an evaluation of the integral in (2.13) is also possible for α(0,1) and α<0. In fact, there is a correspondence between shifts in the integration contour and the put–call parity. Integrals with α<0 give rise to put option prices instead of calls. For an extended discussion of this, we refer the reader to Lee (2004) or Boyarchenko and Levendorskii (2011), in which conformal deformation of the integration contour is exploited in order to achieve improved numerical accuracy.

In Lee (2004) and in our calculations, the parameters equal η+=39.7840, η-=20.2648 and K=5.9311.

Table 1 presents the specific parameters and compares the bound for the VG model with the results obtained by Lee (2004). Based on Table 1, we note that for the VG model presented in Madan et al (1998), we can achieve comparable or better error bounds when compared with the study by Lee (2004).

Table 1: The error bound for European call/put options in the VG model for select examples. [Reference result ¯* from Lee (2004).]
    80 90 100 110 120
12τ=1, α -16.9 -13.8 21.6 29.10 36.3
N=32 a 3.33 6.45 18.1 9.77 3.52
  ωmax 229 229 363 363 424
  ¯ 3.35×10-4 0.00334 0.00562 3.97×10-4 7.33×10-6
  ¯* 6×10-4 0.0032 0.0058 6×10-4 1×10-4
12τ=4, α -13.8 -13.8 22.1 23.7 29.10
N=8 a 6.11 6.11 17.9 15.2 8.75
  ωmax 62.4 42.4 84.9 126 126
  ¯ 3.99×10-4 0.00312 0.00398 3.57×10-4 1.33×10-5
  ¯* 1.3×10-3 0.0057 0.0055 9×10-4 1×10-4

To evaluate the bound, we perform the integration of (3.3) and (3.24) by relying on the Clenshaw–Curtis quadrature method provided in the SciPy package. To supplement Table 1 for a wide range of n, we present the magnitude of the bound compared with the true error in Figure 1.

 The true error and the error bound for evaluating at-the-money options for the VG model test case
Figure 1: The true error and the error bound for evaluating at-the-money options for the VG model test case.

In Figure 1, we see that the choice of numerical parameters for the Fourier inversion has a strong influence on the error of the numerical method. One does not in general have access to the true solution. Thus, the parameters need to be optimized with respect to the bound. Recall that =(α,Δω,a,n) and ¯=¯(α,Δω,n) denote the true and estimated errors, respectively. Keeping the number of quadrature points n fixed, we let (α1,Δω1,a1) and (α2,Δω2) denote the minimizers of the estimated and true errors, respectively:

  (α1,Δω1,a1) =arginf¯,   (4.3)
  (α2,Δω2) =arginf.   (4.4)

Further, we let 1 and 2 denote the true error as a function of the parameters minimizing the estimated error and true error, respectively:

  1 =(α1,Δω1),   (4.5)
  2 =(α2,Δω2).   (4.6)

In Figure 1, we see that the true error increases by approximately an order of magnitude when optimizing to the bound instead of to the true error, which translates into a twofold difference in the number of quadrature points needed for a given tolerance. The difference between 1 and the bound is approximately another order of magnitude, which necessitates another twofold number of quadrature points compared with the theoretical minimum.

The true error ℰ for the two VG test cases presented in Table 1
Figure 2: The true error for the two VG test cases presented in Table 1, and the bound-minimizing configurations (white circle) (α2,Δω2) for the examples. (a) n=8. (b) n=32.

In Figure 2, we present the true error of the Fourier method for the two test cases in Table 1.11The reference value for computing the true error was obtained using the numerical methods with n and Δω, such that the level of accuracy is of the order 10-10. We note that while minimizing error bounds will produce suboptimal results, the numerical parameters that minimize the bound are a good approximation of the true optimal parameters. This, of course, is a consequence of the error bound having qualitatively similar behavior to the true error, especially as one gets further away from the true optimal parameters.

Remark 4.1.

In practice, the Hardy norm in coefficient M reduces to evaluating an L1 norm along the two boundaries of the strip of width 2a. We find that, for practical purposes, the performance of the Clenshaw–Curtis quadrature of the QUADPACK library (provided by the SciPy library) is more than adequate, enabling the evaluation of the bound in a fraction of a second.

For example, the evaluations of the bounds in Table 1 take only around 0.3 seconds on a mid 2014 Macbook Pro equipped with a 2.6 GHz Intel Core i5 processor; this is without attempting to optimize or parallelize the implementation, and while checking for input sanity factors such as the evaluation of the characteristic function in a domain that is a subset in the permitted strip.

We believe that through optimizing routines, skipping sanity checks for inputs and using lower-level computation routines this can be optimized even further, guaranteeing a fast performance even when numerous evaluations are needed.

Remark 4.2.

Like many other authors, we note the exceptional guaranteed accuracy of the FT method with only dozens of quadrature points. This is partially a result of the regularity of the European option price. Numerous Fourier-based methods have been developed for pricing path-dependent options. One might, for the sake of generality, be tempted to use these methods for European options as well, correcting for the lack of early exercise opportunities. This can be done, certainly; but due to weakened regularity, the required number of quadrature points is easily in the thousands, even when no rigorous bound for the error is required.

We raise one point of comparison, the European option pricing example in Jackson et al (2008, Table 2), which indicates a number of quadrature points for pricing the option in the range of thousands. With the method introduced, to guarantee ¯10-3, even with no optimization, n=64 turns out to be sufficient.

4.2 Call options under Kou dynamics

To contrast with the pure-jump process presented above, we also test the performance of the bound for the Kou model. We present the relevant results in Table 2. This model differs from the first example in terms of its being dissipative as well as in terms of regularity, in the sense that the maximal width of the domain Aa is, in the case at hand, considerably narrower. The Lévy measure in the Kou model is given by


with p+q=1. For the characterization given in Toivanen (2007), the values are set as


From the expression of the characteristic exponent (see Kou and Wang 2004)


it is straightforward to see that


This range is considerably narrower than that considered earlier. When transforming the option prices in strike space, the relevant expressions for option prices and the error bounds contain a factor exponential in k. The practical implication of this is that, for deep out-of-the-money calls, it is often beneficial to exploit the put–call parity and compute deep in-the-money calls. However, in the case at hand, the strip width does not permit such a luxury. As a consequence, the parameters that minimize the bound are near-identical over a wide range of moneyness, suggesting that we use the FFT algorithm to evaluate the option prices at once for a range of strikes.

Table 2: Numerical performance of the bound for the Kou model, with the test case in Toivanen (2007) (see also d’Halluin et al 2005), with the number of quadrature points set to n=32. [The point of comparison ¯* refers to the corresponding bound computed with the method described in Lee (2004, Chapters 6.1–6.4). In the ¯, the cutoff error has been evaluated using a computationally more intensive Clenshaw–Curtis quadrature instead of an asymptotic argument with an exponentially decaying upper bound for the option price.]
  80 90 100 110 120
¯ 2.67×10-4 3.49×10-4 4.43×10-4 5.52×10-4 6.77×10-4
α -1.57 -1.57 -1.57 -1.57 -1.57
ωmax 22.9 22.8 22.6 22.5 22.4
¯* 0.34 0.26 0.21 0.17 0.13
¯ 6.87×10-4 1.90×10-3 2.82×10-3 2.72×10-3 2.29×10-3

4.3 Binary option in the Merton model

For the particular case of the Merton model, the Lévy measure is given by


and the characteristic exponent is correspondingly given by


We may employ a fast, semi-closed-form evaluation of the relevant integrals instead of resorting to quadrature methods. We choose the Merton model as an example of bounding the error of the numerical method for such a model. The parameters are adopted from the estimated parameters for the Standard & Poor’s (S&P) 500 index from Andersen and Andreasen (2000):


In Figure 3, we present the bound and true error for the Merton model to demonstrate the bound on another dissipative model. The option presented is a binary option with finite support on [95,105]; no damping was needed or used. We note that, as in the case of the pure-jump module presented above, our bound reproduces the qualitative behavior of the true error. The configuration that results from optimizing the bound is a good approximation of the true error. Such behavior is consistent across the range of n of the most practical relevance.

The true error ℰ1 and the bound ℰ¯ for the dissipative Merton model for a range of quadrature points n
Figure 3: The true error 1 and the bound ¯ for the dissipative Merton model for a range of quadrature points n, along with the bound-minimizing configurations, contrasted with the true error.

4.4 Call options

In Section 3.4, explicit expressions to bound are provided. To evaluate these bounds, it is necessary to compute gα^L and gα^LAa. According to Remark 3.9, once we compute these values, we could use them for any model, provided that they satisfy the conditions considered there.

The payoff of perhaps the most practical relevance is that of a call option. Consider g, defined by


for which the selection of a damping parameter α>0 is necessary to have the damped payoff in L1() as well as to ensure the existence of a Fourier transformation. In this case, we have

  g^α(ω) =S0exp((1-α+iω)x)-exp(k+(iω-α)x)dx   (4.7)
    =S0exp((1-α+iω)k)(1+iω-α)(iω-α)   (4.8)


  |gα^(ω)|2=S02e2(1-α)k(α2+ω2)((1-α)2+ω2).   (4.9)

It is easy to see that the previous expression decreases as |ω| increases. This yields

  g^αL=|g^α(0)|=S0e(1-α)kα2-α   (4.10)
  gα^L[ς,)=|g^α(ς)|.   (4.11)

The maximization of |g^α| in the strip Aa of the complex plane is more subtle. Denoting g^α(η,ρ)=g^α(η+iρ), we look for critical points that satisfy η|g^α|=0. This gives

  4η3+2η(4ρα+2α2-2ρ-2α+ρ2+1)=0.   (4.12)

For fixed ρ, |g^α| has a vanishing derivative with respect to η at a maximum of three points. Of the three roots of the derivative, only the one characterized by η=0 is a local maximum; this gives us that, for call options,

  g^αLAa=supρ[-a,a]|g^α(0,ρ)|.   (4.13)

Now, observe that |g^α(0,ρ)| is a differentiable real function of ρ, whose derivative is given by the following polynomial of second degree:

  p(ρ)k(ρ+α-2ρα-α2-ρ2)-2α-2ρ+1.   (4.14)

We conclude that

  g^αLAa=maxρB{|g^α(0,ρ)|},   (4.15)

where B is the set of no more than four elements consisting of a, -a and the real roots of p that fall in (-a,a).

Remark 4.3.

So far, we have assumed the number of quadrature points n to be constant. In real-life applications, however, this is often not the case. Typically, the user will choose a minimal n that is sufficient to guarantee an error that lies within a predefined error tolerance.

In such a case, we propose the following very simplistic scheme for optimizing numerical parameters and choosing the appropriate n to satisfy an error smaller than ϵ.

  1. (1)

    Select n=n0 and optimize to find the relevant configuration.

  2. (2)

    See if 𝒬+<ϵ; if not, increase n by choosing it from a predetermined increasing sequence n=nj, and repeat the procedure.

Especially when using FFT algorithms to evaluate Fourier transforms, we propose nj=2jn0. We further note that, typically, the optimal configuration for the optimizing configuration for nj+1 quadrature points does not differ too dramatically from the configuration that optimizes bounds for nj.

5 Conclusion

We have presented a decomposition of the error committed in the numerical evaluation of the inverse Fourier transform, which is needed in asset pricing for exponential Lévy models, into truncation and quadrature errors. For a wide class of exponential Lévy models, we have presented an L bound for the error.

This error bound differs from the earlier work of Lee (2004) in the sense that it does not rely on the asymptotic behavior of the option payoff at extreme strikes or option prices. This enables the pricing of a wide variety of nonstandard payoff functions, such as those in Suh and Zapatero (2008). The bound, however, does not take into account path-dependent options. We argue that the methods which allow the evaluation of American, Bermudan or knockoff options are considerably more cumbersome and produce significantly larger errors; so, in implementations where performance is important, such as calibration, using American option pricing methods for European options is not justified.

The bound also provides a general framework in which the truncation error is evaluated using a quadrature method; this remains invariant, regardless of the asymptotic behavior of the option price function. The structure of the bound allows for a modular implementation that decomposes the error components arising from the dynamics of the system and the payoff into a product form for a large class of models, including all dissipative models. In select examples, we also demonstrate performances that are comparable or superior to the relevant points of comparison.

We have focused on the minimization of the bound as a proxy for minimizing numerical error. By doing this, one obtains, for a given parameterization of a model, a rigorous L bound for the error committed in solving the European option price. We have shown that the bound reproduces the qualitative behavior of the actual error. This supports the argument for selecting numerical parameters in a way that minimizes the bound, providing evidence that this selection will, besides guaranteeing numerical precision, be close to the actual minimizing configuration that is not often achievable at an acceptable computational cost.

The bound can be used in the primitive setting of establishing a strict error bound for the numerical estimation of option prices for a given set of physical and numerical parameters, or as a part of a numerical scheme, whereby the end user wishes to estimate an option price either on a single point or in a domain up to a predetermined error tolerance.

In future, the error bounds presented could be used in efforts requiring multiple evaluations of Fourier transformations. Examples of such applications include multi-dimensional Fourier transformations, possibly in sparse tensor grids, as well as time-stepping algorithms for American and Bermudan options. Such applications are sensitive toward the error bound being used, as any numerical scheme will need to be run multiple times, either in high dimension or for multiple time steps (or both).

Declaration of interest

The authors report no conflicts of interest. The authors alone are responsible for the content and writing of the paper.


Häppölä and Crocce are members and Tempone is the director of the KAUST Strategic Research Initiative, Uncertainty Quantification Center. We wish to thank an anonymous referee for their critique, which significantly improved the quality of the manuscript.


Almendral, A., and Oosterlee, C. W. (2005). Numerical valuation of options with jumps in the underlying. Applied Numerical Mathematics 53(1), 1–18 (

Andersen, L., and Andreasen, J. (2000). Jump-diffusion processes: volatility smile fitting and numerical methods for option pricing. Review of Derivatives Research 4(3), 231–262 (

Applebaum, D. (2004). Lévy Processes and Stochastic Calculus. Cambridge Studies in Advanced Mathematics, Volume 93. Cambridge University Press.

Barndorff-Nielsen, O. E. (1997). Normal inverse Gaussian distributions and stochastic volatility modelling. Scandinavian Journal of Statistics 24(1), 1–13 (

Black, F., and Scholes, M. (1973). The pricing of options and corporate liabilities. Journal of Political Economy 81(3), 637–654.

Boyarchenko, S., and Levendorskii, S. (2002). Barrier options and touch-and-out options under regular Lévy processes of exponential type. Annals of Applied Probability 12(4), 1261–1298.

Boyarchenko, S., and Levendorskii, S. (2011). New efficient versions of Fourier transform method in applications to option pricing. Working Paper 1846633, Social Science Research Network.

Briani, M., Chioma, C. L., and Natalini, R. (2004). Convergence of numerical schemes for viscosity solutions to integro-differential degenerate parabolic problems arising in financial theory. Numerische Mathematik 98(4), 607–646.

Carr, P., Geman, H., Madan, D. B., and Yor, M. (2002). The fine structure of asset returns: an empirical investigation. Journal of Business 75(2), 305–333 (

Carr, P., Geman, H., Madan, D. B., and Yor, M. (2003). Stochastic volatility for Lévy processes. Mathematical Finance 13(3), 345–382 (

Cont, R., and Tankov, P. (2004). Financial Modelling with Jump Processes. Chapman & Hall/CRC Financial Mathematics Series. Chapman & Hall/CRC, Boca Raton, FL.

De Innocentis, M., and Levendorskii, S. (2014). Pricing discrete barrier options and credit default swaps under Lévy processes. Quantitative Finance 14(8), 1337–1365 (

d’Halluin, Y., Forsyth, P. A., and Labahn, G. (2004). A penalty method for American options with jump diffusion processes. Numerische Mathematik 97(2), 321–352 (

d’Halluin, Y., Forsyth, P. A., and Vetzal, K. R. (2005). Robust numerical methods for contingent claims under jump diffusion processes. IMA Journal of Numerical Analysis 25(1), 87–112 (

Dotsis, G., Psychoyios, D., and Skiadopoulos, G. (2007). An empirical comparison of continuous-time models of implied volatility indices. Journal of Banking and Finance 31(12), 3584–3603 (

Dupire, B. (1997). Pricing and hedging with smiles. In Mathematics of Derivative Securities, Dempster, S. R., and Pliska, M. A. H. (eds). Cambridge University Press.

Eberlein, E. (2001). Application of generalized hyperbolic Lévy motions to finance. In Lévy Processes, pp. 319–336. Springer (

Feng, L., and Linetsky, V. (2008). Pricing discretely monitored barrier options and defaultable bonds in Levy process models: a fast Hilbert transform approach. Mathematical Finance 18(3), 337–384 (

Hurd, T. R., and Zhou, Z. (2010). A Fourier transform method for spread option pricing. SIAM Journal on Financial Mathematics 1(1), 142–157 (

Jackson, K. R., Jaimungal, S., and Surkov, V. (2008). Fourier space time-stepping for option pricing with Lévy models. The Journal of Computational Finance 12(2), 1–29 (

Kiessling, J., and Tempone, R. (2011). Diffusion approximation of Lévy processes with a view towards finance. Monte Carlo Methods and Applications 17(1), 11–45 (

Kou, S. G. (2002). A jump-diffusion model for option pricing. Management Science 48(8), 1086–1101 (

Kou, S. G., and Wang, H. (2004). Option pricing under a double exponential jump diffusion model. Management Science 50(9), 1178–1192 (

Lee, R. W. (2004). Option pricing by transform methods: extensions, unification and error control. The Journal of Computational Finance 7(3), 51–86 (

Lord, R., Fang, F., Bervoets, F., and Oosterlee, C. W. (2008). A fast and accurate FFT-based method for pricing early-exercise options under Lévy processes. SIAM Journal on Scientific Computing 30(4), 1678–1705 (

Madan, D. B., and Seneta, E. (1990). The variance gamma (VG) model for share market returns. Journal of Business 63(4), 511–524.

Madan, D. B., Carr, P. P., and Chang, E. C. (1998). The variance gamma process and option pricing. European Finance Review 2(1), 79–105 (

Mattner, L. (2001). Complex differentiation under the integral. Nieuw Archief voor Wiskunde 2, 32–35.

Merton, R. C. (1976). Option pricing when underlying stock returns are discontinuous. Journal of Financial Economics 3(12), 125–144 (

Nualart, D., and Schoutens, W. (2001). Backward stochastic differential equations and Feynman–Kac formula for Lévy processes, with applications in finance. Bernoulli 7(5), 761–776 (

Raible, S. (2000). Lévy processes in finance: theory, numerics, and empirical facts. PhD Thesis, Universität Freiburg im Breisgau.

Reed, M., and Simon, B. (1975). Methods of Modern Mathematical Physics, Volume 2: Fourier Analysis, Self-Adjointness. Academic Press.

Rydberg, T. H. (1997). The normal inverse Gaussian Lévy process: simulation and approximation. Communications in Statistics: Stochastic Models 13(4), 887–910 (

Schmelzle, M. (2010). Option pricing formulae using Fourier transform: theory and application. Unpublished manuscript.

Stenger, F. (1993). Numerical Methods Based on Sinc and Analytic Functions. Springer.

Suh, S., and Zapatero, F. (2008). A class of quadratic options for exchange rate stabilization. Journal of Economic Dynamics and Control 32(11), 3478–3501 (

Tankov, P. (2004). Financial Modelling with Jump Processes. CRC Press.

Toivanen, J. (2007). Numerical valuation of options under Kou’s model. PAMM 7(1), 1024001–1024002 (