Chebyshev Interpolation for Parametric Option Pricing

Recurrent tasks such as pricing, calibration and risk assessment need to be executed accurately and in real-time. Simultaneously we observe an increase in model sophistication on the one hand and growing demands on the quality of risk management on the other. To address the resulting computational challenges, it is natural to exploit the recurrent nature of these tasks. We concentrate on Parametric Option Pricing (POP) and show that polynomial interpolation in the parameter space promises to reduce run-times while maintaining accuracy. The attractive properties of Chebyshev interpolation and its tensorized extension enable us to identify criteria for (sub)exponential convergence and explicit error bounds. We show that these results apply to a variety of European (basket) options and affine asset models. Numerical experiments confirm our findings. Exploring the potential of the method further, we empirically investigate the efficiency of the Chebyshev method for multivariate and path-dependent options.


Introduction
The development of fast and accurate computational methods for parametric models is one of the central issues in computational finance. Financial institutions dedicated to trading or assessment of financial derivatives have to cope with the daily tasks of computing numerous characteristic financial quantities. Examples of interest are prices, sensitivities and risk measures for products on different models and for varying parameter constellations. With regard to the ever growing market activities, more and more of these evaluations need to be delivered in realtime. In addition we face a constantly rising model sophistication since the original work of Black and Scholes (1973) and Merton (1973). From the early nineties on, stochastic volatility and Lévy models as well as models based on further classes of stochastic processes have been developed that reflect the observed market data in a more appropriate way. For asset models see e.g. Heston (1993), Eberlein, Keller and Prause (1998), Duffie, Filipović and Schachermayer (2003), Cuchiero, Keller-Ressel and Teichmann (2015). In the case of fixed income models see e.g. Eberlein andÖzkan (2005), Keller-Ressel, Papapantoleon and Teichmann (2013), Filipović, Larsson and Trolle (2014). The aftermath of the financial crisis [2007][2008][2009], moreover, has lead to a new generation of more complex models, for instance by incorporating more risk factors. The usefulness of a pricing model critically depends on how well it captures the relevant aspects of market reality in its numerical implementation. Exploiting new ways to deal with the rising computational complexity therefore supports the evolution of pricing models and touches a core concern of present mathematical finance. A large body of computational tasks in finance needs to be repeatedly performed in real-time for a varying set of parameters. Prominent examples are option pricing and hedging of different option sensitivities, e.g. delta and vega, that also need to be calculated in real-time. In particular for optimization routines arising in model calibration, large parameter sets come into play. Further examples arise in the context of risk control and assessment, such as for quantification and monitoring of risk measures. The following question serves as a starting point of our investigations: How to systemically exploit the recurrent nature of parametric computational problems in finance with the approached objective to gain efficiency? Looking for answers to this question, we focus on Parametric Option Pricing (POP) in the sequel.
In the present literature on computational methods in finance, complexity reduction for parametric problems has largely been addressed by applying Fourier techniques following the seminal works of Carr and Madan (1999) and Raible (2000). See also the monograph Boyarchenko and Levendorskiȋ (2002). Research in this area concentrates on adopting fast Fourier transform (FFT) methods and variants for option pricing. Lee (2004) accurately describes pricing European options with FFT. Further developments are for instance provided by Lord, Fang, Bervoets and Oosterlee (2008) for early exercise options and by Feng and Linetsky (2008) and Kudryavtsev and Levendorskiȋ (2009) for barrier options. Another path to efficiently handle large parameter sets that has been pursued in finance relies on reduced basis methods. These are techniques to solve parametrized partial differential equations. Sachs and Schu (2010), Cont, Lantos and Pironneau (2011), Pironneau (2011) and Haasdonk, Salomon and Wohlmuth (2012) and Burkovska, Haasdonk, Salomon and Wohlmuth (2015) applied this approach to price European, American plain vanilla options and European baskets. FFT methods on the one hand can be highly beneficial when the prices are required in a large number of Fourier variables, e.g. for a large set of strikes of European plain vanillas. On the other hand numerical experiments have shown a promising gain in efficiency of reduced basis methods when an accurate PDE solver is readily available. In essence all these approaches reveal an immense potential of complexity reduction by targeting parameter dependence. Hereto, they exploit the functional architecture of the underlying pricing technique for varying parameters.
Financial institutions have to deal simultaneously with a diversity of models, a multitude of option types, and-as a consequence-a wide variety of underlying pricing techniques. It is therefore worthwhile to explore the possibility of a generic complexity reduction method that is independent of the specific pricing technique. To do so, we focus on the set of option prices and the set of parameters of interest, disregard on purpose the pricing technology and view the option price as a function of the parameters. The core idea is now to introduce interpolation of option prices in the parameter space as a complexity reduction technique for POP.
The resulting procedure naturally splits into two phases: Pre-computation and real-time evaluation. The first one is also called offline phase while the second is also called online phase. In the pre-computation phase the prices are computed for some fixed parameter configurations, namely the interpolation nodes. Here, any appropriate pricing method-for instance based on Fourier, PDE or even Monte Carlo techniques-can be chosen. The real-time evaluation phase then consists of the evaluation of the interpolation. Provided that the evaluation of the interpolation is faster than the benchmark tool, the scheme permits a gain in efficiency in all cases where accuracy can be maintained. Then, we distinguish two use cases: • In comparison to the benchmark pricing routine, the fast evaluation of the interpolation will eventually outweigh the expensive pre-computation phase, if pricing is a task repeatedly employed. Optimization procedures are an obvious instance where this feature becomes advantageous.
• Even if the number of price computations is limited, we can still benefit from the separation of the procedure into its two phases. In this way, e.g., idle times in the financial industry can be put to good use by preparing the interpolation for whenever real-time pricing is needed during business activities.
The question arising at this stage is: Under what circumstances can we hope to find an interpolation method that delivers both reliable results and a considerable gain in efficiency?
One could now be tempted to proceed in a naive manner and first define an equidistant grid and then interpolate piecewise linearly in the parameter space. Numerical experiments for Black&Scholes call prices as function of the volatility, for instance, would then provide convincing evidence that the number of nodes needed for a given accuracy is considerably high. Increasing the polynomial degree might lead to better results. However, convergence might not be guaranteed. Runge (1901) showed that polynomial interpolation on equispaced grids may diverge-even for analytic functions. Second, the evaluation of the polynomial interpolants may be numerically problematic, as it is shown in Runge (1901) that "the interpolation problem for polynomial interpolation on an equidistant grid is exponentially illconditioned", a formulation we borrow from Trefethen (2011). For these reasons we abstain from polynomial interpolation with equidistant grids. Rather we take a step back and ask: Which methods for the interpolation of prices as functions of model and payoff parameters are numerically promising in terms of convergence, stability and implementational ease?
Regarding this research question, we need to take into consideration both the set of interpolation methods as such and the specific features of the functions we investigate. It is well-known that the efficiency of interpolation methods critically depends on the degree of regularity of the approximated function. For the core problem of our study-European (basket) options-we investigate the regularity of the option prices as functions of the parameters. We find that these functions are indeed analytic for a large set of option types, models and parameters. Taking the perspective of approximation theory, this inspires the hope to find suitable interpolation methods.
Empirically, we observe that parameters of interest often range within bounded intervals. One interpolation method that is highly effective for analytic functions on bounded intervals is Chebyshev interpolation. This intensively studied method enjoys excellent numerical properties-in stark contrast to polynomial interpolation on equally spaced nodal points. The interpolation nodes are known beforehand, implementation is straightforward and the method is numerically stable. For univariate functions that are several times differentiable, the method converges polynomially and for univariate analytic functions convergence is exponential. In a remarkable monograph, Trefethen (2013) gives a comprehensive review of Chebyshev interpolation. Its appealing theoretical properties are indeed of practical use as the software tool Chebfun 1 demonstrates. In this implemention Platte and Trefethen (2008) aim "to combine the feel of symbolics with the speed of numerics". Therefore Chebyshev interpolation is our method of choice. 2 Exploring the potential of interpolation methods for more than one single free parameter, we choose a tensorized version of Chebyshev interpolation: For parameters p ∈ R D , where D ∈ N denotes the dimensionality of the parameter space, the price Price p is approximated by tensorized Chebyshev polynomials T j with pre-computed coefficients c j , j ∈ J, as follows, Chebyshev interpolation is a standard numerical method that has proven to be extremely useful for applications in such diverse fields as physics, engineering, statistics and economics. Nevertheless, for pricing tasks in mathematical finance Chebyshev interpolation still seems to be rarely used and its potential is yet to be unfolded. Pistorius and Stolte (2012) use Chebyshev interpolation of Black&Scholes prices in the volatility as an intermediate step to derive a pricing methodology for a time-changed model. Independently from us, Pachon (2016) recently proposed Chebyshev interpolation as a quadrature rule for the computation of option prices with a Fourier type representation, which is comparable to the cosine method.
Our main results are the following: • Theorem 3.2 provides accessible sufficient conditions on options and models that guarantee an asymptotic error decay of order O ̺ − D √ N in the total number N of interpolation nodes where ̺ > 1 is given by the domain of analyticity and D is the number of varying parameters.
• More specific conditions for parametric European options in Lévy models are provided in Corollary 3.6, while Corollary 3.9 provides the framework for parametric basket options in affine models.
These results establish an error analysis that is based on the domain of analyticity of the prices as functions of the parameters. Observing that typical payoff functions are not smooth, we cannot expect an exponential error decay for interpolation with arbitrarily small maturities. Small maturities thus serve as an example that domains of analyticity need to be carefully studied.
• The investigations in Sections 3.2-4.2 show that for a large set of relevant (basket) options, models and free parameters a domain of analyticity can indeed be identified. This gives examples of relevant financial applications where (sub)exponential error decay is guaranteed.
To numerically validate the theoretical results we compare prices obtained by Chebyshev interpolation to benchmark prices by Fourier techniques.
• Numerical experiments in affine models confirm the theoretical error decay empirically for European call options ( Figure 5.3) and digital down&out options (Figure 5.4). For the considered model examples of Black&Scholes, Merton, CGMY and Heston we observe L ∞ -error levels of order 10 −10 using not more than N = (25 + 1) 2 Chebyshev interpolation nodes when D = 2 parameters are varied.
Numerical results show that already a small number of nodes leads to high accuracy. This motivates us to further explore the potential of the Chebyshev method for multivariate options. Here we deliberately go beyond the scope of our theoretical results and consider additional features like path-dependency.
• For multivariate basket and path-dependent options in the Black&Scholes, Heston and Merton model we use Monte Carlo as reference method. In all of our settings in Section 5.2 Chebyshev interpolation achieves an accuracy that is similar to the accuracy of the Monte Carlo simulation itself (10 −3 ) for D = 2.
In addition we present empirical results demonstrating the efficiency of the Chebyshev method.
• The gain in efficiency in comparison to Fourier techniques is first validated for bivariate options of European type in Section 5.3.1.
• Secondly, the explicit gains in efficiency in comparison to Monte Carlo methods are shown in Section 5.3.2 taking multivariate lookback options in the Heston model as examples.
The remainder of the article is organized as follows. In Section 2 we introduce Chebyshev interpolation in detail and present the general convergence results. Section 3 establishes a convergence analysis of Chebyshev interpolation for POP. We formulate analyticity conditions for the payoff profiles and models that guarantee (sub)exponential convergence of the method. These conditions are verified in Section 4 for different option types, models and free parameters. The numerical experiments in Section 5 confirm these findings using Fourier techniques. Pricing basket options, the gain in efficiency is numerically investigated. Experiments based on Monte Carlo and finite differences moreover suggest to further explore the potential of the approach beyond the scope of the theoretical investigations from the previous sections. The resulting conclusion and outlook are presented in Section 6. Finally, the appendix provides the proof of the multivariate convergence result.

