Optimal damping with a hierarchical adaptive quadrature for efficient Fourier pricing of multi-asset options in Lévy models

Christian Bayer; Chiheb Ben Hammouda; Antonis Papapantoleon

1 Introduction

Pricing multi-asset options, such as basket and rainbow options, is an interesting and challenging problem in quantitative finance because prices cannot be analytically computed in most cases; thus, efficient numerical methods are required. Moreover, despite the popularity of the Black–Scholes model, where the stock dynamics follow the geometric Brownian motion (GBM), Lévy models such as the variance Gamma (VG) (Madan and Seneta 1990) and normal inverse Gaussian (NIG) (Barndorff-Nielsen 1997b) have shown a better fit to empirical market behavior (Cont and Tankov 2003; Schoutens 2003) by accounting for market jumps in prices, semiheavy tails and high leptokurtosis.

Under the no-arbitrage assumption, option prices are given as expectations under an (equivalent) martingale measure and approximated using numerical integration methods. In this context, the prevalent numerical method is the Monte Carlo method (Glasserman 2003), which has a convergence rate insensitive to the input space dimensionality and payoff regularity, except for multilevel Monte Carlo methods (Bayer et al 2020b), where Lipschitz continuity is necessary to obtain optimal convergence rates. However, the convergence may be very slow, and we may not be able to exploit the available regularity structure to achieve better convergence rates. An alternative stream of research has emerged on multi-asset option pricing based on continuous-time Markov chain approximation (Mijatović and Pistorius 2013; Xi et al 2019; Kirkby et al 2020), which has been shown to outperform the Monte Carlo approach for options in two and three dimensions. However, pricing in more than three dimensions still poses a great numerical challenge for these approaches. Another class of methods relies on deterministic quadrature techniques whose performance is highly dependent on the input space dimension and integrand regularity. Some studies (see, for example, Bayer et al 2020a, 2023) have combined adaptivity, sparsification techniques and hierarchical representations (the Brownian bridge and the Richardson extrapolation) with quadrature methods to treat the high dimensionality effectively. Moreover, financial payoffs usually have low regularity; therefore, analytic and numerical smoothing techniques were introduced for better convergence (Bayer et al 2018, 2020a; Ben Hammouda 2020; Bayer et al 2023). All the aforementioned improvements were performed in the physical space.

In this work, we propose what we believe to be a novel approach for pricing European multidimensionalbasket and rainbow options under multivariate GBM and Lévy models.¹¹ 1 Rainbow options (Margrabe 1978) are appealing to investors because they offer reduced risk exposure to the market at a cheap cost by betting more on the performance of an individual stock in a group of stocks than the overall performance of a portfolio of stocks when considering basket options (see, for example, Guillaume 2008). Unlike the abovementioned approaches, we recover the high regularity of the integrand by mapping the problem from the physical space to the frequency space, when the Fourier transforms of the payoff and density are well defined and known explicitly. Moreover, when designing our method, we effectively treat two key aspects affecting the numerical complexity: the choice of the damping parameters that ensure integrability and control the regularity class of the integrand; and the high dimensionality of the integration problem. Based on the extension of the one-dimensional (1D) Fourier valuation formula (Raible 2000; Lewis 2001) to the multivariate case, we first smooth the Fourier integrand via an optimal choice of the damping parameters based on a proposed optimization rule. Then we use adaptive sparse grid quadrature (ASGQ) based on sparsification and dimension-adaptivity techniques to accelerate the numerical quadrature convergence in high dimensions.

Fourier-based pricing methods (Carr and Madan 1999; Raible 2000; Lewis 2001; Duffie et al 2003; Fang and Oosterlee 2009; Lord et al 2008; Leentvaar and Oosterlee 2008; Kwok et al 2012; Baschetti et al 2022) map the original problem to the frequency space and obtain the solution in the physical space using the Fourier inversion theorem. The approximation of the resulting integral is performed numerically using direct integration (DI) methods or the fast Fourier transform (FFT). The common ingredient for these approaches is the explicit knowledge of the characteristic function (ie, the Fourier transform of the probability density function) corresponding to the price dynamics.

There are three main popular Fourier valuation approaches. In the first approach, originally proposed by Carr and Madan (1999) (see also Lee 2004; Carr and Wu 2004), a Fourier transform is applied in the log-strike variable, $k$ . Hence, for fixed maturity $T$ , the whole curve of option prices, $C(T,\cdot)$ , is computed. To ensure the existence of the Fourier transform, we must multiply the pricing function by a damping factor with respect to the strike parameter. This method is appropriate for 1D problems; however, extending it to the multi-asset option pricing context is difficult. The strike price is not defined for all stocks, whereas the multivariate density depends on all the underlyings. Moreover, the derivations must be performed separately for each payoff and stock dynamics.

The second approach (Fang and Oosterlee 2009; Ruijter and Oosterlee 2012; Zhang and Oosterlee 2013), called the COS method, relies on the Fourier cosine series expansion of the density function, and relates the cosine coefficients to the characteristic function. Although the COS method has been shown to be efficient at handling 1D and 2D problems, it is still challenging to generalize this class of methods to the multidimensional setting, for a number of reasons. First, in general, the cosine series coefficients of the payoff function are not known analytically; therefore, they need to be recovered numerically by evaluating high-dimensional integrals using the Clenshaw–Curtis quadrature or the discrete cosine transform, as suggested by Ruijter and Oosterlee (2012). Second, even though this approach does not introduce damping parameters to ensure integrability, the truncation parameters of the integration domain must be determined. Junike and Pankrashkin (2022) showed that there exist cases where the method fails to converge to the correct price if these parameters are chosen based on the rule of thumb for cumulants suggested by Fang and Oosterlee (2009) and Ruijter and Oosterlee (2012). To circumvent this issue, Junike and Pankrashkin propose an alternative truncation heuristic for 1D cases, but a practical choice in a high-dimensional setting remains a challenging open problem. To avoid determining an a priori truncation range (but at a higher computational cost), Colldeforns-Papiol et al (2017) replaced the Fourier cosine expansion by expressing the density as a finite combination of Shannon wavelet scaling functions, allowing for adaptive estimation of the truncation range. Finally, the number of Fourier cosine series coefficients required for the density expansion grows exponentially with the number of underlyling assets, as pointed out by Chau and Oosterlee (2019).

Given a characteristic function and Fourier transform of the payoff function, the third approach (Raible 2000; Lewis 2001; Hurd and Zhou 2010; Eberlein et al 2010) uses a highly modular pricing framework. This method separates the underlying stochastic process from the derivative payoff using the Plancherel–Parseval theorem and uses the generalized inverse Fourier transform to recover the option price. In addition, it introduces damping parameters with respect to the stock price variables to ensure integrability, which shifts the integration contour to a line parallel with the real axis in the complex space in order to avoid singularities. This technique is easier to extend to the multivariate case than the aforementioned approaches of Carr and Madan (1999), Lee (2004) and Carr and Wu (2004) and Fang and Oosterlee (2009), Ruijter and Oosterlee (2012) and Zhang and Oosterlee (2013).

Regarding the use of the Parseval-based Fourier valuation approach (as in the works of Raible (2000), Lewis (2001), Hurd and Zhou (2010) and Eberlein et al (2010)), there has, to the best of our knowledge, been no precise analysis of the effect of the damping parameters on the convergence speed of quadrature methods or guidance on choosing them to improve the numerical performance, particularly in the multivariate setting. Previous works have made arbitrary choices for the damping parameter, and only Lord and Kahl (2007) and Kahl and Lord (2012) studied the damping parameter selection for the first type of Fourier valuation approach (as in the works of Carr and Madan (1999), Lee (2004) and Carr and Wu (2004)) in the 1D setting to obtain robust integration behavior.

In this paper, when pricing basket and rainbow options under the multivariate GBM and Lévy models, and based on the extension of the 1D Fourier valuation formula (Raible 2000; Lewis 2001) to the multivariate case, we demonstrate that the choice of damping parameters has a significant impact on the speed of convergence of the numerical quadrature. In addition, motivated by error estimates based on contour integration tools, we propose a general framework for the optimal choice of the damping parameters, which can be tailored and extended to various pricing models, resulting in a smoother integrand and improving the efficiency of the numerical quadrature. Based on this proposed rule, the vector of the optimal damping parameters can be obtained by solving a simple optimization problem. Moreover, we demonstrate the consistent advantage of the optimal damping rule through numerical examples with different dimensions and parameter constellations.

The numerical evaluation of the resulting inverse Fourier integral can be performed using the FFT algorithm (Carr and Madan 1999; Crisóstomo 2018), which can be faster than DI methods because it exploits periodicities and symmetries. However, it cannot satisfy the requirement for matching the pricing algorithm to the structure of the market data and must be assisted by interpolation and extrapolation methods for the smile surface. This is in contrast to DI methods, which allow for flexible strikes.²² 2 We refer the reader to Zhu (2009, Chapter 4) and Lord and Kahl (2007) for further comparisons of FFT and DI. A further downside of the FFT method is that it has an additional truncation error and requires the determination of the upper and lower truncation parameters of the integral. This task is nontrivial for multidimensional integrals because the speed of decay to zero of the integrand depends on the damping factors, which are unknown a priori, creating dependence between the truncation and damping parameters. In this work, we opt for the DI approach combined with an unbounded quadrature (the Gaussian quadrature rule) to evaluate the option price. This approach can be efficiently vectorized, allowing for a faster calibration procedure. The optimal choice of the damping and truncation parameters for FFT when pricing multi-asset options remains an open problem for future work.

Through an extensive numerical study on basket and rainbow options under the multivariate GBM, VG and NIG models, we demonstrate the advantages of the dimension-adaptive quadrature and our rule for choosing the damping parameters on the numerical complexity of the quadrature methods for approximating the Fourier valuation integrals. Moreover, we illustrate cases where our approach outperforms the COS method. Finally, we show that our approach achieves substantial computational gains over the Monte Carlo method for different dimensions and parameter constellations.

The remainder of the paper is structured as follows. In Section 2, we introduce the proposed pricing framework in the Fourier space and the multivariate valuation formula. In Section 3, we explain our methodology: Section 3.1 states our heuristic rule for choosing the damping parameters and provides the motivation behind it; Section 3.2 details the optimal damping rule; and Section 3.3 presents the different hierarchical deterministic quadrature methods used for numerically evaluating the inverse Fourier integrals of interest. Finally, in Section 4, we report and analyze the obtained results, illustrating the advantages of the proposed approach and highlighting the considerable computational gains achieved over the COS and Monte Carlo methods. In Section 5, we state our conclusions.

2 Problem setting and pricing framework

In Section 2.1, we introduce the general Fourier valuation framework for the multi-asset options that we consider in this work. Then, in Section 2.2, we present specific details on the type of options and models investigated in this study.

2.1 The multivariate Fourier pricing valuation formula

We aim to efficiently price European multi-asset options (eg, basket/rainbow options) where the asset’s dynamics follow a certain multivariate stochastic model (eg, the Lévy model). We extend the 1D representation (Lewis 2001) to derive the pricing valuation formula in the Fourier space for the multivariate setting. Before proving the general valuation formula stated in Proposition 2.3, we introduce some necessary notation, definitions and assumptions.

Notation and definitions.

•
We define $\smash{\bm{X}_{t}:=(X_{t}^{1},\dots,X_{t}^{d})}$ , a $d$ -dimensional vector of log asset prices at time $t$ , $0\leq t\leq T$ , where $T$ is the maturity time of the option and $X^{i}_{t}:=\log(S^{i}_{t})$ , $i=1,\dots,d$ , $\{S^{i}_{t}\}_{i=1}^{d}$ are the prices of the underlying assets at time $t$ . The dynamics of $\bm{X}_{t}$ follow a multivariate stochastic model with parameters denoted by the vector $\bm{\varTheta}_{m}$ . We denote by $\rho_{\bm{X}_{t}}(\cdot)$ the corresponding risk-neutral conditional transition probability density function.
•
For $\bm{z}\in\mathbb{C}^{d}$ , $\smash{\varPhi_{\bm{X}_{T}}(\bm{z})}:=\smash{\mathbb{E}_{\rho_{\bm{X}_{T}}}[% \mathrm{e}^{\mathrm{i}\langle\bm{z},\bm{X}_{T}\rangle}]}$ denotes the joint characteristic function of $\bm{X}_{T}$ extended to the complex plane, and $\langle\cdot,\cdot\rangle$ denotes the inner product in $\mathbb{R}^{d}$ extended to $\mathbb{C}^{d}$ ; ie, for $\bm{y}$ , $\bm{z}\in\mathbb{C}^{d}$ , $\smash{\langle\bm{y},\bm{z}\rangle=\sum_{j=1}^{d}y_{j}z_{j}}$ .
•
For $\bm{x}\in\mathbb{R}^{d}$ , $P(\bm{x})$ denotes the payoff function. For $\bm{z}\in\mathbb{C}^{d}$ ,

$\hat{P}(\bm{z}):=\int_{\mathbb{R}^{d}}\mathrm{e}^{-\mathrm{i}\langle\bm{z},\bm% {x}\rangle}P(\bm{x})\mathrm{d}\bm{x}$

is the Fourier transform of $P(\cdot)$ , and $P_{\bm{R}}(\bm{x}):=\mathrm{e}^{\langle\bm{R},\bm{x}\rangle}P(\bm{x})$ is the dampened payoff function.
•
Let $\bm{\varTheta}_{p}=(K,T,r)$ be the vector of payoff and market parameters, with $K$ being the strike price and $r$ the deterministic interest rate.
•
We denote by $\mathrm{i}$ the unit imaginary number and by $\Re[\cdot]$ and $\Im[\cdot]$ the real and imaginary parts of a complex number, respectively.
•
Let $L_{\mathrm{bc}}^{1}(\mathbb{R}^{d})$ denote the space of bounded and continuous functions in $L^{1}(\mathbb{R}^{d})$ .
•
We denote by $\bm{I}_{d}$ the $d\times d$ identity matrix.

Assumption 2.1 (Assumptions on the payoff).

(1)
The payoff function, $\bm{x}\mapsto P(\bm{x})$ , is continuous for all $\bm{x}\in\mathbb{R}^{d}$ .
(2)
There exists $\bm{R}\in\delta_{P}:=\{\bm{R}\in\mathbb{R}^{d}\mid\bm{x}\mapsto P_{\bm{R}}(\bm% {x})\in L_{\mathrm{bc}}^{1}(\mathbb{R}^{d}),\ \bm{u}\mapsto\hat{P}(\bm{u}+% \mathrm{i}\bm{R})\in L^{1}(\mathbb{R}^{d})\}$ , where $\delta_{P}$ is the strip of regularity (analyticity) of the payoff’s Fourier transform.

Assumption 2.2.

(Assumption on the model and the corresponding characteristic function)

(1)
There exists $\bm{R}\in\delta_{X}:=\{\bm{R}\in\mathbb{R}^{d}\mid\bm{u}\mapsto\varPhi_{\bm{X}% _{T}}(\bm{u}+\mathrm{i}\bm{R})\text{~{}exists and~{}}|\varPhi_{\bm{X}_{T}}(\bm% {u}+\mathrm{i}\bm{R})|<\infty~{}\forall\bm{u}\in\mathbb{R}^{d}\}$ , where $\delta_{X}$ is the strip of regularity (analyticity) of the extended characteristic function.

Proposition 2.3.

(Multivariate Fourier pricing valuation formula) We use the definitions above and suppose that Assumptions 2.1 and 2.2 hold and that $\delta_{V}=\delta_{X}\cap\delta_{P}\neq\emptyset$ . Then, for $\bm{R}\in\delta_{V}$ , the option value is given by

$V(\bm{\varTheta}_{m},\bm{\varTheta}_{p})=(2\pi)^{-d}\mathrm{e}^{-rT}\Re\bigg{[% }\int_{\mathbb{R}^{d}}\varPhi_{\bm{X}_{T}}(\bm{u}+\mathrm{i}\bm{R})\hat{P}(\bm% {u}+\mathrm{i}\bm{R})\mathrm{d}\bm{u}\bigg{]}.$

(2.1)

Proof.

Given Assumption 2.1 and using the Fourier inversion theorem (see Hörmander 2015, Chapter 7), we have

$P_{\bm{R}}(\bm{x})=(2\pi)^{-d}\int_{\mathbb{R}^{d}}\mathrm{e}^{\mathrm{i}% \langle\bm{u},\bm{x}\rangle}\hat{P}_{\bm{R}}(\bm{u})\mathrm{d}\bm{u},\quad\bm{% R}\in\delta_{P},~{}\bm{x}\in\mathbb{R}^{d}.$

(2.2)

Moreover, we have that

	$\displaystyle\hat{P}_{\bm{R}}(\bm{u})$	$\displaystyle=\int_{\mathbb{R}^{d}}\mathrm{e}^{-\mathrm{i}\langle\bm{u},\bm{x}% \rangle}\mathrm{e}^{\langle\bm{R},\bm{x}\rangle}P(\bm{x})\mathrm{d}\bm{x}$
		$\displaystyle=\int_{\mathbb{R}^{d}}\mathrm{e}^{-\mathrm{i}\langle\bm{u}+% \mathrm{i}\bm{R},\bm{x}\rangle}P(\bm{x})\mathrm{d}\bm{x}=\hat{P}(\bm{u}+% \mathrm{i}\bm{R}),\quad\bm{u}\in\mathbb{R}^{d},~{}\bm{R}\in\delta_{P},$		(2.3)

where $\hat{P}(\bm{u}+\mathrm{i}\bm{R})$ , sometimes called the generalized Fourier transform (Titchmarsh 1948) or the Fourier–Laplace transform (Hörmander 2015), is a holomorphic extension of the Fourier transform to the horizontal strips (tubes) $\bm{z}\in\mathbb{R}^{d}+\mathrm{i}\delta_{P}\subset\mathbb{C}^{d}$ in the complex domain (Carlsson and Wittsten 2017). Using (2.3), we can write (2.2) as

$P(\bm{x})=\Re\bigg{[}(2\pi)^{-d}\mathrm{e}^{-\langle\bm{R},\bm{x}\rangle}\int_% {\mathbb{R}^{d}}\mathrm{e}^{\mathrm{i}\langle\bm{u},\bm{x}\rangle}\hat{P}(\bm{% u}+\mathrm{i}\bm{R})\mathrm{d}\bm{u}\bigg{]},\quad\bm{x}\in\mathbb{R}^{d},~{}% \bm{R}\in\delta_{P}.$

(2.4)

