Magic Points in Finance: Empirical Integration for Parametric Option Pricing

We propose an oﬄine-online procedure for Fourier transform based option pricing. The method supports the acceleration of such essential tasks of mathematical ﬁnance as model calibration, real-time pricing, and, more generally, risk assessment and parameter risk estimation. We adapt the empirical magic point interpolation method of Barrault, Nguyen, Maday and Patera (2004) to parametric Fourier pricing. In the oﬄine phase, a quadrature rule is tailored to the family of integrands of the parametric pricing problem. In the online phase, the quadrature rule then yields fast and accurate approximations of the option prices. Under analyticity assumptions the pricing error decays exponentially. Numerical experiments in one dimension conﬁrm our theoretical ﬁndings and show a signiﬁcant gain in eﬃciency, even for examples beyond the scope of the theoretical results.


Introduction
Most of the option pricing methods based on Fourier transforms aim at the evaluation of individual option prices.For this use case, existing pricing tools have achieved impressive performance.For real-time applications and those involving repeated evaluations particularly fast run-times are crucial.Therefore, Fast Fourier transforms (FFT) have become highly popular to reduce computational complexity when prices are required simultaneously for a large set of different strikes, following the seminal works of Carr and Madan (1999) and Raible (2000).See also the monograph Boyarchenko and Levendorski ȋ (2002b).In this paper we shift the focus from the pricing problem for one or several strikes to the full parametric option pricing problem, considering all parameters such as strike, maturity, and model parameters.
For a given model and option type, various applications require the evaluation of Fourier pricing routines repeatedly for different parameter constellations.We mention three of those applications: Firstly, during calibration of financial models with Fourier methods, the optimization relies on multiple evaluations of a Fourier integral for varying model and option parameters.Secondly, each (intra-)day recalibration leads to new model parameters for which several option prices and their sensitivities have to be computed in real-time.Thirdly, a variety of other relevant financial quantities that need to be repeatedly evaluated for different parameter constellations can be expressed via Fourier transforms analogously to the Fourier representation of option prices.As one example we refer to Armenti et al. (2015) who propose a method to quantify risk allocation.At the heart of their algorithm lies an optimization routine that requires the repeated evaluation of parametric expectations that can be expressed as Fourier integrals.They give numerical examples, where standard Fourier methods turn out to be too slow for practical use.
In all of these cases, Fourier pricing routines that are fast and accurate for a whole set of relevant parameter constellations are required.Often, numerical experiments on the performance of pricing routines are presented for some fixed parameter constellations, see for instance von Sydow et al. (2015).Such comparisons clearly present highly valuable insight in different methods.It has, however, to be recognized that a method that performs highly efficient for some specific parameters needs not to be as efficient for the whole set of relevant parameters.This phenomenon has already been observed and covered by experiments, see de Innocentis and Levendorski ȋ (2014).In this article we present a new Fourier pricing method that is designed to be highly efficient for a whole range of parameters of interest.
We look at this Parametric Option Pricing (POP) through the lense of offlineonline schemes.The main idea is to achieve fast and accurate real-time pricing, founded on a pre-computation step.The architecture of such methods decomposes into two separate phases.In the so-called offline phase, the algorithm parses the complexity of the parametric pricing problem and extracts a structure bearing all of the important information on the whole problem as such.This is the computationally demanding part.Ideally, the offline phase is only performed once for a selected model class and option type.Thus the offline phase can be seen as part of the implementation of the pricing method.Intuitively, it represents the learning phase of the algorithm.In the so-called online phase, real-time pricing is performed.This second part benefits from the pre-computation and thus yields the desired fast and accurate pricing results.
Two types of offline-online schemes have been proposed for POP in the literature: In Sachs and Schu (2010), Cont et al. (2011), Pironneau (2011) and Haasdonk et al. (2012) offline-online decomposition has been adopted to solve parametric partial differential equations for option pricing.In contrast, Gaß et al. (2015) present polynomial interpolation of the option price in the parameter space.
In this paper, we tailor an offline-online scheme to Fourier pricing.Over the last fifteen years, Fourier based option pricing has been applied successfully in both academia and practice.Pioneered by Stein and Stein (1991) and Heston (1993) for Brownian models, researchers have exploited the flexibility of the approach to create fast and efficient pricing algorithms for a large class of models and option types.(Fast) Fourier pricing of European options in Lévy and the large class of affine models has first been developed by Carr and Madan (1999), Raible (2000) and Duffie et al. (2000).There is also a large and growing literature on Fourier methods to price path dependent options, see e.g.Boyarchenko and Levendorski ȋ (2002c), Feng and Linetsky (2008), Kudryavtsev and Levendorski ȋ (2009), Zhylyevskyy (2010), Fang and Oosterlee (2011), Levendorskiy and Xie (2012), Feng and Lin (2013) and Zeng and Kwok (2014), and see Eberlein et al. (2010) for a general framework and analysis.
The main contributions of this article are to -propose an empirical quadrature rule to efficiently evaluate Fourier integrals for option pricing, -find exponential convergence of the pricing error when the option price satisfies strict analyticity assumptions, -empirically observe exponential convergence of the method, even for examples beyond the scope of the theoretical results, -empirically compare the efficiency of our method to the cosine method of Fang and Oosterlee (2008).
Our numerical and theoretical results show that the offline-online decomposition can be used to find a quadrature rule that offers very satisfying results in terms of accuracy and efficiency.As a further advantage, the offline-online scheme is build in such a way that the resulting quadrature rule satisfies a pre-specified accuracy for the whole parameter domain of interest.Our univariate results lay the cornerstone for further research and show the potential for extensions of the method especially with regard to higher dimensions.
Parametric integrals also naturally arise in many other disciplines of applied mathematics.We refer to Gaß and Glau (2015) for this more general focus.
To achieve our goals, we apply the Empirical Magic Point Interpolation method developed by Barrault et al. (2004) in the context of parametric nonlinear partial differential equations.While they enforce an affine decomposition of parametric operators, we decompose parametric Fourier integrands.Sketching our idea, the starting point is the Fourier representation of the parametric option price, with generalized Fourier transform f K of the payoff function f K and the generalized Fourier transform ϕ T,q of the modelling random variable X q T .We follow the iterative Empirical Interpolation procedure outlined in Maday et al. (2009), that we describe in detail in section 3 below.For M ∈ N the method recursively gives magic points z * 1 , . . ., z * M ∈ Ω, basis functions q 1 , . . ., q M and θ M m := M j=1 (B M ) −1 jm q j and B M jm := q m (z * j ).The resulting price approximation is of the form The algorithm naturally decomposes into two phases.An offline phase, where the just mentioned quantities are constructed, and an online phase where real-time pricing is performed.More precisely, the two phases are described as follows.
Offline phase: For a given parameter space, -identify the magic points z * 1 , . . ., z * M ∈ Ω, and -precompute the integrals Ω θ M m (z) dz for all m ≤ M .
Online phase: For an arbitrary parameter constellation (K, T, q), -evaluate the Fourier integrands f K (−z * m )ϕ T,q (z * m ) for all m ≤ M and -assemble the sum in (1).
In the cases we consider, the number of summands M ranges in the dozens for a high accuracy already.Thus the evaluation of prices by ( 1) is fast and accurate.The following features of our problems at hand are key for the efficiency of the online phase: Typically the mapping (K, T, q, z) → f K (−z)ϕ T,q (z), i.e. the parametric integrand in (1) (i) is explicitly available, and (ii) enjoys desirable analyticity properties.
Thanks to (i), the evaluation of a single summand in (1) is effortless, and thanks to (ii), a few summands already yield high accuracy.In our numerical experiments for option pricing in univariate models, we achieve average absolute pricing accuracies ranging from 10 −6 to 10 −10 for 40 to 50 magic points, depending on the model used.This article is organized as follows: In the next section we revisit the framework for Fourier pricing in detail.In section 3 we adapt the Empirical Magic Point Interpolation method to Fourier pricing and describe the resulting algorithm that we call MagicFT.Based on Theorem 2.4 in Maday et al. (2009), we present in section 4 exponential convergence results under suitable analyticity conditions along with explicit error bounds.We investigate these analyticity properties for different payoff profiles and models in section 5.In section 6 we implement the algorithm and perform an empirical convergence study.In several case studies we investigate the MagicFT approximation for several models individually.Moreover we compare the MagicFT method to the popular cosine method of Fang and Oosterlee (2008).We conclude with an appendix that highlights essential features of Empirical Magic Point Interpolation and, for the sake of self-contained readability, presents detailed proofs of the convergence results.

