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 is modeled by the stochastic process , defined by , where is assumed to be a Lévy process whose jump measure satisfies
Assuming the risk-neutral dynamic for , the price at time of a European option with payoff and maturity time is given by
where is the short rate that we assume to be constant and is the time to maturity. Extensions to nonconstant deterministic short rates are straightforward.
The infinitesimal generator of a Lévy process is given by (see Applebaum 2004)
where is the characteristic triple of the Lévy process. The risk-neutral assumption on implies
and fixes the drift term (see Kiessling and Tempone 2011) of the Lévy process to
Thus, the infinitesimal generator of may be written under the risk-neutral assumption as
Consider as the reward function in log prices (ie, defined by ). Now, take to be defined as
Then, solves the following PIDE:
Observe that and are related by
Consider a damped version of defined by ; we see that .
There are different conventions for the Fourier transform. Here, we consider the operator such that
defined for functions for which the previous integral is convergent. We also use as a shorthand notation of . To recover the original function , we define the inverse Fourier transform as
We have that .
Applying to , we get . Observe also that the Fourier transform applied to gives , where is the characteristic exponent of the process , which satisfies . The explicit expression for is
From the previous considerations, it can be concluded that
Now, , so satisfies the following ODE:
Solving the previous ODE explicitly, we obtain
Observe that the first factor on the right-hand side of the above equation is (ie, ), where denotes the characteristic function of the random variable :
Now, we employ the inverse Fourier transformation to obtain the value function:
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:
Bounding and consequently minimizing the error in the approximation of by
is the main focus of this paper and will be addressed in the following section.
Although we are mainly concerned with option pricing when the payoff function can be damped in order to guarantee regularity in the 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
so the Delta and Gamma of the option equal
Because the expressions involve partial derivatives with respect to only , the results in this work are applicable to the computation of and through a modification of the payoff function:
When the Fourier space payoff function manifests exponential decay, the introduction of a coefficient that is polynomial in does not change the regularity of 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 ’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 simultaneously
The fast Fourier transform (FFT) algorithm provides an efficient way of computing (2.15) for an equidistantly spaced mesh of values for 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 , treat as a constant and write
Using this convention, the time dependence is given by
contrasted with the -space solution
We note that for a call option payoff to be in , we demand that in (2.23) is positive. By omitting the exponential factors that contain the and dependence in (2.23) and (2.24), respectively, one can get from (2.23) to (2.24) using the mapping . Thanks to this, much of the analysis regarding the -space transformation generalizes in a straightforward manner to the -space transform.
3 Error bound
The aim of this section is to compute a bound of the error when approximating the option price by , which is defined in (2.15). Considering
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
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 ;
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 , with , the strip of width 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 is a diffusive process, or there are “enough small jumps”.
Assume that, for ,
the characteristic function of the random variable has an analytic extension to the set
the Fourier transform of is analytic in the strip ,
there exists a continuous function such that for all and for all .
Then, the quadrature error is bounded by
where is given by
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), is the family of functions that are analytic in , such that
Lemma 3.2 (Stenger 1993, Theorem 3.2.1).
Let ; then, define
Proof of Theorem 3.1.
First, observe that H1 and H2 imply that the function is analytic in . H3 allows us to use dominated convergence theorem to prove that is finite and coincides with . 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.
If , and are such that
and then H1 in Theorem 3.1 is fulfilled.
Denoting by the characteristic function of , we want to prove that is analytic in . Considering that , the only nontrivial part of the proof is to verify that
is analytic in , where 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 , ensures the analyticity of , provided that satisfies the following: is -measurable for all , is holomorphic for all and is locally bounded. In our case, we consider the measure space to be , with the Borel -algebra and the Lebesgue measure, and . It is clear that is Borel measurable and is holomorphic. It remains for us to verify that
is locally bounded. To this end, we assume that (and, since , ) and split the integration domain into and to prove that both integrals are uniformly bounded.
Regarding the integral in , we observe that
For the integral in , observe that, denoting , we have for every , for every and for , , . From these observations, we get that the Maclaurin polynomial of degree of is null for every . We can bound by the remainder term, which, in our region of interest, is bounded by ; thus, we obtain
which is finite by the hypothesis on . This finishes the proof. ∎
If, for all , the function is in , 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 or there exists such that 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 , define as
Observe that , and, by our assumptions on the jump measure , is finite. Further, if is such that
then . To see this, note that (3.14) implies the existence of , such that
If , observe that
where, for the first inequality, it was taken into account that and that the integral is increasing with . By combining the two previous infima and considering , we get that .
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, for all .
Assume that and are such that (3.9) holds, and either or for some . Then, the quadrature error is bounded by
Further, if , we have
Considering , defined by
we have that
However, for ,
For the factor involving the characteristic exponent, we have
Now, observe that
If , we bound by 1, getting
Assume that . Using that for it holds that , we can bound the first term of the integral in the following manner:
Taking the previous considerations and integrating in with respect to , we obtain (3.5).
Finally, observing that and bounding it by 0, the bound (3.16) is obtained by evaluating the integral. ∎
In the case of call options, hypothesis H2 implies a dependence between the strip-width parameter and the damping parameter . We have that the damped payoff of the call option is in if and only if ; hence, the appropriate choice of strip-width parameter is given by . 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 . In such a case, we require .
In the case of binary options with payoffs that have finite support (), we can set any (ie, no damping is needed at all, and even if such damping is chosen, it has no effect on the appropriate choice of ).
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 satisfies
for every natural number , then we have that
Further, if is a nonincreasing concave integrable function, we get
When and either or , then the function in (3.23) can be chosen as
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
where is an upper bound of defined in (3.3), and is nonincreasing, integrable and satisfies (3.23). Both and may depend on the parameters of the model and the artificial parameters, but they are independent of and . 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 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 . For this reason, we aim to find the parameters that minimize the bound for a fixed . The following result shows that the bound obtained, as a function of , has a unique local minimum, which is the global minimum.
Fix , , and , and consider the bound as a function of . There exists an optimal such that is decreasing in and increasing in ; thus, a global minimum of is attained at .
Further, the optimal is either the only point at which , with defined in (3.27), changes sign, or if .
Let us simplify the notation by calling , and . We want to prove the existence of such that is decreasing for and increasing for . We have
The first term is differentiable with respect to and goes to 0 if . This allows us to express it as an integral of its derivative. We can then express 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 and is positive if is large enough. This can be denoted by
Taking into account that is integrable, we can compute the limit of the integrand in , obtaining
Let us prove that is increasing with for all , which renders also increasing with . The derivative of with respect to 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 if . ∎
3.4 Explicit error bounds
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.
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 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 becomes nontrivial. As an example, for the particular case of the Merton model, we have that any finite value of 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 and . Thus, the integrals can be expressed in terms of the cumulative normal distribution :
Now we consider the case of pure-jump processes (ie, ) that satisfy the condition for some . In this case, the integrals are expressible in terms of the incomplete gamma function . First, let us define the auxiliary integral:
for . Using this, the integrals become
We have that , defined as the minimum of the right-hand sides of the two previous equations, is
a bound for the integral. Bearing this in mind, we have
provided that .
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 and for , 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)):
By Proposition 3.3, we get that
which, combined with the requirement that (see Remark 3.6), implies
for calls and puts, respectively. We note that an evaluation of the integral in (2.13) is also possible for and . In fact, there is a correspondence between shifts in the integration contour and the put–call parity. Integrals with 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 , and .
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).
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 , we present the magnitude of the bound compared with the true error in Figure 1.
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 and denote the true and estimated errors, respectively. Keeping the number of quadrature points fixed, we let and denote the minimizers of the estimated and true errors, respectively:
Further, we let and denote the true error as a function of the parameters minimizing the estimated error and true error, respectively:
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 and the bound is approximately another order of magnitude, which necessitates another twofold number of quadrature points compared with the theoretical minimum.
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 and , such that the level of accuracy is of the order . 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.
In practice, the Hardy norm in coefficient reduces to evaluating an norm along the two boundaries of the strip of width . 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 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.
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 , even with no optimization, 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 is, in the case at hand, considerably narrower. The Lévy measure in the Kou model is given by
with . 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 . 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.
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 ; 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 of the most practical relevance.
4.4 Call options
In Section 3.4, explicit expressions to bound are provided. To evaluate these bounds, it is necessary to compute and . 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 , defined by
for which the selection of a damping parameter is necessary to have the damped payoff in as well as to ensure the existence of a Fourier transformation. In this case, we have
It is easy to see that the previous expression decreases as increases. This yields
The maximization of in the strip of the complex plane is more subtle. Denoting , we look for critical points that satisfy . This gives
For fixed , 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 is a local maximum; this gives us that, for call options,
Now, observe that is a differentiable real function of , whose derivative is given by the following polynomial of second degree:
We conclude that
where is the set of no more than four elements consisting of , and the real roots of that fall in .
So far, we have assumed the number of quadrature points to be constant. In real-life applications, however, this is often not the case. Typically, the user will choose a minimal 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 to satisfy an error smaller than .
Select and optimize to find the relevant configuration.
See if ; if not, increase by choosing it from a predetermined increasing sequence , and repeat the procedure.
Especially when using FFT algorithms to evaluate Fourier transforms, we propose . We further note that, typically, the optimal configuration for the optimizing configuration for quadrature points does not differ too dramatically from the configuration that optimizes bounds for .
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 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 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 (http://doi.org/fw2zjg).
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 (http://doi.org/bbqjx2).
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 (http://doi.org/fjwqd4).
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 (http://doi.org/fqvqsv).
Carr, P., Geman, H., Madan, D. B., and Yor, M. (2003). Stochastic volatility for Lévy processes. Mathematical Finance 13(3), 345–382 (http://doi.org/d8pqm7).
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 (http://doi.org/btnk).
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 (http://doi.org/d5sn5v).
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 (http://doi.org/dkpm7g).
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 (http://doi.org/d288t3).
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 (http://doi.org/dgp3sh).
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 (http://doi.org/fwbmkb).
Hurd, T. R., and Zhou, Z. (2010). A Fourier transform method for spread option pricing. SIAM Journal on Financial Mathematics 1(1), 142–157 (http://doi.org/dqvtxv).
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 (http://doi.org/btnm).
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 (http://doi.org/fbxt9z).
Kou, S. G. (2002). A jump-diffusion model for option pricing. Management Science 48(8), 1086–1101 (http://doi.org/fg58sf).
Kou, S. G., and Wang, H. (2004). Option pricing under a double exponential jump diffusion model. Management Science 50(9), 1178–1192 (http://doi.org/bqk52m).
Lee, R. W. (2004). Option pricing by transform methods: extensions, unification and error control. The Journal of Computational Finance 7(3), 51–86 (http://doi.org/btnn).
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 (http://doi.org/bd4r4h).
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 (http://doi.org/d4hkvp).
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 (http://doi.org/d659w5).
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 (http://doi.org/cs6qm4).
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 (http://doi.org/fmsknc).
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 (http://doi.org/bd5qf3).
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 (http://doi.org/cdssm5).