Chebyshev Polynomial Interpolation
In this section we introduce the notation for Chebyshev interpolation. Following Trefethen (2013), the one-dimensional version is shown. Then we present the multivariate extension and convergence results. Consider an option price with a single varying parameter An interpolation of Price p with Chebyshev polynomials of degree N is of the form

Multivariate Chebyshev Interpolation
The Chebyshev polynomial interpolation (2.2)-(2.4) has a tensor based extension to the multivariate case, see e.g. Sauter and Schwab (2004). In order to obtain a nice notation, consider the interpolation of prices For a more general hyperrectangular parameter space P = [p 1 , p 1 ]×. . .×[p D , p D ], the appropriate linear transformations need to be performed. Let N := (N 1 , . . . , N D ) with N i ∈ N 0 for i = 1, . . . , D. The interpolation with D i=1 (N i + 1) summands is given by where the summation index j is a multiindex ranging over J := {(j 1 , . . . , j D ) ∈ N D 0 : j i ≤ N i for i = 1, . . . , D}, i.e. (2.7) The basis functions T j for j = (j 1 , . . . , j D ) ∈ J are defined by The coefficients c j for j = (j 1 , . . . , j D ) ∈ J are given by where ′′ indicates that the first and last summand are halved and the Chebyshev nodes p k for multiindex k = (k 1 , . . . , k D ) ∈ J are given by (2.10) p k = (p k 1 , . . . , p k D ) with the univariate Chebyshev nodes p k i = cos π k i

Convergence of Multivariate Chebyshev Interpolation
In the univariate case, it is well known that the error of approximation with Chebyshev polynomials decays polynomially for differentiable functions and exponentially for analytic functions.
In the multivariate case we will extend a convergence result from Sauter and Schwab (2004). We consider parametric option prices of form Price p for p ∈ P (2.11) with P ⊂ R D of hyperrectangular structure, i.e. P = [p 1 , p 1 ]×. . .×[p D , p D ] with real p i ≤ p i for all i = 1, . . . , D. We define the D-variate and transformed analogon of a Bernstein ellipse around the hyperrectangle P with parameter vector ̺ ∈ (1, ∞) D as Theorem 2.2. Let P ∋ p → Price p be a real valued function that has an analytic extension to some generalized Bernstein ellipse B(P, ̺) for some parameter vector ̺ ∈ (1, ∞) D and sup p∈B(P,̺) |Price p | ≤ V . Then The proof of the theorem is provided in Appendix A. Remark 2.4. In particular, under the assumptions required by Theorem 2.2 with N = D i=1 (N i + 1) denoting the total number of nodes, Corollary 2.3 shows that the error decay is of (sub)exponential order O ̺ − D √ N for some ̺ > 1.
In the setting of Theorem 2.2 additionally the derivatives of Price p are approximated by the according derivatives of the Chebyshev interpolation. The onedimensional result is shown in Tadmor (1986) and a multivariate result is derived in Canuto and Quarteroni (1982) for functions in Sobolev spaces. These results allow us to obtain the Chebyshev approximation of derivatives with no additional cost. To state the according convergence results we follow Canuto and Quarteroni (1982) and introduce the weighted Sobolev spaces for σ ∈ N by (2.14) |∂ α φ(p)| 2 ω(p) dp, wherein α = (α 1 , . . . , α D ) ∈ N D 0 is a multiindex and ∂ α = ∂ α 1 · · · ∂ α D and the weight function ω on P given by Then we apply the result of Canuto and Quarteroni (1982, Theorem 3.1) in the following corollary.
Corollary 2.5. Under the assumptions of Theorem 2.2 for any D 2 < σ ∈ N and any σ ≥ µ ∈ N 0 there exists a constant C > 0 such that Proof. In our setting we have P ⊂ R D of hyperrectangular structure. Under the assumptions of Theorem 2.2 it follows that p → Price p ∈ W σ 2 (P) and therewith p → Price p ∈ W σ,ω 2 (P). Before we apply Canuto and Quarteroni (1982, Theorem 3.1), which assumes P = [−1, 1] D , we investigate how the linear transformation τ P , as introduced in the proof of Theorem 2.2, influences the derivatives. Let p → Price p be a function on P. We set h(p) = Price p • τ P (p). Further, let I N ( h)(p) be the Chebyshev interpolation of h(p) on [−1, 1] D . Then, it directly follows First, let us assume D = 1, i.e. P = [p, p], and let α ∈ N 0 . For the partial derivatives it holds