Parametric Option Pricing with Fourier Transform
We compute option prices of the form Here we use payoff and model parameters In order to pass to the pricing formula in terms of Fourier transforms, we impose the following exponential moment condition which allows us to define for every (T, q) ∈ T × Q the extension of the characteristic function of X q T to the complex domain We further introduce the following integrability condition x → e η,x f K (x), ξ → ϕ T,q (ξ + iη) ∈ L 1 (R d ) for all (K, T, q) ∈ P. (Int) Furthermore, we denote the generalized Fourier transform of f K .The Fourier representation of option prices traces back to the pioneering works of Carr and Madan (1999) and Raible (2000).The following version is an immediate consequence of Theorem 3.2 in Eberlein et al. (2010).
Proposition 2.1 (Fourier pricing).Let η ∈ R d such that (Exp) and (Int) are satisfied.Then for every (K, T, q) ∈ P, Typically, that is for the most common option types, the generalized Fourier transform of f K is of the form for every z ∈ R d + iη with some constant c ∈ R and a function F : R d + iη → C.
Then the option prices (5) are indeed parametric Fourier integrals of the form As a first step in the numerical evaluation of (7) we employ an elementary symmetry by virtue of the identity f (−ξ) = f (ξ) for all ξ ∈ R d and all real-valued integrable functions f , and obtain which reduces the numerical effort by half.
In a second step we restrict the domain of integration to a compact set Ω ⊂ R d .The resulting error is determined by the decay of the integrand and will be further analyzed in appendix C. From now on we set Ω := Ω 1 × . . .× Ω d with bounded open intervals Ω 1 ⊂ R + + iη 1 and Ω j ⊂ R + iη j for j = 2, . . ., d.

