Numerical smoothing with hierarchical adaptive sparse grids and quasi-Monte Carlo methods for efficient option pricing

Abstract

When approximating the expectations of a functional of a solution to a stochastic differential equation, the numerical performance of deterministic quadrature methods, such as sparse grid quadrature and quasi-Monte Carlo (QMC) methods, may critically depend on the regularity of the integrand. To overcome this issue and improve the regularity structure of the problem, we consider cases in which analytic smoothing (bias-free mollification) cannot be performed and introduce a novel numerical smoothing approach by combining a root-finding method with a one-dimensional numerical integration with respect to a single well-chosen variable. We prove that, under appropriate conditions, the resulting function of the remaining variables is highly smooth, potentially affording the improved efficiency of adaptive sparse grid quadrature (ASGQ) and QMC methods, particularly when combined with hierarchical transformations (i.e. the Brownian bridge and Richardson extrapolation on the weak error). This approach facilitates the effective treatment of high dimensionality. Our study is motivated by option pricing problems, focusing on dynamics where the discretization of the asset price is necessary. Based on our analysis and numerical experiments, we demonstrate the advantages of combining numerical smoothing with the ASGQ and QMC methods over these methods without smoothing and the Monte Carlo approach. Finally, our approach is generic and can be applied to solve a broad class of problems, particularly approximating distribution functions, computing financial Greeks, and estimating risk quantities.

Keywords:

• Adaptive sparse grid quadrature
• Quasi-Monte Carlo
• Numerical smoothing
• Brownian bridge
• Richardson extrapolation
• Option pricing
• Monte Carlo
• Distribution functions
• Greeks
• Risk estimation

2010 Mathematics Subject Classification:

• 65C05
• 65D30
• 65D32
• 65Y20
• 91G20
• 91G60

Acknowledgments

The authors are incredibly grateful to the anonymous referees for their valuable comments and suggestions that greatly contributed to shaping the final version of the paper.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Notes

1 We consider discontinuities either in the gradients (kinks) or in the function (jumps).

2 The Fourier transform of the density function is available and inexpensive to compute.

3 We present the monotonicity condition for an increasing function without loss of generality. However, the assumption still holds for a decreasing function, which may be the case when considering a spread option.

4 We assume that $\{W^{\left(j\right)}{\}}_{j=1}^d$ are uncorrelated and the correlation terms are included in the diffusion terms $b_{ij}$.

5 The notation $(\cdot ,\cdot )$ denotes the scalar product operator.

6 For ease of presentation, we set the drift term in (12(12) $\begin{array}{r}\mathrm{d}{X}_t^{\left(j\right)}={\sigma }^{\left(j\right)}{X}_t^{\left(j\right)}\mathrm{d}{W}_t^{\left(j\right)},1\le j\le d,\end{array}$(12) ) to 0.

7 Without loss of generality, the correlated Brownian bridge can be obtained via simple matrix multiplication.

8 The formulation of our method is generic; for instance the mapping ${\psi }^j$ may be based on Haar basis functions as in (25(25) $\begin{array}{r}{W}_t^N:=Z_{-1}{\mathrm{\Psi }}_{-1}\left(t\right)+\sum _{n=0}^N\sum _{k=0}^{2^n-1}{Z}_{n,k}{\mathrm{\Psi }}_{n,k}\left(t\right).\end{array}$(25) ) instead of the Brownian bridges. Moreover, a different scheme for the mapping Φ may be considered instead of the Euler–Maruyama scheme used in this work.

9 Note that ${\mathbf{1}}_{1\times d}$ denotes the row vector with dimension d, where all its coordinates are 1.

10 The locations may differ depending on the considered payoff function; for instance, many payoffs in quantitative finance have kinks at the strike price.

11 Of course, the points ${\zeta }_k$ must be selected in a systematic manner depending on ${\overline{y}}_1^{\ast }$.

12 The points ${\zeta }_{\cdot ,\cdot }^{\mathrm{L}\mathrm{a}\mathrm{g}}$ and ${\zeta }_{\cdot ,\cdot }^{\mathrm{L}\mathrm{e}\mathrm{g}}$ must be selected systematically depending on $\{{\overline{y}}_i^{\ast }{\}}_{i=1}^R$.