Then, using (2.4), Assumption 2.2 and Fubini’s theorem, we obtain that

	$\displaystyle V(\bm{\varTheta}_{m},\bm{\varTheta}_{p})$	$\displaystyle=\mathrm{e}^{-rT}\mathbb{E}_{\rho_{\bm{X}_{T}}}[P(\bm{X}_{T})]$
		$\displaystyle=(2\pi)^{-d}\mathrm{e}^{-rT}\mathbb{E}_{\rho_{\bm{X}_{T}}}\bigg{[% }\Re\bigg{[}\mathrm{e}^{-\langle\bm{R},\bm{X}_{T}\rangle}\int_{\mathbb{R}^{d}}% \mathrm{e}^{\mathrm{i}\langle\bm{u},\bm{X}_{T}\rangle}\hat{P}(\bm{u}+\mathrm{i% }\bm{R})\mathrm{d}\bm{u}\bigg{]}\bigg{]},\quad\bm{R}\in\delta_{P},$
		$\displaystyle=(2\pi)^{-d}\mathrm{e}^{-rT}\Re\bigg{[}\int_{\mathbb{R}^{d}}% \mathbb{E}_{\rho_{\bm{X}_{T}}}[\mathrm{e}^{\mathrm{i}\langle\bm{u}+\mathrm{i}% \bm{R},\bm{X}_{T}\rangle}]\hat{P}(\bm{u}+\mathrm{i}\bm{R})\mathrm{d}\bm{u}% \bigg{]},\quad\bm{R}\in\delta_{V}:=\delta_{P}\cap\delta_{X},$
		$\displaystyle=(2\pi)^{-d}\mathrm{e}^{-rT}\Re\bigg{[}\int_{\mathbb{R}^{d}}% \varPhi_{\bm{X}_{T}}(\bm{u}+\mathrm{i}\bm{R})\hat{P}(\bm{u}+\mathrm{i}\bm{R})% \mathrm{d}\bm{u}\bigg{]}.$

The application of Fubini’s theorem is justified by taking $\bm{R}\in\delta_{P}$ to enforce $\hat{P}(\bm{u}+\mathrm{i}\bm{R})\in L^{1}(\mathbb{R}^{d})$ and taking $\bm{R}\in\delta_{X}$ to ensure that $\varPhi_{\bm{X}_{T}}(\mathrm{i}\bm{R})$ exists and is bounded. ∎

In what follows, we define the integrand of interest from (2.1) as

$g(\bm{u};\bm{R},\bm{\varTheta}_{m},\bm{\varTheta}_{p}):=(2\pi)^{-d}\mathrm{e}^% {-rT}\Re[\varPhi_{\bm{X}_{T}}(\bm{u}+\mathrm{i}\bm{R})\hat{P}(\bm{u}+\mathrm{i% }\bm{R})],\quad\bm{u}\in\mathbb{R}^{d},~{}\bm{R}\in\delta_{V}.$

(2.5)

Remark 2.4.

(Connection to the valuation formula of Eberlein et al (2010)) The notation we used for the definition of the Fourier transform and the payoff dampening is the same as that given by Hurd and Zhou (2010) but with a different set of assumptions. The valuation formula of Eberlein et al (2010, Theorem 3.2) can easily be recovered from (2.1) by considering $-\bm{R}$ instead of $\bm{R}$ and $\smash{\hat{P}(-\bm{z})}$ instead of $\smash{\hat{P}(\bm{z})}$ , and by using the relation $\smash{\varPhi_{\bm{X}_{T}}(\bm{u})=M_{\bm{X}_{T}}(\mathrm{i}\bm{u})}$ , where $\smash{M_{\bm{X}_{T}}(\cdot)}$ denotes the moment generating function of $\bm{X}_{T}$ .

Remark 2.5 (The case of discontinuous payoffs).

In general, the regularity assumptions on the payoff, such as its continuity, can be compensated by more regularity assumptions on the model. However, in the particular case of European options and the Lévy models considered, the continuity condition in Assumption 2.1 can be dropped because the processes considered possess a Lebesgue density. We refer the reader to Eberlein et al (2010) for more details.

2.2 Payoffs and multivariate asset models

2.2.1 Payoffs and their Fourier transforms

In this paper, we focus on two specific examples of payoffs: the basket put,³³ 3 To simplify the presentation, we consider the unweighted basket put payoff. The generalization to the weighted case as presented in Section 4 can be done straightforwardly by considering $X_{t}^{i}=\log(S^{i}_{t}/w_{i})$ , $i=1,\dots,d$ . given for $K>0$ by

$P(\bm{X}_{T})=\max\bigg{(}K-\sum_{i=1}^{d}w_{i}\mathrm{e}^{X_{T}^{i}},0\bigg{)% },\quad w_{i}>0\ \forall i=1\dots,d,$

(2.6)

and the call on min, given for $K>0$ by,

$P(\bm{X}_{T})=\max(\min(\mathrm{e}^{X_{T}^{1}},\dots,\mathrm{e}^{X_{T}^{d}})-K% ,0).$

(2.7)

The Fourier transforms for the payoffs in (2.6), (2.7) are given by

	$\displaystyle\hat{P}(\bm{z})$	$\displaystyle=K^{1-\mathrm{i}\sum_{j=1}^{d}z_{j}}\frac{\prod_{j=1}^{d}% \varGamma(-\mathrm{i}z_{j})}{\varGamma(-\mathrm{i}\sum_{j=1}^{d}z_{j}+2)},% \quad\bm{z}\in\mathbb{C}^{d},$
	$\displaystyle\Im[z_{j}]$	$\displaystyle>0\quad\forall j\in\{1,\dots,d\},$		(2.8)

for $w_{i}=1$ , $i=1,\dots,d$ , and

	$\displaystyle\hat{P}(\bm{z})$	$\displaystyle=\frac{K^{1-\mathrm{i}\sum_{j=1}^{d}z_{j}}}{(\mathrm{i}(\sum_{j=1% }^{d}z_{j})-1)\prod_{j=1}^{d}(\mathrm{i}z_{j})},\quad\bm{z}\in\mathbb{C}^{d},$
	$\displaystyle\Im[z_{j}]$	$\displaystyle<0\quad\forall j\in\{1,\dots,d\},\qquad\sum_{j=1}^{d}\Im[z_{j}]<-1,$		(2.9)

where $\smash{\varGamma(z)=\int_{0}^{+\infty}\mathrm{e}^{-t}t^{z-1}\mathrm{d}t}$ is the complex Gamma function defined for $z\in\mathbb{C}$ , with $\Re[z]>0$ . The regularity strips, $\delta_{P}$ , are shown in Table 1. The considered payoffs satisfy Assumption 2.1. We refer the reader to Eberlein et al (2010), Hubalek and Kallsen (2005) and Hurd and Zhou (2010) for further details on the derivation.

Table 1: Strip of regularity,

\delta_{P}

, of the payoff transforms.

Payoff	$\delta_{P}$
Basket put	$\displaystyle\{\bm{R}\in\mathbb{R}^{d},\ R_{i}>0~{}\forall i\in\{1,\dots,d\}\}$
Call on min	$\displaystyle\bigg{\{}\bm{R}\in\mathbb{R}^{d},\ R_{i}<0~{}\forall i\in\{1,% \dots,d\},~{}\sum_{i=1}^{d}R_{i}<-1\bigg{\}}$

2.2.2 Multivariate models and the corresponding characteristic functions

For the asset dynamics, we study the three models given by Examples 2.6–2.8.

Example 2.6 (Multivariate GBM).

The joint risk-neutral dynamics of the stock prices are modeled as follows:

$S_{i}(t)=S_{i}(0)\exp[(r-\tfrac{1}{2}\sigma_{i}^{2})t+\sigma_{i}W_{i}(t)],% \quad i=1,\dots,d,$

(2.10)

where $\sigma_{1},\dots,\sigma_{d}>0$ and $\{W_{1}(t),\dots,W_{d}(t),~{}t\geq 0\}$ are correlated standard Brownian motions with correlation matrix $\bm{C}\in\mathbb{R}^{d\times d}$ having components $(\bm{C})_{i,j}=\rho_{i,j}$ , where $-1\leq\rho_{i,j}\leq 1$ denotes the correlation between $W_{i}$ and $W_{j}$ . Moreover, $\smash{\bm{\varSigma}\in\mathbb{R}^{d\times d}}$ denotes the covariance matrix of the log returns, $\smash{\{\log(S_{i}(t)/S_{i}(0))\}_{i=1}^{d}}$ , with $\bm{\varSigma}_{ij}=\rho_{i,j}\sigma_{i}\sigma_{j}$ .

Example 2.7 (Multivariate VG (Luciano and Schoutens 2006)).

The joint risk-neutral dynamics of the stock prices are modeled as follows:

$S_{i}(t)=S_{i}(0)\exp[(r+\mu_{\mathrm{VG},i})t+\theta_{i}G(t)+\sigma_{i}\sqrt{% G(t)}W_{i}(t)],\quad i=1,\dots,d,$

(2.11)

where $\{W_{1}(t),\dots,W_{d}(t)\}$ are independent standard Brownian motions, and $\{G(t)\mid t\geq 0\}$ is a common Gamma process with parameters $(t/\nu,1/\nu)$ that is independent of all Brownian motions. In addition, $\theta_{i}\in\mathbb{R}$ and $\sigma_{i}>0$ , $1\leq i\leq d$ . The covariance matrix $\smash{\bm{\varSigma}\in\mathbb{R}^{d\times d}}$ satisfies $\varSigma_{i,j}=\sigma_{i}^{2}$ for $i=j$ , and equals 0 otherwise. Finally, $\smash{\bm{\mu}_{\mathrm{VG}}:=(\mu_{\mathrm{VG},1},\dots,\mu_{\mathrm{VG},d})}$ are the martingale correction terms that ensure $\{\mathrm{e}^{-rt}S_{i}(t)\mid t\geq 0\}$ is a martingale and are given by

$\mu_{\mathrm{VG},i}=\frac{1}{\nu}\log(1-\tfrac{1}{2}\sigma_{i}^{2}\nu-\theta_{% i}\nu),\quad i=1,\dots,d.$

(2.12)

Example 2.8.

(Multivariate NIG (Barndorff-Nielsen 1997a, 1977)) The joint risk-neutral dynamics of the stock prices are modeled as follows:

$S_{i}(t)=S_{i}(0)\exp\{(r+\mu_{\mathrm{NIG},i})t+\beta_{i}\mathrm{IG}(t)+\sqrt% {\mathrm{IG}(t)}W_{i}(t)\},\quad i=1,\dots,d,$

(2.13)

where $\{W_{1}(t),\dots,W_{d}(t)\}$ are independent standard Brownian motions, $\{\mathrm{IG}(t)\mid t\geq 0\}$ is a common inverse Gaussian process with parameters $\smash{(\delta^{2}t^{2},\alpha^{2}-\bm{\beta}^{\mathrm{T}}\bm{\Delta}\bm{\beta% })}$ and it is independent of all Brownian motions. In addition, $\smash{\alpha\in\mathbb{R}_{+}}$ , $\smash{\bm{\beta}\in\mathbb{R}^{d}}$ , $\smash{\alpha^{2}>\bm{\beta}^{\mathrm{T}}\bm{\Delta}\bm{\beta}}$ , $\delta>0$ and $\smash{\bm{\Delta}\in\mathbb{R}^{d\times d}}$ is a symmetric positive definite matrix with a unit determinant. Further, $\smash{\{\mu_{\mathrm{NIG},i}\}_{i=1}^{d}}$ are the martingale correction terms that ensure that $\smash{\{\mathrm{e}^{-rt}S_{i}(t)\mid t\geq 0\}}$ is a martingale, given by

$\mu_{\mathrm{NIG},i}=-\delta(\sqrt{\alpha^{2}-\beta_{i}^{2}}-\sqrt{\alpha^{2}-% (\beta_{i}+1)^{2}}),\quad i=1,\dots,d.$

(2.14)

Table 2 shows the expression of the characteristic function for each model in Examples 2.6–2.8: the regularity strips, $\delta_{X}$ , expressed in Table 3. The characteristic functions of the models studied satisfy Assumption 2.2. We refer the reader to Eberlein et al (2010) for further details.

Table 2: The expression of the characteristic function,

\varPhi_{\bm{X}_{T}}(\cdot)

, for the different pricing models. [Here,

\bm{1}_{\mathbb{R}^{d}}

is the

d

-dimensional unit vector.]

Model	$\varPhi_{\bm{X}_{T}}(\bm{z})$
GBM	$\displaystyle\exp(\mathrm{i}\langle\bm{z},\bm{X}_{0}\rangle)$
	$\displaystyle\quad{}\times\exp\bigg{(}\mathrm{i}\langle\bm{z},r\bm{1}_{\mathbb% {R}^{d}}-\frac{1}{2}\operatorname{diag}(\bm{\varSigma})\rangle T-\frac{T}{2}% \langle\bm{z},\bm{\varSigma}\bm{z}\rangle\bigg{)}$ ,	$\bm{z}\in\mathbb{C}^{d},\quad\Im[\bm{z}]\in\delta_{X}$
VG	$\displaystyle\exp(\mathrm{i}\langle\bm{z},\bm{X}_{0}\rangle)$
	$\displaystyle\quad{}\times\exp(\mathrm{i}\langle\bm{z},r\bm{1}_{\mathbb{R}^{d}% }+\bm{\mu}_{\mathrm{VG}}\rangle T)\bigg{(}1-\mathrm{i}\nu\langle\bm{\theta},% \bm{z}\rangle+\frac{1}{2}\nu\langle\bm{z},\bm{\varSigma}\bm{z}\rangle\bigg{)}^% {\!-T/\nu}$ ,	$\bm{z}\in\mathbb{C}^{d},\quad\Im[\bm{z}]\in\delta_{X}$
NIG	$\displaystyle\exp(\mathrm{i}\langle\bm{z},\bm{X}_{0}\rangle)$
	$\displaystyle\quad{}\times\exp\bigg{(}\mathrm{i}\langle\bm{z},r\bm{1}_{\mathbb% {R}^{d}}+\bm{\mu}_{\mathrm{NIG}}\rangle T+\delta T(\sqrt{\alpha^{2}-\langle\bm% {\beta},\bm{\Delta}\bm{\beta}\rangle}-\sqrt{\alpha^{2}-\langle\bm{\beta}+% \mathrm{i}\bm{z},\bm{\Delta}(\bm{\beta}+\mathrm{i}\bm{z})\rangle})\bigg{)}$ ,	$\bm{z}\in\mathbb{C}^{d},\quad\Im[\bm{z}]\in\delta_{X}$

Table 3: Strip of analyticity,

\delta_{X}

, of the characteristic functions for the different pricing models.

Model	$\delta_{X}$
GBM	$\mathbb{R}^{d}$
VG	$\{\bm{R}\in\mathbb{R}^{d},(1+\nu\langle\bm{\theta},\bm{R}\rangle-\frac{1}{2}% \nu\langle\bm{R},\bm{\varSigma}\bm{R}\rangle)>0\}$
NIG	$\{\bm{R}\in\mathbb{R}^{d},(\alpha^{2}-\langle(\bm{\beta}-\bm{R}),\bm{\Delta}(% \bm{\beta}-\bm{R})\rangle)>0\}$

Remark 2.9 (About the strip of regularity).

Unlike in the 1D case, in the multivariate setting, the choice of the vector of damping parameters $\bm{R}$ , which satisfies the analyticity condition in Table 3, is nontrivial, requiring numerical approximations. Moreover, to obtain more intuition for the multivariate NIG model with $\bm{\Delta}=\bm{I}_{d}$ , the strip of regularity $\smash{\delta_{X}^{\mathrm{NIG}}}$ is an open ball centered at $\bm{\beta}$ with radius $\alpha$ . This fact further complicates the arbitrary choice for damping parameters when the integrand is anisotropic because we must first identify the spherical boundary to determine the admissible combinations of values for the damping parameters enforcing the integrability.

Remark 2.10.

(Efficient vectorized implementation for model calibration) It is often more convenient to work with the scaled versions of the payoff, for example,

	$\displaystyle P(\bm{X}_{T})$	$\displaystyle=\max\bigg{(}K-\sum_{i=1}^{d}w_{i}\mathrm{e}^{X_{T}^{i}},0\bigg{)}$
		$\displaystyle=K\max\bigg{(}1-\sum_{i=1}^{d}\mathrm{e}^{{X_{T}^{\prime}}^{i}},0% \bigg{)},\quad X^{\prime i}_{T}=\log\bigg{(}\frac{S_{0}^{i}}{w_{i}K}\bigg{)},$

so the strike variable, $K$ , is taken out of the integral in (2.1). Moreover, the considered models are stochastic processes with independent increments, which allows us to factorize the multivariate characteristic function as $\smash{\varPhi_{\bm{X}_{T}}(\bm{z})=\mathrm{e}^{\mathrm{i}\langle\bm{z},\bm{X}% _{0}\rangle}\phi_{\bm{X}_{T}}(\bm{z})}$ for $\bm{z}\in\mathbb{C}^{d}$ , such that $\phi_{\bm{X}_{T}}(\bm{z})$ is independent of $\bm{X}_{0}$ .

The advantage in eliminating the dependence of the characteristic function and the payoff from the strike-dependent terms is that we can now evaluate them once in the Fourier domain for multiple values of strike $K$ . This allows for an efficient vectorized implementation that is particularly useful for practitioners interested in model calibration.

3 Methodology

3.1 Motivation for and characterization of the optimal damping rule

This section proposes, and outlines the motivation for, a rule for the optimal choice of the damping parameters $\bm{R}$ that can accelerate the convergence of the numerical quadrature in the Fourier space when approximating (2.1) for pricing multi-asset options under the considered pricing models for various parameters. The main idea is to establish a connection between the damping parameter values, integrand properties and quadrature error.

3.1.1 Motivation for the damping rule

Before considering the integral of interest from (2.1), we provide the general motivation for the rule through a simple 1D integration example for a real-valued function $f$ with respect to a weight function $\lambda(\cdot)$ over the support interval $[a,b]$ (finite, half-infinite or doubly infinite):

$I[f]:=\int_{a}^{b}f(x)\lambda(x)\mathrm{d}x\approx\sum_{k=1}^{N}w_{k}f(x_{k}):% =Q_{N}[f],$

(3.1)

where the quadrature estimator, $\smash{Q_{N}[f]}$ , is characterized by the nodes $\smash{\{x_{k}\}_{k=1}^{N}}$ , which are the roots of the appropriate orthogonal polynomial, $\smash{\pi_{k}(x)}$ , and $\smash{\{w_{k}\}_{k=1}^{N}}$ are the appropriate quadrature weights. Moreover, $\smash{\mathcal{E}_{Q_{N}}[f]}$ denotes the quadrature error (remainder), defined as $\smash{\mathcal{E}_{Q_{N}}[f]:=I[f]-Q_{N}[f]}$ .