Magic Point Interpolation for Integration
We present the Empirical Magic Point Interpolation method for parametric integration to approximate parametric integrals of the form with the parametric integrands h p (z) = h (K,T,q) (z) := ℜ f K (−z)ϕ T,q (z) (10) for every p = (K, T, q) in a given parameter set P. With P we associate the set of all parametric integrands.Let us point out that the following iterative procedure is defined for a more general set of parametric integrands that are not required to be of the form (10).
Before we closely follow Maday et al. (2009) to describe the interpolation method let us state our basic assumptions that ensure the well-definedness of the iterative procedure.
Assumptions 3.1.Let (Ω, .∞ ) and (P, .∞ ) be compact, P × Ω ∋ (p, z) → h p (z) bounded and p → h p be sequentially continuous, i.e. for every sequence p i → p we have h p i − h p ∞ → 0.Moreover, U is nontrivial in the sense that the set contains elements other than the function that is constantly zero.
For M ∈ N define a mapping I M from U to a tensor specified by and the Magic Point Integration with M points by where we denote by (B M ) −1 jm the entry in the jth line and mth column of the inverse of matrix B M .By definition, B M is a lower triangular matrix with unity diagonal and is thus invertible, confer also section A.1 in the appendix.The magic points z * 1 , . . ., z * M ∈ Ω and the basis functions q 1 , . . ., q M are recursively defined in the following way: In the first step, let Note that thanks to Assumption 3.1, these operations are well-defined.Then, recursively, as long as there are at least M linearly independent functions in U , u M is chosen according to a greedy procedure: The algorithm chooses u M as the function in the set U which is worst represented by the approximation with the previously identified M − 1 magic points and basis functions, Since every u ∈ U is a parametric function, u = h p for some p ∈ P, it can be identified by the associated parameter p.We call p * M ∈ P identifying u M in (16) the M th magic parameter.In the same spirit, let and we call z * M the M th magic point.The M th basis function is the residual, normed to 1, when evaluated at the new magic point z * M , Note the well-definedness of the operations in the iterative step thanks to Assumption 3.1 and the fact that the denominator in ( 18) is only zero, if all functions in U are perfectly represented by the interpolation I M −1 , in which case they span a linear space of dimension M − 1 or less and the procedure would have stopped already.
We may take three different perspectives on the approach: (i) Magic Point Integration is a quadrature rule for integrating parametric functions, where the interpolation nodes are chosen in a precomputation phase according to the set of integrands at hand.
(ii) Consider that for m = 1, . . ., M the functions θ M m are linear combinations of snapshot integrands h p * m with coefficients β m 1 , . . ., β m M and hence This means that Magic Point Integration is an interpolation method for parametric integrals in the parameter space.Thus, taking the error stemming from truncating the integration domain to Ω into account, equation ( 19) induces an approximation of the option price by a linear combination of snapshot prices, β m j Price K j ,T j ,q j .(20) (iii) In view of the representation of the option prices Price K,T,q as parametric Fourier integrals in (7), we use the Magic Point Integration algorithm to approximate parametric Fourier transforms that we call MagicFT as introduced in Gaß and Glau (2015).
From perspective (i), Magic Point Integration for parametric option pricing is an alternative to standard quadrature rules.Standard integration routines suffer from the curse of dimensionality of the integration domain.In contrast, under suitable analyticity conditions, the approximation error of Magic Point Integration decays exponentially in M , independently of the dimension of Ω, if the parameter space is one-dimensional, see Theorem 4.2 below.
Taking the point of view (ii), Magic Point Integration for parametric option pricing can be compared to a benchmark method for parametric option pricing by interpolation.Standard interpolation methods in the parameter suffer from the curse of dimensionality of the parameter space.In contrast, under suitable analyticity conditions on the integrands, the approximation error of Magic Point Integration decays exponentially in M , independently of the dimension of P, if the integration domain is one-dimensional, see Theorem 4.3 below.

Convergence Analysis of Magic Point Integration
To show the virtue of the method in its full generality, we review a general convergence result for Magic Point Interpolation derived in Maday et al. (2009).This result relates the convergence of Magic Point Interpolation to the best linear n-term approximation that is formally expressed by the Kolmogorov n-width.For a real or complex normed linear space X , • and U ⊂ X , the Kolmogorov n-width is given by where E(X , n) is the set of all n dimensional subspaces of X .We denote by L ∞ (Ω, C), • ∞ the Banach space of functions mapping from Ω ⊂ C d to C that are bounded in the supremum norm.
Then for arbitrary ε > 0 and C := c 4 e α +ε we have for all u ∈ U that The proposition directly follows from Theorem 2.4 in Maday et al. (2009), where a slightly different version that does not explicitly use the Kolmogorov n-width is provided.In order to keep our presentation self-contained and as transparent as possible, we present a detailed proof of the Proposition in Appendix A, where we also highlight essential features of the iterative Magic Point Interpolation procedure.

Exponential Convergence of Magic Point Integration for Parametric Option Pricing
In order to formulate our analyticity assumptions, we define the Bernstein ellipse B([−1, 1], ̺) with parameter ̺ > 1 as the open region in the complex plane bounded by the ellipse with foci ±1 and semiminor and semimajor axis lengths summing up to ̺ with the origin as the center and semimajor axis on the real axis.Moreover, we define for b < b ∈ R the generalized Bernstein ellipse by where the transform τ for every x ∈ C. For an arbitrary set X ⊂ R, we define the generalized Bernstein ellipse by ( 24) In order to estimate the error resulting from the Magic Point Interpolation method, we formulate two analyticity conditions.Condition (B1) is tailored to the case of univariate integration domains and (B2) to the case of univariate parameter spaces.
(B1) The function (p, z) → h p (z) is continuous on P × Ω and there exist functions and H 1 (p, z) has an extension H 1 : P × B(Ω, ̺) → C such that, for all fixed p ∈ P the mapping z → H 1 (p, z) is analytic in the interior of the generalized Bernstein ellipse B(Ω, ̺).
(B2) The function (p, z) → h p (z) is continuous on P × Ω and there exist functions and H 1 (p, z) has an extension H 1 : B(P, ̺) × Ω → C such that, for all fixed z ∈ Ω the mapping p → H 1 (p, z) is analytic in the interior of the generalized Bernstein ellipse B(P, ̺).

Parametric European Options, Generalized Moments and Other Univariate Integrals
In the generic situation where option prices have to be evaluated for a large set of different parameter constellations, a parametric integral of form ( 9) for a high dimensional parameter space and a univariate integration domain needs to be computed.This comprises many well-known examples such as prices of European and exotic options and sensitivities of these prices as expressed by the Greeks for different option and model parameters.Also risk measures like VaR and ES and other generalized moments or parametric univariate integrals fall into the scope of this paragraph.
Theorem 4.2.Let Ω ⊂ R and P ⊂ R D be compact.Fix some η ∈ R, some ̺ > 4 and assume that integrability conditions (Exp) and (Int) as well as analyticity condition (B1) are satisfied.Then for all p ∈ P and M ∈ N, where The proof is provided in Gaß and Glau (2015).In view of a self contained presentation we present the proof in detail in appendix B.