Repeating this step iteratively yields
Therewith, the error on [−1, 1] is scaled with a factor 2 α (p−p) α . Extending this to the D-variate case with, this analogously results with α = (α 1 , . . . , α D ) ∈ N D 0 is a multiindex and From Theorem 3.1 in Canuto and Quarteroni (1982) the assertion follows directly for h(·) on P = [−1, 1] D , i.e. for any D 2 < σ ∈ N and any σ ≥ µ ∈ N 0 there exists a constantC > 0 such that For arbitrary P the constant from 2.16 has to be multiplied with the according factor resulting from the linear transformation τ P .
The result in Corollary 2.5 is given in terms of weighted Sobolev norms. In the following remark we connect the approximation error in the weighted Sobolev norm to the C l (P) norm, where C l (P) is the Banach space of all functions u in P such that u and ∂ α u with |α| ≤ l are uniformly continuous in P and the norm Corollary 2.6. Under the assumptions of Theorem 2.2 for any D 2 < σ ∈ N and any σ ≥ µ ∈ N 0 and l ∈ N 0 with µ − l > D 2 , there exists a constantC(σ) > 0 depending on σ, such that Proof. In the setting of Corollary 2.5, we start with the estimation of the approximation error in the weighted Sobolev norms, On P it holds that w(p) ≥ 1 and therewith we can deduce for the Sobolev norm with a constant weight of 1, Price (·) − I N (Price (·) )(·) W µ,ω 2 (P) ≥ Price (·) − I N (Price (·) )(·) W µ,1 2 (P) . With W µ 2 (P) the usual Sobolev space, |∂ α φ(p)| 2 dp, Corollary 6.2 from Wloka (1987) directly yields that for any l with µ − l > D 2 there exists a constantC such that In formula (2.18) we have derived a lower bound for the left hand side of expression (2.17). Next, we will find an upper bound for the right hand side of (2.17). From the definition of the weighted Sobolev norm, see (2.15), it follows Here, we apply P ω(p) dp = π D and that there exists a constant α 2 (σ) depending on σ such that Finally, using (2.18) and (2.19) in (2.17) yields an estimate of the approximation error in the C l (P) norm, Collecting all constants inC(σ), we achieve

Exponential Convergence of Chebyshev Interpolation for POP
In this section we embed the multivariate Chebyshev interpolation into the option pricing framework. We provide sufficient conditions under which option prices depend analytically on the parameters. Analytic properties of option prices can be conveniently studied in terms of Fourier transforms. First, Fourier representations of option prices are explicitly available for a large class of both option types and asset models. Second, Fourier transformation unveils the analytic properties of both the payoff structure and the distribution of the underlying stochastic quantity in a beautiful way. By contrast, if option prices are represented as expectations, their analyticity in the parameters is hidden. For example the function K → (S T −K) + is not even differentiable, whereas the Fourier transform of the dampened call payoff function evidently is analytic in the strike, compare Table 4.1 on page 19.

Conditions for Exponential Convergence
Let us first introduce a general option pricing framework. We consider option prices of the form where f p 1 is a parametrized family of measurable payoff functions f p 1 : R d → R + with payoff parameters p 1 ∈ P 1 and X p 2 is a family of R d -valued random variables with model parameters p 2 ∈ P 2 . The parameter set Typically we are given a parametrized R d -valued driving stochastic process H p ′ with S p ′ being the vector of asset price processes modeled as an exponential of H p ′ , i.e.
and X p 2 is an F T -measurable R d -valued random variable, possibly depending on the history of the d driving processes, i.e. p 2 = (T, p ′ ) and We now focus on the case that the price (3.1) is given in terms of Fourier transforms. This enables us to provide sufficient conditions under which the parametrized prices have an analytic extension to an appropriate generalized Bernstein ellipse. For most relevant options, the payoff profile f p 1 is not integrable and its Fourier transform over the real axis is not well defined. Instead, there exists an exponential dampening factor η ∈ R d such that e η,· f p 1 ∈ L 1 (R d ). We therefore introduce exponential weights in our set of conditions and denote the Fourier transform of and we denote the Fourier transform of e η,· f ∈ L 1 (R d ) byf (·−iη). The exponential weight of the payoff will be compensated by exponentially weighting the distribution of X p 2 and that weight will reappear in the argument of ϕ p 2 , the characteristic function of X p 2 .
The proof of the theorem builds on the following proposition that can be derived from Eberlein, Glau and Papapantoleon (2010, Theorem 3.2), observing that the Proposition 3.3. Assume conditions (A1), (A3) and that the mapping z → f p 1 (−z− iη)ϕ p 2 (z + iη) belongs to L 1 (R d ) for every p = (p 1 , p 2 ) ∈ P, then the option prices (3.1) have for every p = (p 1 , p 2 ) ∈ P the Fourier transform based representation We are now in a position to prove Theorem 3.2.
Proof of Theorem 3.2. In view of Theorem 2.2, it suffices to show that the mapping p → Price p has an analytic extension to the generalized Bernstein ellipse B(P, ̺). Thanks to Proposition 3.3, we have the following Fourier representation of the option prices, Due to assumptions (A2) and (A4) the mapping has an analytic extension to B(P, ̺). Let γ be a contour of a compact triangle in the interior of B([p i , p i ], ̺ i ) for arbitrary i = 1, . . . , D. Then thanks to assumption (A2) and (A4) we may apply Fubini's theorem and obtain Moreover, thanks to assumptions (A2) and (A4), dominated convergence shows continuity of p → Price p in B(P, ̺) which yields the analyticity of p → Price p in B(P, ̺) thanks to a version of Morera's theorem provided in Jänich (2004, Satz 8).
Similar to Corollary 2.5 in the setting of Theorem 3.2 additionally the according derivatives are approximated as well by the Chebyshev interpolation. A very interesting application of this result in finance is the computation of sensitivities like delta or vega of an option price for risk assessment purposes. Theorem 3.2 together with Corollary 2.5 yield the following corollary.
Corollary 3.4. Under the assumptions of Theorem 3.2, for all l ∈ N, µ and σ with σ > D 2 , 0 ≤ µ ≤ σ and µ − l > D 2 there exist a constant C, such that , where the spaces and norms are defined in Section 2.2.
Remark 3.5. The upper bound in condition (A2) is tailored to the application to payoff profiles with varying option parameters, compare Section 4.1 and particularly Table 4.1, below. Examples of (time-inhomogeneous) Lévy processes satisfying the upper bound in (A4) for P 2 = [T , T ] ⊂ (0, ∞) are provided in Glau (2016), where also its connection to the parabolicity of the Kolmogorov equation of the process is investigated. In a fix parameter setting, condition (A3) for α ∈ (1, 2] in Glau (2016) translates to our upper bound in (A4). This condition is satisfied for a large variety of models.

Convergence Results for Selected Option Prices
In the previous subsection, Conditions 3.1 and Theorem 3.2 introduced a framework in which the Chebyshev approximation achieves (sub)exponential error decay. Now we relate this abstract framework to two specific settings.