The analysis of the quadrature error can be performed through two representations. The first relies on estimates based on high-order derivatives for a smooth function $f$ (Gautschi 2004; Davis and Rabinowitz 2007; Xiang 2012). These error representations are of limited practical value because high-order derivatives are usually challenging to estimate and control. For this reason, to derive our rule, we opt for the second form of quadrature error representation, valid for functions that can be extended holomorphically into the complex plane, which corresponds to the case in (2.5).

There are a number of methods for estimating the error $\smash{\mathcal{E}_{Q_{N}}[f]}$ when $f$ is holomorphic, including

(i)
methods of contour integration (Takahasi and Mori 1971; Donaldson and Elliott 1972),
(ii)
methods based on Hilbert space norm estimates (Davis and Rabinowitz 1954; Donaldson 1973) that consider $\smash{\mathcal{E}_{Q_{N}}}$ as a linear functional on $f$ , and
(iii)
methods based on approximation theory (Babuška et al 2007; Trefethen 2008).

Regardless of the approach taken, the results are often comparable because the error bounds involve the supremum norm of $f$ .

We focus on error estimates based on contour integration tools to showcase these error bounds. This approach uses Cauchy’s theorem in the theory of complex variables to express the value of an analytic function at some point $z$ by means of a contour integral (Cauchy integral) extended over a simple closed curve (or open arc) in the complex plane encircling the point $z$ .

Theorem 3.1.

Assuming that the function $f$ can be extended analytically into a sizable region of the complex plane, which contains the interval $[a,b]$ with no singularities, then the error integral in the approximation (3.1) can be expressed as

$\mathcal{E}_{Q_{N}}[f]=\frac{1}{2\pi\mathrm{i}}\oint_{\mathcal{C}}K_{N}(z)f(z)% \mathrm{d}z,$

(3.2)

where

$K_{N}(z)=\frac{H_{N}(z)}{\pi_{N}(z)},\quad H_{N}(z)=\int_{a}^{b}\lambda(x)% \frac{\pi_{N}(z)}{z-x}\mathrm{d}x,$

(3.3)

and $\mathcal{C}$ is a contour containing the interval $[a,b]$ , within which $f(z)$ has no singularities.⁴⁴ 4 The two most common choices of $\mathcal{C}$ are $\mathcal{C}=\mathcal{C}_{r}$ , the circle $|z|=r$ , $r>1$ , and $C=\mathcal{C}_{\rho}$ , the ellipse with focuses at $a$ and $b$ , where the sum of its semiaxes is equal to $\rho,\rho>1$ . Circles can be used only if the analyticity domain is sufficiently large, and ellipses have the advantage of shrinking to the interval $[a,b]$ when $\rho\rightarrow 1$ , making them suitable for dealing with functions that are analytic on the segment $[a,b]$ .

Proof.

We refer the reader to Donaldson and Elliott (1972) and Gautschi (2004) for a proof. ∎

In the finite case, the contour $\mathcal{C}$ is closed and (3.3) represents an analytic function in the connected domain $\mathbb{C}\setminus[a,b]$ , while we may take $\mathcal{C}$ to lie on the upper and lower edges of the real axis in the infinite case for large $|x|$ . Discussions on choosing adequate contours are found in the work of Elliott (1970), Donaldson and Elliott (1972) and Donaldson (1973). Moreover, precise estimates of $H_{N}(z)$ were derived by Donaldson and Elliott (1972) and Elliott and Tuan (1974).

As $f(\cdot)$ has no singularities within $\mathcal{C}$ , using Theorem 3.1, we obtain

$|\mathcal{E}_{Q_{N}}[f]|\leq\frac{1}{2\pi}\sup_{z\in\mathcal{C}}|f(z)|\oint_{% \mathcal{C}}|K_{N}(z)|\,|\mathrm{d}z|,$

(3.4)

where the quantity $\smash{\oint_{\mathcal{C}}|K_{N}(z)|\,|\mathrm{d}z|}$ depends on the quadrature rule only. We expect that when the size of the contour increases, $\smash{\oint_{\mathcal{C}}|K_{N}(z)|\,|\mathrm{d}z|}$ decreases, whereas $\sup_{z\in\mathcal{C}}|f(z)|$ increases by the maximum modulus theorem. The optimal choice of the contour $\mathcal{C}$ is the one that minimizes the right-hand side of (3.4).

Extending the error bound (3.4) to a multidimensional setting can be performed straightforwardly in a recursive way and by using tensorization tools (see Appendix A (online) for an illustration). Moreover, the term $\sup_{z\in\mathcal{C}}|f(z)|$ in the upper bound of (3.4) is independent of the quadrature method.

3.2 Characterization of the optimal damping rule

Motivated by the error bound (3.4), we propose a rule for choosing the damping parameters that improves the numerical convergence of the designed numerical quadrature method (see Section 3.3) when approximating (2.1). The rule consists in solving the following constrained optimization problem:

$\bm{R}^{\ast}:=\bm{R}^{\ast}(\bm{\varTheta}_{m},\bm{\varTheta}_{p})=% \operatorname*{arg~{}min}_{\bm{R}\in\delta_{V}}\sup_{\bm{u}\in\mathbb{R}^{d}}|% g(\bm{u};\bm{R},\bm{\varTheta}_{m},\bm{\varTheta}_{p})|,$

(3.5)

where we use the notation from Section 2, $g(\cdot)$ is defined in (2.5) and $\bm{R}^{\ast}:=(R_{1}^{\ast},\dots,R_{d}^{\ast})$ denotes the vector of optimal damping parameters.

We can use Proposition 3.2 to reduce (3.5) to a simpler optimization problem. This is done by finding the vector of the damping parameters, $\smash{\bm{R}\in\delta_{V}}$ , that minimize the peak of the integrand in (2.5) at the origin point $\smash{\bm{u}=\bm{0}_{\mathbb{R}^{d}}}$ .

Proposition 3.2.

For $g$ defined by (2.5) and $\bm{R}\in\delta_{V}$ , we have

$\bm{R}^{\ast}=\operatorname*{arg~{}min}_{\bm{R}\in\delta_{V}}\sup_{\bm{u}\in% \mathbb{R}^{d}}|g(\bm{u};\bm{R},\bm{\varTheta}_{m},\bm{\varTheta}_{p})|=% \operatorname*{arg~{}min}_{\bm{R}\in\delta_{V}}g(\bm{0}_{\mathbb{R}^{d}};\bm{R% },\bm{\varTheta}_{m},\bm{\varTheta}_{p}).$

(3.6)

Proof.

Let $\smash{f\colon\mathbb{R}^{d}\mapsto\mathbb{R}_{+}}$ be an arbitrary real-valued nonnegative function, and let $\bm{R}\in\mathbb{R}^{d}$ such that the dampened function $\bm{x}\mapsto f_{\bm{R}}(\bm{x})\in L^{1}(\mathbb{R}^{d})$ , then we have that its Fourier transform satisfies

	$\displaystyle\|\hat{f}_{\bm{R}}(\bm{u})\|$	$\displaystyle\leq\int_{\mathbb{R}^{d}}\|\mathrm{e}^{\mathrm{i}\langle\bm{u},\bm% {x}\rangle}\|\mathrm{e}^{\langle\bm{R},\bm{x}\rangle}f(\bm{x})\mathrm{d}\bm{x}$
		$\displaystyle=\int_{\mathbb{R}^{d}}\mathrm{e}^{\langle\bm{R},\bm{x}\rangle}f(% \bm{x})\mathrm{d}\bm{x}=\hat{f}_{\bm{R}}(\bm{0}),\quad\bm{u}\in\mathbb{R}^{d}.$		(3.7)

Moreover, by (2.3) we have that $\smash{\hat{f}_{\bm{R}}(\bm{u})=\hat{f}(\bm{u}+\mathrm{i}\bm{R})}$ ; hence, (3.7) implies that

$|\hat{f}(\bm{u}+\mathrm{i}\bm{R})|\leq\hat{f}(\mathrm{i}\bm{R})\quad\forall\bm% {u}\in\mathbb{R}^{d}.$

(3.8)

Equation (3.8) is known as the ridge property of Fourier transforms (see Lukacs 1970). Since both the payoff and the probability density functions are real-valued and nonnegative, their Fourier transforms satisfy the ridge property (ie, $|\hat{P}(\bm{u}+\mathrm{i}\bm{R})|\leq\hat{P}(\mathrm{i}\bm{R})$ for all $\bm{u}\in\mathbb{R}^{d}$ , $\bm{R}\in\delta_{P}$ and $|\varPhi_{\bm{X}_{T}}(\bm{u}+\mathrm{i}\bm{R})|\leq\varPhi_{\bm{X}_{T}}(% \mathrm{i}\bm{R})$ for all $\bm{u}\in\mathbb{R}^{d}$ , $\bm{R}\in\delta_{X}$ ); hence, the integrand (2.1) can be bounded by

	$\displaystyle\|g(\bm{u};\bm{R},\bm{\varTheta}_{m},\bm{\varTheta}_{p})\|$	$\displaystyle\leq(2\pi)^{-d}\mathrm{e}^{-rT}\|\varPhi_{\bm{X}_{T}}(\mathrm{i}% \bm{R})\|\,\|\hat{P}(\mathrm{i}\bm{R})\|$
		$\displaystyle=\|g(\bm{0}_{\mathbb{R}^{d}};\bm{R})\|\quad\forall\bm{u}\in\mathbb{% R}^{d},~{}\bm{R}\in\delta_{V}.$		(3.9)

∎

Equation (3.6) cannot be solved analytically, especially in high dimensions; therefore, we solve it numerically, approximating $\bm{R}^{\ast}$ by $\smash{\bar{\bm{R}}=(\bar{R}_{1},\dots,\bar{R}_{d})}$ . In this context, we use the interior point method (Byrd et al 2000, 1999) with an accuracy of order $10^{-6}$ ; other algorithms such as the limited-memory Broyden–Fletcher–Goldfarb–Shanno algorithm (L-BFGS-B) were tested and work effectively.

One-dimensional illustration showing the shape of the integrand with respect to the damping parameter R and calE-sub-R convergence with respect to N, using Gauss--Laguerre quadrature for the European put option under different pricing models (S-sub-0=100, K=100, r=0%, T=1). The left-hand column shows the shape of the integrand with respect to the damping parameter R. The right-hand column shows calE-sub-R convergence with respect to N. (a) GBM, where sigma=0.4. (b) VG, where sigma=0.4, theta=-0.3, nu=0.257 — Figure 1: One-dimensional illustration showing the shape of the integrand with respect to the damping parameter $R$ and $\mathcal{E}_{R}$ convergence with respect to $N$ , using Gauss–Laguerre quadrature for the European put option under different pricing models ( $S_{\text{0}}=\text{100}$ , $K=\text{100}$ , $r=\text{0}\%$ , $T=\text{1}$ ). The left-hand column shows the shape of the integrand with respect to the damping parameter $R$ . The right-hand column shows $\mathcal{E}_{R}$ convergence with respect to $N$ . (a) GBM, where $\sigma=\text{0.4}$ . (b) VG, where $\sigma=\text{0.4}$ , $\theta=-\text{0.3}$ , $\nu=\text{0.257}$ . (c) NIG, where $\alpha=\text{15}$ , $\beta=-\text{3}$ , $\delta=\text{0.2}$ . The relative quadrature error $\mathcal{E}_{R}$ is defined as $\mathcal{E}_{R}=|Q_{N}[g]-\text{RV}|/\text{RV}$ , where $Q_{N}$ is the quadrature estimator of (2.1) based on the Gauss–Laguerre rule.

The effect of the damping parameters on the shape of the integrand in the case of 2D-basket put options under the VG model with parameters bold-sigma=(0.4,0.4), bold-theta=(-0.3,-0.3), nu=0.257. (a) bold-R=(0.2,0.2). (b) bold-R=(1,1). (c) bold-R=(2,2). (d) bold-R=(3,3). — Figure 2: The effect of the damping parameters on the shape of the integrand in the case of 2D-basket put options under the VG model with parameters $\bm{\sigma}=(\text{0.4},\text{0.4})$ , $\bm{\theta}=(-\text{0.3},-\text{0.3})$ , $\nu=\text{0.257}$ . (a) $\bm{R}=(\text{0.2},\text{0.2})$ . (b) $\bm{R}=(\text{1},\text{1})$ . (c) $\bm{R}=(\text{2},\text{2})$ . (d) $\bm{R}=(\text{3},\text{3})$ .

The numerical investigation through different models and parameters confirms that the damping parameters have a considerable effect on the properties of the integrand, particularly its peak, tail-heaviness and oscillatory behavior (Figure 1 shows the single put option under the different models and Figure 2 shows the 2D-basket put under VG). We observe that the damping parameters that produce the lowest peak of the integrand around the origin are associated with a faster convergence of the relative quadrature error than other damping parameters. Moreover, we observe that highly peaked integrands are more likely to oscillate, implying a deteriorated convergence of the numerical quadrature. Regardless of the quadrature method used (see Section 3.3 for an explanation of these methods), this observation is consistent for several parameter constellations under the three tested pricing dynamics (GBM, VG and NIG), and for different dimensions of the basket put and rainbow options. The computational advantage of the optimal damping rule on the error convergence for the multi-asset basket put and call on min options is illustrated for different models in Section 4.1.2.

Remark 3.3 (The case of the isotropic integrand).

The $d$ -dimensional optimization problem (3.6) is simplified further to a 1D problem when the integrand is isotropic.

Remark 3.4 (Improving the damping parameters rule).

Other rules for choosing the damping parameters can be investigated to improve the numerical convergence of quadrature methods. We can account for additional features, such as the distance of the damping parameters to the poles, which affects the choice of the integration contour in (3.2), or controlling the regularity of the integrand via high-order derivative estimates. However, we expect such rules to be more complicated and computationally expensive (eg, the evaluation of the gradient of the integrand). Investigating other rules remains open for future work.

3.3 Numerical evaluation of the inverse Fourier integrals using hierarchical deterministic quadrature methods

We aim to efficiently approximate (2.1) using a tensorization of quadrature formulas over $\mathbb{R}^{d}$ . When using Fourier transforms for option pricing, the standard numerical approach truncates and discretizes the integration domain and uses an FFT based on bounded equispaced quadrature formulas, such as the trapezoidal rule or Simpson’s rule. The FFT is restricted to the use of a uniform quadrature mesh, in contrast to Gaussian quadrature rules, which have higher polynomial exactness (the $N$ -point Gaussian quadrature rule is exact up to polynomials of degree $2N-1$ ). Moreover, using FFT requires the truncation range to be prespecified. This option is efficient in the 1D setting, as the estimation of the truncation intervals based, for instance, on the cumulants was widely covered in the literature. It remains affordable even though the additional cost might be high due to the inappropriate choice of truncation parameters. However, this is not the case in the multidimensional setting because determining the truncation parameters becomes more challenging. Moreover, the truncation errors nontrivially depend on the damping parameter values. Choosing larger than necessary truncation domains leads to a more significant increase in the computational effort for higher dimensions. Finally, FFT-based approaches need to be followed by interpolation techniques in order to obtain the option values at the desired strikes grid, which may lead to a loss of accuracy. This is in contrast to DI methods, which can be efficiently vectorized with respect to the desired strikes grid.

For these reasons, we choose the DI approach with Gaussian quadrature rules. Moreover, our numerical investigation (see Appendix B (online)) suggests that the Gauss–Laguerre quadrature exhibits faster convergence than the Gauss–Hermite rule. Therefore, we use Laguerre quadrature on semi-infinite domains after applying the necessary transformations.

Before defining the multivariate quadrature estimators, we first introduce the notation in the univariate setting (for more details, see Chen (2018)). Let $\beta$ denote a nonnegative integer, referred to as the “discretization level”, and let $m\colon\mathbb{N}\rightarrow\mathbb{N}$ represent a strictly increasing function with $m(0)=0$ and $m(1)=1$ , called a “level-to-nodes function”. At each level $\beta$ , we consider a set of $m(\beta)$ distinct quadrature points

$\mathcal{H}^{m(\beta)}=\{x_{\beta}^{1},x_{\beta}^{2},\dots,x_{\beta}^{m(\beta)% }\}\subset\mathbb{R}$

and a set of quadrature weights

$\bm{\omega}^{m(\beta)}=\{\omega_{\beta}^{1},\omega_{\beta}^{2},\dots,\omega_{% \beta}^{m(\beta)}\}.$

We also let $C^{0}(\mathbb{R})$ be the space of real-valued continuous functions over $\mathbb{R}$ . We define the univariate quadrature operator applied to a function $f\in C^{0}(\mathbb{R})$ as

$Q^{m(\beta)}\colon C^{0}(\mathbb{R})\rightarrow\mathbb{R},\qquad Q^{m(\beta)}[% f]:=\sum_{j=1}^{m(\beta)}f(x_{\beta}^{j})\omega_{\beta}^{j}.$