Basket Options, Multivariate Generalized Moments and Other Multivariate Integrals
The following result is for example well suited for the error analysis of Magic Point Integration for basket options for a single free parameter.In particular, this is interesting for real-time pricing of basket options with either varying strikes or varying maturities in a fixed calibrated asset model.Moreover, the paragraph applies to the computation of generalized moments such as covariances, and general multivariate integrals with a single varying parameter in the integrand.
Theorem 4.3.Let Ω ⊂ R d and P ⊂ R be compact.Fix some η ∈ R d , some ̺ > 4 and assume that integrability conditions (Exp) and (Int) as well as analyticity condition (B2) are satisfied.Then for all p ∈ P and M ∈ N, where The proof is provided in Gaß and Glau (2015).Compared to the proof of Theorem 4.2 in appendix B, the only difference is that now the analyticity properties of H 1 with respect to the parameters p are exploited to derive an estimate for the best n-term approximation of U .
The implementation of Magic Point Interpolation inevitably involves additional problem simplifications and approximations in order to perform the necessary optimizations.In particular, instead of the whole parameter space a training set is fixed in advance.In this context, the results from Theorem 4.2 and 4.3 are only statements for the training set of functions.Rigorous a priori error bounds for integrals corresponding to parameters outside of the training set can be straightforwardly derived from the a priori error bounds for the Magic Point Interpolation method from Eftang et al. (2010).

Examples of Univariate Payoff Profiles
Table 1 presents a selection of payoff profiles f K for option parameter K as function of the logarithm of the underlying asset.We state the range of possible weight values η such that x → e ηx f K (x) ∈ L 1 (R) and the respective generalized Fourier transform exists.

Type
Payoff Weight Fourier transform Table 1: Typical payoff profiles for single stock options and the respective generalized Fourier transform.
Examining the generalized Fourier transforms of the payoff profiles f K in Table 1, we realize that all of them admit a factorization in the spirit of condition (B1) as for some c ∈ R. While all of the payoff profiles f K of Table 1 either are not differentiable or even discontinuous, the mapping z → K iz+c is a holomorphic function and thus perfectly fits the requirements of Theorem 4.2.

Example of a Multivariate Payoff Profile
The payoff profile of a call option on the minimum of d assets is defined as .
A similar decomposition as in ( 27) in the univariate case can directly be read off.
While the mapping K → f K (x) displays a kink, the mapping is analytic in the half space {K ∈ C | ℜ(K) > 0}.This perfectly qualifies the call option on the minimum of d assets for the convergence result provided in Theorem 4.3.

Examples of Asset Models
We present a selection of asset models that we use for pricing options in the numerical experiments in section 6 below.The MagicFT algorithm, as we apply it, operates on Fourier integrands that consist of the generalized Fourier transform of the option profile, f K , as well as the Fourier transform of the process that drives the underlying asset at maturity, ϕ T,q .Theoretically, Theorem 4.2 requires the analytic property from the characteristic function ϕ T,q of the model in the sense of condition (B1).Yet, for some models fulfilling this requirement means strongly restricting the parameter space.This would leave us with parameter spaces that are too limited for practical purposes.Empirically, however, we observe that condition (B1) may be replaced by a much weaker condition while still maintaining exponential convergence.The existence of a shared strip of analyticity S R (η) of width R ∈ (0, ∞) d given by where all ξ → ϕ T,q (ξ), T ∈ T , q ∈ Q, are analytic on, grants exponential convergence of the algorithm, already.Enforcing such a shared strip means imposing conditions on the model parameter space Q, too.Yet these restrictions turn out to be rather mild compared to the stronger condition (B1) of Theorem 4.2.
In the following model presentations we denote by Q the parameter space that the model as such is defined on.From this we derive admissible parameter sets Q such that condition (B1) is satisfied.If this is not possible, they are chosen to guarantee the existence of a shared strip of analyticity according to (30).Throughout the following model introductions, constant r > 0 denotes the risk-free interest rate.

Multivariate Black-Scholes Model
The d-variate Black-Scholes model is driven by a d-variate Brownian motion.The parameter space of the model solely consists of values determining the underlying covariance matrix σ ∈ R d×d , which is symmetric and positive definite.For a concise representation of the parameter space, we define Q as with the function σ : R d(d+1)/2 → R d×d defined by ( 32) σ(q) ij = q (max{i,j}−1) max{i,j}/2+min{i,j} , i, j ∈ {1, . . ., d}.
By construction, σ(q) is symmetric.The characteristic function of the process X q T , T ∈ T , q ∈ Q, driving log-returns in the model is then given by ( 33) For each q ∈ Q given by ( 31), the characteristic function of the d-variate Black-Scholes model is analytic in z on the whole of C d .We thus may choose the parameter set Q for the MagicFT algorithm according to the following remark.
Remark 5.1 (Q for the multivariate Black-Scholes model).Let σ i ≤ σ i ∈ R + for all i ∈ {1, . . ., d(d + 1)/2}.Define with the function σ given by (32).With the parameter set Q defined as above and compact T ⊂ R + , the characteristic function of the univariate Black-Scholes model satisfies condition (B1) of Theorem 4.2.

Univariate Merton Jump Diffusion Model
In the univariate case, the Merton Jump Diffusion model by Merton (1976) naturally extends the Black-Scholes model to a jump diffusion setting.The logarithm of the asset price process is composed of a Brownian part with variance σ 2 > 0 and a compound Poisson jump part consisting of normally N (α, β 2 ) distributed jumps arriving at a rate λ > 0. The model parameter space is thus given by (36) and the characteristic function of X q T with T ∈ T , q ∈ Q computes to As in the univariate Black-Scholes model, for each q ∈ Q and T > 0, the characteristic function ϕ T,q of the Merton model is holomorphic.
With the parameter set Q defined as above and compact T ⊂ R + , the characteristic function of the Merton model satisfies condition (B1) of Theorem 4.2.