European Options in Univariate Lévy Models
Let r be the deterministic and constant interest rate. We consider the parametrized family of asset prices, with t ≥ 0. For fixed π = (σ, b) let L π be a Lévy process with characteristics (σ 2 , b, F ) and |x|>1 |x|F (dx) < ∞. Due to the Theorem of Lévy-Khintchine we have the following representation of the characteristic function of the parametrized Lévy process We assume L π is defined under a risk neutral measure. Therefore, for every π ∈ Π we assume E(e L π t ) < ∞ for some and equivalently all t > 0 and the drift condition to ensure that the discounted asset price process is a martingale. Denoting the pure jump part of the exponent of the characteristic function, this reads as The fair value at time t = 0 of a European option with payoff function f K for In order to guarantee (sub)exponential convergence of the Chebyshev interpolation, in the following corollary, we translate the exponential moment condition (A3), the analyticity condition and the upper bound in (A4) to conditions on the cumulant generating function ψ π .
Corollary 3.6. Let Conditions (A1) and (A2) be satisfied for weight η ∈ R and ̺ 1 > 1 and set and there exist constants C 1 , C 2 > 0 and β < 2 such that If additionally one of the following conditions is satisfied, then there exist constants C > 0 and ̺ > 1 such that Proof. In view of Theorem 3.2 and Corollary 2.3, under the hypothesis of this corollary, it suffices to verify Conditions (A3) and (A4). Thanks to the exponential moment condition (3.11), Sato (1999, Theorem 25.17) implies that the validity of the Lévy-Khintchine formula (3.6) extends to the complex strip U : for every z ∈ U. (3.14) In particular, E(e −ηL π t ) = e tψ π (iη) < ∞ for every π ∈ Π, which yields (A3). Denote s 0 := log(S 0 ). The assumptions of the corollary yield for every z ∈ U the analyticity of (s 0 , t) → ϕ p 2 =(e s 0 ,t) (z) = S iz 0 e tψ π (z) on the generalized Bernstein ellipse B(P 2 , ̺ 2 ) for some parameter ̺ 2 ∈ (1, ∞) 2 , i.e. the validity of the analyticity Condition in (A4). Next we consider An elementary manipulation shows for every z ∈ R, which in combination with the assumption on the imaginary part ofψ yields the upper bound in Condition (A4) under assumption (ii). In view of condition (A4) also follows under assumption (i).
Let us notice that for the Merton model, Condition (3.11) is satisfied for every η ∈ R and Condition (3.12) is satisfied with β = 1. Examples of Lévy jump diffusion modes, i.e. examples satisfying condition (i) of Corollary 3.6 are e.g. the Black&Scholes and the Merton model. Examples of pure jump Lévy models satisfying condition (ii) of Corollary 3.6 are provided in Glau (2016), compare also Remark 3.5 and Table 4.2 below.
For simplicity and in accordance with our numerical experiments, we fixed the model parameters π. The analysis can be extended to varying model parameters π.
Remark 3.7. Assuming (3.11), we observe that the mappings Remark 3.8. Corollary 3.6 allows us to determine explicit error bounds for call options in Lévy models. The fair price at t = 0 of a call option with strike K and maturity T with deterministic interest rate r ≥ 0 is under a risk-neutral probability measure. Noticing that it suffices to interpolate the function This effectively reduces the dimensionality D of the interpolation problem by one.
Let us fix some η < −1 and let P In particular, there exists a constant C > 0 such that Additionally, when fixing the maturity T , letting ζ =

Basket Options in Affine Models
Let X π ′ be a parametric family of affine processes with state space D ⊂ R d for π ′ ∈ Π ′ such that for every π ′ ∈ Π ′ there exists a complex-valued function ν π ′ and a C d -valued function φ π ′ such that for every t ≥ 0, z ∈ R d and x ∈ D. Under mild regularity conditions, the functions ν π ′ and φ π ′ are determined as solutions to generalized Riccati equations. We refer to Duffie, Filipović and Schachermayer (2003) for a detailed exposition. The rich class of affine processes comprises the class of Lévy processes, for which ν π ′ (t, iz) = tψ π ′ (z) with ψ π ′ given as some exponent in the Lévy-Khintchine formula (3.6) and φ π ′ (t, iz) ≡ 0. Moreover, many popular stochastic volatility models such as the Heston model as well as stochastic volatility models with jumps, e.g. the model of Barndorff-Nielsen and Shephard (2001) and time-changed Lévy models, see Carr, Geman, Madan and Yor (2003) and Kallsen (2006), are driven by affine processes.
Consider option prices of the form where f K is a parametrized family of measurable payoff functions f K : R d → R + for K ∈ P 1 . Corollary 3.9. Under the conditions (A1)-(A3) for weight η ∈ R d , ̺ ∈ (1, ∞) D and P = P 1 × P 2 ⊂ R D of hyperrectangular structure assume (i) for every parameter p 2 = (t, x, π ′ ) ∈ P 2 ⊂ R D−m that the validity of the affine property (3.21) extends to z = R + iη, i.e. for every z ∈ R + iη, Then there exist constants C > 0, ̺ > 1 such that Proof. Thanks to Theorem 3.2 and Corollary 2.3 and in view of the assumed validity of Conditions (A1)-(A3), it suffices to verify Condition (A4). While assumptions (i) and (ii) together yield the analyticity condition in (A4), part (iii) provides the upper bound in (A4).
The existence of exponential moments of affine processes has been investigated by Keller-Ressel and Mayerhofer (2015), where also criteria are provided under which formula (3.21) and the related generalized Riccati system can be extended to complex exponential moments z ∈ C d . The question has already been treated for important special cases, which allow for more explicit conditions. Filipović and Mayerhofer (2009) consider affine diffusions and Cheridito and Wugalter (2012) investigate affine processes with killing when the jump measures possess exponential moments of all orders.

Analyticity Properties of Selected Payoff Profiles
We first list in Table 4.1 a selection of payoff profiles f K for option parameter K as function of the logarithm of the underlying. As we have seen in Proposition 3.3, we can represent option prices under certain conditions by their Fourier transform. Therefore, the table provides the Fourier transform f K of the respective option payoff, as well. All examples in Table 4.1 are special cases by k := log(K) of the following lemma that has a straightforward proof.

Type
Payoff Weight Fourier transform f K holomor- Proof. Both assertions are immediate consequences of for all z ∈ R d and all k ∈ C.
Example 4.2 (Call on the minimum of several assets). The payoff profile of a call option on the minimum of d assets is defined as For some payoff profiles it is not possible to transform them to integrable functions by a multiplication with an exponential dampening factor. Thus we split them into summands which can be appropriately dampened. Therefore, we decompose the unity, 1 = 2 d j=1 1 O j (x) a.e. with the distinct orthants O j of R d . More precisely, for j = 1, . . . , 2 d let ζ j := (ζ j 1 , . . . , ζ j d ) with ζ j i ∈ {−1, 1} for the 2 d different possible configurations and let For each j = 1, . . . , 2 d we choose for the call, respectively put, profile (4.3) η j := −(1 + ε)ζ j , respectively η j := εζ j for some ε > 0 for each j = 1, . . . , 2 d .
The following proposition enables us to prove analytic dependence of prices of basket put and call options on the strike. Then for each k ∈ R the payoff profile f k of a basket on a call, respectively put, is of the form f k (x) = e x 1 + . . .
Moreover, for every j = 1, . . . , 2 d and every k ∈ R we have that x → e η j ,x f k j (x) ∈ L 1 (R d ) and for every z ∈ R d the mapping has a holomorphic extension.
Proof. We show the claim for the call option. The put option case is proved analogously. For every j = 1, . . . , 2 d , under the assumptions of the proposition it is elementary to show that x → e η j ,x f k j (x) ∈ L 1 (R d ). Moreover, letting u := z − iη j with z ∈ R d we have From the geometry of A(k) and the holomorphicity of the integrand in (4.6) we can deduce that the integral in (4.6) is holomorphic in k and thus obtain the assertion of the proposition.

Analyticity Properties of Asset Models
In this section, we shortly introduce a range of asset models and present analyticity properties of their Fourier transforms. For some models and some parameters, the domain of analyticity is immediately observable. For some non-trivial cases, we state the domain briefly. Throughout the section, T > 0 denotes the time to maturity of the option while r > 0 refers to the constant risk-free interest rate.