In our case, in (2.1), we have a multivariate integration problem for $g$ (see (2.5)) over $\mathbb{R}^{d}$ . Accordingly, for a multi-index $\smash{\bm{\beta}=(\beta_{i})_{i=1}^{d}\in\mathbb{N}^{d}}$ , the $d$ -dimensional quadrature operator applied to $g$ is defined as

	$\displaystyle Q_{d}^{m(\bm{\beta})}\colon C^{0}(\mathbb{R}^{d})\rightarrow% \mathbb{R},\quad Q_{d}^{m(\bm{\beta})}=\bigotimes_{i=1}^{d}Q^{m(\beta_{i})},$
	$\displaystyle Q_{d}^{m(\bm{\beta})}[g]:=\sum_{j_{1}=1}^{m(\beta_{1})}\cdots% \sum_{j_{d}=1}^{m(\beta_{d})}\omega_{\beta_{1}}^{j_{1}}\cdots\omega_{\beta_{1}% }^{j_{d}}g(x_{\beta_{1}}^{j_{1}},\dots,x_{\beta_{d}}^{j_{d}}):=\sum_{j=1}^{\#% \mathcal{T}^{m(\bm{\beta})}}g(\hat{x}_{j})\bar{\omega}_{j},$

where $\smash{\hat{x}_{j}\in\mathcal{T}^{m(\bm{\beta})}:=\prod_{i=1}^{d}\mathcal{H}^{% m(\beta_{i})}}$ (with cardinality $\smash{\#\mathcal{T}^{m(\bm{\beta})}=\prod_{i=1}^{d}m(\beta_{i})}$ and $\smash{m(\beta_{i})=N_{i}}$ quadrature points in the dimension of $x_{i}$ ),⁵⁵ 5 The $n$ th quadrature operator acts only on the $n$ th variable of $g$ . and $\bar{\omega}_{j}$ is a product of the weights of the univariate quadrature rule. To simplify the notation, we replace $\smash{Q_{d}^{m(\bm{\beta})}}$ with $\smash{Q_{d}^{\bm{\beta}}}$ .

We define the set of differences $\Delta Q_{d}^{\bm{\beta}}$ for indexes $i\in\{1,\dots,d\}$ as follows:

$\Delta_{i}Q_{d}^{\bm{\beta}}:=\begin{cases}Q_{d}^{\bm{\beta}}-Q_{d}^{\bm{\beta% }^{\prime}}&\text{with~{}}\bm{\beta}^{\prime}=\bm{\beta}-\bm{e}_{i},\ \beta_{i% }>0,\\ \\ Q_{d}^{\bm{\beta}}&\text{otherwise},\end{cases}$

(3.10)

where $\bm{e}_{i}$ denotes the $i$ th $d$ -dimensional unit vector. Then, using the telescopic property, the quadrature estimator, defined with respect to a choice of the set of multi-indexes $\mathcal{4}\subset\mathbb{N}^{d}$ , is expressed by⁶⁶ 6 For instance, when $d=2$ , $\Delta Q_{2}^{\bm{\beta}}=\Delta_{2}\Delta_{1}Q_{2}^{(\beta_{1},\beta_{2})}=Q_% {2}^{(\beta_{1},\beta_{2})}-Q_{2}^{(\beta_{1},\beta_{2}-1)}-Q_{2}^{(\beta_{1}-% 1,\beta_{2})}+Q_{2}^{(\beta_{1}-1,\beta_{2}-1)}$ . ${}^{,}$ ⁷⁷ 7 To ensure the validity of the telescoping sum expansion, the index set $\mathcal{4}$ must satisfy the admissibility condition (ie, $\bm{\beta}\in\mathcal{4}$ , $\bm{\alpha}\leq\bm{\beta}\Rightarrow\bm{\alpha}\in\mathcal{4}$ , where $\bm{\alpha}\leq\bm{\beta}$ is defined as $\alpha_{i}\leq\beta_{i}$ , $i=1,\dots,d$ ).

$Q_{d}^{\mathcal{4}}=\sum_{\bm{\beta}\in\mathcal{4}}\Delta Q_{d}^{\bm{\beta}}% \quad\text{with }\Delta Q_{d}^{\bm{\beta}}=\bigg{(}\bigotimes_{i=1}^{d}\Delta_% {i}\bigg{)}Q_{d}^{\bm{\beta}},$

(3.11)

and the quadrature error can be written as

$\mathcal{E}_{Q}=|Q^{\infty}_{d}[g]-Q_{d}^{\mathcal{4}}[g]|\leq\sum_{\bm{\beta}% \in\mathbb{N}^{d}\backslash\{\mathcal{4}\}}|\Delta Q_{d}^{\bm{\beta}}[g]|,$

(3.12)

where

$Q_{d}^{\infty}:=\sum_{\beta_{1}=0}^{\infty}\cdots\sum_{\beta_{d}=0}^{\infty}% \Delta Q_{d}^{(\beta_{1},\dots,\beta_{d})}=\sum_{\bm{\beta}\in\mathbb{N}^{d}}% \Delta Q_{d}^{\bm{\beta}}.$

(3.13)

In (3.11), the strategy chosen for construction of the index set $\mathcal{4}$ and the hierarchy of quadrature points determined by $m(\cdot)$ defines the different hierarchical quadrature methods. Table 4 shows the details of the methods considered in this paper.

Table 4: Construction details for the quadrature methods. [Here,

l\in\mathbb{N}

represents a given level and

\bar{T}\in\mathbb{R}

is a threshold value.]

Quadrature
method	$m(\cdot)$	$\mathcal{4}$
TP	$m(\beta)=\beta$	$\displaystyle\mathcal{4}^{\text{TP}}(l)=\bigg{\{}\bm{\beta}\in\mathbb{N}^{d}% \colon\max_{\text{1}\leq i\leq d}(\beta_{i}-\text{1})\leq l\bigg{\}}$
SM sparse grids	$m(\beta)=\text{2}^{\beta-\text{1}}+\text{1}$ ,	$\displaystyle\mathcal{4}^{\text{SM}}(l)=\bigg{\{}\bm{\beta}\in\mathbb{N}^{d}% \colon\sum_{\text{1}\leq i\leq d}(\beta_{i}-\text{1})\leq l\bigg{\}}$
	$\beta>\text{1}$ , $m(\text{1})=\text{1}$
ASGQ	$m(\beta)=\text{2}^{\beta-\text{1}}+\text{1}$ ,	$\displaystyle\mathcal{4}^{\text{ASGQ}}=\{\bm{\beta}\in\mathbb{N}_{+}^{d}\colon P% _{\bm{\beta}}\geq\bar{T}\}$
	$\beta>\text{1}$ , $m(\text{1})=\text{1}$	(see (3.14) and (3.15))

In many situations, the tensor product (TP) estimator can become rapidly unaffordable because the number of function evaluations increases exponentially with the problem dimensionality, known as the “curse of dimensionality”. We use Smolyak (SM) and ASGQ methods based on sparsification and dimension-adaptivity techniques to overcome this issue. For both TP and SM methods, the construction of the index set is performed a priori. However, ASGQ allows for the a posteriori and adaptive construction of the index set $\mathcal{4}$ by greedily exploiting the mixed regularity of the integrand during the actual computation of the quantity of interest. The construction of $\mathcal{4}^{\mathrm{ASGQ}}$ is performed through profit thresholding, where new indexes are selected iteratively based on the error versus cost–profit rule, with a hierarchical surplus defined by

$P_{\bm{\beta}}=\frac{|\Delta E_{\bm{\beta}}|}{\Delta\mathcal{W}_{\bm{\beta}}},$

(3.14)

where $\smash{\Delta\mathcal{W}_{\bm{\beta}}}$ is the work contribution (ie, the computational cost required to add $\smash{\Delta Q_{d}^{\bm{\beta}}}$ to $Q_{d}^{\mathcal{4}^{\mathrm{ASGQ}}}$ ) and $\Delta E_{\bm{\beta}}$ is the error contribution (ie, a measure of how much the quadrature error would decrease once $\smash{\Delta Q_{d}^{\bm{\beta}}}$ has been added to $\smash{Q_{d}^{\mathcal{4}^{\mathrm{ASGQ}}}}$ ):

	$\displaystyle\Delta E_{\bm{\beta}}$	$\displaystyle=\|Q_{d}^{\mathcal{4}^{\mathrm{ASGQ}}\cup\{\bm{\beta}\}}[g]-Q_{d}^% {\mathcal{4}^{\mathrm{ASGQ}}}[g]\|,$		(3.15)
	$\displaystyle\Delta\mathcal{W}_{\bm{\beta}}$	$\displaystyle=\operatorname{Work}[Q_{d}^{\mathcal{4}^{\mathrm{ASGQ}}\cup\{\bm{% \beta}\}}[g]]-\operatorname{Work}[Q_{d}^{\mathcal{4}^{\mathrm{ASGQ}}}[g]].$

The convergence speed for all quadrature methods in this work is determined by the behavior of the quadrature error defined in (3.12). In this context, given the model and option parameters, the convergence rate depends on the damping parameter values, which control the regularity of the integrand $g$ in the Fourier space.

We let $\smash{N:=\prod_{i=1}^{d}m(\beta_{i})}$ denote the total number of quadrature points used by each method. For the TP method, we have the following (Davis and Rabinowitz 2007):

$\mathcal{E}^{\mathrm{TP}}_{Q}(N;\bm{R})=\mathcal{O}(N^{-r_{t}/d})$

(3.16)

for functions with bounded total derivatives up to order $r_{t}:=r_{t}(\bm{R})$ . When using SM sparse grids (not adaptive), we obtain (Smolyak 1963; Wasilkowski and Wozniakowski 1995; Gerstner and Griebel 1998; Barthelmann et al 2000)

$\mathcal{E}^{\mathrm{SM}}_{Q}(N;\bm{R})=\mathcal{O}(N^{-r_{m}}(\log N)^{(d-1)(% r_{m}+1)})$

(3.17)

for functions with bounded mixed partial derivatives up to order $r_{m}:=r_{m}(\bm{R})$ . Moreover, it was observed by Gerstner and Griebel (2003) that the convergence is even spectral for analytic functions ( $r_{m}\rightarrow+\infty$ ). For the ASGQ method, we get (Chen 2018)

$\mathcal{E}^{\mathrm{ASGQ}}_{Q}(N;\bm{R})=\mathcal{O}(N^{-r_{w}/2}),$

(3.18)

where $r_{w}$ is related to the degree of weighted mixed regularity of the integrand.

In (3.16)–(3.18), we emphasize the dependence of the convergence rates on the damping parameters $\bm{R}$ , which is valid only in this context because these parameters control the regularity of the integrand in the Fourier space. Moreover, our optimized choice of $\bm{R}$ using (3.5) is used not only to increase the order of boundedness of the derivatives of $g$ but also to reduce the bounds on these derivatives.

4 Numerical experiments and results

In this section, we present the results of different numerical experiments conducted for pricing multi-asset European equally weighted basket put ( $w_{i}=1/d$ , $i=1,\dots,d$ ) and call on min options. These examples are tested under the multivariate GBM, VG and NIG models with various parameter constellations for different dimensions $d\in\{2,4,6\}$ . The tested model parameters are justified based on the literature on model calibration (Kirkby 2015; van de Wiel 2015; Choi 2018; Bayer et al 2018; Aguilar 2020; Healy 2021). The detailed illustrated examples are presented in Tables 5–7. To compare the methods in this work, we consider relative errors normalized by the reference prices. When using quadrature methods, the error is the relative quadrature error defined as

$\mathcal{E}_{R}=\frac{|Q_{d}^{\mathcal{4}}[g]-\text{RV}|}{\text{RV}},$

where RV denotes the reference value, and the 95% relative statistical error of the Monte Carlo method is estimated by virtue of the central limit theorem (CLT) as

$\mathcal{E}_{R}\approx\frac{C_{\alpha}\times\sigma_{M}}{\text{RV}\times\sqrt{M% }},$

(4.1)

where $C_{\alpha}=1.96$ for the 95% confidence level, $M$ is the number of Monte Carlo samples and $\sigma_{M}$ is the standard deviation of the quantity of interest.

Table 5: Examples of multi-asset options under the multivariate GBM model. [In all examples,

S_{\text{0}}^{i}=\text{100}

,

i=\text{1},\dots,d

,

T=\text{1}

,

r=\text{0}

. Reference values are computed with Monte Carlo using

\text{10}^{\text{9}}

samples, with the 95% statistical error estimates reported between parentheses.

\bar{R}

is rounded to one decimal place.]

Example	Option	Parameters	value	parameters $\bar{\bm{R}}$
			Reference	Optimal damping
1	2D-basket put	$\bm{\sigma}=(\text{0.4},\text{0.4})$ ,	11.4474	$(\text{2.5},\text{2.5})$
		$\bm{C}=I_{\text{2}}$ , $K=\text{100}$	$(\text{8.0}\text{e}^{-\text{4}})$
2	2D-basket put	$\bm{\sigma}=(\text{0.4},\text{0.8})$ ,	17.831	$(\text{2.1},\text{1.2})$
		$\bm{C}=I_{\text{2}}$ , $K=\text{100}$	$(\text{1.2}\text{e}^{-\text{3}})$
3	2D-call on min	$\bm{\sigma}=(\text{0.4},\text{0.4})$ ,	0 3.4603	$(-\text{3.4},-\text{3.4})$
		$\bm{C}=I_{\text{2}}$ , $K=\text{100}$	$(\text{6.0}\text{e}^{-\text{4}})$
4	2D-call on min	$\bm{\sigma}=(\text{0.4},\text{0.8})$ ,	0 3.7411	$(-\text{3.6},-\text{1.8})$
		$\bm{C}=I_{\text{2}}$ , $K=\text{100}$	$(\text{8.2}\text{e}^{-\text{4}})$
5	4D-basket put	$\bm{\sigma}=(\text{0.4},\text{0.4},\text{0.4},\text{0.4})$ ,	0 8.193	$(\text{2.1},\text{2.1},\text{2.1},\text{2.1})$
		$\bm{C}=I_{\text{4}}$ , $K=\text{100}$	$(\text{6.0}\text{e}^{-\text{4}})$
6	4D-basket put	$\bm{\sigma}=(\text{0.2},\text{0.4},\text{0.6},\text{0.8})$ ,	11.3014	$(\text{2.4},\text{1.9},\text{1.5},\text{1.2})$
		$\bm{C}=I_{\text{4}}$ , $K=\text{100}$	$(\text{8.0}\text{e}^{-\text{4}})$
7	4D-call on min	$\bm{\sigma}=(\text{0.4},\text{0.4},\text{0.4},\text{0.4})$ ,	0 0.317	$(-\text{3.1},-\text{3.1},-\text{3.1},-\text{3.1})$
		$\bm{C}=I_{\text{4}}$ , $K=\text{100}$	$(\text{2.0}\text{e}^{-\text{4}})$
8	4D-call on min	$\bm{\sigma}=(\text{0.2},\text{0.4},\text{0.6},\text{0.8})$ ,	0 0.2382	$(-\text{6.4},-\text{3.1},-\text{2.1},-\text{1.6})$
		$\bm{C}=I_{\text{4}}$ , $K=\text{100}$	$(\text{1.0}\text{e}^{-\text{4}})$
9	6D-basket put	$\bm{\sigma}=(\text{0.4},\text{0.4},\text{0.4},\text{0.4},\text{0.4},\text{0.4})$ ,	0 0.0041	$(\text{2.0},\text{2.0},\text{2.0},\text{2.0},\text{2.0},\text{2.0})$
		$\bm{C}=I_{\text{6}}$ , $K=\text{60}$	$(\text{8.8}\text{e}^{-\text{6}})$
10	6D-basket put	$\bm{\sigma}=(\text{0.2},\text{0.3},\text{0.4},\text{0.5},\text{0.6},\text{0.7})$ ,	0 0.012702	$(\text{2.3},\text{2.1},\text{1.9},\text{1.7},\text{1.5},\text{1.3})$
		$\bm{C}=I_{\text{6}}$ , $K=\text{60}$	$(\text{1.8}\text{e}^{-\text{5}})$
11	6D-call on min	$\bm{\sigma}=(\text{0.4},\text{0.4},\text{0.4},\text{0.4},\text{0.4},\text{0.4})$ ,	0 0.038	$(-\text{3.0},-\text{3.0},-\text{3.0},-\text{3.0},-\text{3.0},-\text{3.0})$
		$\bm{C}=I_{\text{6}}$ , $K=\text{100}$	$(\text{4.4}\text{e}^{-\text{5}})$
12	6D-call on min	$\bm{\sigma}=(\text{0.2},\text{0.3},\text{0.4},\text{0.5},\text{0.6},\text{0.7})$ ,	0 0.0301	$(-\text{6.0},-\text{3.9},-\text{3.0},-\text{2.4},-\text{2.0},-\text{1.8})$
		$\bm{C}=I_{\text{6}}$ , $K=\text{100}$	$(\text{3.7}\text{e}^{-\text{5}})$

Table 6: Examples of multi-asset options under the multivariate VG model. [In all examples,

S_{\text{0}}^{i}=\text{100},i=\text{1},\dots,d,T=\text{1},r=\text{0}

. Reference values are computed with Monte Carlo using

\text{10}^{\text{9}}

samples, with the 95% statistical error estimates reported between parentheses.

\bar{R}

is rounded to one decimal place.]

Example	Option	Parameters	value	parameters $\bar{\bm{R}}$
			Reference	Optimal damping
13	2D-basket put	$\bm{\sigma}=(\text{0.4},\text{0.4})$ ,	11.7589	$(\text{1.7},\text{1.7})$
		$\bm{\theta}=(-\text{0.3},-\text{0.3})$ , $\nu=\text{0.257}$ , $K=\text{100}$	$(\text{1.0}\text{e}^{-\text{3}})$
14	2D-basket put	$\bm{\sigma}=(\text{0.4},\text{0.8})$ ,	17.6688	$(\text{1.7},\text{1.0})$
		$\bm{\theta}=(-\text{0.3},\text{0})$ , $\nu=\text{0.257}$ , $K=\text{100}$	$(\text{1.2}\text{e}^{-\text{3}})$
15	2D-call on min	$\bm{\sigma}=(\text{0.4},\text{0.4})$ ,	0 3.9601	$(-\text{3.5},-\text{3.5})$
		$\bm{\theta}=(-\text{0.3},-\text{0.3})$ , $\nu=\text{0.257}$ , $K=\text{100}$	$(\text{7.0}\text{e}^{-\text{4}})$
16	2D-call on min	$\bm{\sigma}=(\text{0.4},\text{0.8})$ ,	0 3.3422	$(-\text{4.0},-\text{3.5})$
		$\bm{\theta}=(-\text{0.3},\text{0})$ , $\nu=\text{0.257}$ , $K=\text{100}$	$(\text{8.0}\text{e}^{-\text{4}})$
17	4D-basket put	$\bm{\sigma}=(\text{0.4},\text{0.4},\text{0.4},\text{0.4})$ ,	0 8.9441	$(\text{1.2},\text{1.2},\text{1.2},\text{1.2})$
		$\bm{\theta}=(-\text{0.3},-\text{0.3},-\text{0.3},-\text{0.3})$ , $\nu=\text{0.257}$ , $K=\text{100}$	$(\text{8.0}\text{e}^{-\text{4}})$
18	4D-basket put	$\bm{\sigma}=(\text{0.2},\text{0.4},\text{0.6},\text{0.8})$ ,	11.2277	$(\text{1.6},\text{1.4},\text{1.1},\text{0.9})$
		$\bm{\theta}=(-\text{0.3},-\text{0.2},-\text{0.1},\text{0})$ , $\nu=\text{0.257}$ , $K=\text{100}$	$(\text{8.0}\text{e}^{-\text{4}})$
19	4D-call on min	$\bm{\sigma}=(\text{0.4},\text{0.4},\text{0.4},\text{0.4})$ ,	0 0.6137	$(-\text{3.2},-\text{3.2},-\text{3.2},-\text{3.2})$
		$\bm{\theta}=(-\text{0.3},-\text{0.3},-\text{0.3},-\text{0.3})$ , $\nu=\text{0.257}$ , $K=\text{100}$	$(\text{2.0}\text{e}^{-\text{4}})$
20	4D-call on min	$\bm{\sigma}=(\text{0.2},\text{0.4},\text{0.6},\text{0.8})$ ,	0 0.2384	$(-\text{6.6},-\text{3.0},-\text{2.0},-\text{1.5})$
		$\bm{\theta}=(-\text{0.3},-\text{0.2},-\text{0.1},\text{0})$ , $\nu=\text{0.257}$ , $K=\text{100}$	$(\text{1.0}\text{e}^{-\text{4}})$
21	6D-basket put	$\bm{\sigma}=(\text{0.4},\text{0.4},\text{0.4},\text{0.4},\text{0.4},\text{0.4})$ ,	0 0.1691	$(\text{1.1},\text{1.1},\text{1.1},\text{1.1},\text{1.1},\text{1.1})$
		$\bm{\theta}=-(\text{0.3},\text{0.3},\text{0.3},\text{0.3},\text{0.3},\text{0.3})$ , $\nu=\text{0.257}$ , $K=\text{60}$	$(\text{1.0}\text{e}^{-\text{6}})$
22	6D-basket put	$\bm{\sigma}=(\text{0.2},\text{0.3},\text{0.4},\text{0.5},\text{0.6},\text{0.7})$ ,	0 0.04634	$(\text{2.1},\text{1.9},\text{1.7},\text{1.6},\text{1.4},\text{1.2})$
		$\bm{\theta}=(-\text{0.3},-\text{0.2},-\text{0.1},\text{0},\text{0.1},\text{0.2})$ , $\nu=\text{0.257}$ , $K=\text{60}$	$(\text{5.0}\text{e}^{-\text{5}})$
23	6D-call on min	$\bm{\sigma}=(\text{0.4},\text{0.4},\text{0.4},\text{0.4},\text{0.4},\text{0.4})$ ,	0 0.16248	$(-\text{3.1},-\text{3.1},-\text{3.1},-\text{3.1},-\text{3.1},-\text{3.1})$
		$\bm{\theta}=-(\text{0.3},\text{0.3},\text{0.3},\text{0.3},\text{0.3},\text{0.3})$ , $\nu=\text{0.257}$ , $K=\text{100}$	$(\text{1.0}\text{e}^{-\text{4}})$
24	6D-call on min	$\bm{\sigma}=(\text{0.2},\text{0.3},\text{0.4},\text{0.5},\text{0.6},\text{0.7})$ ,	0 0.02269	$(-\text{6.5},-\text{3.7},-\text{2.6},-\text{2.0},-\text{1.7},-\text{1.4})$
		$\bm{\theta}=(-\text{0.3},-\text{0.2},-\text{0.1},\text{0},\text{0.1},\text{0.2})$ , $\nu=\text{0.257}$ , $K=\text{100}$	$(\text{4.0}\text{e}^{-\text{5}})$

Table 7: Examples of multi-asset options under the multivariate NIG model. [In all examples,

S_{\text{0}}^{i}=\text{100},i=\text{1},\dots,d,T=\text{1},r=\text{0}

. Reference values are computed with Monte Carlo using

\text{10}^{\text{9}}

samples, with the 95% statistical error estimates reported between parentheses.

\bar{R}

is rounded to one decimal place.]

Example	Option	Parameters	value	parameters $\bar{\bm{R}}$
			Reference	Optimal damping
25	2D-basket put	$\bm{\beta}=(-\text{3},-\text{3})$ ,	3.3199	$(\text{6.1},\text{6.1})$
		$\alpha=\text{15}$ , $\delta=\text{0.2},\bm{\Delta}=\bm{I}_{\text{2}},K=\text{100}$	$(\text{3.0}\text{e}^{-\text{4}})$
26	2D-basket put	$\bm{\beta}=(-\text{3},\text{0})$ ,	3.8978	$(\text{4.6},\text{4.8})$
		$\alpha=\text{10}$ , $\delta=\text{0.2},\bm{\Delta}=\bm{I}_{\text{2}},K=\text{100}$	$(\text{4.0}\text{e}^{-\text{4}})$
27	2D-call on min	$\bm{\beta}=(-\text{3},-\text{3})$ ,	1.2635	$(-\text{9.9},\text{-9.9})$
		$\alpha=\text{15}$ , $\delta=\text{0.2},\bm{\Delta}=\bm{I}_{\text{2}},K=\text{100}$	$(\text{2.0}\text{e}^{-\text{4}})$
28	2D-call on min	$\bm{\beta}=(-\text{3},\text{0})$ ,	1.4476	$(-\text{7.5},-\text{6.8})$
		$\alpha=\text{10}$ , $\delta=\text{0.2},\bm{\Delta}=\bm{I}_{\text{2}},K=\text{100}$	$(\text{2.0}\text{e}^{-\text{4}})$
29	4D-basket put	$\bm{\beta}=(-\text{3},-\text{3},-\text{3},-\text{3})$ ,	2.554	$(\text{4.0},\text{4.0},\text{4.0},\text{4.0})$
		$\alpha=\text{15}$ , $\delta=\text{0.4},\bm{\Delta}=\bm{I}_{\text{4}},K=\text{100}$	$(\text{3.0}\text{e}^{-\text{4}})$
30	4D-basket put	$\bm{\beta}=(-\text{3},-\text{2},-\text{1},\text{0})$ ,	3.307	$(\text{4.0},\text{4.2},\text{4.2},\text{4.2})$
		$\alpha=\text{15}$ , $\delta=\text{0.4},\bm{\Delta}=\bm{I}_{\text{4}},K=\text{100}$	$(\text{3.0}\text{e}^{-\text{4}})$
31	4D-call on min	$\bm{\beta}=(-\text{3},-\text{3},-\text{3},-\text{3})$ ,	0.17374	$(-\text{8.8},-\text{8.8},-\text{8.8},-\text{8.8})$
		$\alpha=\text{15}$ , $\delta=\text{0.4},\bm{\Delta}=\bm{I}_{\text{4}},K=\text{100}$	$(\text{5.0}\text{e}^{-\text{5}})$
32	4D-call on min	$\bm{\beta}=(-\text{3},-\text{2},-\text{1},\text{0})$ ,	0.20327	$(-\text{6.5},-\text{6.4},-\text{6.3},-\text{6.2})$
		$\alpha=\text{15}$ , $\delta=\text{0.4},\bm{\Delta}=\bm{I}_{\text{4}},K=\text{100}$	$(\text{7.0}\text{e}^{-\text{5}})$
33	6D-basket put	$\bm{\beta}=(-\text{3},-\text{3},-\text{3},-\text{3},-\text{3},-\text{3})$ ,	0.01039	$(\text{3.1},\text{3.1},\text{3.1},\text{3.1},\text{3.1},\text{3.1})$
		$\alpha=\text{15},\delta=\text{0.2},\bm{\Delta}=\bm{I}_{\text{6}},K=\text{80}$	$(\text{2.0}\text{e}^{-\text{5}})$
34	6D-basket put	$\bm{\beta}=(-\text{3},-\text{2},-\text{1},\text{0},\text{1},\text{2})$ ,	$\text{4.39}\text{e}^{-\text{4}}$	$(\text{4.5},\text{4.6},\text{4.7},\text{4.8},\text{4.8},\text{4.9})$
		$\alpha=\text{15},\delta=\text{0.2},\bm{\Delta}=\bm{I}_{\text{6}},K=\text{80}$	$(\text{3.0}\text{e}^{-\text{6}})$
35	6D-call on min	$\bm{\beta}=(-\text{3},-\text{3},-\text{3},-\text{3},-\text{3},-\text{3})$ ,	$\text{6.034}\text{e}^{-\text{5}}$	$(-\text{4.0},-\text{4.0},-\text{4.0},-\text{4.0},-\text{4.0},-\text{4.0})$
		$\alpha=\text{15},\delta=\text{0.2},\bm{\Delta}=\bm{I}_{\text{6}},K=\text{110}$	$(\text{4.0}\text{e}^{-\text{6}})$
36	6D-call on min	$\bm{\beta}=(-\text{3},-\text{2},-\text{1},\text{0},\text{1},\text{2})$ ,	$\text{1.572}\text{e}^{-\text{4}}$	$(-\text{3.2},-\text{3.2},-\text{3.1},-\text{3.2},-\text{3.2},-\text{3.2})$
		$\alpha=\text{15}$ , $\delta=\text{0.2},\bm{\Delta}=\bm{I}_{\text{6}},K=\text{110}$	$(\text{2.0}\text{e}^{-\text{6}})$

The numerical results were obtained using a cluster machine with a central processing unit (CPU) with 72 cores, a 2.1 GHz clock speed and a memory per node of 256 GB. Further, the computer code was written in MATLAB (MATLAB 2021). The ASGQ implementation was based on the Sparse Grids Matlab Kit.⁸⁸ 8 URL: https://sites.google.com/view/sparse-grids-kit. (For more details on the implementation, we refer the reader to Piazzola and Tamellini (2022).)

In Section 4.1.1, we demonstrate, through various tested examples, the importance of sparsification and adaptivity in the Fourier space for accelerating quadrature convergence. Moreover, in Section 4.1.2, we reveal the importance of the choice of the damping parameters on the numerical complexity of the used quadrature methods. In Section 4.2, we compare our approach against one of the state-of-the-art Fourier pricing methods (namely, the COS method (Fang and Oosterlee 2009; von Sydow et al 2015)) for 1D and 2D cases, and we show the advantage of our approach when the damping parameters are tuned appropriately. Finally, in Section 4.3, we show that our approach achieves substantial computational gains over the Monte Carlo method for different dimensions and parameter constellations to meet a certain relative error tolerance (of practical interest) that we set to be sufficiently small.

4.1 Combining the optimal damping heuristic rule with hierarchical deterministic quadrature methods

4.1.1 Effect of sparsification and dimension adaptivity

In this section, we analyze the effect of dimension adaptivity and sparsification on the acceleration of the convergence of the relative quadrature error, $\mathcal{E}_{R}$ . We elaborate on the comparison between the TP, SM and ASGQ methods when optimal damping parameters are used. Table 10 summarizes these findings. In the numerical experiments, ASGQ consistently outperformed SM. Moreover, for the 2D options, the performance of the ASGQ and TP methods is model-dependent, with ASGQ being the best method for options under the GBM model. For $d=4$ , for options under the GBM and VG models, ASGQ performs better than TP, which is not the case for options under the NIG model. As for 6D options, ASGQ performs better than TP in most cases. These observations confirm that the effect of adaptivity and sparsification becomes more important as the dimension of the option increases. For illustrative purposes, Figures 3–5 compare ASGQ and TP for 4D options with anisotropic parameter sets under different pricing models when optimal damping parameters are used. Figure 3(a) reveals that, for the 4D-basket put option under the GBM model, the ASGQ method achieves an $\mathcal{E}_{\mathcal{R}}$ below 1% using 13.3% of the work of the TP quadrature. Moreover, Figure 4(a) indicates that, for the 4D-basket put option under the VG model, the ASGQ method achieves an $\mathcal{E}_{\mathcal{R}}$ below 0.1% using 25% of the work of the TP quadrature. In contrast, for the 4D-basket put option under the NIG model, Figure 5(a) reveals that the TP quadrature attains $\mathcal{E}_{\mathcal{R}}$ below 0.1% using 10% of the work of the ASGQ.

GBM: convergence of the relative quadrature error, calE-sub-R, with respect to N for TP, SM and ASGQ methods for European four-asset options, when the optimal damping parameters, bar-bold-R, are used. (a) Example 6 in Table .... (b) Example 8 in Table .... — Figure 3: GBM: convergence of the relative quadrature error, $\mathcal{E}_{R}$ , with respect to $N$ for TP, SM and ASGQ methods for European four-asset options, when the optimal damping parameters, $\bar{\bm{R}}$ , are used. (a) Example 6 in Table 5. (b) Example 8 in Table 5.

VG: convergence of the relative quadrature error, calE-sub-R, with respect to N for TP, SM and ASGQ methods for European four-asset options, when the optimal damping parameters, bar-bold-R, are used. (a) Example 18 in Table .... (b) Example 20 in Table .... — Figure 4: VG: convergence of the relative quadrature error, $\mathcal{E}_{R}$ , with respect to $N$ for TP, SM and ASGQ methods for European four-asset options, when the optimal damping parameters, $\bar{\bm{R}}$ , are used. (a) Example 18 in Table 6. (b) Example 20 in Table 6.

NIG: convergence of the relative quadrature error, calE-sub-R, with respect to N for TP, SM and ASGQ methods for European four-asset options, when the optimal damping parameters, bar-bold-R, are used. (a) Example 30 in Table .... (b) Example 32 in Table .... — Figure 5: NIG: convergence of the relative quadrature error, $\mathcal{E}_{R}$ , with respect to $N$ for TP, SM and ASGQ methods for European four-asset options, when the optimal damping parameters, $\bar{\bm{R}}$ , are used. (a) Example 30 in Table 7. (b) Example 32 in Table 7.

4.1.2 Effect of the optimal damping rule

GBM: convergence of the relative quadrature error, calE-sub-R, with respect to N for the ASGQ method for different damping parameter values. (a) Example 6 in Table .... (b) Example 8 in Table .... — Figure 6: GBM: convergence of the relative quadrature error, $\mathcal{E}_{R}$ , with respect to $N$ for the ASGQ method for different damping parameter values. (a) Example 6 in Table 5. (b) Example 8 in Table 5.

In this section, we present the computational benefit of using the optimal damping rule proposed in Section 3.1 on the convergence speed of the relative quadrature error of various methods when pricing the multi-asset European basket and rainbow options. Figures 6–8 illustrate that the optimal damping parameters lead to substantially better error convergence behavior. For instance, Figure 6(a) shows that, for the 4D-basket put option under the GBM model, ASGQ achieves an $\mathcal{E}_{\mathcal{R}}$ below 0.1% using around $N=1500$ quadrature points when using optimal damping parameters, compared with around $N=5000$ points to achieve a similar accuracy for damping parameters shifted by $+1$ in each direction with respect to the optimal values. When using damping parameters shifted by $+2$ in each direction with respect to the optimal values, we do not reach $\mathcal{E}_{\mathcal{R}}=10\%$ , even using $N=5000$ quadrature points. Similarly, for the 4D-call on min option under the VG model, Figure 7(b) shows that ASGQ achieves an $\mathcal{E}_{\mathcal{R}}$ below 0.1% using around $N=500$ quadrature points when using the optimal damping parameters. In contrast, ASGQ cannot achieve $\mathcal{E}_{\mathcal{R}}$ below 1% when using damping parameters shifted by $-1$ in each direction with respect to the optimal values with the same number of quadrature points. Finally, for the 4D-basket put option under the NIG model, Figure 8(a) shows that, when using the optimal damping parameters, the TP quadrature crosses $\mathcal{E}_{\mathcal{R}}=0.1\%$ using 22% of the work it would have used with damping parameters shifted by $-2$ in each direction with respect to the optimal values.

VG: convergence of the relative quadrature error, calE-sub-R, with respect to N for the ASGQ method for different damping parameter values. (a) Example 18 in Table .... (b) Example 20 in Table .... — Figure 7: VG: convergence of the relative quadrature error, $\mathcal{E}_{R}$ , with respect to $N$ for the ASGQ method for different damping parameter values. (a) Example 18 in Table 6. (b) Example 20 in Table 6.

NIG: convergence of the relative quadrature error, calE-sub-R, with respect to N for the TP method for different damping parameter values. (a) Example 30 in Table .... (b) Example 32 in Table .... — Figure 8: NIG: convergence of the relative quadrature error, $\mathcal{E}_{R}$ , with respect to $N$ for the TP method for different damping parameter values. (a) Example 30 in Table 7. (b) Example 32 in Table 7.

In summary, in all experiments, small shifts in both directions with respect to the optimal damping parameters lead to worse error convergence behavior, suggesting that the region of optimality of the damping parameters is tight and that our rule is sufficient to obtain optimal quadrature convergence behavior, irrespective of the method. Moreover, arbitrary choices of damping parameters may lead to extremely poor convergence of the quadrature, as illustrated by the purple curves in Figures 6(a), 6(b), 7(a) and 8(b). All the damping parameters compared belong to the strip of regularity of the integrand $\delta_{V}$ defined in Section 2. Finally, although we provide only a limited number of plots to illustrate these findings, consistent results were observed for the different models and damping parameters.

4.2 Comparison of our approach with the COS method

This section presents an empirical comparison of our proposed approach with the COS method (Fang and Oosterlee 2009). Our aim is to compare the performance of the two methods when they are both appropriately tuned, which seems to be a gap in the benchmarking literature (von Sydow et al 2015; Crisóstomo 2018) due to the absence of guidance on the suitable choice of damping parameters. We use two metrics for the comparison: the CPU time required to achieve a predefined relative error, $\smash{\mathcal{E}_{R}=\mathrm{e}^{-3}}$ (see Tables 8 and 9), and the number of times the characteristic function, $N_{\mathrm{CF}}$ , is evaluated to reach this accuracy (see Figures 9 and 10). The second metric is particularly interesting when the costly part of the approximation formula is the evaluation of the characteristic function, and it has the merit of being independent of the implementation and the used computer characteristics.

In this work, the numerical comparison of our approach with the COS method is restricted for options with up to two underlyings because the implementation of the COS method in higher dimensions is not available to us. Since the CPU time needed to achieve a certain accuracy is highly dependent on the way the methods are implemented, we did not use the sparse grids kit (Piazzola and Tamellini 2022) for the comparison. To the extent possible, both of the methods were implemented in a similar manner so as to give reproducible results. To aid a fair comparison, we compare the optimal damping rule with isotropic TP quadrature and the Gauss–Laguerre rule (ODTPQ) with the isotropic version of the COS method. The reported CPU times in Tables 8 and 9 are given in seconds and are computed from the average over $10^{3}$ replications of each experiment. Moreover, the ODTPQ CPU time in Tables 8 and 9 includes both the cost of the quadrature and the cost of the optimization to obtain the damping parameters, $\bar{\bm{R}}$ . Reference values are computed using Monte Carlo simulation with $M=10^{9}$ samples.

4.2.1 Implementation details of the COS method

The 2D-COS formula to approximate the option value is given by (Ruijter and Oosterlee 2012)

	$\displaystyle V(\bm{\varTheta}_{m},\bm{\varTheta}_{p})$	$\displaystyle\approx\mathrm{e}^{-rT}\frac{(b-a)^{2}}{4}\sum_{k_{1}=0}^{N_{% \mathrm{COS}}}\sum_{k_{2}=0}^{N_{\mathrm{COS}}}\frac{1}{2}\bigg{(}\Re\bigg{\{}% \phi_{\bm{X}_{T}}\bigg{(}\frac{k_{1}\pi}{b-a},\frac{k_{2}\pi}{b-a}\bigg{)}\exp% \bigg{(}\mathrm{i}k_{1}\pi\frac{X_{0}^{1}-a}{b-a}+\mathrm{i}k_{2}\pi\frac{X_{0% }^{2}-a}{b-a}\bigg{)}\bigg{\}}$
		$\displaystyle\qquad+\Re\bigg{\{}\phi_{\bm{X}_{T}}\bigg{(}\frac{k_{1}\pi}{b-a},% -\frac{k_{2}\pi}{b-a}\bigg{)}\exp\bigg{(}\mathrm{i}k_{1}\pi\frac{X_{0}^{1}-a}{% b-a}-\mathrm{i}k_{2}\pi\frac{X_{0}^{2}-a}{b-a}\bigg{)}\bigg{\}}\bigg{)}P_{k_{1% },k_{2}}(T),$		(4.2)

where

$\varPhi_{\bm{X}_{T}}(\bm{u})=\mathrm{e}^{\mathrm{i}\langle\bm{u},\bm{X}_{0}% \rangle}\phi_{\bm{X}_{T}}(\bm{u})$

is the characteristic function (see Table 2), and $\smash{P_{k_{1},k_{2}}}$ are the Fourier cosine coefficients of the payoff function $P(\cdot)$ . We use the isotropic version where the number of Fourier modes is the same in each dimension, $N_{1}=N_{2}=N_{\mathrm{COS}}$ . For the truncation range, we use the domain $[a,b]^{d}$ , as suggested by Ruijter and Oosterlee (2012, Section 5.2), which is given by

$\left.\begin{aligned} \displaystyle a&\displaystyle=\min_{1\leq i\leq 2}\Big{% \{}X_{0}^{i}+c^{i}_{1}-L\sqrt{c^{i}_{2}+\sqrt{c^{i}_{4}}}\Big{\}},\\ \\ \displaystyle b&\displaystyle=\max_{1\leq i\leq 2}\Big{\{}X_{0}^{i}+c^{i}_{1}+% L\sqrt{c^{i}_{2}+\sqrt{c^{i}_{4}}}\Big{\}},\end{aligned}\right\}$

(4.3)

where $\smash{c_{n}^{i}}$ is the $n$ th cumulant of the random variable $\smash{X_{T}^{i}}$ and $L=10$ . For the table of cumulants used in this work, we refer the reader to Oosterlee and Grzelak (2019). In the case of 2D-basket put and 2D-call on min options, we approximate $\smash{\{P_{k_{1},k_{2}}\}_{k_{1},k_{2}=0}^{N_{\mathrm{COS}}-1}}$ numerically using a discrete cosine transform (the dct2 function in MATLAB). The number of terms used in each spatial dimension of the DCT approximation is denoted by $Q$ (see Ruijter and Oosterlee (2012, Section 3.2.1) for more details). To the best of our knowledge, there is no rule for the choice of $Q$ . Ruijter and Oosterlee (2012) state only that the number of terms used in the DCT to approximate the payoff cosine coefficients must satisfy $Q\geq\max(N_{1},N_{2})=N_{\mathrm{COS}}$ . We observe that the choice $Q=N_{\mathrm{COS}}$ may result in oscillatory behavior in the error convergence. For this reason, we use $Q=1000$ to ensure that the payoff cosine coefficients are approximated accurately, as shown by Ruijter and Oosterlee (2012). We note that $\smash{N_{\mathrm{CF}}=2^{d-1}N_{\mathrm{COS}}^{d}}$ in $d$ dimensions.

4.2.2 Numerical experiments

Convergence of the relative error with respect to N-sub-CF, the number of characteristic function evaluations of both COS and ODTPQ under the VG model. (a) 2D-call on min. (b) 2D-basket put. The model used, payoff and damping parameters are given in Tables ... and .... — Figure 9: Convergence of the relative error with respect to $N_{\text{CF}}$ , the number of characteristic function evaluations of both COS and ODTPQ under the VG model. (a) 2D-call on min. (b) 2D-basket put. The model used, payoff and damping parameters are given in Tables 8 and 9.

Convergence of the relative error with respect to N-sub-CF, the number of characteristic function evaluations of both COS and ODTPQ under the NIG model. (a) 2D-call on min. (b) 2D-basket put. The model used, payoff and damping parameters are given in Tables ... and .... — Figure 10: Convergence of the relative error with respect to $N_{\text{CF}}$ , the number of characteristic function evaluations of both COS and ODTPQ under the NIG model. (a) 2D-call on min. (b) 2D-basket put. The model used, payoff and damping parameters are given in Tables 8 and 9.

Figures 9 and 10 show that if the damping parameters are chosen appropriately based on the proposed rule in (3.5), our approach achieves a desired relative error with significantly fewer characteristic function evaluations, $N_{\mathrm{CF}}$ , for both 2D-call on min and 2D-basket put options tested under VG and NIG. For instance, Figure 9(a) shows that the COS method achieves an $\mathcal{E}_{R}=\mathrm{e}^{-3}$ using $N_{\mathrm{CF}}=25\,538$ , whereas ODTPQ reaches $\mathcal{E}_{R}=6\mathrm{e}^{-4}$ using only $N_{\mathrm{CF}}=108$ .

Table 8 demonstrates that the ODTPQ approach achieves the relative error, $\mathcal{E}_{R}=\mathrm{e}^{-3}$ , approximately 13 to 29.5 times faster than the COS method for the tested 2D-call on min option under GBM, VG and NIG. Moreover, Table 9 shows that ODTPQ reaches the relative tolerance, $\mathcal{E}_{R}=\mathrm{e}^{-3}$ , approximately 2.5 to 6.5 times faster than the COS method for 2D-basket put options under GBM, VG and NIG.

Table 8: 2D-call on min: CPU time in seconds of ODTPQ and COS to achieve the relative error

\mathcal{E}_{R}=\text{1.0}\text{e}^{-\text{3}}

with

\bm{S}_{\text{0}}=(\text{100},\text{100})

,

K=\text{100}

,

r=\text{0}

,

T=\text{1}

.

			ODTPQ	COS
Model	Parameters	$\bar{\bm{R}}$	CPU time	CPU time
GBM	$\bm{\sigma}=(\text{0.2},\text{0.8})$ , $\bm{C}=\bm{I}_{d}$	$-(\text{7.18},\text{1.65})$	$\text{2.4}\text{e}^{-\text{3}}$	$\text{4.2}\text{e}^{-\text{2}}$
VG	$\bm{\sigma}=(\text{0.2},\text{0.8})$ ,	$-(\text{7.38},\text{1.79})$	$\text{2.2}\text{e}^{-\text{3}}$	$\text{6.5}\text{e}^{-\text{2}}$
	$\bm{\theta}=(-\text{0.3},-\text{0.1})$ , $\nu=\text{0.5}$
NIG	$\bm{\beta}=(-\text{3},-\text{3})$ ,	$-(\text{6.88},\text{6.88})$	$\text{1.8}\text{e}^{-\text{3}}$	$\text{2.36}\text{e}^{-\text{2}}$
	$\bm{\Delta}=\bm{I}_{d}$ , $\alpha=\text{15}$ , $\delta=\text{0.5}$

Table 9: 2D-basket put: CPU time in seconds of ODTPQ and COS to achieve the relative error

\mathcal{E}_{R}=\text{1.0}\text{e}^{-\text{3}}

with

\bm{S}_{\text{0}}=(\text{100},\text{100})

,

K=\text{100}

,

r=\text{0}

,

T=\text{1}

.

			ODTPQ	COS
Model	Parameters	$\bar{\bm{R}}$	CPU time	CPU time
GBM	$\bm{\sigma}=(\text{0.2},\text{0.8})$ , $\bm{C}=\bm{I}_{d}$	$(\text{3.05},\text{1.36})$	$\text{5.4}\text{e}^{-\text{3}}$	$\text{2.7}\text{e}^{-\text{2}}$
VG	$\bm{\sigma}=(\text{0.2},\text{0.8})$ ,	$(\text{1.81},\text{0.89})$	$\text{9.1}\text{e}^{-\text{3}}$	$\text{2.5}\text{e}^{-\text{2}}$
	$\bm{\theta}=(-\text{0.3},-\text{0.1})$ , $\nu=\text{0.5}$
NIG	$\bm{\beta}=(-\text{3},-\text{3})$ ,	$(\text{4.5},\text{4.5})$	$\text{3.6}\text{e}^{-\text{3}}$	$\text{2.4}\text{e}^{-\text{2}}$
	$\bm{\Delta}=\bm{I}_{d}$ , $\alpha=\text{15}$ , $\delta=\text{0.5}$

Remark 4.1 (About the COS method in multiple dimensions).

For insights on the performance of the COS method in more than two dimensions, we refer the reader to Junike and Stier (2023), where the numerical experiments indicate that the Monte Carlo method outperforms the COS method for cash-or-nothing put options with more than three underlyings if the target error tolerance is of order $\mathrm{e}^{-3}$ .

4.3 Computational comparison of quadrature methods with optimal damping and Monte Carlo

This section compares the Monte Carlo method and our proposed approach based on the best quadrature method in the Fourier space combined with the optimal damping parameters in terms of errors and computation time. The comparison is performed for all option examples in Tables 5–7. We fix a sufficiently small relative error tolerance in the price estimates, and then we compare the time taken by different methods to compute it. This is done in the following way.

(1)
We find the least number of quadrature points to reach a predefined relative quadrature error.
(2)
We estimate, using the CLT formula (4.1), the required number of Monte Carlo samples to achieve the same relative error achieved by the quadrature method.
(3)
We compare the CPU times of both methods, including the cost of numerical optimization of (3.6) preceding the numerical quadrature for the Fourier approach. The Monte Carlo CPU time is obtained through an average of 10 runs.

The results presented in Table 10 show that our approach significantly outperforms the Monte Carlo method for all the tested options with various models, parameter sets and dimensions. In particular, for all tested 2D and 4D options, the proposed approach requires less than 20% (even less than 1% in most cases) of the Monte Carlo work to achieve a total relative error below 0.1%. In general, these gains degrade for the tested 6D options. For example 21 in Table 6, this approach requires around 43% of the work of Monte Carlo in order to achieve a total relative error below 1%. The magnitude of the CPU gain varies depending on different factors, such as the model and payoff parameters affecting the integrand differently in physical space (related to the Monte Carlo estimator variance), and the integrand regularity in Fourier space (related to the quadrature error for quadrature methods). Finally, numerical experiments suggest that the advantage of employing ASGQ over TP is more pronounced when pricing options with more than two dimensions, except for the NIG model, where TP quadrature performs exceptionally well even in the 4D case. Nevertheless, our empirical results demonstrate that in the case of six dimensions or more, the use of ASGQ is preferable to TP for all pricing models.

Table 10: Comparison of the Fourier approach combined with the optimal damping rule and the best quadrature (quad) method with the Gauss–Laguerre rule and the Monte Carlo method for the European basket and rainbow options under the multivariate GBM, VG and NIG pricing dynamics for various dimensions. [Tables 5–7 present the selected parameter sets for each pricing model, the reference values with their corresponding statistical errors and the optimal damping parameters. CPU times are given in seconds.

M

is the number of Monte Carlo samples.

N

is the number of quadrature points. MC denotes Monte Carlo. CPU time ratio (Quad/MC) values are given in percent.]

	Best		MC		Quad		CPU
Example	quad	$\mathcal{E}_{\bm{R}}$	CPU time	$𝑴$	CPU time	$𝑵$	time ratio
0 1 in Table 5	ASGQ	$\text{7.00}\text{e}^{-\text{4}}$	00 7.36	$\text{1.2}\times\text{10}^{\text{7}}$	0.63	0 33	8.5
0 2 in Table 5	ASGQ	$\text{3.70}\text{e}^{-\text{4}}$	0 20.7	$\text{3.3}\times\text{10}^{\text{7}}$	0.65	0 67	3.14
13 in Table 6	TP	$\text{2.90}\text{e}^{-\text{4}}$	0 44	$\text{8.8}\times\text{10}^{\text{7}}$	0.25	0 64	0.57
14 in Table 6	TP	$\text{1.80}\text{e}^{-\text{4}}$	0 70.9	$\text{1.4}\times\text{10}^{\text{8}}$	0.23	0 64	0.32
25 in Table 7	TP	$\text{2.90}\text{e}^{-\text{4}}$	0 75.3	$\text{1.1}\times\text{10}^{\text{8}}$	0.2	0 36	0.26
26 in Table 7	TP	$\text{5.86}\text{e}^{-\text{4}}$	0 17.2	$\text{2.6}\times\text{10}^{\text{7}}$	0.2	0 25	1.16
0 3 in Table 5	ASGQ	$\text{7.00}\text{e}^{-\text{4}}$	0 47.3	$\text{7.6}\times\text{10}^{\text{7}}$	0.6	0 37	1.26
0 4 in Table 5	ASGQ	$\text{5.80}\text{e}^{-\text{4}}$	102	$\text{1.4}\times\text{10}^{\text{8}}$	0.63	0 37	0.62
15 in Table 6	ASGQ	$\text{8.26}\text{e}^{-\text{4}}$	0 19.5	$\text{4.1}\times\text{10}^{\text{7}}$	0.54	0 25	2.77
16 in Table 6	TP	$\text{5.37}\text{e}^{-\text{4}}$	0 87.1	$\text{1.4}\times\text{10}^{\text{8}}$	0.16	0 49	0.18
26 in Table 7	TP	$\text{6.70}\text{e}^{-\text{4}}$	0 35.8	$\text{5.3}\times\text{10}^{\text{7}}$	0.22	100	0.61
27 in Table 7	TP	$\text{6.46}\text{e}^{-\text{4}}$	0 42.2	$\text{6.5}\times\text{10}^{\text{7}}$	0.22	0 64	0.52
0 5 in Table 5	ASGQ	$\text{2.46}\text{e}^{-\text{4}}$	0 207	$\text{1.00}\times\text{10}^{\text{8}}$	7.8	0 5 257	0 3.77
0 6 in Table 5	ASGQ	$\text{8.12}\text{e}^{-\text{4}}$	00 14.5	$\text{7.90}\times\text{10}^{\text{6}}$	2.73	0 1 433	18.83
17 in Table 6	ASGQ	$\text{2.58}\text{e}^{-\text{4}}$	0 106.3	$\text{1.23}\times\text{10}^{\text{8}}$	5	0 3 013	0 4.7
18 in Table 6	ASGQ	$\text{3.58}\text{e}^{-\text{4}}$	00 38.7	$\text{4.50}\times\text{10}^{\text{7}}$	2	0 1 109	0 5.17
27 in Table 7	TP	$\text{4.57}\text{e}^{-\text{4}}$	00 50.2	$\text{4.70}\times\text{10}^{\text{7}}$	0.5	00 256	0 1
28 in Table 7	TP	$\text{4.10}\text{e}^{-\text{4}}$	00 49.4	$\text{4.80}\times\text{10}^{\text{7}}$	0.52	00 256	0 1
0 7 in Table 5	ASGQ	$\text{5.70}\text{e}^{-\text{4}}$	1 147	$\text{7.00}\times\text{10}^{\text{8}}$	1	00 435	0 0.09
0 8 in Table 5	ASGQ	$\text{5.50}\text{e}^{-\text{4}}$	1 580	$\text{9.60}\times\text{10}^{\text{8}}$	0.95	00 654	0 0.06
19 in Table 6	ASGQ	$\text{5.90}\text{e}^{-\text{4}}$	0 220	$\text{3.00}\times\text{10}^{\text{8}}$	1.25	00 567	0 0.57
20 in Table 6	ASGQ	$\text{8.90}\text{e}^{-\text{4}}$	0 249	$\text{3.30}\times\text{10}^{\text{8}}$	1.4	00 862	0 0.56
29 in Table 7	TP	$\text{7.20}\text{e}^{-\text{4}}$	0 193.5	$\text{2.00}\times\text{10}^{\text{8}}$	8.7	20 736	0 4.5
30 in Table 7	TP	$\text{4.20}\text{e}^{-\text{4}}$	0 716	$\text{7.80}\times\text{10}^{\text{8}}$	0.8	0 2 401	0 0.11
0 9 in Table 5	ASGQ	$\text{2.90}\text{e}^{-\text{2}}$	00 18.53	$\text{5.5}\times\text{10}^{\text{6}}$	0 2	0 318	11
10 in Table 5	ASGQ	$\text{3.30}\text{e}^{-\text{3}}$	0 548	$\text{1.5}\times\text{10}^{\text{8}}$	0 2.1	0 340	0 0.38
21 in Table 6	ASGQ	$\text{7.80}\text{e}^{-\text{3}}$	000 5.4	$\text{4.7}\times\text{10}^{\text{6}}$	0 2.3	0 453	42.6
22 in Table 6	ASGQ	$\text{5.40}\text{e}^{-\text{3}}$	00 31.5	$\text{2.5}\times\text{10}^{\text{7}}$	0 3.5	0 566	11
31 in Table 7	ASGQ	$\text{1.47}\text{e}^{-\text{2}}$	00 14.2	$\text{1.0}\times\text{10}^{\text{7}}$	0 3.4	0 616	24
32 in Table 7	TP	$\text{3.75}\text{e}^{-\text{2}}$	00 33.5	$\text{2.5}\times\text{10}^{\text{7}}$	11.7	4 096	35
11 in Table 5	ASGQ	$\text{1.40}\text{e}^{-\text{3}}$	2 635	$\text{6.9}\times\text{10}^{\text{8}}$	0 6	3 070	0 0.23
12 in Table 5	ASGQ	$\text{1.70}\text{e}^{-\text{3}}$	2 110	$\text{5.3}\times\text{10}^{\text{8}}$	0 4.5	1 642	0 0.21
23 in Table 6	ASGQ	$\text{2.00}\text{e}^{-\text{3}}$	00 85	$\text{6.8}\times\text{10}^{\text{7}}$	19.5	7 401	23
24 in Table 6	ASGQ	$\text{2.60}\text{e}^{-\text{3}}$	0 360	$\text{2.8}\times\text{10}^{\text{8}}$	0 4.6	1 671	0 1.28
33 in Table 7	ASGQ	$\text{5.70}\text{e}^{-\text{2}}$	00 85.5	$\text{6.3}\times\text{10}^{\text{7}}$	0 1	0 105	0 1.17
34 in Table 7	ASGQ	$\text{3.79}\text{e}^{-\text{2}}$	0 108	$\text{7.5}\times\text{10}^{\text{7}}$	0 1.4	0 340	0 1.3

5 Conclusions and future work

In this work, we introduced what we believe to be a novel heuristic rule for an optimal choice of the damping parameters that arise in Lewis’s valuation formula (Lewis 2001) to ensure analyticity of the integrand and to achieve highly efficient pricing of European-style multi-asset options. The rule was motivated analytically by a quadrature error bound controlled by the infimum norm of the integrand and was shown to be very effective through numerical examples. We also proposed the use of ASGQ to alleviate the curse of dimensionality and profit from the high smoothness of the integrand in the Fourier space. Numerical experiments showed that our Fourier-based approach combining the optimal damping rule and the ASGQ under the Gauss–Laguerre measure consistently outperformed the Monte Carlo method for the pricing of multivariate basket and rainbow options under various pricing models (namely, GBM, VG and NIG) and for different parameter constellations and dimensions in these models.

Declaration of interest

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

Acknowledgements

Christian Bayer gratefully acknowledges support from the German Research Foundation (DFG) via the Cluster of Excellence MATH $+$ (Project AA4-2). This paper is based on work supported by the King Abdullah University of Science and Technology (KAUST) Office of Sponsored Research (OSR) under award OSR-2019-CRG8-4033 and by the Alexander von Humboldt Foundation. Antonis Papapantoleon gratefully acknowledges the financial support from the Hellenic Foundation for Research and Innovation (grant HFRI-FM17-2152).

References

Aguilar, J.-P. (2020). Some pricing tools for the variance Gamma model. International Journal of Theoretical and Applied Finance 23(4), Paper 2050025 (https://doi.org/10.1142/S0219024920500259).
Babuška, I., Nobile, F., and Tempone, R. (2007). A stochastic collocation method for elliptic partial differential equations with random input data. SIAM Journal on Numerical Analysis 45(3), 1005–1034 (https://doi.org/10.1137/050645142).
Barndorff-Nielsen, O. (1977). Exponentially decreasing distributions for the logarithm of particle size. Proceedings of the Royal Society A 353(1674), 401–419 (https://doi.org/10.1098/rspa.1977.0041).
Barndorff-Nielsen, O. E. (1997a). Normal inverse Gaussian distributions and stochastic volatility modelling. Scandinavian Journal of Statistics 24(1), 1–13 (https://doi.org/10.1111/1467-9469.00045).
Barndorff-Nielsen, O. E. (1997b). Processes of normal inverse Gaussian type. Finance and Stochastics 2(1), 41–68 (https://doi.org/10.1007/s007800050032).
Barthelmann, V., Novak, E., and Ritter, K. (2000). High dimensional polynomial interpolation on sparse grids. Advances in Computational Mathematics 12(4), 273–288 (https://doi.org/10.1023/A:1018977404843).
Baschetti, F., Bormetti, G., Romagnoli, S., and Rossi, P. (2022). The SINC way: a fast and accurate approach to Fourier pricing. Quantitative Finance 22(3), 427–446 (https://doi.org/10.1080/14697688.2021.1965192).
Bayer, C., Siebenmorgen, M., and Tempone, R. (2018). Smoothing the payoff for efficient computation of basket option prices. Quantitative Finance 18(3), 491–505 (https://doi.org/10.1080/14697688.2017.1308003).
Bayer, C., Ben Hammouda, C., and Tempone, R. (2020a). Hierarchical adaptive sparse grids and quasi-Monte Carlo for option pricing under the rough Bergomi model. Quantitative Finance 20(9), 1457–1473 (https://doi.org/10.1080/14697688.2020.1744700).
Bayer, C., Ben Hammouda, C., and Tempone, R. (2020b). Multilevel Monte Carlo with numerical smoothing for robust and efficient computation of probabilities and densities. Preprint, arXiv (https://doi.org/10.48550/arXiv.2003.05708).
Bayer, C., Ben Hammouda, C., and Tempone, R. (2023). Numerical smoothing with hierarchical adaptive sparse grids and quasi-Monte Carlo methods for efficient option pricing. Quantitative Finance 23(2), 209–227 (https://doi.org/10.1080/14697688.2022.2135455).
Ben Hammouda, C. (2020). Hierarchical approximation methods for option pricing and stochastic reaction networks. PhD Thesis, King Abdullah University of Science and Technology, Thuwal, Kingdom of Saudi Arabia.
Byrd, R. H., Hribar, M. E., and Nocedal, J. (1999). An interior point algorithm for large-scale nonlinear programming. SIAM Journal on Optimization 9(4), 877–900 (https://doi.org/10.1137/S1052623497325107).
Byrd, R. H., Gilbert, J. C., and Nocedal, J. (2000). A trust region method based on interior point techniques for nonlinear programming. Mathematical Programming 89(1), 149–185 (https://doi.org/10.1007/PL00011391).
Carlsson, M., and Wittsten, J. (2017). A note on holomorphic functions and the Fourier–Laplace transform. Mathematica Scandinavica 120(2), 225–248 (https://doi.org/10.7146/math.scand.a-25612).
Carr, P., and Madan, D. B. (1999). Option valuation using the fast Fourier transform. The Journal of Computational Finance 2(4), 61–73 (https://doi.org/10.21314/JCF.1999.043).
Carr, P., and Wu, L. (2004). Time-changed Lévy processes and option pricing. Journal of Financial Economics 71(1), 113–141 (https://doi.org/10.1016/S0304-405X(03)00171-5).
Chau, K. W., and Oosterlee, C. W. (2019). Exploration of a cosine expansion lattice scheme. Preprint, arXiv (https://doi.org/10.48550/arXiv.1907.02758).
Chen, P. (2018). Sparse quadrature for high-dimensional integration with Gaussian measure. ESAIM: Mathematical Modelling and Numerical Analysis 52(2), 631–657 (https://doi.org/10.1051/m2an/2018012).
Choi, J. (2018). Sum of all Black–Scholes–Merton models: an efficient pricing method for spread, basket, and Asian options. Journal of Futures Markets 38(6), 627–644 (https://doi.org/10.1002/fut.21909).
Colldeforns-Papiol, G., Ortiz-Gracia, L., and Oosterlee, C. W. (2017). Two-dimensional Shannon wavelet inverse Fourier technique for pricing European options. Applied Numerical Mathematics 117, 115–138 (https://doi.org/10.1016/j.apnum.2017.03.002).
Cont, R., and Tankov, P. (2003). Financial Modelling with Jump Processes. Chapman & Hall/CRC, Boca Raton, FL.
Crisóstomo, R. (2018). Speed and biases of Fourier-based pricing choices: a numerical analysis. International Journal of Computer Mathematics 95(8), 1565–1582 (https://doi.org/10.1080/00207160.2017.1322691).
Davis, P., and Rabinowitz, P. (1954). On the estimation of quadrature errors for analytic functions. Mathematical Tables and Other Aids to Computation 8(48), 193–203 (https://doi.org/10.2307/2002092).
Davis, P. J., and Rabinowitz, P. (2007). Methods of Numerical Integration. Courier Corporation, North Chelmsford, MA.
Donaldson, J. D. (1973). Estimates of upper bounds for quadrature errors. SIAM Journal on Numerical Analysis 10(1), 13–22 (https://doi.org/10.1137/0710003).
Donaldson, J. D., and Elliott, D. (1972). A unified approach to quadrature rules with asymptotic estimates of their remainders. SIAM Journal on Numerical Analysis 9(4), 573–602 (https://doi.org/10.1137/0709051).
Duffie, D., Filipović, D., and Schachermayer, W. (2003). Affine processes and applications in finance. Annals of Applied Probability 13(3), 984–1053 (https://doi.org/10.1214/aoap/1060202833).
Eberlein, E., Glau, K., and Papapantoleon, A. (2010). Analysis of Fourier transform valuation formulas and applications. Applied Mathematical Finance 17(3), 211–240 (https://doi.org/10.1080/13504860903326669).
Elliott, D. (1970). Uniform asymptotic expansions of the classical orthogonal polynomials and some associated functions. Technical Report 21, Mathematics Department, The University of Tasmania Hobart, Tasmania.
Elliott, D., and Tuan, P. D. (1974). Asymptotic estimates of Fourier coefficients. SIAM Journal on Mathematical Analysis 5(1), 1–10 (https://doi.org/10.1137/0505001).
Fang, F., and Oosterlee, C. W. (2009). A novel pricing method for European options based on Fourier-cosine series expansions. SIAM Journal on Scientific Computing 31(2), 826–848 (https://doi.org/10.1137/080718061).
Gautschi, W. (2004). Orthogonal Polynomials: Computation and Approximation. Oxford University Press (https://doi.org/10.1093/oso/9780198506720.001.0001).
Gerstner, T., and Griebel, M. (1998). Numerical integration using sparse grids. Numerical Algorithms 18(3), 209–232 (https://doi.org/10.1023/A:1019129717644).
Gerstner, T., and Griebel, M. (2003). Dimension-adaptive tensor-product quadrature. Computing 71, 65–87 (https://doi.org/10.1007/s00607-003-0015-5).
Glasserman, P. (2003). Monte Carlo Methods in Financial Engineering. Stochastic Modelling and Applied Probability, Volume 53. Springer (https://doi.org/10.1007/978-0-387-21617-1).
Guillaume, T. (2008). Making the best of best-of. Review of Derivatives Research 11(1), 1–39 (https://doi.org/10.1007/s11147-008-9022-1).
Healy, J. (2021). The pricing of vanilla options with cash dividends as a classic vanilla basket option problem. Preprint, arXiv (https://doi.org/10.48550/arXiv.2106.12971).
Hörmander, L. (2015). The Analysis of Linear Partial Differential Operators. I. Distribution Theory and Fourier Analysis. Classics in Mathematics. Springer (https://doi.org/10.1007/978-3-642-61497-2).
Hubalek, F., and Kallsen, J. (2005). Variance-optimal hedging and Markowitz-efficient portfolios for multivariate processes with stationary independent increments with and without constraints. Working Paper, TU München, Germany.
Hurd, T. R., and Zhou, Z. (2010). A Fourier transform method for spread option pricing. SIAM Journal on Financial Mathematics 1(1), 142–157 (https://doi.org/10.1137/090750421).
Junike, G., and Pankrashkin, K. (2022). Precise option pricing by the COS method: how to choose the truncation range. Applied Mathematics and Computation 421, Paper 126935 (https://doi.org/10.1016/j.amc.2022.126935).
Junike, G., and Stier, H. (2023). The multidimensional COS method for option pricing. Preprint, arXiv (https://doi.org/10.48550/arXiv.2307.12843).
Kahl, C., and Lord, R. (2012). Fourier inversion methods in finance. In Handbook of Computational Finance, Duan, J.-C., Härdle, W. K., and Gentle, J. E. (eds). Springer.
Kirkby, J. L. (2015). Efficient option pricing by frame duality with the fast Fourier transform. SIAM Journal on Financial Mathematics 6(1), 713–747 (https://doi.org/10.1137/140989480).
Kirkby, J. L., Nguyen, D. H., and Nguyen, D. (2020). A general continuous time Markov chain approximation for multi-asset option pricing with systems of correlated diffusions. Applied Mathematics and Computation 386, Paper 125472 (https://doi.org/10.1016/j.amc.2020.125472).
Kwok, Y. K., Leung, K. S., and Wong, H. Y. (2012). Efficient options pricing using the fast Fourier transform. In Handbook of Computational Finance, Duan, J.-C., Härdle, W. K., and Gentle, J. E. (eds), pp. 579–604. Springer (https://doi.org/10.1007/978-3-642-17254-0_21).
Lee, R. W. (2004). Option pricing by transform methods: extensions, unification and error control. The Journal of Computational Finance 7(3), 51–86 (https://doi.org/10.21314/JCF.2004.121).
Leentvaar, C. C. W., and Oosterlee, C. W. (2008). Multi-asset option pricing using a parallel Fourier-based technique. The Journal of Computational Finance 12(1), 1–26 (https://doi.org/10.21314/JCF.2008.184).
Lewis, A. L. (2001). A simple option formula for general jump-diffusion and other exponential Lévy processes. Working Paper, Social Science Research Network (https://doi.org/10.2139/ssrn.282110).
Lord, R., and Kahl, C. (2007). Optimal Fourier inversion in semi-analytical option pricing. Discussion Paper 2006-066/2, Tinbergen Institute, Amsterdam (https://doi.org/10.2139/ssrn.921336).
Lord, R., Fang, F., Bervoets, F., and Oosterlee, C. W. (2008). A fast and accurate FFT-based method for pricing early-exercise options under Lévy processes. SIAM Journal on Scientific Computing 30(4), 1678–1705 (https://doi.org/10.1137/070683878).
Luciano, E., and Schoutens, W. (2006). A multivariate jump-driven financial asset model. Quantitative Finance 6(5), 385–402 (https://doi.org/10.1080/14697680600806275).
Lukacs, E. (1970). Characteristic Functions. Charles Griffin & Co.
Madan, D. B., and Seneta, E. (1990). The variance Gamma (VG) model for share market returns. Journal of Business 63(4), 511–524 (https://doi.org/10.1086/296519).
Margrabe, W. (1978). The value of an option to exchange one asset for another. Journal of Finance 33(1), 177–186 (https://doi.org/10.1111/j.1540-6261.1978.tb03397.x).
MATLAB (2021). MATLAB, Version 9.11. The MathWorks Inc., Natick, MA.
Mijatović, A., and Pistorius, M. (2013). Continuously monitored barrier options under Markov processes. Mathematical Finance 23(1), 1–38 (https://doi.org/10.1111/j.1467-9965.2011.00486.x).
Oosterlee, C. W., and Grzelak, L. A. (2019). Mathematical Modeling and Computation in Finance: With Exercises and Python and MATLAB Computer Codes. World Scientific (https://doi.org/10.1142/q0236).
Piazzola, C., and Tamellini, L. (2022). The Sparse Grids Matlab kit: a Matlab implementation of sparse grids for high-dimensional function approximation and uncertainty quantification. Preprint, arXiv:2203.09314v1 [cs.MS] (https://doi.org/10.48550/arXiv.2203.09314).
Raible, S. (2000). Lévy processes in finance: theory, numerics, and empirical facts. PhD Thesis, University of Freiburg, Germany. URL: https://d-nb.info/961285192/34.
Ruijter, M. J., and Oosterlee, C. W. (2012). Two-dimensional Fourier cosine series expansion method for pricing financial options. SIAM Journal on Scientific Computing 34(5), B642–B671 (https://doi.org/10.1137/120862053).
Schoutens, W. (2003). Lévy Processes in Finance: Pricing Financial Derivatives. Wiley (https://doi.org/10.1002/0470870230).
Smolyak, S. A. (1963). Quadrature and interpolation formulas for tensor products of certain classes of functions. Doklady Akademii Nauk 148(5), 1042–1045.
Takahasi, H., and Mori, M. (1971). Estimation of errors in the numerical quadrature of analytic functions. Applicable Analysis 1(3), 201–229 (https://doi.org/10.1080/00036817108839015).
Titchmarsh, E. C. (1948). Introduction to the Theory of Fourier Integrals, 2nd edn. Clarendon Press, Oxford.
Trefethen, L. N. (2008). Is Gauss quadrature better than Clenshaw–Curtis? SIAM Review 50(1), 67–87 (https://doi.org/10.1137/060659831).
van de Wiel, D. P. (2015). Valuation of insurance products using a normal inverse Gaussian distribution. Master’s Thesis, Tilburg University, Tilburg, Netherlands. URL: https://arno.uvt.nl/show.cgi?fid=146193.
von Sydow, L., Höök, L. J., Larsson, E., Lindström, E., Milovanović, S., Persson, J., Shcherbakov, V., Shpolyanskiy, Y., Sirén, S., Toivanen, J., Waldén, J., Wiktorsson, M., Levesley, J., Li, J., Oosterlee, C. W., Ruijter, M. J., Toropov, A., and Zhao, Y. (2015). BENCHOP: the BENCHmarking project in option pricing. International Journal of Computer Mathematics 92(12), 2361–2379 (https://doi.org/10.1080/00207160.2015.1072172).
Wasilkowski, G. W., and Wozniakowski, H. (1995). Explicit cost bounds of algorithms for multivariate tensor product problems. Journal of Complexity 11(1), 1–56 (https://doi.org/10.1006/jcom.1995.1001).
Xi, Y., Ding, K., and Ning, N. (2019). Simultaneous two-dimensional continuous-time Markov chain approximation of two-dimensional fully coupled Markov diffusion processes. Working Paper, Social Science Research Network (https://doi.org/10.2139/ssrn.3461115).
Xiang, S. (2012). Asymptotics on Laguerre or Hermite polynomial expansions and their applications in Gauss quadrature. Journal of Mathematical Analysis and Applications 393(2), 434–444 (https://doi.org/10.1016/j.jmaa.2012.03.056).
Zhang, B., and Oosterlee, C. W. (2013). Efficient pricing of European-style Asian options under exponential Lévy processes based on Fourier cosine expansions. SIAM Journal on Financial Mathematics 4(1), 399–426 (https://doi.org/10.1137/110853339).
Zhu, J. (2009). Applications of Fourier Transform to Smile Modeling: Theory and Implementation. Springer (https://doi.org/10.1007/978-3-642-01808-4).

[bib.bib1] Aguilar, J.-P. (2020). Some pricing tools for the variance Gamma model. International Journal of Theoretical and Applied Finance 23(4), Paper 2050025 (https://doi.org/10.1142/S0219024920500259).

[bib.bib2] Babuška, I., Nobile, F., and Tempone, R. (2007). A stochastic collocation method for elliptic partial differential equations with random input data. SIAM Journal on Numerical Analysis 45(3), 1005–1034 (https://doi.org/10.1137/050645142).

[bib.bib3] Barndorff-Nielsen, O. (1977). Exponentially decreasing distributions for the logarithm of particle size. Proceedings of the Royal Society A 353(1674), 401–419 (https://doi.org/10.1098/rspa.1977.0041).

[bib.bib4] Barndorff-Nielsen, O. E. (1997a). Normal inverse Gaussian distributions and stochastic volatility modelling. Scandinavian Journal of Statistics 24(1), 1–13 (https://doi.org/10.1111/1467-9469.00045).

[bib.bib5] Barndorff-Nielsen, O. E. (1997b). Processes of normal inverse Gaussian type. Finance and Stochastics 2(1), 41–68 (https://doi.org/10.1007/s007800050032).

[bib.bib6] Barthelmann, V., Novak, E., and Ritter, K. (2000). High dimensional polynomial interpolation on sparse grids. Advances in Computational Mathematics 12(4), 273–288 (https://doi.org/10.1023/A:1018977404843).

[bib.bib7] Baschetti, F., Bormetti, G., Romagnoli, S., and Rossi, P. (2022). The SINC way: a fast and accurate approach to Fourier pricing. Quantitative Finance 22(3), 427–446 (https://doi.org/10.1080/14697688.2021.1965192).

[bib.bib8] Bayer, C., Siebenmorgen, M., and Tempone, R. (2018). Smoothing the payoff for efficient computation of basket option prices. Quantitative Finance 18(3), 491–505 (https://doi.org/10.1080/14697688.2017.1308003).

[bib.bib9] Bayer, C., Ben Hammouda, C., and Tempone, R. (2020a). Hierarchical adaptive sparse grids and quasi-Monte Carlo for option pricing under the rough Bergomi model. Quantitative Finance 20(9), 1457–1473 (https://doi.org/10.1080/14697688.2020.1744700).

[bib.bib10] Bayer, C., Ben Hammouda, C., and Tempone, R. (2020b). Multilevel Monte Carlo with numerical smoothing for robust and efficient computation of probabilities and densities. Preprint, arXiv (https://doi.org/10.48550/arXiv.2003.05708).

[bib.bib11] Bayer, C., Ben Hammouda, C., and Tempone, R. (2023). Numerical smoothing with hierarchical adaptive sparse grids and quasi-Monte Carlo methods for efficient option pricing. Quantitative Finance 23(2), 209–227 (https://doi.org/10.1080/14697688.2022.2135455).

[bib.bib12] Ben Hammouda, C. (2020). Hierarchical approximation methods for option pricing and stochastic reaction networks. PhD Thesis, King Abdullah University of Science and Technology, Thuwal, Kingdom of Saudi Arabia.

[bib.bib13] Byrd, R. H., Hribar, M. E., and Nocedal, J. (1999). An interior point algorithm for large-scale nonlinear programming. SIAM Journal on Optimization 9(4), 877–900 (https://doi.org/10.1137/S1052623497325107).

[bib.bib14] Byrd, R. H., Gilbert, J. C., and Nocedal, J. (2000). A trust region method based on interior point techniques for nonlinear programming. Mathematical Programming 89(1), 149–185 (https://doi.org/10.1007/PL00011391).

[bib.bib15] Carlsson, M., and Wittsten, J. (2017). A note on holomorphic functions and the Fourier–Laplace transform. Mathematica Scandinavica 120(2), 225–248 (https://doi.org/10.7146/math.scand.a-25612).

[bib.bib16] Carr, P., and Madan, D. B. (1999). Option valuation using the fast Fourier transform. The Journal of Computational Finance 2(4), 61–73 (https://doi.org/10.21314/JCF.1999.043).

[bib.bib17] Carr, P., and Wu, L. (2004). Time-changed Lévy processes and option pricing. Journal of Financial Economics 71(1), 113–141 (https://doi.org/10.1016/S0304-405X(03)00171-5).

[bib.bib18] Chau, K. W., and Oosterlee, C. W. (2019). Exploration of a cosine expansion lattice scheme. Preprint, arXiv (https://doi.org/10.48550/arXiv.1907.02758).

[bib.bib19] Chen, P. (2018). Sparse quadrature for high-dimensional integration with Gaussian measure. ESAIM: Mathematical Modelling and Numerical Analysis 52(2), 631–657 (https://doi.org/10.1051/m2an/2018012).

[bib.bib20] Choi, J. (2018). Sum of all Black–Scholes–Merton models: an efficient pricing method for spread, basket, and Asian options. Journal of Futures Markets 38(6), 627–644 (https://doi.org/10.1002/fut.21909).

[bib.bib21] Colldeforns-Papiol, G., Ortiz-Gracia, L., and Oosterlee, C. W. (2017). Two-dimensional Shannon wavelet inverse Fourier technique for pricing European options. Applied Numerical Mathematics 117, 115–138 (https://doi.org/10.1016/j.apnum.2017.03.002).

[bib.bib22] Cont, R., and Tankov, P. (2003). Financial Modelling with Jump Processes. Chapman & Hall/CRC, Boca Raton, FL.

[bib.bib23] Crisóstomo, R. (2018). Speed and biases of Fourier-based pricing choices: a numerical analysis. International Journal of Computer Mathematics 95(8), 1565–1582 (https://doi.org/10.1080/00207160.2017.1322691).

[bib.bib24] Davis, P., and Rabinowitz, P. (1954). On the estimation of quadrature errors for analytic functions. Mathematical Tables and Other Aids to Computation 8(48), 193–203 (https://doi.org/10.2307/2002092).

[bib.bib25] Davis, P. J., and Rabinowitz, P. (2007). Methods of Numerical Integration. Courier Corporation, North Chelmsford, MA.

[bib.bib26] Donaldson, J. D. (1973). Estimates of upper bounds for quadrature errors. SIAM Journal on Numerical Analysis 10(1), 13–22 (https://doi.org/10.1137/0710003).

[bib.bib27] Donaldson, J. D., and Elliott, D. (1972). A unified approach to quadrature rules with asymptotic estimates of their remainders. SIAM Journal on Numerical Analysis 9(4), 573–602 (https://doi.org/10.1137/0709051).

[bib.bib28] Duffie, D., Filipović, D., and Schachermayer, W. (2003). Affine processes and applications in finance. Annals of Applied Probability 13(3), 984–1053 (https://doi.org/10.1214/aoap/1060202833).

[bib.bib29] Eberlein, E., Glau, K., and Papapantoleon, A. (2010). Analysis of Fourier transform valuation formulas and applications. Applied Mathematical Finance 17(3), 211–240 (https://doi.org/10.1080/13504860903326669).

[bib.bib30] Elliott, D. (1970). Uniform asymptotic expansions of the classical orthogonal polynomials and some associated functions. Technical Report 21, Mathematics Department, The University of Tasmania Hobart, Tasmania.

[bib.bib31] Elliott, D., and Tuan, P. D. (1974). Asymptotic estimates of Fourier coefficients. SIAM Journal on Mathematical Analysis 5(1), 1–10 (https://doi.org/10.1137/0505001).

[bib.bib32] Fang, F., and Oosterlee, C. W. (2009). A novel pricing method for European options based on Fourier-cosine series expansions. SIAM Journal on Scientific Computing 31(2), 826–848 (https://doi.org/10.1137/080718061).

[bib.bib33] Gautschi, W. (2004). Orthogonal Polynomials: Computation and Approximation. Oxford University Press (https://doi.org/10.1093/oso/9780198506720.001.0001).

[bib.bib34] Gerstner, T., and Griebel, M. (1998). Numerical integration using sparse grids. Numerical Algorithms 18(3), 209–232 (https://doi.org/10.1023/A:1019129717644).

[bib.bib35] Gerstner, T., and Griebel, M. (2003). Dimension-adaptive tensor-product quadrature. Computing 71, 65–87 (https://doi.org/10.1007/s00607-003-0015-5).

[bib.bib36] Glasserman, P. (2003). Monte Carlo Methods in Financial Engineering. Stochastic Modelling and Applied Probability, Volume 53. Springer (https://doi.org/10.1007/978-0-387-21617-1).

[bib.bib37] Guillaume, T. (2008). Making the best of best-of. Review of Derivatives Research 11(1), 1–39 (https://doi.org/10.1007/s11147-008-9022-1).

[bib.bib38] Healy, J. (2021). The pricing of vanilla options with cash dividends as a classic vanilla basket option problem. Preprint, arXiv (https://doi.org/10.48550/arXiv.2106.12971).

[bib.bib39] Hörmander, L. (2015). The Analysis of Linear Partial Differential Operators. I. Distribution Theory and Fourier Analysis. Classics in Mathematics. Springer (https://doi.org/10.1007/978-3-642-61497-2).

[bib.bib40] Hubalek, F., and Kallsen, J. (2005). Variance-optimal hedging and Markowitz-efficient portfolios for multivariate processes with stationary independent increments with and without constraints. Working Paper, TU München, Germany.

[bib.bib41] Hurd, T. R., and Zhou, Z. (2010). A Fourier transform method for spread option pricing. SIAM Journal on Financial Mathematics 1(1), 142–157 (https://doi.org/10.1137/090750421).

[bib.bib42] Junike, G., and Pankrashkin, K. (2022). Precise option pricing by the COS method: how to choose the truncation range. Applied Mathematics and Computation 421, Paper 126935 (https://doi.org/10.1016/j.amc.2022.126935).

[bib.bib43] Junike, G., and Stier, H. (2023). The multidimensional COS method for option pricing. Preprint, arXiv (https://doi.org/10.48550/arXiv.2307.12843).

[bib.bib44] Kahl, C., and Lord, R. (2012). Fourier inversion methods in finance. In Handbook of Computational Finance, Duan, J.-C., Härdle, W. K., and Gentle, J. E. (eds). Springer.

[bib.bib45] Kirkby, J. L. (2015). Efficient option pricing by frame duality with the fast Fourier transform. SIAM Journal on Financial Mathematics 6(1), 713–747 (https://doi.org/10.1137/140989480).

[bib.bib46] Kirkby, J. L., Nguyen, D. H., and Nguyen, D. (2020). A general continuous time Markov chain approximation for multi-asset option pricing with systems of correlated diffusions. Applied Mathematics and Computation 386, Paper 125472 (https://doi.org/10.1016/j.amc.2020.125472).

[bib.bib47] Kwok, Y. K., Leung, K. S., and Wong, H. Y. (2012). Efficient options pricing using the fast Fourier transform. In Handbook of Computational Finance, Duan, J.-C., Härdle, W. K., and Gentle, J. E. (eds), pp. 579–604. Springer (https://doi.org/10.1007/978-3-642-17254-0_21).

[bib.bib48] Lee, R. W. (2004). Option pricing by transform methods: extensions, unification and error control. The Journal of Computational Finance 7(3), 51–86 (https://doi.org/10.21314/JCF.2004.121).

[bib.bib49] Leentvaar, C. C. W., and Oosterlee, C. W. (2008). Multi-asset option pricing using a parallel Fourier-based technique. The Journal of Computational Finance 12(1), 1–26 (https://doi.org/10.21314/JCF.2008.184).

[bib.bib50] Lewis, A. L. (2001). A simple option formula for general jump-diffusion and other exponential Lévy processes. Working Paper, Social Science Research Network (https://doi.org/10.2139/ssrn.282110).

[bib.bib51] Lord, R., and Kahl, C. (2007). Optimal Fourier inversion in semi-analytical option pricing. Discussion Paper 2006-066/2, Tinbergen Institute, Amsterdam (https://doi.org/10.2139/ssrn.921336).

[bib.bib52] Lord, R., Fang, F., Bervoets, F., and Oosterlee, C. W. (2008). A fast and accurate FFT-based method for pricing early-exercise options under Lévy processes. SIAM Journal on Scientific Computing 30(4), 1678–1705 (https://doi.org/10.1137/070683878).

[bib.bib53] Luciano, E., and Schoutens, W. (2006). A multivariate jump-driven financial asset model. Quantitative Finance 6(5), 385–402 (https://doi.org/10.1080/14697680600806275).

[bib.bib54] Lukacs, E. (1970). Characteristic Functions. Charles Griffin & Co.

[bib.bib55] Madan, D. B., and Seneta, E. (1990). The variance Gamma (VG) model for share market returns. Journal of Business 63(4), 511–524 (https://doi.org/10.1086/296519).

[bib.bib56] Margrabe, W. (1978). The value of an option to exchange one asset for another. Journal of Finance 33(1), 177–186 (https://doi.org/10.1111/j.1540-6261.1978.tb03397.x).

[bib.bib57] MATLAB (2021). MATLAB, Version 9.11. The MathWorks Inc., Natick, MA.

[bib.bib58] Mijatović, A., and Pistorius, M. (2013). Continuously monitored barrier options under Markov processes. Mathematical Finance 23(1), 1–38 (https://doi.org/10.1111/j.1467-9965.2011.00486.x).

[bib.bib59] Oosterlee, C. W., and Grzelak, L. A. (2019). Mathematical Modeling and Computation in Finance: With Exercises and Python and MATLAB Computer Codes. World Scientific (https://doi.org/10.1142/q0236).

[bib.bib60] Piazzola, C., and Tamellini, L. (2022). The Sparse Grids Matlab kit: a Matlab implementation of sparse grids for high-dimensional function approximation and uncertainty quantification. Preprint, arXiv:2203.09314v1 [cs.MS] (https://doi.org/10.48550/arXiv.2203.09314).

[bib.bib61] Raible, S. (2000). Lévy processes in finance: theory, numerics, and empirical facts. PhD Thesis, University of Freiburg, Germany. URL: https://d-nb.info/961285192/34.

[bib.bib62] Ruijter, M. J., and Oosterlee, C. W. (2012). Two-dimensional Fourier cosine series expansion method for pricing financial options. SIAM Journal on Scientific Computing 34(5), B642–B671 (https://doi.org/10.1137/120862053).

[bib.bib63] Schoutens, W. (2003). Lévy Processes in Finance: Pricing Financial Derivatives. Wiley (https://doi.org/10.1002/0470870230).

[bib.bib64] Smolyak, S. A. (1963). Quadrature and interpolation formulas for tensor products of certain classes of functions. Doklady Akademii Nauk 148(5), 1042–1045.

[bib.bib65] Takahasi, H., and Mori, M. (1971). Estimation of errors in the numerical quadrature of analytic functions. Applicable Analysis 1(3), 201–229 (https://doi.org/10.1080/00036817108839015).

[bib.bib66] Titchmarsh, E. C. (1948). Introduction to the Theory of Fourier Integrals, 2nd edn. Clarendon Press, Oxford.

[bib.bib67] Trefethen, L. N. (2008). Is Gauss quadrature better than Clenshaw–Curtis? SIAM Review 50(1), 67–87 (https://doi.org/10.1137/060659831).

[bib.bib68] van de Wiel, D. P. (2015). Valuation of insurance products using a normal inverse Gaussian distribution. Master’s Thesis, Tilburg University, Tilburg, Netherlands. URL: https://arno.uvt.nl/show.cgi?fid=146193.

[bib.bib69] von Sydow, L., Höök, L. J., Larsson, E., Lindström, E., Milovanović, S., Persson, J., Shcherbakov, V., Shpolyanskiy, Y., Sirén, S., Toivanen, J., Waldén, J., Wiktorsson, M., Levesley, J., Li, J., Oosterlee, C. W., Ruijter, M. J., Toropov, A., and Zhao, Y. (2015). BENCHOP: the BENCHmarking project in option pricing. International Journal of Computer Mathematics 92(12), 2361–2379 (https://doi.org/10.1080/00207160.2015.1072172).

[bib.bib70] Wasilkowski, G. W., and Wozniakowski, H. (1995). Explicit cost bounds of algorithms for multivariate tensor product problems. Journal of Complexity 11(1), 1–56 (https://doi.org/10.1006/jcom.1995.1001).

[bib.bib71] Xi, Y., Ding, K., and Ning, N. (2019). Simultaneous two-dimensional continuous-time Markov chain approximation of two-dimensional fully coupled Markov diffusion processes. Working Paper, Social Science Research Network (https://doi.org/10.2139/ssrn.3461115).

[bib.bib72] Xiang, S. (2012). Asymptotics on Laguerre or Hermite polynomial expansions and their applications in Gauss quadrature. Journal of Mathematical Analysis and Applications 393(2), 434–444 (https://doi.org/10.1016/j.jmaa.2012.03.056).

[bib.bib73] Zhang, B., and Oosterlee, C. W. (2013). Efficient pricing of European-style Asian options under exponential Lévy processes based on Fourier cosine expansions. SIAM Journal on Financial Mathematics 4(1), 399–426 (https://doi.org/10.1137/110853339).

[bib.bib74] Zhu, J. (2009). Applications of Fourier Transform to Smile Modeling: Theory and Implementation. Springer (https://doi.org/10.1007/978-3-642-01808-4).

	$\displaystyle\|g(\bm{u};\bm{R},\bm{\varTheta}_{m},\bm{\varTheta}_{p})\|$	$\displaystyle\leq(2\pi)^{-d}\mathrm{e}^{-rT}\|\varPhi_{\bm{X}_{T}}(\mathrm{i}% \bm{R})\|\,\|\hat{P}(\mathrm{i}\bm{R})\|$
		$\displaystyle=\|g(\bm{0}_{\mathbb{R}^{d}};\bm{R})\|\quad\forall\bm{u}\in\mathbb{% R}^{d},~{}\bm{R}\in\delta_{V}.$		(3.9)

Journal of Computational Finance

1 Introduction

2 Problem setting and pricing framework

2.1 The multivariate Fourier pricing valuation formula

Notation and definitions.

Assumption 2.1 (Assumptions on the payoff).

Assumption 2.2.

Proposition 2.3.

Proof.

Remark 2.4.

Remark 2.5 (The case of discontinuous payoffs).

2.2 Payoffs and multivariate asset models

2.2.1 Payoffs and their Fourier transforms

2.2.2 Multivariate models and the corresponding characteristic functions

Example 2.6 (Multivariate GBM).

Example 2.7 (Multivariate VG (Luciano and Schoutens 2006)).

Example 2.8.

Remark 2.9 (About the strip of regularity).

Remark 2.10.

3 Methodology

3.1 Motivation for and characterization of the optimal damping rule

3.1.1 Motivation for the damping rule

Theorem 3.1.

Proof.

3.2 Characterization of the optimal damping rule

Proposition 3.2.

Proof.

Remark 3.3 (The case of the isotropic integrand).

Remark 3.4 (Improving the damping parameters rule).

3.3 Numerical evaluation of the inverse Fourier integrals using hierarchical deterministic quadrature methods

4 Numerical experiments and results

4.1 Combining the optimal damping heuristic rule with hierarchical deterministic quadrature methods

4.1.1 Effect of sparsification and dimension adaptivity

4.1.2 Effect of the optimal damping rule

4.2 Comparison of our approach with the COS method

4.2.1 Implementation details of the COS method

4.2.2 Numerical experiments

Remark 4.1 (About the COS method in multiple dimensions).

4.3 Computational comparison of quadrature methods with optimal damping and Monte Carlo

5 Conclusions and future work

Declaration of interest

Acknowledgements

References

You are currently on corporate access.