Univariate CGMY Model
Another well-known Lévy model that we consider is the univariate CGMY model by Carr, Geman, Madan and Yor (2002).This class is also known as Koponen and KoBoL in the literature, see e.g.Boyarchenko and Levendorski ȋ (2002a) and as tempered stable processes.With the model parameter space given by ( 40) for all z ∈ R, where Γ(•) denotes the Gamma function.For no-arbitrage pricing we set the drift b ∈ R to The condition (M − 1) Y ∈ R in (40) guarantees b ∈ R. Contrary to Black-Scholes' and Merton's model, the domain in C that the characteristic function of the CGMY model is analytic on does not exist independently of its parametrization.Consequently, Theorem 4.2 does not apply to pricing in the CGMY model unless the parameter set that the algorithm may choose from is unreasonably restricted.Yet, empirically we maintain exponential convergence in the CGMY model case when Q and η are chosen such that all ξ → ϕ T,q (ξ), T ∈ T , q ∈ Q, share a common strip of analyticity S R (η) as introduced in (30) depending on η ∈ R and R > 0, the desired strip width.In the following, we derive conditions which guarantee the existence of such a strip.The result of our analysis will consist in a combined suggestion for the weight value η that complies with the restriction posed by the option choice as outlined by Table 1 and a set of restrictions on the parameter space.These restrictions guarantee a shared strip of analyticity as described above achieving a certain prescribed width R > 0.
Strip of analyticity for CGMY Before we are able to derive conditions on the parameter space that originate a shared strip of analyticity, let us first determine the strip of maximal width R > 0 that an individually parameterized characteristic function of the CGMY model ϕ T,q , T ∈ T , q ∈ Q, is analytic on.This strip in C is derived by analyzing the characteristic function ϕ T,q , T ∈ T , q ∈ Q, of the CGMY process on the domain of integration in (5) of Proposition 2.1 for different weight values.Let η ∈ R and consider the characteristic function ϕ T,q on the line The values of η for which ϕ T,q is analytic on the associated line (43) determine the width of the strip of analyticity of ϕ T,q .For these values of η ∈ R, both mappings need to be analytic on R. By (43), we have For analyticity of these two quantities on R we need to ensure that both hold.Inequalities (44) and (45) yield bounds η − , η + given by These two bounds span the strip of analyticity S R (η) for an individually parametrized characteristic function of the CGMY model, wherein η = (η + + η − )/2 = (G − M )/2 and diameter 2R = G + M , as shown in Figure 1.
Now we can translate these findings to conditions on the model parameter set to derive a compact set Q ⊂ Q and a value for η ∈ R that ensure a common strip of analyticity S R (η) for all mappings ξ → ϕ T,q (ξ), T ∈ T , q ∈ Q.From our considerations during the derivation above and in particular by ( 46) we conclude that such a Q and η need to satisfy (47) max We limit the rest of this analysis to the case of a call option where we necessarily have

R
Figure 1: For fixed parametrization q ∈ Q, the hatched area visualizes the strip of analyticity of the characteristic function of the CGMY process at T ∈ T , X q T .Its bounds are determined by G ≥ 0 and M ≥ 0.
by Table 1.With G ≥ 0 due to the model parametrization (40), the second inequality in (47) trivially holds automatically.Combining ( 47) and ( 48) thus yields condition ( 49) max A strip width of R > 0 consequently follows if the final strip condition (M + 1)/2 yields a strip of analyticity S R (η) with diameter 2R that all of the mappings ξ → ϕ T,q (ξ), T ∈ T , q ∈ Q, share.We collect and summarize these results in the following remark.
All ϕ T,q , T ∈ T , q ∈ Q, share a common strip of analyticity S R (η) with While the characteristic function of the CGMY model parametrized by Q of (52) in general does not satisfy condition (B1) of Theorem 4.2, empirically we still observe exponential convergence of the MagicFT algorithm.
Additionally, to avoid forcing the algorithm to support unrealistic parameter constellations, impose the following additional plausibility restriction.

Univariate Normal Inverse Gaussian Model
Another Lévy model we present is the univariate Normal Inverse Gaussian (NIG) model.The parameterization consists of δ, α > 0, β ∈ R, with α 2 > β 2 .The model parameter set Q is thus given by ( 54) The characteristic function of X q T for this model is given by for T ∈ T , q ∈ Q, wherein the no-arbitrage condition requires The second condition in (54), As in the CGMY model, the analyticity condition (B1) posed by Theorem 4.2 is not satisfied by all realistic parameter choices q ∈ Q.We therefore, analogously to the CGMY case, derive a common strip of analyticity.
(57) All ϕ T,q , T ∈ T , q ∈ Q, share a common strip of analyticity S R (η) with While Q of (57) in general does not satisfy (B1) of Theorem 4.2, empirically we still observe exponential convergence of the MagicFT algorithm.
Remark 5.6 (Plausibility constraint on Q in the univariate NIG model).Let q ∈ Q of (54).The implied variance σ 2 NIG of a univariate NIG process at t = 1, X q 1 , is given by , confer Prause (1999).To keep volatilities supported by the MagicFT algorithm within reasonable bounds 0 < σ − < σ + add the final restriction for all q ∈ Q of (57).

The univariate Heston Model
The models introduced above are all Lévy models.We now introduce the model by Heston (1993) that does not fall into this class but is an affine stochastic volatility model, instead.In the univariate Heston model, the asset price process (S q t ) t≥0 follows the stochastic differential equation dS q=(v 0 ,κ,θ,σ,ρ) t with the two Brownian motions W 1 , W 2 correlated by ρ ∈ [−1, 1] and with q ∈ Q defined by The Feller condition σ 2 ≤ 2κθ in Q of (62) ensures an almost surely non-negative volatility process (v t ) t≥0 .With T ∈ T , q ∈ Q, the characteristic function ϕ T,q of the log-asset price process (log(S t /S 0 )) t≥0 at T is given by for all z ∈ R, with supporting functions defined by confer Schoutens et al. (2004).We simply choose Q ⊂ Q to be a bounded subset of the parameter space.
Remark 5.7 (Q for the univariate Heston model).Choose bounds for the initial value of the volatility process, 0 < v 0 ≤ v 0 , for its speed of mean reversion, 0 < κ ≤ κ, the long-term volatility mean, 0 < θ ≤ θ, and the volatility of the volatility process itself, 0 < σ ≤ σ, and a domain for the correlation parameter, Despite the fact that Q defined above in general might not satisfy condition (B1) of Theorem 4.2, we still observe exponential convergence of the MagicFT algorithm.
For an analysis of the strip of analyticity in the Heston model, see Levendorskiy (2012).