Multivariate Black&Scholes Model
In the multivariate Black&Scholes model, the characteristic function of L π T , the multivariate analogon of (3.5), with π = (b, σ) is given by with covariance matrix σ ∈ R d×d such that det(σ) > 0 and drift b ∈ R d adhering to the drift condition In the Black&Scholes model, analyticity in the parameters is immediately confirmed, i.e. p 2 → ϕ p 2 (z), given in (4.7), is holomorphic for every z ∈ R d . The admissible parameter domain, however, is restricted to parameter constellations such that σ is a covariance matrix.
Remark 4.4. Let η ∈ R d be the chosen weight in Conditions 3.1 and let the open set U be given by Then for every z ∈ R d , (T, b, σ) → ϕ p 2 =(T,b,σ( σ)) (z + iη) given by (4.7) is analytic on U . Note that U does not depend on η.

Univariate Merton Jump Diffusion Model
As second Lévy model we state the Merton Jump Diffusion model by Merton (1976). The logarithm of the asset price process follows a jump diffusion with volatility σ 2 enriched by jumps arriving at a rate λ > 0 with normally N (α, β 2 ) distributed jump sizes. Here, the characteristic function of L π T in (3.5) with π = (b, σ, α, β, λ) is given by (4.10) where σ > 0, α ∈ R and β ≥ 0. Since the characteristic function in the Merton Jump Diffusion model is composed of analytic functions, it is itself analytic on the whole parameter domain.

Univariate CGMY Model
We consider the CGMY model by Carr, Geman, Madan and Yor (2002), which is a special case of the extended Koponen's family of Boyarchenko and Levendorskiȋ (2000), with characteristic function (3.6) of L π T in (3.5) with π = (b, C, G, M, Y ) where C > 0, G > 0, M > 0, 0 < Y < 2 and Γ(·) denotes the Gamma function. The drift b ∈ R is set to to comply with the drift condition (3.7) for martingale pricing.
Then for every z ∈ R, (T, b, C, G, M, Y ) → ϕ p 2 =(T,b,C,G,M,Y ) (z + iη) for the characteristic function ϕ p 2 of the CGMY model (4.12) is analytic on U (η). Table 4.2 displays for selected Lévy models conditions on the weight η ∈ R and the index α ∈ (1, 2] that guarantee (A3) and the upper bound in (A4) for a fixed model parameter constellation.

Class
(A3) and the upper bound in (A4) hold for η such that and α such that For the CGMY model with parameters C, G, M > 0 and Y < 2 it is obvious that (3.11) is statisfied for every η ∈ R with G > η and M > −η. Moreover, equation (4.12) in Glau (2015) shows that for Y ≤ 1, (3.12) is satisfied with β = 1.

Numerical Experiments
We apply the Chebyshev interpolation method to parametric option pricing considering a variety of option types in different well known option pricing models. Moreover, we conduct both an error analysis as well as a convergence study. The first focuses on the accuracy that can be achieved with a reasonable number of Chebyshev interpolation points. The latter confirms the theoretical order of convergence derived in Section 3, when the number of Chebyshev points increases. Finally, we study the gain in efficiency for selected multivariate options. We measure the numerical accuracy of the Chebyshev method by comparing derived prices with prices coming from a reference method. We employ the reference method not only for computing reference prices but also for computing prices at Chebyshev nodes Price p (k 1 ,...,k D ) with (k 1 , . . . , k D ) ∈ J during the precomputation phase of the Chebyshev coefficients c j , j ∈ J, in (2.9). Thereby, a comparability between Chebyshev prices and reference prices is maintained.
We implemented the Chebyshev method for applications with two parameters. To that extent we pick two free parameters p i 1 , p i 2 out of (3.2), 1 ≤ i 1 < i 2 ≤ D, in each model setup and fix all other parameters at reasonable constant values. We then evaluate option prices for different products on a discrete parameter grid (5.1) Once the prices have been derived on P, we compute the discrete L ∞ (P) and L 2 (P) error measures, , to interpret the accuracy of our implementation and of the Chebyshev method as such.

European Options
We consider a plain vanilla European call option on one asset as well as a European digital down&out option, first. Both derivatives have been introduced in Table 4.1. For these products we investigate the performance of the Chebyshev interpolation method for the Heston model and the Lévy models of Black&Scholes, Merton and the CGMY model. We keep the strike parameter constant, p 1 = k = log(K) for K = 1. As previously discussed in Remark 3.8, this is no restriction of generality. We also disregard interest rates, setting r = 0. For the three Lévy models we vary the maturity T (in years) as well as the option moneyness S 0 /K. All three models fall within the scope of Corollary 3.8, where the error is analyzed explicitly. Thus we expect Model fixed parameters free parameters p 1 p 2 p 1 p 2 (sub)exponential convergence for the three Lévy models. For the Heston model we vary S 0 /K and let v 0 as one of the model parameters float. Due to the analyticity of the Fourier transform of the payoff in S 0 /K and of the characteristic function of the process in v 0 , compare Section 4.2.4, we expect this convergence for the Heston model as well. A detailed overview of the realistically chosen parametrization is given by Table 5.1. For numerical integration in Fourier pricing we use Matlab's quadgk routine over the interval [0, ∞) with absolute precision bound of ε < 10 −14 .
The first question we address concerns the achievable accuracy with a fixed number of Chebyshev polynomials. We set N 1 = N 2 = 10 and precompute the Chebyshev coefficients as defined in (2.9) with D = 2 using the parametrization of Table 5.1 for the models therein. We evaluate the resulting polynomial over a parameter grid of dimension D = 2 and compute the approximate European option prices in each node. As a comparison, we also compute the respective Fourier price via numerical integration of the accordingly parametrized integrand in (3.4). Figure 5.1 shows the results for the European call option. The accuracy achieved by N = N 1 = N 2 = 10 shows a significant spread over the four different models but reaches very satisfying error levels from 10 −7 to 10 −10 . Increasing the number of Chebyshev points further improves the accuracy. Since at its core the implementation of the Chebyshev method consists of summing up matrices, this refinement comes at virtually no additional cost.
In the same parametrization setting we price a European digital down&out option. While a call payoff profile is not differentiable but at least continuous, the digital payoff function is not even continuous, compare Table 4.1. This reduction in smoothness of the payoff function reduces the accuracy of the interpolation p → Price (p) as well. We compare the Chebyshev interpolation again with  We compare the Chebyshev interpolation with N = N1 = N2 = 10 to classic Fourier pricing by numerical integration. The parametrization of the models and the option has been chosen according to Table 5.1. We observe that the achieved accuracy varies significantly between the models. The order of accuracy of Black&Scholes prices is identical to the accuracy for the CGMY case. The Merton model falls behind by one order of magnitude, while the Heston model shows a very strong fit between prices and their approximation by the Chebyshev method. N 1 = N 2 = 10 to classic Fourier pricing by numerical integration. The parametrization of the models and the option has been chosen again according to Table 5.1, where K = 1 now denotes the parameter of the digital option. Figure 5.2 shows the results. Comparing the results of the call option pricing in Figure 5.1 to the results of the digital option pricing in Figure 5.2 we observe that the accuracy achieved by N 1 = N 2 = 10 is reduced by a factor 10 1 to 10 2 . This demonstrates how little the reduced smoothness of the payoff profile effects the accuracy of the Chebyshev method in these cases. Furthermore we conduct an empirical convergence study for this very same setting of option and model parametrization. For an increasing degree N = N 1 = N 2 , the Chebyshev polynomial is set up and prices over a parameter grid of structure (5.19) are computed. Again, Fourier pricing serves as a comparison. For each N ∈ {1, . . . , 30}, the error measures ε L ∞ and ε L 2 , defined by (5.2) on the discrete parameter grid P defined in (5.19), are evaluated. We observe exponential convergence for all four models. Figure 5.3 shows the decay for the European call option while Figure 5.4 does so for the European digital down&out option.
Defining ζ = T +T T −T and setting ̺ = ζ + ζ 2 − 1, the theoretical convergence analysis predicts a slope of the error decays in both Figure 5.3 as well as Figure 5.4 of at