13 The cardinality of ${\mathcal{T}}^{m\left(\text{β}\right)}$ is $\mathrm{\#}{\mathcal{T}}^{m\left(\text{β}\right)}=\prod _{n=1}^{dN-1}m\left({\text{β}}_n\right)$ with $m\left({\text{β}}_n\right)$ quadrature points along the nth dimension, and $m:\mathbb{N}\to \mathbb{N}$ is a strictly increasing function with $m\left(0\right)=0$ and $m\left(1\right)=1$.

14 We offer a note of caution regarding the convergence as $N\to \mathrm{\infty }$. Although the sequence of random processes $X_{\cdot }^N$ converges to the solution of (26(26) $\begin{array}{r}\mathrm{d}{X}_t=b\left(X_t\right)\mathrm{d}{W}_t,X_0=x\in \mathbb{R}.\end{array}$(26) ) (under the usual assumptions on b), this is not true in any sense for deterministic functions.

15 As an example, let us assume that $X_T^N(y,{\mathbf{z}}^N)=\cos \left(y\right)+z_{n,k}$. Then, (30(30) $\begin{array}{rl}\frac{\mathrm{\partial }{H}_{\delta }^N\left({\mathbf{z}}^N\right)}{\mathrm{\partial }{z}_{n,k}}& ={\int }_{\mathbb{R}}\frac{\mathrm{\partial }{g}_{\delta }\left(X_T^N\right(y,{\mathbf{z}}^N\left)\right)}{\mathrm{\partial }y}\\ & \times {\left(\frac{\mathrm{\partial }{X}_T^N}{\mathrm{\partial }y}\right(y,{\mathbf{z}}^N\left)\right)}^{-1}\frac{\mathrm{\partial }{X}_T^N}{\mathrm{\partial }{z}_{n,k}}(y,{\mathbf{z}}^N)\frac{1}{\sqrt{2\pi }}{\mathrm{e}}^{-\frac{y^2}{2}}\mathrm{d}y.\end{array}$(30) ) is generally not integrable.

16 It is probably difficult to argue that a deterministic constant C may exist in Assumption 3.2.

17 For the parts of the domain separated by the discontinuity location, the derivatives of G w.r.t. $y_1$ are bounded up to order s.

18 A similar argument holds for the part of G located on the right of the discontinuity.

19 We refer to Chen (2018) and Ernst et al. (2018) for a clear characterization of p.

20 ${\mathcal{W}}_{dN-1,\mathit{\gamma }}$ is equipped with the (unanchored) norm $\Vert f{\Vert }_{{\mathcal{W}}_{dN-1,\mathit{\gamma }}}^2=\sum _{\mathit{\alpha }\subseteq \{1:dN-1\}}\frac{1}{{\gamma }_{\alpha }}{\int }_{[0,1{]}^{|\mathit{\alpha }|}}\left({\int }_{[0,1{]}^{dN-1-|\mathit{\alpha }|}}\frac{{\mathrm{\partial }}^{|\mathit{\alpha }|}}{\mathrm{\partial }{\mathbf{y}}_{\mathit{\alpha }}}f\right(\mathbf{y})\mathrm{d}{\mathbf{y}}_{-\mathit{\alpha }}{)}^{2}d{\mathbf{y}}_{\mathit{\alpha }}$, where $\mathbf{y}:=(y_j{)}_{j\in \mathit{\alpha }}$ and ${\mathbf{y}}_{-\mathit{\alpha }}:=(y_j{)}_{j\in \{1:dN-1\}\setminus \mathit{\alpha }}$.

21 The payoff g is expressed by $g\left(\mathbf{x}\right)=max(\sum _{j=1}^d{c}_j{x}^{\left(j\right)}-K,0)$, where $\{c_j{\}}_{j=1}^d$ denote the weights of the basket.

22 The dimension of the integration problem is N for the GBM examples and $2N$ for the Heston examples.

23 When $F^N(Z_{-1},{\mathbf{Z}}^N)=\mathcal{O}\left(c\right)$ for some deterministic constant c, this property is retained when integrating out one of the rdvs, i.e. we still achieve ${\int }_{\mathbb{R}}{F}^N(y,{\mathbf{Z}}^N)\frac{1}{\sqrt{2\pi }}{\mathrm{e}}^{-\frac{y^2}{2}}\mathrm{d}y=\mathcal{O}\left(c\right).$

Additional information

Funding

C. Bayer gratefully acknowledges support from the German Research Foundation (DFG) via the Cluster of Excellence MATH+ (project AA4-2) and the individual grant BA5484/1. This publication is based on work supported by the King Abdullah University of Science and Technology (KAUST) Office of Sponsored Research (OSR) under Award No. OSR-2019-CRG8-4033 and the Alexander von Humboldt Foundation.