Numerical Experiments
In the previous sections we introduced the MagicFT algorithm for option pricing and presented several asset models and option types.We also proved theoretical claims for option pricing with the MagicFT algorithm.In this section we numerically validate these theoretical claims and provide empirical indication that the scope of the algorithm extends to a much wider class of pricing applications than suggested by the theorems earlier.

Implementation
The implementation of the algorithm in Matlab introduces some simplifications as suggested by e.g.Remark 3 in Maday et al. (2009).A theoretical argumentation for the discretization approach described in the following can be found in Eftang et al. (2010) and Maday et al. (2016).The continuous parameter space P is thus replaced by a discrete parameter cloud randomly sampled.Each magic parameter that the algorithm selects is a member of this discrete set.Consequently, the set U that the algorithm is trained on is replaced by a discrete set, as well.Additionally, we take Ω to be a discrete set with a finite number of points in each spacial dimension.Each function u ∈ U is then represented by its evaluation on this discrete Ω and is thus replaced by a finite-dimensional vector, numerically.The optimization steps from ( 15)-( 17) thus reduce to a search on finite sets.When all h p * m ∈ U for m = 1, . . ., M are identified, they are integrated using Matlab's quadgk routine (with an absolute tolerance requirement of 10 −14 , a relative tolerance requirement of 10 −12 , a maximum number of intervals of 200, 000) and linearly assembled to derive the quantities Ω θ M m (z) dz for m = 1, . . ., M .

Model fixed parameters free parameters
Table 2: In the numerical experiments, we price European call options as an example.
Various models have been selected.In the implementation, the Fourier integrands that the algorithm constructs the basis functions q m with are parametrized according to the intervals above.For each model investigated, U consists of a pool of |U | = 4000 Fourier integrands.