Basket and Path-dependent Options
In this section we use the Chebyshev method to price basket and path-dependent options. First, we apply the method to interpolate Monte-Carlo estimates of prices of financial products and check the resulting accuracy. To this aim we exemplarily choose basket, barrier and lookback options in 5-dimensional Black&Scholes, Heston and Merton models. Second, we combine the Chebyshev method with a Crank-Nicolson finite difference solver with Brennan Schwartz approximation, see Brennan and Schwartz (1977), for pricing a univariate American put option in the Black&Scholes model.
In our Monte-Carlo simulation we use 10 6 sample paths, antithetic variates as variance reduction technique and 400 time steps per year. The error of the Monte-Carlo method cannot be computed directly. We thus turn to statistical error analysis and use the well-known 95% confidence bounds to determine the accuracy. Following the assumption of a normally distributed Monte-Carlo estimator with mean equal to the estimator's value and variance equal to the empirical variance of the payoff on the Monte-Carlo samples these bounds are derived. The confidence bounds then yield a range around the mean that includes the true price with 95% probability. We pick two free parameters p i 1 , p i 2 out of (3.2), 1 ≤ i 1 < i 2 ≤ D, in each model setup and fix all other parameters at reasonable constant values. In this section we define the discrete parameter grid P ⊆ and call P test grid. On this test grid the largest confidence bound is 0.025 an on average lees than 0.013. For the finite difference method we find that the absolute error between numerical approximation and option price is below 0.005 on all computed parameter tuples in P. This error bound was computed by comparing each approximation to the limit of the sequence of finite difference approximations with increasing grid size. In our calculations we work with a grid size in time as well as in space (log-moneyness) of 50 · max{1, T } and compared the result to the resulting prices using grid sizes of 1000 · max{1, T }. This grid size has been determined as sufficient for the limit due to hardly any changes compared to grid sizes of 500 · max{1, T }.
Here, our main concern is the accuracy of the Chebyshev interpolation when we vary for each option the parameters strike and maturity analogously to the previous section. For N ∈ {5, 10, 30}, we precompute the Chebyshev coefficients as defined in (2.9) with D = 2 where always N 1 = N 2 = N . An overview of fixed and free parameters in our model selection is given in Table 5.2. For computational simplicity in the Monte-Carlo simulation, we assume uncorrelated underlyings.
Let us briefly define the payoffs of the multivariate basket and path-dependent options. The payoff profile of a basket option for d underlyings is given as We denote S t = (S 1 t , . . . , S d t ), S j T := min 0≤t≤T S j t and S j T := max 0≤t≤T S j t . A lookback option for d underlyings is defined as As a multivariate barrier option on d underlyings we define the payoff For an American put option the payoff is the same as for a European put,  but the option holder has the right to exercise the option at any time t up to maturity T . We now turn to the results of our numerical experiments. In order to evaluate the accuracy of the Chebyshev interpolation we look for the worst case error ε L ∞ . The absolute error of the Chebychev interpolation method can be directly computed by comparing the interpolated option prices with those obtained by the reference numerical algorithm i.e. either the Monte-Carlo or the Finite Difference method. 4.298 · 10 −3 6.6358 2.309 · 10 −2 6.6315 Table 5.5: Interpolation of exotic options with Chebyshev interpolation. N = 30 and d = 5 in all cases. In addition to the L ∞ errors the table displays the Monte-Carlo (MC) prices, the Monte-Carlo confidence bounds and the Chebyshev Interpolation (CI) prices for those parameters at which the L ∞ error is realized.

Model Option
Since the Chebychev interpolation matches the reference method on the Chebychev nodes, we will use the out-of-sample test grid as in (5.3). Table 5.3 shows the numerical results for the basket and path-dependent options for N = 5, Table 5.4 for N = 10 and Table 5.5 for N = 30. In addition to the L ∞ errors the tables display the Monte-Carlo (MC) prices, the Monte-Carlo confidence bounds and the Chebyshev Interpolation (CI) prices for those parameters at which the L ∞ error is realized.
The results show that for N = 30 the accuracy is for all selected options at a level of 10 −3 . We see that the Chebyshev interpolation error is dominated by the Monte-Carlo confidence bounds to a degree which renders it negligible in a comparison between the two. For basket and barrier options the L ∞ error already reaches satisfying levels of order 10 −3 at N = 10 already. Again, the Chebyshev approximation falls within the confidence bounds of the Monte-Carlo approximation. Thus, Chebyshev interpolation with only 121 = (10+1) 2 nodes suffices for mimicking the Monte Carlo pricing results. This statement does not hold for lookback options, where the L ∞ error still differs noticeably when comparing N = 10 to N = 30. As can be seen from Table 5.3 Chebyshev interpolation with N = 5 may yield unreliable pricing results. For lookback options in the Heston model we even observe negative prices in individual cases. Chebyshev pricing of American options in the N ε L ∞ FD price CI price 5 3.731 · 10 −3 1.9261 1.9224 10 1.636 · 10 −3 12.0730 12.0746 30 3.075 · 10 −3 6.3317 6.3286 Table 5.6: Interpolation of one-dimensional American puts with Chebyshev interpolation in the Black&Scholes model. In addition to the L ∞ errors the table displays the Finite Differences (FD) prices and the Chebyshev Interpolation (CI) prices for those parameters at which the L ∞ error is realized.
Black&Scholes model is even more accurate as illustrated in Table 5.6. Here, already for N = 5 the accuracy of the reference method is achieved. We conclude that the Chebyshev interpolation is highly promising for the evaluation of multivariate basket and path-dependent options. Yet the accuracy of the interpolation critically depends on the accuracy of the reference method at the nodal points which motivates further analysis that we perform in the subsequent subsection.

Interaction of Approximation Errors at Nodal Points and Interpolation Errors
The Chebyshev method is most promising for use cases, where computationally intensive pricing methods are required. Then, for computing the prices at the Chebyshev nodes in order to set up the interpolation, the issue of distorted prices at the Chebyshev nodes and their consequences rises naturally. The observed noisy prices at the Chebyshev nodes are where ε p (k 1 ,...,k D ) is the approximation error introduced by the underlying numerical technique at the Chebyshev nodes. Due to linearity, the resulting interpolation is of the form with the error function with the coefficients c ε j for j = (j 1 , . . . , j D ) ∈ J given by If ε p (k 1 ,...,k D ) ≤ ε for all Chebyshev nodes p (k 1 ,...,k D ) , we obtain since the Chebyshev polynomials are bounded by 1. This yields the following remark.
In this example, the accuracy of the reference method has to reach a level of 10 −5 to guarantee an overall error of order 10 −3 . This demonstrates a trade-off between increasing N 1 and N 2 compared to the accuracy of the reference method. The error bound above is rather conservative. Our experiments from the previous section suggest that this bound highly overestimates the errors empirically observed. However, the presented error bound from Remark 5.1 can guarantee a desired accuracy by determining an adequate number of Chebyshev nodes and the corresponding accuracy of the reference method used at the Chebyshev nodes. For practical implementation we suggest the following procedure. For a prescribed accuracy the N i , i = 1, . . . , D, can be determined from the first term in (5.8) by choosing N i , i = 1, . . . , D, as small as possible such that the prescribed accuracy is attained. Accordingly, the accuracy that the reference method needs to achieve is bounded by the second term. A very accurate reference method in combination with small N i , i = 1, . . . , D, promise best results. With this rule of thumb in mind the experiments of Section 5.3.2 below have been conducted.

Study of the gain of efficiency
In the previous section we investigated the accuracy of the Chebyshev polynomial interpolation method using Fourier, Monte-Carlo and Finite Difference as reference pricing methods. Finally, we investigate the gain in efficiency achieved by the method in comparison with Fourier pricing as well as in comparison to Monte Carlo pricing. In Section 5.3.1 we compute the results on a standard PC with an Intel i5 CPU, 2.50 GHz with cache size 3 MB. In Section 5.3.2 we used a PC with Intel Xeon CPU with 3.10 GHz with 20 MB SmartCache. All codes are written in Matlab R2014a.

Comparison to Fourier pricing
Here, we compare the method to Fourier pricing. We choose the pricing problem of a call option on the minimum of two assets as an example. Today's values of the underlying two assets are fixed at (5.9) S 1 0 = 1, Modeling the future development of the underlyings, (S j t ) t≥0 , j ∈ {1, 2}, we consider two bivariate models, separately. First, the two assets will be driven by the bivariate Black&Scholes model of Section 4.2.1. The bivariate Black&Scholes model is parametrized by a covariance matrix σ ∈ R 2×2 that we choose to be given by (5.10) σ 11 = 0.2 2 , σ 12 = 0.01, In a second efficiency study, asset movements follow the more involved bivariate Heston model in the version of Section 4.2.4 above for which we choose the parametrization In both cases we neglect interest rates, thus setting r = 0. The benchmark method, that is Fourier pricing, is evaluated using Matlab's quad2d routine. We prescribe an absolute and relative accuracy of at least 10 −8 from the integration result and integrate the Fourier integrand over the domain Ω = [−50, 50] × [0, 50], with a maximum number of 4000 function evaluations. We use Fourier integration with the same accuracy specifications at the Chebyshev nodes to set up the Chebyshev method for pricing based on strike K and maturity T as the two free parameters taking values in the intervals For a fair comparison, the number of Chebyshev polynomials is chosen such that Chebyshev interpolation prices yield an accuracy that matches the accuracy of the benchmark method resulting in (5.13) N BS Cheby = 11 and N Heston Cheby = 23, for the bivariate Black&Scholes model and the bivariate Heston model, respectively. Figure 5.5 illustrates the absolute errors over the whole K × T domain of interest between Fourier pricing and the Chebyshev method for both models, with the Chebyshev interpolator being based on N BS Cheby +1 polynomials in the Black&Scholes model case and N Heston Cheby + 1 polynomials in the Heston model case. In order to set up an efficiency study, we will compute sets of prices with increasing number of parameter tupels. To this end, when the offline phase of the Chebyshev method has been completed we compute 98 pricing surfaces, that is for Here, Chebyshev interpolation is based on N Heston Cheby + 1 = 24 Chebyshev polynomials. We achieve an absolute accuracy of order 10 −8 in both cases, thus matching the accuracy that the benchmark method Fourier pricing reaches.
each M ∈ {3, . . . , 100} we compute prices for all parameter tuples from Θ M defined by The computation time consumed by the Chebyshev offline phase is stored. Also, for each M ∈ {3, . . . , 100}, run-times for deriving all |Θ M | = M 2 prices are measured and stored for both routines, the Fourier pricing method and the Chebyshev interpolation algorithm. Figure 5.6 depicts these run-time measurements and Table 5.7 provides a second perspective.
In the Black&Scholes model case, the offline phase required T BS offline = 8 seconds for deriving option prices at all (N BS Cheby + 1) 2 = 144 Chebyshev nodes. The Heston model required T Heston offline = 101 seconds for the (N Heston Cheby + 1) 2 = 576 supporting prices. Taking this initial investment into account deems pricing with the Chebyshev method rather costly when only few option prices are derived after the offline phase has been completed. Yet, as Figure 5.6 shows, the increase in pricing speed that is achieved once the Chebyshev algorithm has been set up eventually outpaces Fourier pricing. From our experiments we conclude that the Chebyshev method already outruns its benchmark Fourier pricing in terms of total run-times when the number of prices to be computed exceeds (N BS Cheby + 1) 2 or (N Heston Cheby + 1) 2 , respectively. Additionally, Table 5.7 highlights that in both models for a total number of 50 2 parameter tuples, the Chebyshev method exhibits a significant decrease in (total) pricing run-times. For the maximal number investigated, i.e. 100 2 parameter tuples, pricing in the Black&Scholes model results in 95% of run-time savings and 90% runtime savings in the Heston model case in our implementation. This results from the fact the online phase in the Chebyshev method consist of computationally cheap polynomial evaluations and elementary assembling.   Figure 5.6. With increasing number of computed prices, the Chebyshev algorithm increasingly profits from the initial investment of the offline phase.

Comparison to Monte-Carlo pricing
In this section we choose a multivariate lookback option in the Heston model, based on 5 underlyings, as example. For the efficiency study we first vary one parameter, then we vary two. Variation of one model parameter For the multivariate lookback option in the Heston model the following parameters are fixed with j = 1, . . . , 5 as (5.15) S j 0 = 100, r = 0.005, K = 100, T = 1, As free parameter in the Chebyshev interpolation we pick the volatility of the volatility coefficient σ = σ j , j = 1, . . . , 5, The benchmark method is the Monte-Carlo pricing, again with 10 6 sample paths, antithetic variates as variance reduction technique and 400 time steps per year. We refer to this setting as benchmark setting. Following the discussion from Section 5.2.1, for the evaluation of the prices at the nodal points we guarantee a smallε by the Monte-Carlo method we enrich the Monte-Carlo setting to 5 · 10 6 sample paths, antithetic variates and 400 time steps per year. In Table 5 (5.17) The highest observed error on the test grid is at a level of 10 −2 . On the same test grid the benchmark Monte-Carlo setting has a worst case confidence bound of 1.644 · 10 −2 and by comparing the benchmark Monte-Carlo prices to the enriched Monte-Carlo prices on this test grid, the maximal absolute error is 7.361 · 10 −3 . Therefore, we conclude that the Monte-Carlo benchmark setting and the presented Chebyshev interpolation method have a roughly comparable accuracy and on the basis of this accuracy study we now turn to the comparison of run-times.
We compare run-times of the Chebyshev interpolation with N Heston Cheby = 6, in which the offline phase is based on the enriched Monte-Carlo setting, to the run-times of the Monte-Carlo benchmark setting described above. Table 5.9 provides the respective results. The results for M = 1 have been empirically measured, all others have been extrapolated from that since for each parameter set, the same amount of computation time would have to be invested. The table indicates that from M = 50 onwards interpolation by Chebyshev is faster. In Figure 5.7 we present additionally for each M = 1, . . . , 100 the run-times of the Chebyshev interpolation method, including the offline phase, compared to the Monte-Carlo method. Here we observe that for M = 35 both lines intersect and for Varying ε L ∞ MC price MC conf. bound CI price σ 9.970·10 −3 18.6607 4.592 · 10 −3 18.6707 Table 5.8: Interpolation of multivariate lookback options with Chebyshev interpolation for N = 6 based on an enriched Monte-Carlo setting with 5 · 10 6 sample paths, antithetic variates and 400 time steps per year. In addition to the L ∞ error on the test grid we also report the Monte-Carlo (MC) price, the Monte-Carlo confidence bound and the Chebyshev Interpolation (CI) price for those parameters at which the L ∞ error is realized. We observe that the accuracy of the Chebyshev interpolation for N = 6 is roughly in the same range as the accuracy of the benchmark Monte-Carlo setting (worst case confidence bound of 1.644 · 10 −2 and worst case error of 7.361 · 10 −3 ).
M > 35 the Chebyshev interpolation method is faster. Contrary to Section 5.3.1, the intersection of the two lines does not occur at M = N Heston Cheby + 1. This reflects the fact that in the offline phase for the Chebyshev interpolation a Monte-Carlo method with more sample paths has been used. Variation of two model parameters Now, we vary two parameters, we choose ρ j = ρ, j = 1, . . . , 5, and vary fixing all other parameters to the values of setting (5.15). Guaranteeing a roughly comparable accuracy between the Chebyshev interpolation method and the benchmark Monte-Carlo pricing, we use the following test grid P ⊆ [σ min , σ max ]×[ρ min , ρ max ], Varying ε L ∞ MC price MC conf. bound CI price σ, ρ 5.260 · 10 −2 5.239 1.428 · 10 −2 5.292 Table 5.10: Interpolation of multivariate lookback options with Chebyshev interpolation for N = 6 based on an enriched Monte-Carlo setting with 5 · 10 6 sample paths, antithetic variates and 400 time steps per year. In addition to the L ∞ error on the test grid we also report the Monte-Carlo (MC) price, the Monte-Carlo confidence bound and the Chebyshev Interpolation (CI) price for those parameters at which the L ∞ error is realized. We observe that the accuracy of the Chebyshev interpolation N = 6 is roughly in the same range as the accuracy of the benchmark Monte-Carlo setting (worst case confidence bound of 6.783 · 10 −2 and worst case error of 2.791 · 10 −2 ).
In Table 5.10 we present the accuracy results for the Chebyshev interpolation with N Heston Cheby = 6 based on the enriched Monte-Carlo setting. Comparing the benchmark Monte-Carlo setting and the enriched Monte-Carlo setting on this test grid, we observe that the maximal absolute error is 2.791 · 10 −2 and the confidence bounds of the benchmark Monte-Carlo setting do not exceed 6.783 · 10 −2 .
For a run-time comparison, we show for different values of M the run-times necessary to compute the prices for M 2 parameter tupels. Again, the run-times are measured for M = 1 and extrapolated for other values of M . Table 5.11 presents the results. In Figure 5.8 for each M = 1, . . . , 100 the run-times of the Chebyshev interpolation method, including the offline phase, compared to the Monte-Carlo method are presented. We observe that for M = 15 both lines intersect and for M > 15 the Chebyshev method outperforms its benchmark. Contrary to the case of varying only one parameter, the intersection of both lines occurs at a significantly lower value of M due to the fact that for each M pricing for M 2 parameter tupels is performed.
Additionally, Table 5.11 highlights that for a total number of 50 2 parameter tuples, the Chebyshev method exhibits a significant decrease in (total) pricing runtimes. For the maximal number of 100 2 parameter tuples that we investigated, pricing in either model resulted in more than 97% of run-time savings in our implementation. While the computation of 100 2 Heston prices via the Monte-Carlo method consumes up to 39 days, the Chebyshev method computes the very same prices in 23 hours only. Note that only 7 seconds of this time span are consumed by actual pricing during the online phase.