Empirical Convergence
We study the empirical convergence of our implementation of the MagicFT pricing algorithm.A plain vanilla European call option on one asset serves as an example.
We investigate the convergence in several models.For each model we set up a pool U of parametrized Fourier integrands that the algorithm picks from.For each model, the discrete parameter pool is chosen as a uniform sample of magnitude |P| = 4000 from the free parameter ranges enlisted in Table 2.
It is interesting to note that, not necessarily all model parameters have to be considered in the parametric option pricing.For example the parameter Y in the CGMY model reflects the degree of roughness of the paths of the process and can in principle be estimated from the historical stock price data.In contrast to other model parameters, such as those determining variance and skewness, the path behaviour is generally consider be stable over time.We can therefore fix the parameter Y = 1.1 in the CGMY model.Additionally, for the NIG and CGMY model, a shared strip of analyticity of width R = 1/2 is enforced such that for all investigated models, the dampening factor η could be set to η = −1.5.Furthermore, all model restric-tions stated in section 5.3 are respected.Also, implied variances are kept in the interval [0.01 2 , 0. Figure 2 shows the empirically observed error decay during the offline phase of the algorithm for all five considered models in the number of basis functions M .For each model, the quantity max z∈Ω u M (z) − I M −1 (u M )(z) is shown for increasing values of M .The algorithm has been instructed to construct basis functions q m until an error threshold of 10 −10 has been reached in step ( 16) or until M has reached the value 50.This offline phase based on the Matlab implementation described in the previous Section 6.1 takes less than 1 minute of time on a standard laptop computer for each model.
We observe exponential error decay in all considered models.Recall that Theorem 4.2 predicts this behavior only for the Black-Scholes and the Merton model where analyticity of the associated Fourier integrands is parameter independent.For the other two Lévy models, however, the existence of a shared strip of analyticity results in exponential error decay, as well.In case of the Heston model, the issue of analyticity of the Fourier integrands in U has not been investigated here.Still, we observe exponential error decay too.The empirical results depicted in Figure 2 thus indicate that it might be promising to investigate a theoretical result providing exponential error decay beyond the scope of Theorem 4.2.

Out of sample pricing study
In the previous paragraph we studied empirical convergence during the offline phase of the algorithm.More precisely, we investigated for several models how accurately all Fourier integrands in the given pool U could be approximated on their integration interval Ω by the M selected integrands or rather by the basis functions q m , m = 1, . . ., M , constructed thereof.Now we analyze, how the observed accuracy on the level of in sample integrands translates to the accuracy in an out of sample call option pricing exercise.MagicFT pricing error decay Figure 3: Pricing error decay study on 1000 out of sample parameter constellations for different models.In each model, for increasing values of M , the L ∞ error over the randomly drawn parameter sets is evaluated.The parameter sets have been drawn from the intervals given by Table 2.
To this extent we randomly draw 1000 parameter constellations for each model according to the same rules as in the offline phase.For each such sample we compute the respective Fourier price by numerical integration on [0, 65] thus containing the discrete Ω that the MagicFT algorithm has been trained on.We integrate using Matlab's quadgk with absolute tolerance of 10 −12 and 200, 000 integration intervals.Additionally, in each model we approximate all prices associated with the randomly drawn parameters for increasing values of M , evaluate the L ∞ error and study its decay in M as depicted in Figure 3.
We observe exponential rates for all considered models.Curiously, the error decay attains plateau-like shapes, especially for higher values of M .We explain this decay structure by assuming that each plateau is associated with a certain single parameter realization from the random sample that dominates the L ∞ error until a magic parameter close to it or rather the respective basis function contributes to the approximation of the belonging price.Due to such outliers, the order in which the offline phase errors were decaying in Figure 2 has changed.
In Figure 4, we depict evaluations of the absolute as well as the relative pricing Figure 4: Results of the out of sample pricing exercise.For each of the five considered models, 1000 parameter sets have been drawn from the intervals given by Table 2.For each set, the Fourier price as well as the MagicFT price have been calculated.On the left column, all absolute errors are depicted.On the right, the relative errors are shown.errors for all out of sample parameter sets, individually.Here, relative errors have been computed only for prices larger than 10 −3 to exclude numerical noise.In each model, M is set to its final value assigned during the respective model's offline phase and can be read off from Figure 2.
Pricing accuracy in this out of sample pricing exercise reaches very satisfactory levels albeit the achieved accuracies vary between the considered models.For all models, average absolute pricing accuracy reaches levels between avg min ≈ 10 −12 in the Black-Scholes model and avg max ≈ 10 −8 for the Merton model.Average relative pricing accuracy ranges between 10 −12 and 10 −7 .We observe individual outliers for all models.Even occasional mispricing, however, stays within practically acceptable bounds smaller than 10 −5 .The ten worst absolute errors in each model are further addressed in the next section.

Individual Case Studies
We take a closer look into the numerical results for each model individually.For this we are interested in the distribution of magic parameters in each model.

Black-Scholes
During the offline phase of the algorithm for the Black-Scholes model only the option strike K has been fixed, K = 1.The model parameter σ as well as the two other parameters S 0 /K and maturity T were allowed to vary within the bounds assigned by Table 2.In the Black-Scholes case, the individual parameter intervals tensorize meaning that any combination of parameter values respecting the individual bounds can be picked by the algorithm.As the left part in Figure 5 demonstrates for the magic parameter choices for S 0 /K and T , however, rather extreme constellations have been selected.The right part of Figure 5 provides a complete overview over all parameter combinations selected in the offline phase of the algorithm for the Black-Scholes model.With the exception of T and σ combinations, rather extreme parameter pairs have been selected.This special behavior is not surprising, since T and σ always appear together as a product in the Fourier integrands of the Black-Scholes model, compare the definition of the characteristic function in the Black-Scholes model in (33).The even distribution of the (T, σ) parameter pairs thus reflects the even distribution of all individual parameters over their domain, observable on the elements on the main diagonal of the figure.

Comparison with the Cosine method
A wide range of different methods for the efficient evaluation of Fourier integrals for option pricing have been successfully applied.One of these methods is the popular cosine method of Fang and Oosterlee (2008).We use the cosine method as relevant benchmark to our MagicFT method, also since an implementation for the Black-Scholes, Heston and Merton model (among others) by the original authors is publicly available in the BENCHOP project, see von Sydow et al. (2015).Despite their similarities as Fourier pricing routines, both methods differ conceptually since MagicFT is an offline-online scheme.
In order to compare the MagicFT and the cosine method for the Black-Scholes, Heston and Merton model, we use for each model the parameter sample set from Section 6.3.The accuracies of both methods will be measured against the Fourier integral with Ω = [0, 65] as benchmark.
To allow the comparison of both methods, we will only consider parameter constellations for which this benchmark does not exhibit truncation errors above a threshold of ε params > 10 −8 .Therefore 11 of the 1000 parameter tuples have been omitted according to the criterion The cosine method allows to set a specific integration range via the parameter L in equation ( 49) of Fang and Oosterlee (2008).As mentioned in de Innocentis and Levendorski ȋ (2014) this parameter needs to be carefully selected to guarantee convergence of a whole range of pricing parameters.In the preparation of our numerical studies, we have identified that the parameter L = 14, L = 18 and L = 3.1 for the models Black-Scholes, Heston and Merton, respectively, lead to the best possible convergence results.
The accuracy of both numerical methods for varying numbers of nodes M = 1 to M = 50 is shown in Figure 6.Instead of comparing milliseconds in CPU time we use the number of summands M of both methods as a measure for the computational complexity of both pricing routines.We deem this approach justified by the fact that pricing in both methods consists of assembling sums of known coefficients multiplied essentially by the characteristic function of the underlying model, which is available in closed form for all examples considered.Thus Figure 6 can be directly used to infer the efficiency of the numerical methods.The plots show that the errors of the MagicFT method is significantly below the error of the cosine method from M = 15 onwards.In the Black-Scholes case, both methods show a similar rate of convergence, while MagicFT is more accurate in absolute values.In the Heston model, both methods show an exponential rate of convergence.Here, MagicFT exhibits a higher rate and is more accurate.For the Merton model, we observe that the cosine method is not accurate for the whole range of parameters.To give a quantitative example: for M = 50, the error of MagicFT reaches levels of the order of 10 −6 , whereas the cosine method still shows errors of 3.8 • 10 −4 for the Black-Scholes, 1.1 • 10 −2 for the Heston and 0.13 for the Merton model.
As mentioned before, the literature on Fourier-based pricing offers many other reliable and efficient approaches to evaluate Fourier integrals.In particular, Levendorski ȋ (2016) proposes to choose an appropriate deformation of the contours of integration prior to the discretization.This approach leads to very accurate results already for few discretization points and thus is especially attractive.It is interesting to consider a combination of the MagicFT method with the the approach of Levendorski ȋ (2016).Both methods thereby would benefit mutually and we expect further gains in efficiency.To be more precise, we suggest to first choose the contour of integration optimally in regard to analyticity.Then, as a second step the Mag-icFT algorithm can be applied for the resulting parametric integrals.This would allow MagicFT to benefit from the improved region of analyticity of the integrands.

Outlook
The results of the experiments for univariate Fourier integrals indicate that extensions of MagicFT to multivariate option pricing is promising.Firstly, the univariate method benefits from the offline phase when compared to the cosine method.Secondly, the architecture of this offline-online method guarantees a rate of convergence that does not intrinsically suffer from the curse of dimensionality since it directly relates to the best n-term approximation.The core of the offline phase is an optimization procedure, that in contrast will directly be affected by the curse of dimensionality.It is, however, crucial to notice that the offline phase only needs to be performed once for a whole model class and option type.Therefore the run-time of the offline phase can be seen as part of the implementation phase of the pricing routine.For all of our univariate examples the offline phase has required only one minute on a standard laptop.For moderate dimensions we therefore expect that a direct extension of the algorithm will still lead to a practically useful method.
Our implementation of the offline phase of the MagicFT algorithm follows the standard procedure in the literature.Here a global optimization routine that requires the evaluation of the integrands for a large number of samples in the parameter space is used to find the magic points and basis functions.This optimization algorithm exhibits reasonable run-times in our numerical examples.For applications in higher dimensions, refinements of this implementation of the offline phase can be beneficial and are currently being investigated in the empirical integration literature, see Maday et al. (2016).

Conclusion
We have introduced the MagicFT algorithm for parametric option pricing (POP).Analyticity conditions theoretically guarantee an exponential rate of convergence of the method in the number of magic points.The numerical experiments confirm these findings and suggest an exponential rate of convergence even for models and options beyond the scope of our theoretical results.This gives rise to the hope that further valuable theoretical results can be established.
Thanks to its architecture, the method is highly efficient for a pre-specified range of parameter constellations of interest.In contrast to other interpolation methods, there are no generic geometric constraints for the choice of the parameter space.We have compared experimentally the performance of the MagicFT method to the cosine method for a whole range of parameters.This comparison indicates that the MagicFT method is beneficial when the efficiency for a whole range of parameters is crucial.

A Properties of Magic Point Interpolation
For the reader's convenience we state useful features of the algorithm and give a detailed proof of the convergence result Proposition 4.1, which basically coincides with Theorem 2.4 of Maday et al. (2009).

A.1 General features
The Magic Point Interpolation algorithm satisfies some immediate properties, which are identified by Barrault et al. (2004) and Maday et al. (2009) and summarised in the sequel: Exact interpolation at magic points For all functions u ∈ U , the interpolation is exact at the magic points, in the sense that for every m = 1, . . ., M Next, we define the residuals for all z ∈ Ω.Our assumption on the Kolmogorov n-width guarantees the existence of constants c > 0 and α > log(4) such that for every n ∈ N, where E(X , n) is the set of all n dimensional subspaces of X = L ∞ (Ω, C).Thus, for M ∈ N and every c 1 > c there exists a linear subspace U M −1 ⊂ X such that for all q j , j < m, there exists v j ∈ U M −1 such that (74) In the sequel we derive a bound for r M ∞ in the case M > o.For all u ∈ U we conclude that where we used identities ( 70) and (71).Equation ( 68) shows that α m j = α m j (u) is independent of m and thus I m (u) − I m−1 (u) = α m m q m .By equation ( 67) we know that q m is maximal at z * m and together with equation ( 66 Finally, with inequality (75) we conclude and this proves the claim.

B Proof of Theorem 4.2
The following proof is taken from Gaß and Glau (2015).
Proof.In principle, we exploit the analyticity property of H 1 from condition (B1) to estimate the Kolmogorov n-width of the set U .This can conveniently be achieved by inserting an example of an interpolation method that is equipped with exact error bounds.We choose Chebyshev polynomial interpolation for this task.For polynomials of degree N ∈ N, the Chebyshev nodes are given by z k = cos π 2k+1 2N +2 for k = 0, . . ., N and the basis functions are defined as

C Truncation Error in Fourier Pricing
We introduce the following condition that is satisfied for a large class of models and payout profiles: For every N > 0, there exist constants α, C 1 , C 2 , m > 0 such that uniformly for every (K, T, q) ∈ P = K × T × Q, ℜ log(ϕ T,q (ξ + iη)) ≤ −C 1 |ξ| α for all |ξ| > N , (Gård) for all |ξ| > N .(Poly) Virtually for every payoff profile f K the generalized Fourier transform f K (•+iη) exists for some η ∈ R d and decays polynomially, uniformly in a reasonably large set of parameters K. Condition (Gård) already appears in another context, where it implies that the related bilinear form satisfies a so-called Gårding condition with respect to fractional Sobolev spaces of order α.This helps to classify the solution spaces of related weak solutions to associated Kolmogorov equations.For a proof of this implication as well as for numerous examples of classes of (time-inhomogeneous) Lévy processes satisfying the condition we refer to Glau (2016b) for the case η = 0 and to Glau (2016a) for η = 0.The following proposition is immediate.
Lemma C.1 (Truncation error).Assume for every (K, T, q) ∈ P, (Gård) and (Poly) and let Ω ⊂ R + × R d + iη and denote |Ω| the diameter of the largest ball centered in the origin that is contained in Ω.For every β < α there exists a constant c > 0 such that uniformly for every (K, T, q) ∈ P,

Figure 2 :
Figure 2: A study of the empirical order of convergence of the error in step (17) during the offline phase of the MagicFT algorithm.Five different models and European call options are considered.Both the models and the option are parametrized according to Table 2.The convergence result is theoretically backed by Theorem 4.2 for the Black-Scholes and the Merton model.A shared strip of analyticity of the respective Fourier integrands of width R = 1/2 has been enforced for the NIG and CGMY model.

Figure 5 :
Figure 5: Left: Parameter pairs (S 0 /K, T ) selected by the MagicFT algorithm in the offline phase of the Black-Scholes model.Right: All magic parameters selected during the offline phase of the algorithm for the Black-Scholes model (empty blue circles).The filled orange circles denote the ten parameter constellations that resulted in the maximal absolute pricing errors during the out of sample pricing exercise.

Figure 6 :
Figure 6: Efficiency study of the MagicFT and cosine method for the Black-Scholes, Heston and Merton models.The plots show the L ∞ error across 1000 random parameter constellations.