Conclusion and Outlook
This article focuses on applying the Chebyshev method to European option pricing. Within this scope, Sections 2-4 establish theoretical convergence results and numerical case studies in Section 5 confirm these findings. Moreover, in experiments with Fourier pricing an accuracy of 10 −5 is observed with less than ten Chebyshev nodes in each parameter, see Figure 5.3. The financial implications of the high precision achieved with such a small number of interpolation nodes are twofold. First, it shows that there are interesting cases in which we observe an accuracy in the range of machine precision. In this comfortable situation without methodological risk we can ignore the fact that approximations are implemented. Second, compared with other sources of risk, already errors from a much lesser accuracy level can be ignored. If we agree that an accuracy of 10 −4 is satisfactory, already 36-49 interpolation nodes for the approximation of call option prices as reported in Figure 5.3 are sufficient. Additionally, also the numerical experiments for American, barrier and lookback options display promising results. A theoretical error analysis for nonlinear pricing problems is beyond the scope of the present article, while we are convinced that further investigations in this direction are valuable. For instance our analytic approach based on examinations of the Fourier representations can be adopted to barrier options in Lévy models leading to the involvement of Wiener-Hopf factorizations, see Eberlein et al. (2011). In general we expect the regularity analysis to become more challenging in the presence of nonlinearities. For American options the current work of Teichmann (2015) may lead to regularity assertions for American options that are inherited from their corresponding European counterparts.
The theoretical and experimental results of our case studies show that the method can perform considerably well when few parameters are varied. As a consequence, we recommend the interpolation method for this case and also when solely the strike of a plain vanilla option is varied and fast Fourier methods are available. For calibration purposes for example, strikes are not given in a discrete logarithmic scale, which makes an additional interpolation necessary in order to apply FFT. Here, Chebyshev polynomials offer an attractive alternative. In particular, the maturity can be used as supplementary free variable. Moreover, for models with a low number of parameters, another path could be beneficial: Interpolating the objective function of the parameters directly. Then the optimization would boil down to a minimization of a tensorized polynomial, which could be exploited in further research. As may be seen from Armenti et al. (2015), where the present article is applied for the first time, this advantage can also be exploited for other optimization procedures in finance for example in risk allocation.
The multivariate construction of the interpolation and the theoretical error analysis suggest that the empirically observed error behaviour extends to three and more varying parameters, as well, as long as analyticity is provided. More precisely, in case of analyticity, the rate is of order ρ − D √ N for some constant ρ depending on the domain of analyticity and N the total number of interpolation nodes. For multivariate polynomial interpolation, the introduction of sparsity techniques promise higher efficiency, as for instance by compression techniques for tensors as reviewed by Kolda and Bader (2009). We address the issue of curse of dimensionality further in Gaß, Glau and Mair (2015), where we take a different route by replacing the Chebyshev interpolation with an empirical interpolation for Fourier pricing methods.

A Proof of Theorem 2.2
Proof. In Sauter and Schwab (2004, Proof of Lemma 7.3.3) the proof is given for the following error bound: where N is the number of interpolation points in each of the D dimensions, ̺ min := min D i=1 ̺ i and V the bound of f on B(P, ̺) with P = [−1, 1] D . Here, we extend the proof by incorporating the different values of N i , i = 1, . . . , D, as well as expressing the error bound with the different ̺ i , i = 1, . . . , D.
In general we work with a parameter space P of hyperrectangular structure, P = [p 1 , p 1 ] × . . . × [p D , p D ]. With the introduced linear transformation in Section 2.1 we have a transformation τ P : [−1, 1] D → P with .
Applying the error estimation from Sauter and Schwab (2004, Lemma 7.3.
Starting with (A.2), we first focus on ||E|| ̺ , At this step we apply Lemma A.1 and obtain Overall, using for x > 0, µ j ∈ N 0 and j = 1, . . . , D and this leads to From this point on we use the convergence of the geometric series since |̺ −2 j | < 1, j = 1, . . . , D, Recalling (A.2), we have to estimate f ̺ , From π D 2T 0 = 1 it directly follows that 1 2 ̺ = π D 2 2 T 0 2 ̺ = π D and therewith Combining the results leads to The following lemma shows that the Chebyshev interpolation of a polynomial with a degree as most as high as the degree of the interpolating Chebyshev polynomial is exact and furthermore determines an upper bound for interpolating Chebyshev polynomials with a higher degree. Proof. Uniqueness properties of the Chebyshev interpolation directly imply (A.3). The proof of (A.4) is similar to Sauter and Schwab (2004, Proof of Hilfssatz 7.3.1). They use the zeros of the Chebyshev polynomial as interpolation points, whereas we use the extreme points and therefore, we use a different orthogonality property in this proof. We first focus on the one-dimensional case. Recalling (2.2), the Chebyshev interpolation of T µ , µ > N , is given as where x k denotes the k-th extremum of T N . Here, we can apply the following orthogonality (Rivlin (1990, p.54)), For j ≤ N and µ > N this yields the existence of γ ≤ N such that I N (T µ ) = T γ . (A.6) (A.6) follows elementarily from the case that for any µ > N only for one 0 ≤ j ≤ N the orthogonality can lead to a coefficient c j > 0.
Proving the claim, we distinguish several cases. In all of these cases we assume that there exists 0 ≤ j ≤ N such that N k=0 ′′ T µ (x k )T j (x k ) = 0. We will then show that for all other 0 ≤ i ≤ N, i = j it follows N k=0 ′′ T µ (x k )T j (x k ) = 0.
Therewith, (A.6) holds and it directly follows that |T µ − I N (T µ )| ≤ |T µ | + |I N (T µ )| ≤ 1 + 1 = 2. Thus (A.4) holds in the one-dimensional case. The extension to the multivariate case follows analogously by applying the triangle inequality and inserting the one-dimensional result to each tensor component.