On pricing of discrete Asian and Lookback options under the Heston model

Abstract

We propose a new, data-driven approach for efficient pricing of – fixed- and floating-strike – discrete arithmetic Asian and Lookback options when the underlying process is driven by the Heston model dynamics. The method proposed in this article constitutes an extension of Perotti and Grzelak [Fast sampling from time-integrated bridges using deep learning, J. Comput. Math. Data Sci. 5 (2022)], where the problem of sampling from time-integrated stochastic bridges was addressed. The model relies on the Seven-League scheme [S. Liu et al. The seven-league scheme: Deep learning for large time step Monte Carlo simulations of stochastic differential equations, Risks 10 (2022), p. 47], where artificial neural networks are employed to ‘learn’ the distribution of the random variable of interest utilizing stochastic collocation points [L.A. Grzelak et al. The stochastic collocation Monte Carlo sampler: Highly efficient sampling from expensive distributions, Quant. Finance 19 (2019), pp. 339–356]. The method results in a robust procedure for Monte Carlo pricing. Furthermore, semi-analytic formulae for option pricing are provided in a simplified, yet general, framework. The model guarantees high accuracy and a reduction of the computational time up to thousands of times compared to classical Monte Carlo pricing schemes.

Keywords:

• Discrete arithmetic Asian option
• discrete Lookback option
• Heston model
• stochastic collocation (SC)
• artificial neural network (ANN)
• seven-league scheme (7L)

2020 Mathematics Subject Classifications:

• 91G60
• 91G30
• 91-08

1. Introduction

A non-trivial problem in the financial field is the pricing of path-dependent derivatives, as for instance Asian and Lookback options. The payoffs of such derivatives are expressed as functions of the underlying process monitored over the life of the option. The monitoring can either be continuous or discrete. Depending on different specifications of the underlying dynamics, only a few case-specific theoretical formulae exist. For example, under the Black-Scholes and Merton log-normal dynamics, closed-form formulae were derived for continuously-monitored geometric Asian options (see, for instance, Devreese et al. [8]). In the same model framework, Goldman et al. [10], Conze [6] derived an analytic formula for the continuously-monitored Lookback option, using probabilistic arguments such as the reflection principle. However, options whose payoffs are discretely-monitored are less tractable analytically, and so approximations are developed, as for the discrete Lookback option under the lognormal dynamics [15]. Furthermore, stochastic volatility frameworks are even more challenging for the pricing task, where no applicable closed-form theoretical solutions are known.

Whenever an exact theoretical pricing formula is not available, a rich literature on numerical methods and approximations exists. The three main classes of approaches are Monte Carlo (MC) methods (e.g. [16]), partial differential equations (PDEs) techniques (see the extensive work in [27,28]), and Fourier-inversion based techniques (among many relevant works, we report [3,31]).

Monte Carlo methods are by far the most flexible approaches, since they neither require any particular assumption on the underlying dynamics, nor on the targeted payoff. Furthermore, they benefit from a straightforward implementation based on the discretization of the time horizon as, for instance, the well-known Euler-Maruyama scheme. The cost to be paid, however, is typically a significant computational time to get accurate results. PDE approaches are more problem-specific since they require the derivation of the partial differential equation which describes the evolution of the option value over time. Then, the PDE is usually solved using finite difference methods. Fourier-inversion-based techniques exploit the relationship between the probability density function (PDF) and the characteristic function (ChF) to recover the underlying transition density by means of Fast Fourier Transform (FFT). Thanks to the swift algorithm for FFT, such methods produce high-speed numerical evaluation, but they are often problem-specific, depending on the underlying dynamics. A relevant example is [7], where a numerical method is proposed for the pricing of discrete arithmetic Asian options under stochastic volatility model dynamics. In the same group of techniques, we also refer to [19] where a unified framework for the pricing of discrete Asian options is described, allowing for regime-switching jump diffusion underlying dynamics. A further example is [21], where a close approximation of the classic Heston model is proposed – via the CTMC model – which allows for a general SWIFT-based (see [20]) pricing approach with application to exotic derivatives, as discrete Asian options.

In this article, we propose a data-driven¹ extension to MC schemes that allows for efficient pricing of discretely-monitored Asian and Lookback options, without losing the flexibility typical of MC methods. We develop the methodology in the complex framework of the stochastic volatility model of Heston [14], with an extensive application to the case of Feller condition not satisfied. Under this dynamics, we show how to price fixed- or floating-strike discrete Asian and Lookback options. Moreover, the pricing model is applied also for the challenging task of pricing options with both a fixed- and a floating-strike component. We underline that the strengths of the method are its speed and accuracy, coupled with significant flexibility. The procedure is, indeed, independent of the underlying dynamics (it could be applied, for instance, for any stochastic volatility model), and it is not sensitive to the targeted payoff.

Inspired by the works in [22,25], the method relies on the technique of Stochastic Collocation (SC) [12], which is employed to accurately approximate the targeted distribution by means of piecewise polynomials. Artificial neural networks (ANNs) are used for fast recovery of the coefficients which uniquely determine the piecewise approximation. Given these coefficients, the pricing can be performed in a ‘MC fashion’ sampling from the target (approximated) distribution and computing the numerical average of the discounted payoffs. Furthermore, in a simplified setting, we provide a semi-analytic formula that allows to directly price options without the need of sampling from the desired distribution. In both situations (MC and semi-analytic pricing) we report a significant computational speed-up, without affecting the accuracy of the result which remains comparable with the one of highly expensive MC methods.

The remainder of the paper is as follows. In Section 2, we formally define discrete arithmetic Asian and Lookback options, as well as the model framework for the underlying process. Then, in Section 3, the pricing model is described. Two different cases are considered, in increasing order of complexity, to handle efficiently both unconditional sampling Section 3.2 and conditional sampling Section 3.3 for pricing of discrete arithmetic Asian and Lookback options. Section 4 provides theoretical support to the given numerical scheme. The quality of the methodology is also inspected empirically with several numerical experiments, reported in Section 5. Section 6 concludes.

2. Discrete arithmetic Asian and Lookback options

In a generic setting, given the present time $t_0\ge 0$, the payoff at time $T>t_0$ of a discrete arithmetic Asian or Lookback option, with underlying process $S\left(t\right)$, can be written as: (1) $H_{\omega }(T;S)=max(\omega \left(A\right(S)-K_{1}S(T)-K_2),0),$(1) where $A\left(S\right)\equiv A\left(S\right(t);t_0<t\le T)$ is a deterministic function of the underlying process S, the constants $K_1,K_2\ge 0$ control the floating- and fixed-strikes of the option, and $\omega =\pm 1$. Particularly, discrete arithmetic Asian and Lookback options are obtained by setting the quantity $A\left(S\right)$ respectively as follows: (2) $A\left(S\right):=\frac{1}{N}\sum _{n\in I}S\left(t_n\right),A\left(S\right):=\omega \underset{n\in I}{max}\mathrm{\omega S}\left(t_n\right),$(2) with $I=\{1,\dots ,N\}$ a set of indexes, $t_1<t_2<\cdots <t_N=T$, a discrete set of future monitoring dates, and ω as in (1(1) $H_{\omega }(T;S)=max(\omega \left(A\right(S)-K_{1}S(T)-K_2),0),$(1) ).

Note that in both the cases of discrete arithmetic Asian and Lookback options, $A\left(S\right)$ is expressed as a deterministic transformation of the underlying process' path, which is the only requirement to apply the proposed method. Therefore, in the paper, we refer always to the class of discrete arithmetic Asian options, often just called Asian options, for simplicity. However, the theory holds for both classes of products. Actually, the pricing model applies to any product with a path-dependent European-type payoff , requiring only a different definition of $A\left(S\right)$.

2.1. Pricing of arithmetic Asian options and Heston framework

This section focuses on the risk-neutral pricing of arithmetic – fixed- and floating-strike – Asian options, whose payoff is given in (1(1) $H_{\omega }(T;S)=max(\omega \left(A\right(S)-K_{1}S(T)-K_2),0),$(1) ) with $A\left(S\right)$ in (2(2) $A\left(S\right):=\frac{1}{N}\sum _{n\in I}S\left(t_n\right),A\left(S\right):=\omega \underset{n\in I}{max}\mathrm{\omega S}\left(t_n\right),$(2) ). By setting $K_1=0$ or $K_2=0$, we get two special cases: the fixed- or the floating-strike arithmetic Asian option. From Equation (1(1) $H_{\omega }(T;S)=max(\omega \left(A\right(S)-K_{1}S(T)-K_2),0),$(1) ), for $K_1=0$, the simplified payoff of a fixed-strike arithmetic Asian option reads: (3) $H_{\omega }^{\mathtt{fx}}(T;S)=max(\omega \left(A\right(S)-K_2),0),$(3) therefore, with the risk-neutral present value: (4) $V_{\omega }^{\mathtt{fx}}\left(t_0\right)=\frac{M\left(t_0\right)}{M\left(T\right)}{\mathbb{E}}_{t_0}^{\mathbb{Q}}[max(\omega \left(A\right(S)-K_2),0)],$(4) where we assume the money-savings account $M\left(t\right)$ to be defined through the deterministic dynamics $\mathrm{d}M\left(t\right)=\mathrm{rM}\left(t\right)\mathrm{d}t$ with constant interest rate $r\ge 0$. For $K_2=0$, however, the payoff in Equation (1(1) $H_{\omega }(T;S)=max(\omega \left(A\right(S)-K_{1}S(T)-K_2),0),$(1) ), becomes the one of a floating-strike arithmetic Asian option: (5) $H_{\omega }^{\mathtt{fl}}(T;S)=max(\omega \left(A\right(S)-K_{1}S(T\left)\right),0).$(5) The payoff in (5(5) $H_{\omega }^{\mathtt{fl}}(T;S)=max(\omega \left(A\right(S)-K_{1}S(T\left)\right),0).$(5) ) is less tractable when compared to the one in (3(3) $H_{\omega }^{\mathtt{fx}}(T;S)=max(\omega \left(A\right(S)-K_2),0),$(3) ), because of the presence of two dependent stochastic unknowns, namely $S\left(T\right)$ and $A\left(S\right)$. However, a similar representation to the one in Equation (4(4) $V_{\omega }^{\mathtt{fx}}\left(t_0\right)=\frac{M\left(t_0\right)}{M\left(T\right)}{\mathbb{E}}_{t_0}^{\mathbb{Q}}[max(\omega \left(A\right(S)-K_2),0)],$(4) ) can be achieved, allowing for a unique pricing approach in both cases. By a change of measure from the risk-neutral measure $\mathbb{Q}$ to the measure ${\mathbb{Q}}^S$ associated with the numéraire $S\left(t\right)$, i.e. the stock measure, we prove the following proposition.

Proposition 2.1

Pricing of floating-strike arithmetic Asian option under the stock measure

Under the stock measure ${\mathbb{Q}}^S$, the value at time $t_0\ge 0$ of an arithmetic floating-strike Asian option, with maturity $T>t_0$ and future monitoring dates $t_n$, $n\in \{1,\dots ,N\}$, reads: (6) $V_{\omega }^{\mathtt{fl}}\left(t_0\right)=S\left(t_0\right){\mathbb{E}}_{t_0}^S[max(\omega \left(A^{\mathtt{fl}}\right(S)-K_1),0)],$(6) with $A^{\mathtt{fl}}\left(S\right)$, defined as: $A^{\mathtt{fl}}\left(S\right):=A\left(\frac{S(\cdot )}{S\left(T\right)}\right)=\frac{A\left(S\right)}{S\left(T\right)},$where $A\left(S\right)$ is defined in (2(2) $A\left(S\right):=\frac{1}{N}\sum _{n\in I}S\left(t_n\right),A\left(S\right):=\omega \underset{n\in I}{max}\mathrm{\omega S}\left(t_n\right),$(2) ).

Proof.

For a proof, see Appendix A.1.

Both the representations in Equations (4(4) $V_{\omega }^{\mathtt{fx}}\left(t_0\right)=\frac{M\left(t_0\right)}{M\left(T\right)}{\mathbb{E}}_{t_0}^{\mathbb{Q}}[max(\omega \left(A\right(S)-K_2),0)],$(4) ) and (6(6) $V_{\omega }^{\mathtt{fl}}\left(t_0\right)=S\left(t_0\right){\mathbb{E}}_{t_0}^S[max(\omega \left(A^{\mathtt{fl}}\right(S)-K_1),0)],$(6) ) can be treated in the same way. This means that the present value of both a fixed- and a floating-strike Asian option can be computed similarly, as stated in the following proposition (see, e.g. [13]).

Proposition 2.2

Symmetry of fixed- and floating-strike Asian option present value

Let us consider the process $S\left(t\right)$ and the money-savings account $M\left(t\right)$, for $t\ge t_0$. Then, the same representation holds for the value at time $t_0$ of both fixed- or floating-strike Asian options, with maturity $T>t_0$, underlying process $S\left(t\right)$, and future monitoring dates $t_n$, $n\in \{1,\dots ,N\}$. The present value is given by: (7) $V_{\omega }^{\lambda }\left(t_0\right)=\{\begin{array}{l}{\displaystyle \frac{M\left(t_0\right)}{M\left(T\right)}{\mathbb{E}}_{t_0}^{\mathbb{Q}}[max(\omega \left(A\right(S)-K_2),0)],\mathrm{for}\lambda =\mathtt{fx},}\\ S\left(t_0\right){\mathbb{E}}_{t_0}^S[max(\omega \left(A^{\mathtt{fl}}\right(S)-K_1),0)],\mathrm{for}\lambda =\mathtt{fl}.\end{array}$(7)

Proof.

The proof follows by direct comparison between Equations (4(4) $V_{\omega }^{\mathtt{fx}}\left(t_0\right)=\frac{M\left(t_0\right)}{M\left(T\right)}{\mathbb{E}}_{t_0}^{\mathbb{Q}}[max(\omega \left(A\right(S)-K_2),0)],$(4) ) and (6(6) $V_{\omega }^{\mathtt{fl}}\left(t_0\right)=S\left(t_0\right){\mathbb{E}}_{t_0}^S[max(\omega \left(A^{\mathtt{fl}}\right(S)-K_1),0)],$(6) ).

When $K_1\ne 0$ and $K_2\ne 0$, Equation (1(1) $H_{\omega }(T;S)=max(\omega \left(A\right(S)-K_{1}S(T)-K_2),0),$(1) ) is the payoff of a fixed- and floating-strikes arithmetic Asian option. Its present value does not allow any simplified representation, and we write it as the expectation of the discounted payoff under the risk-neutral measure $\mathbb{Q}$: (8) $V_{\omega }\left(t_0\right)=\frac{M\left(t_0\right)}{M\left(T\right)}{\mathbb{E}}_{t_0}^{\mathbb{Q}}[max(\omega \left(A\right(S)-K_{1}S(T)-K_2),0)].$(8) By comparing Equations (7(7) $V_{\omega }^{\lambda }\left(t_0\right)=\{\begin{array}{l}{\displaystyle \frac{M\left(t_0\right)}{M\left(T\right)}{\mathbb{E}}_{t_0}^{\mathbb{Q}}[max(\omega \left(A\right(S)-K_2),0)],\mathrm{for}\lambda =\mathtt{fx},}\\ S\left(t_0\right){\mathbb{E}}_{t_0}^S[max(\omega \left(A^{\mathtt{fl}}\right(S)-K_1),0)],\mathrm{for}\lambda =\mathtt{fl}.\end{array}$(7) ) and (8(8) $V_{\omega }\left(t_0\right)=\frac{M\left(t_0\right)}{M\left(T\right)}{\mathbb{E}}_{t_0}^{\mathbb{Q}}[max(\omega \left(A\right(S)-K_{1}S(T)-K_2),0)].$(8) ) a difference in the two settings is unveiled. Equation (7(7) $V_{\omega }^{\lambda }\left(t_0\right)=\{\begin{array}{l}{\displaystyle \frac{M\left(t_0\right)}{M\left(T\right)}{\mathbb{E}}_{t_0}^{\mathbb{Q}}[max(\omega \left(A\right(S)-K_2),0)],\mathrm{for}\lambda =\mathtt{fx},}\\ S\left(t_0\right){\mathbb{E}}_{t_0}^S[max(\omega \left(A^{\mathtt{fl}}\right(S)-K_1),0)],\mathrm{for}\lambda =\mathtt{fl}.\end{array}$(7) ) is characterized by a unique unknown stochastic quantity $A\left(S\right)$, whereas in (8(8) $V_{\omega }\left(t_0\right)=\frac{M\left(t_0\right)}{M\left(T\right)}{\mathbb{E}}_{t_0}^{\mathbb{Q}}[max(\omega \left(A\right(S)-K_{1}S(T)-K_2),0)].$(8) ) an additional term appears, namely the stock price at final time, $S\left(T\right)$. Furthermore, the two stochastic quantities in (8(8) $V_{\omega }\left(t_0\right)=\frac{M\left(t_0\right)}{M\left(T\right)}{\mathbb{E}}_{t_0}^{\mathbb{Q}}[max(\omega \left(A\right(S)-K_{1}S(T)-K_2),0)].$(8) ) are not independent. This might suggest that different procedures should be employed for the different payoffs. In particular, in a MC setting, to value (7(7) $V_{\omega }^{\lambda }\left(t_0\right)=\{\begin{array}{l}{\displaystyle \frac{M\left(t_0\right)}{M\left(T\right)}{\mathbb{E}}_{t_0}^{\mathbb{Q}}[max(\omega \left(A\right(S)-K_2),0)],\mathrm{for}\lambda =\mathtt{fx},}\\ S\left(t_0\right){\mathbb{E}}_{t_0}^S[max(\omega \left(A^{\mathtt{fl}}\right(S)-K_1),0)],\mathrm{for}\lambda =\mathtt{fl}.\end{array}$(7) ) we only have to sample from the unconditional distribution of $A\left(S\right)$; while in (8(8) $V_{\omega }\left(t_0\right)=\frac{M\left(t_0\right)}{M\left(T\right)}{\mathbb{E}}_{t_0}^{\mathbb{Q}}[max(\omega \left(A\right(S)-K_{1}S(T)-K_2),0)].$(8) ) the MC scheme requires dealing with both the sampling of $S\left(T\right)$ and the conditional sampling of $A\left(S\right)|S\left(T\right)$.

Let us define the stochastic volatility dynamics of Heston for the underlying stochastic process $S\left(t\right)$, with initial value $S\left(t_0\right)=S_0$, through the following system of stochastic differential equations (SDEs): (9) $\begin{array}{rl}\mathrm{d}S\left(t\right)& =rS\left(t\right)\mathrm{d}t+\sqrt{v\left(t\right)}S\left(t\right)\mathrm{d}{W}_x\left(t\right),S\left(t_0\right)=S_0,\end{array}$(9) (10) $\begin{array}{rl}\mathrm{d}v\left(t\right)& =\kappa (\overline{v}-v(t\left)\right)\mathrm{d}t+\gamma \sqrt{v\left(t\right)}\mathrm{d}{W}_v\left(t\right),v\left(t_0\right)=v_0,\end{array}$(10) with $r,\kappa ,\overline{v},v_0\ge 0$, $\gamma >0$ the constant rate, the speed of mean reversion, the long-term mean of the variance process, the initial variance, and the volatility-of-volatility, respectively. $W_x\left(t\right)$ and $W_v\left(t\right)$ are Brownian Motions (BMs) under the risk-neutral measure $\mathbb{Q}$ with correlation coefficient $\rho \in [-1,1]$, i.e. $\mathrm{d}{W}_x\left(t\right)\mathrm{d}{W}_v\left(t\right)=\rho \mathrm{d}t$.² The dynamics in (9(9) $\begin{array}{rl}\mathrm{d}S\left(t\right)& =rS\left(t\right)\mathrm{d}t+\sqrt{v\left(t\right)}S\left(t\right)\mathrm{d}{W}_x\left(t\right),S\left(t_0\right)=S_0,\end{array}$(9) ) and (10(10) $\begin{array}{rl}\mathrm{d}v\left(t\right)& =\kappa (\overline{v}-v(t\left)\right)\mathrm{d}t+\gamma \sqrt{v\left(t\right)}\mathrm{d}{W}_v\left(t\right),v\left(t_0\right)=v_0,\end{array}$(10) ) are defined in the risk-neutral framework. However, Proposition 2.1 entails a different measure framework, whose dynamics still fall in the class of stochastic volatility model of Heston, with adjusted parameters (see Proposition Appendix A.1).

3. Swift numerical pricing using deep learning

This section focuses on the efficient pricing of discrete arithmetic Asian options in a MC setting. The method uses a Stochastic Collocation [12] (SC) based approach to approximate the target distribution. Then, artificial neural networks (ANNs) ‘learn’ the proxy of the desired distribution, allowing for fast recovery [22,25].

3.1. ‘Compressing’ distribution with stochastic collocation

In the framework of MC methods, the idea of SC – based on the probability integral transform³ – is to approximate the relationship between a ‘computationally expensive’ random variable, say A, and a ‘computationally cheap’ one, say ξ. The approximation is then used for sampling. A random variable is ‘expensive’ if its inverse CDF is not known in analytic form, and needs to be computed numerically. With SC, the sampling of A is performed at the cost of sampling ξ (see [12]). Formally, the following mapping is used to generate samples from A: (11) $A\stackrel{\mathrm{d}}{=}{F}_A^{-1}\left(F_{\xi }\right(\xi \left)\right)=:g\left(\xi \right)\approx \tilde{g}\left(\xi \right),$(11) with $F_A$ and $F_{\xi }$ being respectively the CDFs of A and ξ, and the function $\tilde{g}$ a suitable, easily evaluable approximation of g. The reason why we prefer $\tilde{g}$ to g, is that, by definition, every evaluation of g requires the numerical inversion of $F_A$, the CDF of A.

Many possible choices of $\tilde{g}$ exist. In [12,25], $\tilde{g}$ is an $(M-1)$-degree polynomial expressed in Lagrange basis, defined on collocation points (CPs) $\mathit{\xi }:=\{{\xi }_k{\}}_{k=1}^M$ computed as Gauss-Hermite quadrature nodes, i.e.: (12) $\tilde{g}\left(x\right):=\sum _{k=1}^M{a}_k{\ell }_k\left(x\right),{\ell }_k\left(x\right):=\prod _{\begin{array}{c}1\le j\le M\\ j\ne k\end{array}}\frac{x-{\xi }_j}{{\xi }_k-{\xi }_j},k=1,\dots ,M,$(12) where the coefficients $\mathbf{a}:=\{a_k{\}}_{k=1}^M$ of the polynomial in the Lagrange basis representation, called collocation values (CVs), are derived by imposing the system of equations: (13) $g\left({\xi }_k\right)=\tilde{g}\left({\xi }_k\right)=:a_k,k=1,\dots ,M,$(13) which requires only M evaluations of g.

In this work, we define $\tilde{g}$ as a piecewise polynomial. Particularly, we divide the domain $\mathbb{R}$ of the random variable ξ (which for us is standard normally distributed⁴) in three regions. In each region, we define $\tilde{g}$ as a polynomial. In other words, given the partition of the real line: (14) ${\mathrm{\Omega }}_{-}\cup {\mathrm{\Omega }}_M\cup {\mathrm{\Omega }}_{+}:=(-\mathrm{\infty },-\overline{\xi })\cup [-\overline{\xi },\overline{\xi }]\cup (\overline{\xi },+\mathrm{\infty }),$(14) for a certain $\overline{\xi }>0$, $\tilde{g}$ is specified as: (15) $\tilde{g}\left(\xi \right):=g_{-}\left(\xi \right)\cdot {\mathbb{1}}_{{\mathrm{\Omega }}_{-}}\left(\xi \right)+g_M\left(\xi \right)\cdot {\mathbb{1}}_{{\mathrm{\Omega }}_M}\left(\xi \right)+g_{+}\left(\xi \right)\cdot {\mathbb{1}}_{{\mathrm{\Omega }}_{+}}\left(\xi \right),$(15) where ${\mathbb{1}}_{(\cdot )}$ is the indicator function, and $g_{-},g_M,g_{+}$ are suitable polynomials. To ensure high accuracy in the approximation, $g_M$ is defined as a Lagrange polynomial of high-degree M−1. The CPs $\mathit{\xi }$, which identify the Lagrange basis in (12(12) $\tilde{g}\left(x\right):=\sum _{k=1}^M{a}_k{\ell }_k\left(x\right),{\ell }_k\left(x\right):=\prod _{\begin{array}{c}1\le j\le M\\ j\ne k\end{array}}\frac{x-{\xi }_j}{{\xi }_k-{\xi }_j},k=1,\dots ,M,$(12) ), are chosen as Chebyshev nodes in the bounded interval ${\mathrm{\Omega }}_M=[-\overline{\xi },\overline{\xi }]$ [9]. Instead, the CVs $\mathbf{a}$ are defined as in (13(13) $g\left({\xi }_k\right)=\tilde{g}\left({\xi }_k\right)=:a_k,k=1,\dots ,M,$(13) ). The choice of Chebyshev nodes allows for increasing the degree of the interpolation (i.e. the number of CPs and CVs), avoiding the Runge's phenomenon within the interval ${\mathrm{\Omega }}_M$ [26]. However, the behaviour of $g_M$ outside ${\mathrm{\Omega }}_M$ is out of control. We expect the high-degree polynomial $g_M$ to be a poor approximation of g in ${\mathrm{\Omega }}_{-}$ and ${\mathrm{\Omega }}_{+}$. Therefore, we define $g_{-}$ and $g_{+}$ as linear (or at most quadratic) polynomials, with degree $M_{-}-1$ and $M_{+}-1$, built on the extreme CPs of $\mathit{\xi }$.

Summarizing, $g_{-}$, $g_M$ and $g_{+}$ are all defined as Lagrange polynomials: (16) $g_{(\cdot )}\left(x\right):=\sum _{k\in I_{(\cdot )}}{a}_k{\ell }_k\left(x\right),{\ell }_k\left(x\right):=\prod _{\begin{array}{c}j\in I_{(\cdot )}\\ j\ne k\end{array}}\frac{x-{\xi }_j}{{\xi }_k-{\xi }_j},k\in I_{(\cdot )},$(16) where the sets of indexes for $g_{-}$, $g_M$ and $g_{+}$ are $I_{-}=\{1,2\}$, $I_M=\{1,\dots ,M\}$ and $I_{+}=\{M-1,M\}$, respectively (if a quadratic extrapolation is preferred, we get $I_{-}=\{1,2,3\}$, and $I_{+}=\{M-2,M-1,M\}$).

Remark 3.1

‘Compressed’ distributions

The SC technique is a tool to ‘compress’ the information regarding A in (11(11) $A\stackrel{\mathrm{d}}{=}{F}_A^{-1}\left(F_{\xi }\right(\xi \left)\right)=:g\left(\xi \right)\approx \tilde{g}\left(\xi \right),$(11) ), into a small number of coefficients, the CVs $\mathbf{a}$. Indeed, the relationship between A and $\mathbf{a}$ is bijective, provided the distribution of the random variable ξ, and the corresponding CPs $\mathit{\xi }$ (or, equivalently, the Lagrange basis in (12(12) $\tilde{g}\left(x\right):=\sum _{k=1}^M{a}_k{\ell }_k\left(x\right),{\ell }_k\left(x\right):=\prod _{\begin{array}{c}1\le j\le M\\ j\ne k\end{array}}\frac{x-{\xi }_j}{{\xi }_k-{\xi }_j},k=1,\dots ,M,$(12) )), are specified a priori.

3.2. Semi-analytical pricing of fixed- or floating-strike Asian options

Let us first consider the pricing of the fixed- or the floating-strike Asian options. Both the products allow for the same representation in which the only unknown stochastic quantity is $A^{\lambda }\left(S\right)$, $\lambda \in \{\mathtt{\text{fx}},\mathtt{\text{fl}}\}$, as given in Proposition 2.2. For the sake of simplicity, in the absence of ambiguity, we call $A^{\lambda }\left(S\right)$ just A.

For pricing purposes, we can benefit of the SC technique presented in the previous section, provided we know the map $\tilde{g}$ (or, equivalently, the CVs $\mathbf{a}$), i.e.: (17) $\begin{array}{rl}{V}_{\omega }\left(t_0\right)& =C\mathbb{E}[max(\omega (A-K),0)]\\ & \approx C\mathbb{E}[max(\omega \left(\tilde{g}\right(\xi )-K),0)]=:{\tilde{V}}_{\omega }\left(t_0\right),\end{array}$(17) where ξ is a standard normally distributed random variable, and C is a constant coherent with Proposition 2.2. We note, that ${\tilde{V}}_{\omega }\left(t_0\right)$ is the expectation of (the positive part of) polynomials of a standard normal distribution. Hence, a semi-analytic formula exists, in a similar fashion as the one given in [11].

Proposition 3.2

Semi-analytic pricing formula

Let ${\tilde{V}}_{\omega }\left(t_0\right)$ (and C) be defined as in Equation (17(17) $\begin{array}{rl}{V}_{\omega }\left(t_0\right)& =C\mathbb{E}[max(\omega (A-K),0)]\\ & \approx C\mathbb{E}[max(\omega \left(\tilde{g}\right(\xi )-K),0)]=:{\tilde{V}}_{\omega }\left(t_0\right),\end{array}$(17) ), with $\tilde{g}$ defined in (15(15) $\tilde{g}\left(\xi \right):=g_{-}\left(\xi \right)\cdot {\mathbb{1}}_{{\mathrm{\Omega }}_{-}}\left(\xi \right)+g_M\left(\xi \right)\cdot {\mathbb{1}}_{{\mathrm{\Omega }}_M}\left(\xi \right)+g_{+}\left(\xi \right)\cdot {\mathbb{1}}_{{\mathrm{\Omega }}_{+}}\left(\xi \right),$(15) ). Assume further that ${\mathit{\alpha }}_{-}$, ${\mathit{\alpha }}_M$ and ${\mathit{\alpha }}_{+}$ are the coefficients in the canonical basis of monomials for the three polynomials $g_{-}$, $g_M$ and $g_{+}$ respectively of degree $M_{-}-1$, M−1 and $M_{+}-1$. Then, using the notation $a\vee b=max(a,b)$, the following semi-analytic pricing approximation holds: $\begin{array}{rl}\frac{{\tilde{V}}_{\omega }\left(t_0\right)}{\mathrm{\omega C}}& =\left[\sum _{i=0}^{M_{-\omega }-1}{\omega }^i{\alpha }_{-\omega ,i}{m}_i\right(\omega c_K,-\overline{\xi }\vee \omega c_K\left)\right]\cdot \left(F_{\xi }\right(-\overline{\xi }\vee \omega c_K)-F_{\xi }(\omega c_K\left)\right)\\ & +\left[\sum _{i=0}^{M-1}{\omega }^i{\alpha }_{M,i}{m}_i\right(-\overline{\xi }\vee \omega c_K,\overline{\xi }\vee \omega c_K\left)\right]\cdot \left(F_{\xi }\right(\overline{\xi }\vee \omega c_K)-F_{\xi }(-\overline{\xi }\vee \omega c_K\left)\right)\\ & +\left[\sum _{i=0}^{M_{\omega }-1}{\omega }^i{\alpha }_{\omega ,i}{m}_i\right(\overline{\xi }\vee \omega c_K,+\mathrm{\infty }\left)\right]\cdot (1-F_{\xi }(\overline{\xi }\vee \omega c_K\left)\right)-K\cdot (1-F_{\xi }(\omega c_K\left)\right),\end{array}$where $F_{\xi }$ is the CDF of a standard normal random variable, ξ, $m_i(a,b):=\mathbb{E}[{\xi }^i|a\le \xi \le b]$⁵, $c_K$ satisfies $K=\tilde{g}\left(c_K\right)$, and $\omega =\pm 1$ according to the call/put case.

Proof.

For proof of the previous proposition, see Appendix A.2.

We note also that Proposition 3.2 uses as input in the pricing formula the coefficients in the canonical basis, not the ones in the Lagrange basis, $\mathbf{a}$, in (13(13) $g\left({\xi }_k\right)=\tilde{g}\left({\xi }_k\right)=:a_k,k=1,\dots ,M,$(13) ).

Remark 3.3

Change of basis

Given a Lagrange basis identified by M collocation points $\mathit{\xi }$, any $(M-1)$-degree polynomial $g\left(\xi \right)$ is uniquely determined by the corresponding M coefficients $\mathbf{a}$. A linear transformation exists that connects the M coefficients in the Lagrange basis with the M coefficients $\mathit{\alpha }$ in the canonical basis of monomials. In particular, it holds: $\mathbf{M}\mathit{\alpha }=\mathit{a},$with $\mathbf{M}\equiv \mathbf{M}\left(\mathit{\xi }\right)$ a $M\times M$ Vandermonde matrix with element ${\mathbf{M}}_{k,i}:={\xi }_k^{i-1}$ in position $(k,i)$. The matrix $\mathbf{M}$ admits an inverse; thus, the coefficients $\mathit{\alpha }$ in the canonical basis are the result of matrix-vector multiplication, provided the coefficients $\mathbf{a}$ in the Lagrange basis are known. Moreover, since the matrix $\mathbf{M}$ only depends on $\mathit{\xi }$, its inverse can be computed a priori once the CPs $\mathit{\xi }$ are fixed.

Proposition 3.2 provides a semi-analytic formula for the pricing of fixed- or floating-strike Asian options. Indeed, it requires the inversion of the map $\tilde{g}$ which typically is not available in analytic form. On the other hand, since both the CPs $\mathit{\xi }$ and the CVs $\mathbf{a}$ are known, a proxy of ${\tilde{g}}^{-1}$ is easily achievable by interpolation on the pairs of values $(a_k,{\xi }_k)$, $k=1,\dots ,M$.

The last problem is to recover the CVs $\mathbf{a}$ (which identify $\tilde{g}$) in an accurate and fast way. We recall that, for $k=1,\dots ,M$, each CV $a_k$ is defined in terms of the exact map g and the CP ${\xi }_k$ by the relationship: $a_k:=g\left({\xi }_k\right)=F_A^{-1}\left(F_{\xi }\right({\xi }_k\left)\right).$The presence of $F_A^{-1}$ makes it impossible to directly compute $\mathbf{a}$ efficiently. On the other hand, by definition, the CVs $\mathbf{a}$ are quantiles of the random variable $A\equiv A\left(S\right)$, which depends on the parameters $\mathbf{p}$ of the underlying process S. As a consequence, there must exist some unknown mapping H which links $\mathbf{p}$ to the corresponding $\mathbf{a}$. We approximate such a mapping from synthetic data setting a regression problem, which is solved with an ANN $\tilde{H}$ (in the same fashion as in [22,25]). We have the following mapping: $\mathbf{p}\mapsto \mathbf{a}:=H\left(\mathbf{p}\right)\approx \tilde{H}\left(\mathbf{p}\right),\mathbf{p}\in {\mathrm{\Omega }}_{\mathbf{p}},\mathbf{a}\in {\mathrm{\Omega }}_{\mathbf{a}},$with ${\mathrm{\Omega }}_{\mathbf{p}}$ and ${\mathrm{\Omega }}_{\mathbf{a}}$ being the spaces of the underlying model parameters and of the CVs, respectively, while the ANN $\tilde{H}$ is the result of an optimization process on a synthetic training set⁶: (18) $\mathbf{T}=\left\{\right({\mathbf{p}}_i,{\mathbf{a}}_i):i\in \{1,\dots ,N_{\mathtt{pairs}}\left\}\right\}.$(18) The pricing procedure is summarized in the following algorithm.

Algorithm: Semi-analytic pricing

1. Fix the M collocation points $\mathit{\xi }$.

2. Given the parameters $\mathbf{p}$, approximate the M collocation values, i.e. $\mathbf{a}\approx \tilde{H}\left(\mathbf{p}\right)$.

3. Given $\mathbf{a}$, compute the coefficients ${\mathit{\alpha }}_{-}$, ${\mathit{\alpha }}_M$ and ${\mathit{\alpha }}_{+}$ for $g_{-}$, $g_M$ and $g_{+}$ (see Remark 3.3).

4. Given K and $\mathbf{a}$, compute $c_K$ of Proposition 3.2 interpolating ${\tilde{g}}^{-1}$ on $(\mathbf{a},\mathit{\xi })$.

5. Given the coefficients ${\mathit{\alpha }}_{-}$, ${\mathit{\alpha }}_M$, ${\mathit{\alpha }}_{+}$, and $c_K$, use Proposition 3.2 to compute ${\tilde{V}}_{\omega }\left(t_0\right)$.

3.3. Swift Monte Carlo pricing of fixed- and floating-strikes Asian options

Let us consider the case of an option whose payoff has both a fixed- and a floating-strike. The present value of such a derivative is given by: $V_{\omega }\left(t_0\right)=\frac{M\left(t_0\right)}{M\left(T\right)}{\mathbb{E}}_{t_0}^{\mathbb{Q}}[max(\omega \left(A\right(S)-K_{1}S(T)-K_2),0)],$hence the price of the option is a function of the two dependent quantities $A\left(S\right)$ and $S\left(T\right)$. This means that, even in a MC setting, the dependency between $A\left(S\right)$ and $S\left(T\right)$ has to be fulfilled. Therefore, a different methodology with respect to the one proposed in the previous section needs to be developed.

Due to the availability of efficient and accurate sampling techniques for the underlying process S at a given future time T (we use the COS method [24] enhanced with SC [12], and we call it COS-SC), the main issue is the sampling of the conditional random variable $A\left(S\right)|S\left(T\right)$. This task is addressed in the same fashion as it is done in [25], where ANNs and stochastic collocation are applied for the efficient sampling from time-integral of stochastic bridges⁷, namely ${\int }_{t_0}^{T}S\left(t\right)\mathrm{d}t$ given the value of $S\left(T\right)$. The underlying idea here is the same since the random variable $A\left(S\right)$ is conditional to $S\left(T\right)$. Especially, in the previous sections we pointed out that the distribution of $A\left(S\right)$ has an unknown parametric form which depends on the set of Heston parameters $\mathbf{p}$. Similarly, we expect the distribution of $A|S\left(T\right)=\hat{S}$ to be parametric into the ‘augmented’ set of parameters ${\mathbf{p}}_{\hat{S}}:=\mathbf{p}\cup \left\{\hat{S}\right\}$. Hence, there exists a mapping H which links ${\mathbf{p}}_{\hat{S}}$ with the CVs, ${\mathbf{a}}_{\hat{S}}$, corresponding to the conditional distribution $A\left(S\right)|S\left(T\right)=\hat{S}$. We approximate H by means of a suitable ANN $\text{∼}H$, getting the following mapping scheme: ${\mathbf{p}}_{\hat{S}}\mapsto {\mathbf{a}}_{\hat{S}}:=H\left({\mathbf{p}}_{\hat{S}}\right)\approx \tilde{H}\left({\mathbf{p}}_{\hat{S}}\right),{\mathbf{p}}_{\hat{S}}\in {\mathrm{\Omega }}_{{\mathbf{p}}_S},{\mathbf{a}}_{\hat{S}}\in {\mathrm{\Omega }}_{{\mathbf{a}}_S},$where ${\mathrm{\Omega }}_{{\mathbf{p}}_S}$ and ${\mathrm{\Omega }}_{{\mathbf{a}}_S}$ are respectively the spaces of the underlying model parameters (augmented with $\hat{S}$) and of the CVs (corresponding to the conditional distribution $A\left(S\right)|S\left(T\right)$), and $\tilde{H}$ is the result of a regression problem on a suitable training set $\mathbf{T}$ (see Equation (18(18) $\mathbf{T}=\left\{\right({\mathbf{p}}_i,{\mathbf{a}}_i):i\in \{1,\dots ,N_{\mathtt{pairs}}\left\}\right\}.$(18) )). We propose a first brute force sampling scheme.

Algorithm: Brute force conditional sampling and pricing

1. Fix the M collocation points $\mathit{\xi }$.

2. Given the parameters $\mathbf{p}$, for $j=1,\dots ,N_{\mathtt{\text{paths}}}$, repeat:

   1. generate the sample ${\hat{S}}_j$ from $S\left(T\right)$ (e.g. with COS-SC method [12,24]);

   2. given ${\mathbf{p}}_{{\hat{S}}_j}$, approximate the M conditional CVs, i.e. ${\mathbf{a}}_{{\hat{S}}_j}\approx \tilde{H}\left({\mathbf{p}}_{{\hat{S}}_j}\right)$;

   3. given ${\mathbf{a}}_{{\hat{S}}_j}$, use SC to generate the conditional sample ${\hat{A}}_j$.

3. Given the pairs $({\hat{S}}_j,{\hat{A}}_j)$, for $j=1,\dots ,N_{\mathtt{\text{paths}}}$, and any desired $(K_1,K_2)$, evaluate: $V_{\omega }\left(t_0\right)\approx \frac{1}{N_{\mathtt{paths}}}\frac{M\left(t_0\right)}{M\left(T\right)}\sum _{j=1}^{N_{\mathtt{paths}}}max(\omega ({\hat{A}}_j-K_1{\hat{S}}_j-K_2),0).$

Nonetheless, the brute force sampling proposed above requires $N_{\mathtt{\text{paths}}}$ evaluations of $\tilde{H}$ (see 2(b) in the previous algorithm). This is a massive computational cost, even if a single evaluation of an ANN is high-speed. We can, however, benefit from a further approximation. We compute the CVs using $\tilde{H}$ only at specific reference values for $S\left(T\right)$. Then, the intermediate cases are derived utilizing (linear) interpolation. We choose a set of Q equally-spaced values $\{S^1,S^2,\dots ,S^Q\}$ for $S\left(T\right)$, defined as: (19) $S^q:=S_{min}+\frac{q-1}{Q-1}(S_{max}-S_{min}),q=1,\dots ,Q,$(19) where the boundaries are quantiles corresponding to the probabilities $p_{min},p_{max}\in (0,1)$, i.e. $S_{min}:=F_{S\left(T\right)}^{-1}\left(p_{min}\right)$ and $S_{max}:=F_{S\left(T\right)}^{-1}\left(p_{max}\right)$.

Calling ${\mathbf{p}}^q={\mathbf{p}}_{S^q}$, and ${\mathbf{a}}^q={\mathbf{a}}_{S^q}$, $q=1,\dots ,Q$, we compute the grid $\mathbf{G}$ of reference CPs, with only Q ANN evaluations, namely: (20) $\left[\begin{array}{c}{\mathbf{a}}^1\\ \cdots \\ {\mathbf{a}}^q\\ \cdots \\ {\mathbf{a}}^Q\end{array}\right]=\left[\begin{array}{c}H\left({\mathbf{p}}^1\right)\\ \cdots \\ H\left({\mathbf{p}}^q\right)\\ \cdots \\ H\left({\mathbf{p}}^Q\right)\end{array}\right]\approx \left[\begin{array}{c}\tilde{H}\left({\mathbf{p}}^1\right)\\ \cdots \\ \tilde{H}\left({\mathbf{p}}^q\right)\\ \cdots \\ \tilde{H}\left({\mathbf{p}}^Q\right)\end{array}\right]=:\mathbf{G},$(20) where ${\mathbf{a}}^q$, $H\left({\mathbf{p}}^q\right)$ and $\tilde{H}\left({\mathbf{p}}^q\right)$, $q=1,\dots ,Q$, are row vectors. The interpolation on $\mathbf{G}$ is much faster than the evaluation of $\tilde{H}$. Therefore, the grid-based conditional sampling results more efficient than the brute force one, particularly when sampling a huge number of MC samples.

The algorithm for the grid-based sampling procedure, to be used instead of point 2. in the previous algorithm, is reported here.

Algorithm: Grid-based conditional sampling

2.1.	Fix the boundary probabilities $p_{min},p_{max}\in (0,1)$ and compute the boundary quantiles $S_{min}:=F_{S\left(T\right)}^{-1}\left(p_{min}\right)$ and $S_{max}:=F_{S\left(T\right)}^{-1}\left(p_{max}\right)$ (e.g. with the COS method [24]).
2.2.	Compute the reference values $S^q:=S_{min}+\frac{q-1}{Q-1}(S_{max}-S_{min})$, $q=1,\dots ,Q$.
2.3.	Given the ‘augmented’ parameters ${\mathbf{p}}^q$, evaluate Q times $\tilde{H}$ to compute $\mathbf{G}$ (see (20(20) $\left[\begin{array}{c}{\mathbf{a}}^1\\ \cdots \\ {\mathbf{a}}^q\\ \cdots \\ {\mathbf{a}}^Q\end{array}\right]=\left[\begin{array}{c}H\left({\mathbf{p}}^1\right)\\ \cdots \\ H\left({\mathbf{p}}^q\right)\\ \cdots \\ H\left({\mathbf{p}}^Q\right)\end{array}\right]\approx \left[\begin{array}{c}\tilde{H}\left({\mathbf{p}}^1\right)\\ \cdots \\ \tilde{H}\left({\mathbf{p}}^q\right)\\ \cdots \\ \tilde{H}\left({\mathbf{p}}^Q\right)\end{array}\right]=:\mathbf{G},$(20) )).
2.4.	Given the parameters $\mathbf{p}$ and the grid $\mathbf{G}$, for $j=1,\dots ,N_{\mathtt{\text{paths}}}$, repeat: generate the sample ${\hat{S}}_j$ from $S\left(T\right)$ (e.g. with COS-SC method [12,24]); given ${\hat{S}}_j$, approximate the M conditional CVs, i.e. ${\mathbf{a}}_{{\hat{S}}_j}$, by interpolation in $\mathbf{G}$; given ${\mathbf{a}}_{{\hat{S}}_j}$, use SC to generate the conditional sample ${\hat{A}}_j$.

4. Error analysis

This section is dedicated to the assessment and discussion of the error introduced by the main approximations used in the proposed pricing method. Two primary sources of error are identifiable. The first error is due to the SC technique: in Section 3.1 the exact map g is approximated by means of the piecewise polynomial $\tilde{g}$. The second one is a regression error, which is present in both Sections 3.2 and 3.3. ANNs $\tilde{H}$ are used instead of the exact mappings H. For the error introduced by the SC technique, we bound the ‘$L^2$-distance’, ${ϵ}_{\mathrm{SC}}$, between the exact distribution and its SC proxy showing that $\tilde{g}{\mathbb{1}}_{{\mathrm{\Omega }}_M}=g_M$ in ${\mathrm{\Omega }}_M$ is an analytic function. ${ϵ}_{\mathrm{SC}}$ is used to provide a direct bound on the option price error, ${ϵ}_P$. On the other hand, regarding the approximation of H via $\tilde{H}$ we provide a general convergence result for ReLU-architecture ANN, i.e. ANN with Rectified Linear Units as activation functions.

4.1. Stochastic collocation error using Chebyshev polynomials

Let us consider the error introduced in the methodology using the SC technique Section 3.1, and investigate how this affects the option price. We restrict the analysis to the case of fixed- or floating-strike discrete arithmetic Asian and Lookback options Section 3.2. We define the error ${ϵ}_P$ as the ‘$L^1$-distance’ between the real price $V_{\omega }\left(t_0\right)$ and its approximation ${\tilde{V}}_{\omega }\left(t_0\right)$, i.e.: (21) ${ϵ}_P:=\left|{\tilde{V}}_{\omega }\right(t_0)-V_{\omega }(t_0\left)\right|.$(21) Given the standard normal kernel $\xi \sim \mathcal{N}(0,1)$, we define the SC error as the (squared) $L^2$-norm of $g-\tilde{g}$, i.e.: (22) ${ϵ}_{\mathrm{SC}}:=\mathbb{E}\left[\right(g-\tilde{g}{)}^2\left(\xi \right)].$(22) We decompose ${ϵ}_{\mathrm{SC}}$ accordingly to the piecewise definition of $\tilde{g}$, namely: $\begin{array}{rl}{ϵ}_{\mathrm{SC}}& =\mathbb{E}\left[\right(g{\mathbb{1}}_{{\mathrm{\Omega }}_{-}}-g_{-}{)}^2\left(\xi \right)]+\mathbb{E}\left[\right(g{\mathbb{1}}_{{\mathrm{\Omega }}_M}-g_M{)}^2\left(\xi \right)]+\mathbb{E}\left[\right(g{\mathbb{1}}_{{\mathrm{\Omega }}_{+}}-g_{+}{)}^2\left(\xi \right)]\\ & =:{ϵ}_{-}+{ϵ}_M+{ϵ}_{+}.\end{array}$with the domains ${\mathrm{\Omega }}_{(\cdot )}$ defined in Equation (14(14) ${\mathrm{\Omega }}_{-}\cup {\mathrm{\Omega }}_M\cup {\mathrm{\Omega }}_{+}:=(-\mathrm{\infty },-\overline{\xi })\cup [-\overline{\xi },\overline{\xi }]\cup (\overline{\xi },+\mathrm{\infty }),$(14) ), i.e. for $\overline{\xi }>0$: ${\mathrm{\Omega }}_{-}=(-\mathrm{\infty },-\overline{\xi }),{\mathrm{\Omega }}_M=[-\overline{\xi },\overline{\xi }],{\mathrm{\Omega }}_{+}=(\overline{\xi },+\mathrm{\infty }).$To deal with the ‘extrapolation’ errors ${ϵ}_{-}$ and ${ϵ}_{+}$, we formulate the following assumption.

Assumption 4.1

The functions $(g{\mathbb{1}}_{{\mathrm{\Omega }}_{-}}-g_{-}{)}^2$ and $(g{\mathbb{1}}_{{\mathrm{\Omega }}_{+}}-g_{+}{)}^2$ are $O(\exp x^2/2)$. Equivalently, $g^2{\mathbb{1}}_{{\mathrm{\Omega }}_{-}}$ and $g^2{\mathbb{1}}_{{\mathrm{\Omega }}_{+}}$ are $O(\exp x^2/2)$ (since $g_{{}_{-}}$ and $g_{{}_{+}}$ are polynomials).

Given Assumption 4.1 and the fact that $\xi \sim \mathcal{N}(0,1)$⁸, then the ‘extrapolation’ errors ${ϵ}_{-}$ and ${ϵ}_{+}$ vanish, with exponential rate, as $\overline{\xi }$ tends to infinity, i.e. ${ϵ}_{-}={ϵ}_{-}\left(\overline{\xi }\right)$, ${ϵ}_{+}={ϵ}_{+}\left(\overline{\xi }\right)$, and: (23) $\{\begin{array}{l}{ϵ}_{-}\left(\overline{\xi }\right)\to 0,\\ {ϵ}_{+}\left(\overline{\xi }\right)\to 0,\end{array}\mathrm{for}\overline{\xi }\to +\mathrm{\infty }.$(23) An illustration of the speed of convergence is reported in Figure 1. Figure 1a shows that the growth of $g{\mathbb{1}}_{{\mathrm{\Omega }}_{+}}$ is (much) less than exponential (consistently with Assumption 4.1), whereas Figure 1(b) illustrates the exponential decay of ${ϵ}_{-}$ and ${ϵ}_{+}$ when $\overline{\xi }$ increases.

Figure 1. Left: (slow) growth of $g{\mathbb{1}}_{{\mathrm{\Omega }}_{+}}$ compared to the linear extrapolation $g_{+}$. Right: exponential decays of ${ϵ}_{-}\left(\overline{\xi }\right)$ and ${ϵ}_{+}\left(\overline{\xi }\right)$ in (23(23) $\{\begin{array}{l}{ϵ}_{-}\left(\overline{\xi }\right)\to 0,\\ {ϵ}_{+}\left(\overline{\xi }\right)\to 0,\end{array}\mathrm{for}\overline{\xi }\to +\mathrm{\infty }.$(23) ) when the $\overline{\xi }$ increases. The upper x-axis represents the probability $F_{\xi }\left(\overline{\xi }\right)$.

Therefore, if $\overline{\xi }$ is taken sufficiently big, the error ${\xi }_{\mathrm{SC}}$ in (22(22) ${ϵ}_{\mathrm{SC}}:=\mathbb{E}\left[\right(g-\tilde{g}{)}^2\left(\xi \right)].$(22) ) is mainly driven by the ‘interpolation’ error ${ϵ}_M$, whose estimate is connected to error bounds for Chebyshev polynomial interpolation, and it is the focus of the next part.

Theorem 4.1

Error bound for analytic function [9,29]

Let f be a real function on $[-1,1]$ and $f_M$ be its $(M-1)$-degree polynomial interpolation built on Chebyshev nodes ${\xi }_k:=\cos \frac{k-1}{M-1}\pi$, $k=1,\dots ,M$. If f has an analytic extension in a Bernstein ellipse $\mathcal{B}$ with foci $\pm 1$ and major and minor semiaxis lengths summing up to $\varrho >1$ such that $\underset{\mathcal{B}}{sup}|f|\le \frac{\varrho -1}{4}\overline{C}$ for some constant $\overline{C}>0$, then, for each $M\ge 1$, the following bound holds: $||f-f_M|{|}_{L^{\mathrm{\infty }}\left(\right[-1,1\left]\right)}\le \overline{C}{\varrho }^{1-M}.$

Since $g{\mathbb{1}}_{{\mathrm{\Omega }}_M}$ is approximated by means of the $(M-1)$-degree polynomials $g_M$, built on Chebyshev nodes, to apply Theorem 4.1, we verify the required assumptions, namely the boundedness of $g{\mathbb{1}}_{{\mathrm{\Omega }}_M}$ in ${\mathrm{\Omega }}_M$ and its analyticity.

We recall that: (24) $g=F_A^{-1}\circ F_{\xi },$(24) with $F_A$ and $F_{\xi }$ the CDFs of $A\left(S\right)$ and ξ, respectively. Hence, the boundedness on the compact domain ${\mathrm{\Omega }}_M$ is satisfied because the map g is monotone increasing (as a composition of monotone increasing functions), and defined everywhere in ${\mathrm{\Omega }}_M$.

Furthermore, since the CDF of a standard normal, $F_{\xi }$, is analytic, from (24(24) $g=F_A^{-1}\circ F_{\xi },$(24) ) it follows that g is analytic if $F_A^{-1}$ is analytic. The analyticity of $F_A^{-1}$ is fulfilled if $F_A$ is analytic and $F_A^{\prime }=f_A$ does not vanish in the domain ${\mathrm{\Omega }}_M$. Observe that, by restricting the domain to ${\mathrm{\Omega }}_M$, the latter condition is trivially satisfied because we are ‘far’ from the tails of $A\left(S\right)$ (corresponding to the extrapolation domains ${\mathrm{\Omega }}_{-}$ and ${\mathrm{\Omega }}_{+}$), and $F_A^{\prime }$ do not vanish in other regions than the tails.

On the contrary, proving that $F_A$ is analytic is not trivial because of the lack of an explicit formula for $F_A$. However, it is beneficial to represent $F_A$ through the characteristc function (ChF) of $A\left(S\right)$, ${\varphi }_A$. For that purpose, we use a well-known inversion result.

Theorem 4.2

ChF inversion theorem

Let us denote by F and ϕ the CDF and the ChF of a given real-valued random variable defined on $\mathbb{R}$. Then, it is possible to retrieve F from ϕ according to the inversion formula: $F\left(x\right)-F\left(0\right)=\frac{1}{2\pi }{\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}\varphi \left(u\right)\frac{1-{\mathrm{e}}^{-\mathrm{iux}}}{\mathrm{iu}}\mathrm{d}u,$with the integral being understood as a principal value.

Proof.

For detailed proof, we refer to [17].

Thanks to Theorem 4.2, we have that if ${\varphi }_A$ is analytic, so it is $F_A$ (as long as the integral in (eqn: L2error) is well-defined). Thus, the problem becomes to determine if ${\varphi }_A$ is analytic. We rely on a characterization of entire⁹ ChFs, which can be used in this framework to show that in the cases of – fixed- or floating-strike – discrete arithmetic Asian and Lookback options, the (complex extension of the) function ${\varphi }_A$ is analytic in a certain domain.

Theorem 4.3

Characterization of entire ChFs [4]

Let Y be a real random variable. Then, the complex function $\varphi \left(z\right):=\mathbb{E}\left[{\mathrm{e}}^{\mathrm{izY}}\right]$, $z\in \mathbb{C}$, is entire if and only if the absolute moments of Y exist for any order, i.e. $\mathbb{E}\left[|Y{|}^k\right]<+\mathrm{\infty }$ for any $k\in \mathbb{N}$, and the following limit holds: (25) $\underset{k\to +\mathrm{\infty }}{lim}{\left(\frac{|\mathbb{E}\left[Y^k\right]|}{k!}\right)}^{\frac{1}{k}}=0.$(25)

Proof.

A reference for proof is given in [4].

When dealing with the Heston model, there is no closed-form expression for the moments of the underlying process $S\left(t\right)$, as well as for the moments of its transform $A\left(S\right)$. Nonetheless, a conditional case can be studied and employed to provide a starting point for a convergence result.

Proposition 4.4

Conditional ChF ${\varphi }_{A|\mathbf{V}}$ is entire

Let us define the N-dimensional random vector $\mathbf{V}$, with values in ${\mathrm{\Omega }}_{\mathbf{V}}:={\mathbb{R}}_{+}^N$, as: $\mathbf{V}:=\left[I_v\right(t_1),I_v(t_2),\dots ,I_v(t_N){]}^T,I_v(t_n):={\int }_{t_0}^{t_n}v(\tau )\mathrm{d}\tau ,n=1,\dots ,N.$Let the complex conditional characteristic function ${\varphi }_{A|\mathbf{V}}\left(z\right):=\mathbb{E}\left[{\mathrm{e}}^{\mathrm{izA}}|\mathbf{V}\right]$, $z\in \mathbb{C}$, be the extended ChF of the conditional random variable $A|\mathbf{V}$, with $A\equiv A\left(S\right)$ as given in Equation (2(2) $A\left(S\right):=\frac{1}{N}\sum _{n\in I}S\left(t_n\right),A\left(S\right):=\omega \underset{n\in I}{max}\mathrm{\omega S}\left(t_n\right),$(2) ).

Then, ${\varphi }_{A|\mathbf{V}}\left(z\right)$ is entire.

Proof.

See Appendix A.4.

From now on, using the notation of Proposition 4.4, we consider the following assumption on the tail behaviour of the random vector $(\mathbf{V},A)$ satisfied. Informally, we require that the density of the joint distribution of $(\mathbf{V},A)$ has uniform (w.r.t. $\mathbf{V}$) exponential decay for A going to $+\mathrm{\infty }$.

Assumption 4.2

There exists a point $z^{\ast }\in \mathbb{C}$, $z^{\ast }=x^{\ast }-i{y}^{\ast }$, with $x^{\ast },y^{\ast }\in \mathbb{R}$ and $y^{\ast }>0$, such that: ${\int }_{{\mathrm{\Omega }}_{\mathbf{V}}\times {\mathbb{R}}^{+}}{\mathrm{e}}^{y^{\ast }a}\mathrm{d}{F}_{\mathbf{V},A}(\mathbf{v},a)<+\mathrm{\infty },$with $F_{\mathbf{V},A}(\cdot ,\cdot )$ the joint distribution of the random vector $\mathbf{V}$ and the random variable A.

Thanks to Assumption 4.2, the ChF ${\varphi }_A\left(z\right)$ is well defined for any $z\in {\mathcal{S}}_{y^{\ast }}\subset \mathbb{C}$, with the strip ${\mathcal{S}}_{y^{\ast }}:=\mathbb{R}+i\cdot [-y^{\ast },y^{\ast }]$. Moreover, applying Fubini's Theorem, for any $z\in {\mathcal{S}}_{y^{\ast }}$, we have: (26) ${\varphi }_A\left(z\right)={\int }_{{\mathrm{\Omega }}_{\mathbf{V}}}{\varphi }_{A|\mathbf{V}=\mathbf{v}}\left(z\right)\mathrm{d}{F}_{\mathbf{V}}\left(\mathbf{v}\right).$(26) Thus, we can show that the ChF ${\varphi }_A\left(z\right)$ is analytic in the strip ${\mathcal{S}}_{y^{\ast }}$ (the details are given in Appendix A.2).

Proposition 4.5

ChF ${\varphi }_A$ is analytic

Let ${\varphi }_A\left(z\right):=\mathbb{E}\left[{\mathrm{e}}^{\mathrm{izA}}\right]$, $z\in {\mathcal{S}}_{y^{\ast }}$, with $A\equiv A\left(S\right)$. Then, ${\varphi }_A\left(z\right)$ is analytic in ${\mathcal{S}}_{y^{\ast }}$.

Proof.

A proof is given in Appendix A.2.

Thanks to Proposition 4.5, and consistently with the previous discussion, we conclude that the map g in (24(24) $g=F_A^{-1}\circ F_{\xi },$(24) ), is analytic on the domain ${\mathrm{\Omega }}_M$. Therefore, we can apply Theorem 4.1, which yields the following error estimate: $||g{\mathbb{1}}_{{\mathrm{\Omega }}_M}-g_M|{|}_{L^{\mathrm{\infty }}\left({\mathrm{\Omega }}_M\right)}\le \overline{C}{\varrho }^{1-M},$for certain $\varrho >1$ and $\overline{C}>0$. As a consequence, the following bound for the $L^2$-error ${ϵ}_M$ holds: (27) ${ϵ}_M=\mathbb{E}\left[\right(g{\mathbb{1}}_{{\mathrm{\Omega }}_M}-g_M{)}^2\left(\xi \right)]\le {\overline{C}}^2{\varrho }^{2-2M}.$(27) Furthermore, the exponential convergence is also confirmed numerically, as reported in Figure 2. In Figure 2(a) we can appreciate the improvement in the approximation of $\tilde{g}{\mathbb{1}}_{{\mathrm{\Omega }}_M}$ by means of $g_M$, when M is increased, whereas Figure 2(b) reports the exponential decay of ${ϵ}_M$.

Figure 2. Left: exact map $g{\mathbb{1}}_{{\mathrm{\Omega }}_M}$ (blue) compared to the interpolation $g_M$ in the domain ${\mathrm{\Omega }}_M$, for M = 4, 5. Right: $L^2$-error ${ϵ}_M$ in (27(27) ${ϵ}_M=\mathbb{E}\left[\right(g{\mathbb{1}}_{{\mathrm{\Omega }}_M}-g_M{)}^2\left(\xi \right)]\le {\overline{C}}^2{\varrho }^{2-2M}.$(27) ) exponential decay in the order of the polynomial $g_M$.

Using (27(27) ${ϵ}_M=\mathbb{E}\left[\right(g{\mathbb{1}}_{{\mathrm{\Omega }}_M}-g_M{)}^2\left(\xi \right)]\le {\overline{C}}^2{\varrho }^{2-2M}.$(27) ), the $L^2$-norm of $g-\tilde{g}$, ${ϵ}_{\mathrm{SC}}$ in (22(22) ${ϵ}_{\mathrm{SC}}:=\mathbb{E}\left[\right(g-\tilde{g}{)}^2\left(\xi \right)].$(22) ), is bounded by: ${ϵ}_{\mathrm{SC}}\equiv {ϵ}_{\mathrm{SC}}(\overline{\xi },M)=\mathbb{E}\left[\right(g\left(\xi \right)-\tilde{g}\left(\xi \right){)}^2]\le {ϵ}_{-}\left(\overline{\xi }\right)+{\overline{C}}^2{\varrho }^{2-2M}+{ϵ}_{+}\left(\overline{\xi }\right),$which goes to zero when $\overline{\xi }\in {\mathbb{R}}^{+}$ and $M\in \mathbb{N}$ tend to $+\mathrm{\infty }$. Therefore, for any $ϵ>0$ there exist ${\overline{\xi }}^{\ast }\in {\mathbb{R}}^{+}$ and $M^{\ast }\in \mathbb{N}$ such that: (28) ${ϵ}_{\mathrm{SC}}({\overline{\xi }}^{\ast },M^{\ast })<{ϵ}^2,$(28) and because of the exponential decay, we expect ${\overline{\xi }}^{\ast }$ and $M^{\ast }$ do not need to be taken too big.

Eventually, we can benefit from the bound in (28(28) ${ϵ}_{\mathrm{SC}}({\overline{\xi }}^{\ast },M^{\ast })<{ϵ}^2,$(28) ) to control the pricing error, ${ϵ}_P$ in (21(21) ${ϵ}_P:=\left|{\tilde{V}}_{\omega }\right(t_0)-V_{\omega }(t_0\left)\right|.$(21) ). By employing the well-known inequality $max(a+b,0)\le max(a,0)+max(b,0)$ and the Cauchy-Schwarz inequality, we can write: $\begin{array}{rl}{\tilde{V}}_{\omega }\left(t_0\right)=\mathbb{E}\left[{\left(\omega \left(\tilde{g}\right(\xi )-K)\right)}^{+}\right]& \le \sqrt{\mathbb{E}\left[{\left(\tilde{g}\right(\xi )-g(\xi \left)\right)}^2\right]}+\mathbb{E}\left[{\left(\omega \left(g\right(\xi )-K)\right)}^{+}\right]\\ & \le \sqrt{{ϵ}_{\mathrm{SC}}(\overline{\xi },M)}+V_{\omega }\left(t_0\right),\end{array}$and using the same argument twice (exchanging the roles of g and $\tilde{g}$), we end up with the following bound for the option price error: ${ϵ}_P\le \sqrt{{ϵ}_{\mathrm{SC}}({\overline{\xi }}^{\ast },M^{\ast })}\le ϵ,$with ${\overline{\xi }}^{\ast }$ and $M^{\ast }$ as in (28(28) ${ϵ}_{\mathrm{SC}}({\overline{\xi }}^{\ast },M^{\ast })<{ϵ}^2,$(28) ).

4.2. Artificial neural network regression error

As the final part of the error analysis, we investigate when ANNs are suitable approximating maps. In particular, we focus on ANNs with ReLU-architectures, namely ANNs whose activation units are all Rectified Linear Units defined as $\varphi \left(x\right)=x{\mathbb{1}}_{x>0}\left(x\right)$.

Consider the Sobolev space $\left({\mathcal{W}}^{n,\mathrm{\infty }}\right([0,1{]}^d),||\cdot |{|}_{n,d}^{\mathrm{\infty }})$, with $n,d\in \mathbb{N}\mathrm{\setminus }\left\{0\right\}$, namely the space of functions $C^{n-1}\left(\right[0,1{]}^d)$ whose derivatives up to the $(n-1)$th order are all Lipschitz continuous, equipped with the norm $||\cdot |{|}_{n,d}^{\mathrm{\infty }}$ defined as: (29) $||f|{|}_{n,d}^{\mathrm{\infty }}=\underset{|\mathbf{n}|\le n}{max}\underset{\mathbf{x}\in [0,1{]}^d}{\mathrm{ess}\sup }|D^{\mathbf{n}}f\left(\mathbf{x}\right)|,$(29) with $\mathbf{n}:=(n_1,\dots ,n_d)\in {\mathbb{N}}^d$, $|\mathbf{n}|=\sum _{i=1}^d{n}_i$, and $D^{\mathbf{n}}$ the weak derivative operator. Furthermore, we define the unit ball $B_{n,d}:=\{f\in {\mathcal{W}}^{n,\mathrm{\infty }}([0,1{]}^d):||f|{|}_{n,d}^{\mathrm{\infty }}\le 1\}$. Then, the following approximation result holds:

Theorem 4.6

Convergence for ReLU ANN

For any choice of $d,n\in \mathbb{N}\mathrm{\setminus }\left\{0\right\}$ and $ϵ\in (0,1)$, there exists an architecture $H(\mathrm{x}|\cdot )$ based on ReLU (Rectified Linear Unit) activation functions ϕ, i.e. $\varphi \left(x\right)=x{\mathbb{1}}_{x>0}\left(x\right)$, such that:

(1)	$H(\mathrm{x}|\cdot )$ is able to approximate any function $f\in B_{n,d}$ with an error smaller than ϵ, i.e. there exists a matrix of weights $\mathrm{W}$ such that $||f(\cdot )-H(\cdot |\mathrm{W})|{|}_{\mathit{n}\mathit{,}\mathit{d}}^{\mathrm{\infty }}<ϵ$;
(2)	H has at most $c(\ln 1/ϵ+1)$ layers and at most $c{ϵ}^{-d/n}(\ln 1/ϵ+1)$ weights and neurons, with $c=c(d,n)$ an appropriate constant depending on d and n.

Proof.

A proof is available in [30].

Essentially, Theorem 4.6 states that there always exists a ReLU-architecture (with bounded number of layers and activation units) suitable to approximate at any desired precision functions with a certain level of regularity (determined by $\left({\mathcal{W}}^{n,\mathrm{\infty }}\right([0,1{]}^d),||\cdot |{|}_{n,d}^{\mathrm{\infty }})$).

Remark 4.7

Input scaling

We emphasize that although Theorem 4.6 applies to (a subclass of sufficiently regular) functions with domain the d-dimensional hypercube $[0,1{]}^d$, this is not restrictive. Indeed, as long as the regularity conditions are fulfilled, Theorem 4.6 holds for any function defined on a d-dimensional hyperrectangle since it is always possible to linearly map its domain into the d-dimensional hypercube.

Furthermore, we observe that all the results of convergence for ANN rely on the assumption that the training is performed successfully, and so the final error in the optimization process is negligible. Under this assumption, Theorem 4.6 provides a robust theoretical justification for using ReLU-based ANNs as regressors. The goodness of the result can also be investigated empirically, as shown in the next section (see, for instance, Figure 4).

Figure 3. Illustration of a dense ANN with (from the left) one input layer, three hidden layers, and one output layer. Each white node represents an activation unit.

Figure 4. First experiment: FxA. Left: scatter plot of the real CVs ($a_k$, $k=1,\dots ,21$, with different colours) against the predicted ones. Right: zoom on the ‘worst’ case, namely $a_{10}$.

5. Numerical experiments

In this part of the paper, we detail some numerical experiments. We focus on applying the methodology given in Section 3.2 for the numerical pricing of fixed-strike discrete arithmetic Asian and Lookback options. We address the general case of discrete arithmetic Asian options described in Section 3.3. For each pricing experiment errors and timing results are given. The ground truth benchmarks are computed via MC using the almost exact simulation of the Heston models, detailed in Result Appendix A.3.

All the computations are implemented and run on a MacBook Air (M1, 2020) machine, with chip Apple M1 and 16 GB of RAM. The code is written in Python, and torch is the library used for the design and training of the ANN, as in [25].

5.1. A benchmark from the literature

To assess the quality of the proposed methodology, we compare the method against the benchmarks available in Table 5 from [19]. In the experiment, we consider prices of 5 discrete Asian call options with n = 201 equally spaced monitoring dates from time $t_0=0$ to T = 0.25. The underlying initial value is $S\left(t_0\right)=100$, while the other Heston parameters (Set BM) as well as the target strikes are given in Table 2. To produce those results, a toy model has been trained based on the ranges provided in Table 1. The ANN employed consists of 2 hidden layers each one with 20 hidden units. The results are computed with $N_{\mathtt{\text{paths}}}={10}^6$ Monte Carlo paths (consistent with the benchmark [19]) and are presented in Table 2. For both the benchmark and the SC technique are reported the absolute value (V) and the $95\mathrm{\%}$ confidence interval (95CI) of the option price. For the SA technique, instead, only the absolute value is reported since no sampling is involved and so there is no information on the variance of the estimate. All the results are within the BM $95\mathrm{\%}$ confidence interval, and hence confirm the high accuracy of the proposed method.

5.2. Experiments' specifications

Among the three examples of applications presented, two of them rely on the technique given in Section 3.2, while the third is based on the theory in Section 3.3. The first experiment is the pricing of fixed-strike discrete arithmetic Asian options (FxA) with an underlying stock price process following the Heston dynamics. The second example, instead, is connected to the ‘interest rate world’, and is employed for the pricing of fixed-strike discrete Lookback swaptions (FxL). We assume the underlying swap rate is driven by a displaced Heston model with drift-less dynamics, typically used for interest rates. The last one is an application to the pricing of fixed- and floating-strikes discrete arithmetic Asian options (FxFlA) on a stock price driven by the Heston dynamics. In the first (FxA) and last experiment (FxFlA), $A\left(S\right)$ in (2(2) $A\left(S\right):=\frac{1}{N}\sum _{n\in I}S\left(t_n\right),A\left(S\right):=\omega \underset{n\in I}{max}\mathrm{\omega S}\left(t_n\right),$(2) ) is specified as: (30) $A\left(S\right)=\frac{1}{5}\sum _{n=1}^{5}S\left(t_n\right),t_n:=T-(5-n){\tau }_{\mathtt{A}}\mathtt{,}\mathtt{n}\mathtt{=}\mathtt{1}\mathtt{,}\dots \mathtt{,}\mathtt{5}\mathtt{,}$(30) with ${\tau }_{\mathtt{\text{A}}}=\frac{1}{12}$ as monitoring time lag, and $T>4{\tau }_{\mathtt{\text{A}}}$ as option maturity. Differently, in the second experiment (FxL) $A\left(S\right)$ is given by: (31) $A\left(S\right)=\underset{n}{min}S\left(t_n\right),t_n:=T-(30-n){\tau }_{\mathtt{L}}\mathtt{,}\mathtt{n}\mathtt{=}\mathtt{1}\mathtt{,}\dots \mathtt{,}\mathtt{30}\mathtt{,}$(31) with ${\tau }_{\mathtt{\text{L}}}=\frac{1}{120}$ as monitoring time lag, and $T>29{\tau }_{\mathtt{\text{L}}}$ as option maturity. Observe that, assuming the unit is 1 year, with 12 identical months and 360 days, the choices of ${\tau }_{\mathtt{\text{A}}}$ and ${\tau }_{\mathtt{\text{L}}}$ correspond respectively to 1 month and 3 days of time lag in the monitoring dates.

Table 1. Training set for benchmark replication.

$\mathbf{p}$	$N_{\mathtt{\text{val}}}$	range	met.	MC
r	100	[0.04, 0.06]	LHS	$N_{\mathtt{\text{paths}}}$	${10}^6$
κ	100	[2.90, 3.10]	LHS	$\mathrm{\Delta t}$	$\frac{1}{800}$
γ	100	[0.08, 0.12]	LHS
ρ	100	[−0.11, −0.09]	LHS	SC
$\overline{v}$	100	[0.03, 0.05]	LHS	M	21
$v_0$	100	[0.03, 0.05]	LHS	$\overline{\xi }$	2.46
T	25	[0.25, 0.28]	EQ-SP	$F_{\overline{\xi }}\left(\overline{\xi }\right)$	0.993
$N_{\mathtt{\text{pairs}}}$	$2500(100\times 25)$

Table 2. Results benchmark (BM) replication (see, Table 5 from [19]).

	Set BM
	r	κ	γ	ρ	$\overline{v}$	$v_0$	T
	0.05	3.0	0.1	-0.1	0.04	0.04	0.25
	BM			SC			SA
$K_2$	V	95CI		V	95CI		V
90	10.5439	[10.5329, 10.5550]		10.5459	[10.5349, 10.5570]		10.5486
95	6.0168	[6.0069, 6.0267]		6.0173	[6.0074, 6.0271]		6.0196
100	2.6026	[2.5953, 2.6098]		2.5992	[2.5920, 2.6064]		2.5996
105	0.7902	[0.7862, 0.7943]		0.7867	[0.7827, 0.7907]		0.7865
110	0.1622	[0.1604, 0.1639]		0.1612	[0.1595, 0.1630]		0.1615

5.3. Artificial neural network development

In this section, we provide the details about the generation of the training set, for each experiment, and the consequent training of the ANN used in the pricing model.

5.3.1. Training set generation

The training sets are generated through MC simulations, using the almost exact sampling in Result Appendix A.3. In the first two applications (FxA and FxL) the two training sets are defined as in (18(18) $\mathbf{T}=\left\{\right({\mathbf{p}}_i,{\mathbf{a}}_i):i\in \{1,\dots ,N_{\mathtt{pairs}}\left\}\right\}.$(18) ), and particularly they read: $\begin{array}{rl}{\mathbf{T}}_{\mathtt{FxA}}& =\{\left(\right\{r,\kappa ,\gamma ,\rho ,\overline{v},v_0,T{\}}_i,\{a_1,\dots ,a_{21}{\}}_i):i\in \{1,\dots ,N_{\mathtt{pairs}}^{\mathtt{FxA}}\left\}\right\},\\ {\mathbf{T}}_{\mathtt{FxL}}& =\{\left(\right\{\kappa ,\gamma ,\rho ,\overline{v},v_0,T{\}}_i,\{a_1,\dots ,a_{21}{\}}_i):i\in \{1,\dots ,N_{\mathtt{pairs}}^{\mathtt{FxL}}\left\}\right\}.\end{array}$The Heston parameters, i.e. $\mathbf{p}\mathrm{\setminus }\left\{T\right\}$¹⁰, are sampled using Latin Hypercube Sampling (LHS), to ensure the best filling of the parameter space [22,25]. For each set $\mathbf{p}\mathrm{\setminus }\left\{T\right\}$, $N_{\mathtt{\text{paths}}}$ paths are generated, with a time step of $\mathrm{\Delta t}$ and a time horizon up to $T_{max}$. The underlying process S is monitored at each time T for which there are enough past observations to compute $A\left(S\right)$, i.e.: $\begin{array}{rl}T& \ge 4{\tau }_{\mathtt{A}}+\mathrm{\Delta t},\mathrm{for}\mathtt{FxA}\mathtt{,}\\ T& \ge 29{\tau }_{\mathtt{L}}+\mathrm{\Delta t},\mathrm{for}\mathtt{FxL}\mathtt{.}\end{array}$Consequently, the product between the number of Heston parameters' set and the number of available maturities determines the magnitude of the two training sets (i.e. $N_{\mathtt{\text{pairs}}}^{\mathtt{\text{FxA}}}$ and $N_{\mathtt{\text{pairs}}}^{\mathtt{\text{FxL}}}$).

For each $\mathbf{p}$, the CVs $\mathbf{a}$ corresponding to $A\left(S\right)$ are computed as: $a_k:=F_A^{-1}\left(F_{\xi }\right({\xi }_k\left)\right)\approx Q_A\left(F_{\xi }\right({\xi }_k\left)\right),k\in \{1,\dots ,21\},$where $Q_A$ is the empirical quantile function of $A\left(S\right)$, used as a numerical proxy of $F_A^{-1}$, and ${\xi }_k$ are the CPs computed as Chebyshev nodes: ${\xi }_k:=-\overline{\xi }\cos \frac{k-1}{20}\pi ,k\in \{1,\dots ,21\},$with $\overline{\xi }:=F_{\xi }^{-1}\left(0.993\right)\approx 2.46$. We note that the definition of $\overline{\xi }$ avoids any CV $a_k$ to be ‘deeply’ in the tails of $A\left(S\right)$, which are more sensitive to numerical instability in a MC simulation.

The information about the generation of the two training sets is reported in Table 3. Observe that ${\mathbf{T}}_{\mathtt{\text{FxA}}}$ is richer in elements than ${\mathbf{T}}_{\mathtt{\text{FxL}}}$ because of computational constraints. Indeed, the higher number of monitoring dates of $A\left(S\right)$ in FxL makes the generation time of ${\mathbf{T}}_{\mathtt{\text{FxL}}}$ more than twice the one of ${\mathbf{T}}_{\mathtt{\text{FxA}}}$ (given the same number of pairs).

Table 3. Training sets ${\mathbf{T}}_{\mathtt{\text{FxA}}}$ and ${\mathbf{T}}_{\mathtt{\text{FxL}}}$ generation details.

	FxA			FxL
$\mathbf{p}$	$N_{\mathtt{\text{val}}}$	range	met.	$N_{\mathtt{\text{val}}}$	range	met.	MC
r	700	[0.00, 0.05]	LHS				$N_{\mathtt{\text{paths}}}$	${10}^6$
κ	700	[0.20, 1.10]	LHS	300	[0.80, 1.60]	LHS	$\mathrm{\Delta t}$	$\frac{1}{120}$
γ	700	[0.80, 1.10]	LHS	300	[0.40, 1.00]	LHS
ρ	700	[−0.95, −0.20]	LHS	300	[−0.80, −0.30]	LHS	SC
$\overline{v}$	700	[0.02, 0.15]	LHS	300	[0.10, 0.20]	LHS	M	21
$v_0$	700	[0.02, 0.15]	LHS	300	[0.10, 0.20]	LHS	$\overline{\xi }$	2.46
T	160	[0.34, 1.67]	EQ-SP	170	[0.26, 1.67]	EQ-SP	$F_{\xi }\left(\overline{\xi }\right)$	0.993
$N_{\mathtt{\text{pairs}}}$	$112000(700\times 160)$			$51000(300\times 170)$

Since in the general procedure (see Section 3.3) ANNs are used to learn the conditional distribution $A\left(S\right)|S\left(T\right)$ (not just $A\left(S\right)$!), the third experiment requires a training set which contains also information about the conditioning value, $S\left(T\right)$. We define ${\mathbf{T}}_{\mathtt{\text{FxFlA}}}$ as: ${\mathbf{T}}_{\mathtt{FxFlA}}=\{\left(\right\{r,\kappa ,\gamma ,\rho ,\overline{v},v_0,T,S^q,p^q{\}}_i,\{a_1,\dots ,a_{14}{\}}_i):i\in \{1,\dots ,N_{\mathtt{pairs}}\left\}\right\},$where $p^q$ is the probability corresponding to the quantile $S^q$, given as in (19(19) $S^q:=S_{min}+\frac{q-1}{Q-1}(S_{max}-S_{min}),q=1,\dots ,Q,$(19) ), i.e.: $S^q=S_{min}+\frac{q-1}{14}(S_{max}-S_{min}),q\in \{1,\dots ,15\},$with the $S_{min}=F_{S\left(T\right)}^{-1}\left(p_{min}\right)$ and $S_{max}=F_{S\left(T\right)}^{-1}\left(p_{max}\right)$. Heuristic arguments drove the choice of adding in the input set $\mathbf{p}$ the probability $p^q:=F_{S\left(T\right)}\left(S^q\right)$, i.e. the probability implied by the final value $p^q$. Indeed, the ANN training process results more accurate when both $S^q$ and $p^q$ are included in $\mathbf{p}$. Similarly as before, the sets of Heston parameters are sampled using LHS. For each set, $N_{\mathtt{\text{tot}}}$ paths are generated, with a time step of $\mathrm{\Delta t}$ and a time horizon up to $T_{max}$. The underlying process S is monitored at each time T for which there are enough past observations to compute $A\left(S\right)$, i.e.: $T\ge 4{\tau }_{\mathtt{A}}+\mathrm{\Delta t}.$For any maturity T and any realization $S^q$, the inverse CDF of the conditional random variable $A\left(S\right)|S\left(T\right)=S^q$ is approximated with the empirical quantile function, $Q_{A|S^q}$. The quantile function $Q_{A|S^q}$ is built on the $N_{\mathtt{\text{paths}}}$ ‘closest’ paths to $S^q$, i.e. those $N_{\mathtt{\text{paths}}}$ paths whose final values $S\left(T\right)$ are the closest to $S^q$.

Eventually, for any input set $\mathbf{p}=\{r,\kappa ,\gamma ,\rho ,\overline{v},v_0,T,S^q,p^q\}$, the CVs $\mathbf{a}$ corresponding to $A\left(S\right)|S\left(T\right)=S^q$ are computed as: $a_k:=F_{A|S^q}^{-1}\left(F_{\xi }\right({\xi }_k\left)\right)\approx Q_{A|S^q}\left(F_{\xi }\right({\xi }_k\left)\right),k\in \{1,\dots ,14\},$with ${\xi }_k$ the Chebyshev nodes: ${\xi }_k:=-\overline{\xi }\cos \frac{k-1}{13}\pi ,k\in \{1,\dots ,14\},$and $\overline{\xi }:=F_{\xi }^{-1}\left(0.993\right)\approx 2.46$.

The information about the generation of the training set ${\mathbf{T}}_{\mathtt{\text{FxFlA}}}$ are reported in Table 4.

Table 4. Training set ${\mathbf{T}}_{\mathtt{\text{FxFlA}}}$ generation details.

$\mathbf{p}$	$N_{\mathtt{\text{val}}}$	range	met.	MC
r	180	[0.00, 0.05]	LHS	$N_{\mathtt{\text{tot}}}$	$3\cdot {10}^6$
κ	180	[0.20, 1.10]	LHS	$N_{\mathtt{\text{paths}}}$	$2\cdot {10}^5$
γ	180	[0.80, 1.10]	LHS	$\mathrm{\Delta t}$	$\frac{1}{120}$
ρ	180	[−0.92, −0.28]	LHS
$\overline{v}$	180	[0.03, 0.10]	LHS	SC
$v_0$	180	[0.03, 0.10]	LHS	Q	15
T	160	[0.34, 1.67]	EQ-SP	M	14
$S^q$	15	[0.35, 2.01]	EQ-SP	$\overline{\xi }$	2.46
$p^q$	15	[0.05, 0.85]	IMPL	$F_{\xi }\left(\overline{\xi }\right)$	0.993
$N_{\mathtt{\text{pairs}}}$	$432000(180\times 160\times 15)$

5.3.2. Artificial neural network training

Each training set ${\mathbf{T}}_{(\cdot )}$ store a finite amount of pairs $(\mathbf{p},\mathbf{a})$, in which each $\mathbf{p}$ and each corresponding $\mathbf{a}$ are connected by the mapping H. The artificial neural network $\tilde{H}$ is used to approximate and generalize H for inputs $\mathbf{p}$ not in ${\mathbf{T}}_{(\cdot )}$. The architecture of $\tilde{H}$ was initially chosen accordingly to [22,25], then suitably adjusted by heuristic arguments.

$\tilde{H}$ is a fully connected (or dense) ANN with five layers – one input, one output, and three hidden (HidL), as the one illustrated in Figure 3. Input and output layers have a number of units (neurons) – input size (InS) and output size (OutS) – coherent with the targeted problem (FxA, FxL, or FxFlA). Each hidden layer has the same hidden size (HidS) of 200 neurons, selected as the optimal one among different settings. ReLU (Rectified Linear Unit) is the non-linear activation unit (ActUn) for each neuron, and it is defined as $\varphi \left(x\right):=max(x,0)$ [23]. The loss function (LossFunc) is the Mean Squared Error (MSE) between the actual outputs, $\mathbf{a}$ (available in $\mathbf{T}$), and the ones predicted by the ANN, $\tilde{H}\left(\mathbf{p}\right)$. The optimization process is composed of 3,000 epochs (E). During each epoch, the major fraction (70%) of the $\mathbf{T}$ (the actual training set) is ‘back-propagated’ through the ANN in batches of size 1024 (B). The stochastic gradient-based optimizer (Opt) Adam [18] is employed in the optimization. Particularly, the optimizer updates the ANN weights based on the gradient computed on each random batch (during each epoch). The initial learning rate (InitLR) is ${10}^{-3}$, with a decaying rate (DecR) of 0.1 and a decaying step (DecS) of 1,000 epochs. The details are reported in Table 5.

Furthermore, during the optimization routine, the 20% of $\mathbf{T}$ is used to validate the result (namely, to avoid the overfitting of the training set). Eventually, the remaining 10% of $\mathbf{T}$ is used for testing the quality of the ANN. Figure 4 provides a visual insight into the high accuracy the ANN reaches at the end of the training process. Figure 4(a) shows the scatter plot of the real CVs $a_k$, $k=1,\dots ,21$, against the ones predicted using the ANN, for the experiment FxA; in Figure 4(b) a zoom is applied to the ‘worst’ case, namely the CV $a_{10}$, for which anyway is reached the extremely high $R^2$ score of 0.9994.

Table 5. Artificial neural network and optimization details.

ANN architecture						ANN optimization
	FxA	FxL	FxFlA	HidL	3	E	3,000	InitLR	${10}^{-3}$	Opt	Adam
InS	7	6	9	HidS	200	B	1024	DecR	0.1	LossFunc	MSE
OutS	21	21	14	ActUn	ReLU			DecS	1,000

5.4. Sampling and pricing

Given the trained model from the previous section, we can now focus on the actual sampling and/or pricing of options. In particular, for the first two experiments, we consider the following payoffs: (32) $\begin{array}{rl}& \mathrm{FxA}:max(\omega \left(A\right(S)-K_2),0),\frac{K_2}{S\left(t_0\right)}\in [0.8,1.2],\end{array}$(32) (33) $\begin{array}{rl}& \mathrm{FxL}:max(\omega \left(A\right(S)-K_2),0),\frac{K_2}{S\left(t_0\right)}\in [0.8,1.2],\end{array}$(33) whereas for the third, FxFlA, we have: (34) $max(\omega \left(A\right(S)-K_{1}S(T)-K_2),0),(\frac{K_1}{S\left(t_0\right)},\frac{K_2}{S\left(t_0\right)})\in [0.4,0.6]\times [0.4,0.6],$(34) with $A\left(S\right)$ defined as in (30(30) $A\left(S\right)=\frac{1}{5}\sum _{n=1}^{5}S\left(t_n\right),t_n:=T-(5-n){\tau }_{\mathtt{A}}\mathtt{,}\mathtt{n}\mathtt{=}\mathtt{1}\mathtt{,}\dots \mathtt{,}\mathtt{5}\mathtt{,}$(30) ) for FxA and FxFlA, and as in (31(31) $A\left(S\right)=\underset{n}{min}S\left(t_n\right),t_n:=T-(30-n){\tau }_{\mathtt{L}}\mathtt{,}\mathtt{n}\mathtt{=}\mathtt{1}\mathtt{,}\dots \mathtt{,}\mathtt{30}\mathtt{,}$(31) ) for FxL.

All the results in the following sections are compared to a MC benchmark obtained using the almost exact simulation described in Appendix A.3.

5.4.1. Numerical results for FxA

The procedure described in Section 3.2 is employed to solve the problem of pricing fixed-strike discrete Asian options with payoffs as in (32(32) $\begin{array}{rl}& \mathrm{FxA}:max(\omega \left(A\right(S)-K_2),0),\frac{K_2}{S\left(t_0\right)}\in [0.8,1.2],\end{array}$(32) ), with underlying stock price initial value $S\left(t_0\right)=1$. In this experiment, the ANN is trained on Heston model parameters' ranges, which include the examples proposed in [1] representing some real applications. Furthermore, we note the following aspect.

Remark 5.1

Scaled underlying process and (positive) homogeneity of A

The unit initial value is not restrictive. Indeed, the underlying stock price dynamics in Equations (9(9) $\begin{array}{rl}\mathrm{d}S\left(t\right)& =rS\left(t\right)\mathrm{d}t+\sqrt{v\left(t\right)}S\left(t\right)\mathrm{d}{W}_x\left(t\right),S\left(t_0\right)=S_0,\end{array}$(9) ) and (10(10) $\begin{array}{rl}\mathrm{d}v\left(t\right)& =\kappa (\overline{v}-v(t\left)\right)\mathrm{d}t+\gamma \sqrt{v\left(t\right)}\mathrm{d}{W}_v\left(t\right),v\left(t_0\right)=v_0,\end{array}$(10) ) are independent of $S\left(t_0\right)$, with the initial value only accounting as a multiplicative constant (this can be easily proved by means of Itô's lemma). Moreover, since $A\left(S\right)$ is (positive) homogeneous in S also $A\left(S\right)$ can be easily ‘scaled’ according to the desired initial value. Particularly, given the constant c>0, $\mathrm{cA}\left(S\right)\stackrel{\mathrm{d}}{=}A\left(\mathrm{cS}\right)$.

The methodology is tested on different randomly chosen sets of Heston parameters. We report the details for two specific sets, Set I and Set II, available in Table 6. For Set I, in Figure 5, we compare the population from $A\left(S\right)$ obtained employing SC with the MC benchmark (both with $N_{\mathtt{\text{paths}}}={10}^5$ paths each). Figure 5(a) shows the highly accurate approximation of the exact map $g=F_A^{-1}\circ F_{\xi }$ by means of the piecewise polynomial approximation $\tilde{g}$. As a consequence, both the PDF (see Figure 5(b)) and the CDF (see Figure 5(c)) perfectly match. Moreover, the methodology is employed to value fixed-strike arithmetic Asian options (calls and puts) for two sets of parameters (Set I and Set II) and 50 different strikes $K_2$. The resulting prices are reported in Figure 6(a,c). Figure 6(b,d) display the standard errors for MC and SC on the left y-axis, and the pricing error ${ϵ}_P$ for SC and SA (in units of the corresponding MC standard error) on the right y-axis. The pricing error tends to be more significant the more out of money the option is, due to the smaller SE.

The timing results are reported (in milliseconds) in Table 7, for different choices of $N_{\mathtt{\text{paths}}}$. Both SC and SA times only refer to the online pricing computational time, namely the time required for the pricing procedure, excluding the training sets generation and the ANNs training, which are performed – only once – offline. The semi-analytic formula requires a constant evaluation time, as well as the SC technique (if $N_{\mathtt{\text{paths}}}$ is fixed), whereas the MC simulation is dependent on the parameter T (since we decided to keep the same MC step in every simulation). Therefore, the methodology becomes more convenient the longer the maturity of the option T. The option pricing computational time is reduced by tens of times when using SC to generate the population from A, while it is reduced by hundreds of times if the semi-analytic (SA) formula is employed.

Table 6. Tested Heston parameter sets.

Set	r	κ	γ	ρ	$\overline{v}$	$v_0$	T
I	0.04	0.5	1.0	−0.8	0.08	0.05	1.0
Set	r	κ	γ	ρ	$\overline{v}$	$v_0$	T
II	0.02	1.0	0.9	−0.6	0.10	0.13	1.5
Set	r	κ	γ	ρ	$\overline{v}$	$v_0$	T
III	0.01	0.46	0.99	−0.79	0.09	0.11	1.0

Figure 5. Left: comparison between maps g from MC and $\tilde{g}$ from SC (with linear extrapolation). Center: comparison between MC histogram of $A\left(S\right)$ and the numerical PDF from SC. Right: comparison between MC numerical CDF of $A\left(S\right)$ and the numerical CDF from SC.

Figure 6. Left: fixed-strike discrete Asian option prices for 50 different strikes under the Heston model dynamics. Right: MC and SC standard errors (left y-axis) and pricing errors in units of the corresponding MC standard errors SE, obtained with $N_{\mathtt{\text{paths}}}={10}^5$ for both MC and SC (right y-axis). Up: call with Set I. Down: put with Set II.

Table 7. Timing results for option pricing.

	Fx Asian call (Set I)				Fx Asian put (Set II)				Fx-Fl Asian call (Set III)
	$N_{\mathtt{\text{paths}}}$	time	speed-up		$N_{\mathtt{\text{paths}}}$	time	speed-up		$N_{\mathtt{\text{paths}}}$	time	speed-up
MC	${10}^5$	1951 ms	–	MC	${10}^5$	3134 ms	–	MC	${10}^5$	1996 ms	–
SC	${10}^5$	24 ms	81	SC	${10}^5$	29 ms	107	SC	${10}^5$	69 ms	29
SA	–	15 ms	130	SA	–	16 ms	196	–	–	–	–
MC	$2\times {10}^5$	3905 ms	–	MC	$2\times {10}^5$	6386 ms	–	MC	$5\times {10}^4$	953 ms	–
SC	$2\times {10}^5$	58 ms	68	SC	$2\times {10}^5$	55 ms	116	SC	$5\times {10}^4$	43 ms	22
SA	–	15 ms	263	SA	–	15 ms	440	–	–	–	–

Eventually, the error distribution of 10,000 different option prices (one call and one put with 100 values for $K_2$ each for 50 randomly chosen Heston parameters' sets and maturities) is given in Figure 7(a). The SA prices (assuming a linear extrapolation) are compared with MC benchmarks. The outcome is satisfactory and shows the robustness of the methodology proposed. The error is smaller than three times the MC standard error in more than 90% of the cases when $N_{\mathtt{\text{paths}}}={10}^5$ (red histogram), and in more than 80% of the cases when $N_{\mathtt{\text{paths}}}=2\times {10}^5$ (blue histogram).

5.4.2. Numerical results for FxL

In this section, we use the procedure to efficiently value the pipeline risk typically embedded in mortgages. The pipeline risk (in mortgages) is one of the risks a financial institution is exposed to any time a client buys a mortgage. Indeed, when a client decides to buy a mortgage there is a grace period (from one to three months in The Netherlands), during which (s)he is allowed to pick the most convenient rate, namely the minimum.

Figure 7. Pricing error ${ϵ}_P$ distribution. The error is expressed in units of the corresponding Monte Carlo benchmark standard error (SE), and reported for two different numbers of MC paths. Left: FxA (10,000 values). Center: FxL (10,000 values). Right: FxFlA (54,000 values).

Observe now that a suitable Lookback option on a payer swap, namely a Lookback payer swaption, perfectly replicates the optionality offered to the client. In other words, the ‘cost’ of the pipeline risk is assessed by evaluating a proper Lookback swaption. In particular, we price fixed-strike discrete Lookback swaptions with a monitoring period of 3 months and 3-day frequency (see 31(31) $A\left(S\right)=\underset{n}{min}S\left(t_n\right),t_n:=T-(30-n){\tau }_{\mathtt{L}}\mathtt{,}\mathtt{n}\mathtt{=}\mathtt{1}\mathtt{,}\dots \mathtt{,}\mathtt{30}\mathtt{,}$(31) ).

We assume the underlying swap rate $S\left(t\right)$, $0\le t\le T$, is driven by the dynamics given in Equations (9(9) $\begin{array}{rl}\mathrm{d}S\left(t\right)& =rS\left(t\right)\mathrm{d}t+\sqrt{v\left(t\right)}S\left(t\right)\mathrm{d}{W}_x\left(t\right),S\left(t_0\right)=S_0,\end{array}$(9) ) and (10(10) $\begin{array}{rl}\mathrm{d}v\left(t\right)& =\kappa (\overline{v}-v(t\left)\right)\mathrm{d}t+\gamma \sqrt{v\left(t\right)}\mathrm{d}{W}_v\left(t\right),v\left(t_0\right)=v_0,\end{array}$(10) ) with $S\left(t_0\right)=0.05$ and parallel shifted of $\theta =0.03$. By introducing a shift, we handle also the possible situation of negative rates, which otherwise would require a different model specification.

Remark 5.2

Parallel shift of $S\left(t\right)$ and $A\left(S\right)$

A parallel shift θ does not affect the training set generation. Indeed, since $A\left(S\right)=\underset{n}{min}S\left(t_n\right)$, it holds $A(S-\theta )\stackrel{\mathrm{d}}{=}A\left(S\right)-\theta$. Then, it is enough to sample from $A\left(S\right)$ (built from the paths of $S\left(t\right)$ without shift) and perform the shift afterward, to get the desired distribution.

The timing results from the application of the procedure are comparable to the ones in Section 5.4.1 (see Table 7). Furthermore, in Figure 7(b), we report the pricing error distribution obtained by pricing call and put options for 50 randomly chosen Heston parameters' sets and 100 values for $K_2$. In this experiment, we observe that over 95% of the errors are within three MC SE when $N_{\mathtt{\text{paths}}}={10}^5$, and the percentage is about 90% when $N_{\mathtt{\text{paths}}}=2\times {10}^5$.

5.4.3. Numerical results for FxFlA

The third and last experiment consists in the conditional sampling of $A\left(S\right)|S\left(T\right)$. The samples are then used, together with $S\left(T\right)$, for pricing of fixed- and floating-strikes discrete Asian options.

The procedure is tested on 30 randomly chosen Heston parameters' sets. Both the MC benchmark and the SC procedure are based on populations with $N_{\mathtt{\text{paths}}}={10}^5$ paths. In the SC pricing, the process $S\left(T\right)$ is sampled using the COS method [24] combined with SC (COS-SC) to avoid a huge number of CDF numerical inversions [12], and so increase efficiency. Then, we apply the grid-based algorithm of Section 3.3. We evaluate the ANN at a reduced number of reference quantiles, and we compute the CVs corresponding to each sample of $S\left(T\right)$ by means of linear interpolation. The CVs identify the map $\tilde{g}$, which is employed for the conditional sampling. Figure 8(a) shows the cloud of points (for parameters' Set III in Table 6) of the bivariate distribution $\left(S\right(T),A(S\left)\right)$ generated using the procedure against the MC benchmark, while Figure 8(b) only focuses on the marginal distribution of $A\left(S\right)$. We can appreciate a good matching between the two distributions.

Figure 8. Left: joint distribution of $\left(S\right(T),A(S\left)\right)$, for the Heston set of parameters in Set III (see Table 6). Right: marginal distribution of $A\left(S\right)$, for the same set of Heston parameters.

For each set, we price $30\times 30$ call and put options with equally-spaced strikes $K_1$ and $K_2$ in the ranges of (34(34) $max(\omega \left(A\right(S)-K_{1}S(T)-K_2),0),(\frac{K_1}{S\left(t_0\right)},\frac{K_2}{S\left(t_0\right)})\in [0.4,0.6]\times [0.4,0.6],$(34) ). The results for the particular case of call options with Set III in Table 6 and $K_2=0.5$ are illustrated in Figure 9. The SC option prices and the corresponding MC benchmarks are plotted on the left. On the right, the standard errors for MC and SC are reported (left y-axis), and the absolute pricing error ${ϵ}_P$ is shown in units of the corresponding MC standard error (right y-axis). The timing results are reported in Table 7 for $N_{\mathtt{\text{paths}}}={10}^5$ and $N_{\mathtt{\text{paths}}}=5\times {10}^4$ and they keep into account the computational time for the pricing of all the 900 different call options (according to each combination of $K_1$ and $K_2$). Figure 7(c) displays the pricing error ${ϵ}_P$ distribution for the 30 randomly chosen Heston parameter sets (for each set 900 call and 900 put options are priced for every combination of $K_1$ and $K_2$, so the overall number of data is 54,000). About 92% of the errors are within three MC SE when $N_{\mathtt{\text{paths}}}=5\times {10}^4$. The percentage is around 80% when $N_{\mathtt{\text{paths}}}={10}^5$ are used.

Figure 9. Left: fixed-float-strike discrete Asian call option prices for 30 different $K_1$, $K_2=0.5$, and Heston parameter Set III. Right: MC and SC standard errors (left y-axis) and pricing errors in units of the corresponding MC standard errors SE, obtained with $N_{\mathtt{\text{paths}}}={10}^5$ for both MC and SC (right y-axis).

It might look surprising that the performance of the general procedure is better than the special one, but actually, it is not. Indeed, an important aspect needs to be accounted for. The high correlation between $S\left(T\right)$ and $A\left(S\right)$ makes the task of the ANN easier, in the sense that the distribution of $A\left(S\right)|S\left(T\right)$ typically has a low variance around $S\left(T\right)$. In other words, the ANN has to ‘learn’ only a small correction to get $A\left(S\right)|S\left(T\right)$ from $S\left(T\right)$ ($S\left(T\right)$ is an input for the ANN!), whereas the ANN in the special procedure learns the unconditional distribution of $A\left(S\right)$ with no information on the final value $S\left(T\right)$, and so only on the Heston parameters. The result is that a small lack in accuracy due to a not-perfect training process, or most likely to a not-perfect training set, is less significant in the conditional case rather than in the unconditional.

6. Conclusion

In this work, we presented a robust, data-driven procedure for the pricing of fixed- and floating-strike discrete Asian and Lookback options, under the stochastic volatility model of Heston. The usage of Stochastic Collocation techniques combined with deep artificial neural networks allows the methodology to reach a high level of accuracy while reducing the computational time by tens of times when compared to Monte Carlo benchmarks. Furthermore, we provide a semi-analytic pricing formula for European-type options with a payoff given by piecewise polynomial mapping of a standard normal random variable. Such a result allows to even increase the speed-up up to hundreds of times, without deterioration on the accuracy. An analysis of the error provides theoretical justification for the proposed scheme, and the problem of sampling from both unconditional and conditional distributions is further investigated from a numerical perspective. Finally, the numerical results provide clear evidence of the quality of the method.

Acknowledgments

We would like to thank the two anonymous reviewers whose insightful comments and constructive feedback greatly contributed to the enhancement of the quality of this article. Additionally, we extend our appreciation to Rabobank (the Netherlands) for funding this project.

Disclosure statement

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

Notes

1 The meaning of ‘data-driven’ here is the one given in [22]. The (empirical) distribution of interest is computed for a set of structural parameters and stored. Such synthetic ‘data’ are used to ‘drive’ the training of a suitable model.

2 We remark that in most real applications (with a few exceptions such as some commodities and FX rates) the correlation is negative. Furthermore, the phenomenon of ‘moment explosion’ for certain choices of Heston parameters involving positive correlation is discussed in [2].

3 Given the two random variables X, Y, with CDFs $F_X,F_Y$, it holds: $F_X\left(X\right)\stackrel{\mathrm{d}}{=}{F}_Y\left(Y\right)$.

4 The two main reasons for ξ being standard normal are the availability of such a distribution in most of the computing tools, and the ‘similarity’ between a standard normally r.v. and (the logarithm of) $A\left(S\right)$ (see [12] for more details).

5 A recursive formula for the computation of $m_i(a,b)$ is given in Appendix A.2.

6 The synthetic data are generated via MC simulation, as explained in Section 5.

7 By stochastic bridge, we mean any stochastic process conditional to both its initial and final values.

8 Th PDF $f_{\xi }$ works as a dumping factor in ${ϵ}_{-}=\mathbb{E}\left[\right(g{\mathbb{1}}_{{\mathrm{\Omega }}_{-}}-g_{-}{)}^2\left(\xi \right)]$ and ${ϵ}_{+}=\mathbb{E}\left[\right(g{\mathbb{1}}_{{\mathrm{\Omega }}_{+}}-g_{+}{)}^2\left(\xi \right)]$.

9 Entire functions are complex analytic functions in the whole complex plane $\mathbb{C}$.

10 $\mathbf{p}\mathrm{\setminus }\left\{T\right\}=\{r,\kappa ,\gamma ,\rho ,\overline{v},v_0\}$ and $\mathbf{p}\mathrm{\setminus }\left\{T\right\}=\{\kappa ,\gamma ,\rho ,\overline{v},v_0\}$ for FxA and FxL, respectively.

11 Under the risk-neutral measure, $\mathbb{Q}$. However, the same scheme applies under the underlying process measure ${\mathbb{Q}}^S$, with only a minor difference, i.e. $k_1:=(\mathrm{\rho \kappa }/\gamma +1/2)\mathrm{\Delta t}-\rho /\gamma$.

Appendix

Proofs and lemmas

A.1. Underlying process measure for floating-strike options

Proof of Proposition 2.1

Under the risk-neutral measure $\mathbb{Q}$ the value at time $t_0\ge 0$ of a floating-strike Asian Option, with maturity $T>t_0$, underlying process $S\left(t\right)$, and future monitoring dates $t_n$, $n\in \{1,\dots ,N\}$, is given by: $V_{\omega }^{\mathtt{fl}}\left(t_0\right)={\mathbb{E}}_{t_0}^{\mathbb{Q}}[\frac{M\left(t_0\right)}{M\left(T\right)}max(\omega \left(A\right(S)-K_{1}S(T\left)\right),0)].$We define a Radon-Nikodym derivative to change the measure from the risk-neutral measure $\mathbb{Q}$ to the stock measure ${\mathbb{Q}}^S$, namely the measure associated with the numéraire S: $\frac{\mathrm{d}{\mathbb{Q}}^S}{\mathrm{d}\mathbb{Q}}=\frac{S\left(T\right)}{S\left(t_0\right)}\frac{M\left(t_0\right)}{M\left(T\right)},$which yields the following present value, expressed as an expectation under the measure ${\mathbb{Q}}^S$: $\begin{array}{rl}{V}_{\omega }^{\mathtt{fl}}\left(t_0\right)& ={\mathbb{E}}_{t_0}^S[\frac{M\left(t_0\right)}{M\left(T\right)}max(\omega \left(A\right(S)-K_{1}S(T\left)\right),0)\frac{S\left(t_0\right)}{S\left(T\right)}\frac{M\left(T\right)}{M\left(t_0\right)}]\\ & =S\left(t_0\right){\mathbb{E}}_{t_0}^S[max(\omega (\frac{A\left(S\right)}{S\left(T\right)}-K_1),0)].\end{array}$

Proposition Appendix A.1

The Heston model under the underlying process measure

Using the same notation as in (9(9) $\begin{array}{rl}\mathrm{d}S\left(t\right)& =rS\left(t\right)\mathrm{d}t+\sqrt{v\left(t\right)}S\left(t\right)\mathrm{d}{W}_x\left(t\right),S\left(t_0\right)=S_0,\end{array}$(9) ) and (10(10) $\begin{array}{rl}\mathrm{d}v\left(t\right)& =\kappa (\overline{v}-v(t\left)\right)\mathrm{d}t+\gamma \sqrt{v\left(t\right)}\mathrm{d}{W}_v\left(t\right),v\left(t_0\right)=v_0,\end{array}$(10) ), under the stock $S\left(t\right)$ measure, ${\mathbb{Q}}^S$, the Heston framework yields the following dynamics for the process $S\left(t\right)$: $\begin{array}{rl}\mathrm{d}S\left(t\right)& =(r+v(t\left)\right)S\left(t\right)\mathrm{d}t+\sqrt{v\left(t\right)}S\left(t\right)\mathrm{d}{W}_x^S\left(t\right),S\left(t_0\right)=S_0,\\ \mathrm{d}v\left(t\right)& ={\kappa }^{\ast }({\overline{v}}^{\ast }-v(t\left)\right)\mathrm{d}t+\gamma \sqrt{v\left(t\right)}\mathrm{d}{W}_v^S\left(t\right),v\left(t_0\right)=v_0,\end{array}$with ${\kappa }^{\ast }=\kappa -\mathrm{\gamma \rho }$, ${\overline{v}}^{\ast }=\kappa \overline{v}/{\kappa }^{\ast }$, and $W_x^S\left(t\right)$ and $W_v^S\left(t\right)$ are BMs under the underlying process measure ${\mathbb{Q}}^S$ such that $\mathrm{d}{W}_s^S\left(t\right)\mathrm{d}{W}_v^S\left(t\right)=\rho \mathrm{d}t$.

Proof.

Under the stock measure ${\mathbb{Q}}^S$, implied by the stock $S\left(t\right)$ as numéraire, all the assets discounted with S must be martingales. Particularly, this entails that $M\left(t\right)/S\left(t\right)$ must be a martingale, where $M\left(t\right)$ is the money-savings account defined as $\mathrm{d}M\left(t\right)=\mathrm{rM}\left(t\right)\mathrm{d}t$.

From (9(9) $\begin{array}{rl}\mathrm{d}S\left(t\right)& =rS\left(t\right)\mathrm{d}t+\sqrt{v\left(t\right)}S\left(t\right)\mathrm{d}{W}_x\left(t\right),S\left(t_0\right)=S_0,\end{array}$(9) ) and (10(10) $\begin{array}{rl}\mathrm{d}v\left(t\right)& =\kappa (\overline{v}-v(t\left)\right)\mathrm{d}t+\gamma \sqrt{v\left(t\right)}\mathrm{d}{W}_v\left(t\right),v\left(t_0\right)=v_0,\end{array}$(10) ), using Cholesky decomposition, the Heston model can be expressed in terms of independent Brownian motions, ${\tilde{W}}_x\left(t\right)$ and ${\tilde{W}}_v\left(t\right)$, through the following system of SDEs: $\begin{array}{rl}\mathrm{d}S\left(t\right)& =\mathrm{rS}\left(t\right)\mathrm{d}t+\sqrt{v\left(t\right)}S\left(t\right)\mathrm{d}{\tilde{W}}_x\left(t\right),\\ \mathrm{d}v\left(t\right)& =\kappa (\overline{v}-v(t\left)\right)\mathrm{d}t+\gamma \sqrt{v\left(t\right)}\left[\rho \mathrm{d}{\tilde{W}}_x\right(t)+\sqrt{1-{\rho }^2}\mathrm{d}{\tilde{W}}_v(t\left)\right].\end{array}$After application of Itô's Lemma we find: $\mathrm{d}\frac{M\left(t\right)}{S\left(t\right)}=\frac{1}{S\left(t\right)}\mathrm{rM}\left(t\right)\mathrm{d}t-\frac{M\left(t\right)}{S^2\left(t\right)}\left(\mathrm{rS}\right(t)\mathrm{d}t+\sqrt{v\left(t\right)}S(t\left)\mathrm{d}{\tilde{W}}_x\right(t\left)\right)+\frac{M\left(t\right)}{S^3\left(t\right)}v\left(t\right)S^2\left(t\right)\mathrm{d}t,$which implies the following measure transformation: $\mathrm{d}{\tilde{W}}_x\left(t\right)=\mathrm{d}{\tilde{W}}_x^S\left(t\right)+\sqrt{v\left(t\right)}\mathrm{d}t.$Thus, under the stock measure ${\mathbb{Q}}^S$, the dynamics of $S\left(t\right)$ reads: $\begin{array}{rl}\mathrm{d}S\left(t\right)& =\mathrm{rS}\left(t\right)\mathrm{d}t+\sqrt{v\left(t\right)}S\left(t\right)\left(\mathrm{d}{\tilde{W}}_x^S\right(t)+\sqrt{v\left(t\right)}\mathrm{d}t)\\ & =(r+v(t\left)\right)S\left(t\right)\mathrm{d}t+\sqrt{v\left(t\right)}S\left(t\right)\mathrm{d}{\tilde{W}}_x^S\left(t\right),\end{array}$while for the dynamics of $v\left(t\right)$ we find: $\begin{array}{rl}\mathrm{d}v\left(t\right)& =\kappa (\overline{v}-v(t\left)\right)\mathrm{d}t+\gamma \sqrt{v\left(t\right)}[\rho \left(\mathrm{d}{\tilde{W}}_x^S\right(t)+\sqrt{v\left(t\right)}\mathrm{d}t)+\sqrt{1-{\rho }^2}\mathrm{d}{\tilde{W}}_v(t\left)\right]\\ & =[\kappa (\overline{v}-v(t\left)\right)+\mathrm{\gamma \rho v}(t\left)\right]\mathrm{d}t+\gamma \sqrt{v\left(t\right)}\left[\rho \mathrm{d}{\tilde{W}}_x^S\right(t)+\sqrt{1-{\rho }^2}\mathrm{d}{\tilde{W}}_v(t\left)\right].\end{array}$Setting ${\kappa }^{\ast }:=\kappa -\mathrm{\gamma \rho }$, ${\overline{v}}^{\ast }:=\kappa \overline{v}/{\kappa }^{\ast }$, $W_x^S\left(t\right):={\tilde{W}}_x^S\left(t\right)$, and $W_v^S\left(t\right):=\rho {\tilde{W}}_x^S\left(t\right)+\sqrt{1-{\rho }^2}{\tilde{W}}_v\left(t\right)$ the proof is complete.

A.2. Semi-analytic pricing formula

Result Appendix A.1

Moments of truncated standard normal distribution

Let $\xi \sim \mathcal{N}(0,1)$ and $a,b\in [-\mathrm{\infty },+\mathrm{\infty }]$, a<b. Then, the recursive expression for: $m_i(a,b):=\mathbb{E}[{\xi }^i|a\le \xi \le b],$the ith moment of the truncated standard normal distribution $\xi |a\le \xi \le b$, reads: $m_i(a,b)=(i-1)m_{i-2}(a,b)-\frac{b^{i-1}{f}_{\xi }\left(b\right)-a^{i-1}{f}_{\xi }\left(a\right)}{F_{\xi }\left(b\right)-F_{\xi }\left(a\right)},i\in \mathbb{N}\mathrm{\setminus }\left\{0\right\},$where $m_{-1}(a,b):=0$, $m_0(a,b):=1$, and $f_{\xi }$ and $F_{\xi }$ are the PDF and the CDF of ξ, respectively.

Result Appendix A.2

Expectation of polynomial of truncated normal distribution

Let $p\left(x\right)=\sum _{i=0}^{M-1}{\alpha }_i{x}^i$ be a $(M-1)$-degree polynomial and let $\xi \sim \mathcal{N}(0,1)$, with $f_{\xi }$, $F_{\xi }$ its PDF and CDF, respectively. Then, for any $a,b\in [-\mathrm{\infty },+\mathrm{\infty }]$ with a<b, the following holds: (A1) ${\int }_a^{b}p\left(x\right)f_{\xi }\left(x\right)\mathrm{d}x=\sum _{i=0}^{M-1}{\alpha }_i{m}_i(a,b)\left(F_{\xi }\right(b)-F_{\xi }(a\left)\right),$(A1) with $m_i(a,b)$ as defined in Result Appendix A.1.

Proof.

The proof immediately follows thanks to the following equalities: (A2) $\begin{array}{rl}{\int }_a^{b}p\left(x\right)f_{\xi }\left(x\right)\mathrm{d}x& =\sum _{i=0}^{M-1}{\alpha }_i{\int }_a^b{x}^i{f}_{\xi }\left(x\right)\mathrm{d}x=\sum _{i=0}^{M-1}{\alpha }_i\mathbb{E}\left[{\xi }^i{\mathbb{1}}_{[a,b]}\right(\xi \left)\right]\\ & =\sum _{i=0}^{M-1}{\alpha }_i\mathbb{E}[{\xi }^i|a\le \xi \le b]\cdot \mathbb{P}\left[a\le \xi \le b\right].\end{array}$(A2)

Proof of Proposition 3.2

The approximation $\tilde{g}$ is strictly increasing in the domain of interest. Then, setting $c_K={\tilde{g}}^{-1}\left(K\right)$, we have: (A3) $\begin{array}{rl}\frac{{\tilde{V}}_{\omega }\left(t_0\right)}{C}& ={\int }_{-\mathrm{\infty }}^{+\mathrm{\infty }}max(\omega \left(\tilde{g}\right(x)-K),0)f_{\xi }\left(x\right)\mathrm{d}x\\ & ={\int }_{\omega c_K}^{+\mathrm{\infty }}\omega \left(\tilde{g}\right(\mathrm{\omega y})-K)f_{\xi }\left(\mathrm{\omega y}\right)\mathrm{d}y\\ & =\omega \left({\int }_{\omega c_K}^{+\mathrm{\infty }}\tilde{g}\right(\mathrm{\omega y}\left)f_{\xi }\right(y)\mathrm{d}y-K\mathbb{P}[\xi >\omega c_K\left]\right),\end{array}$(A3) where the first equality holds by definition of expectation, the second one relies on a suitable change of variable (y = −x) and the last one holds thanks to the even symmetry of $f_{\xi }$. We define the integral $I_{\omega }\left(c_K\right)$ as: $I_{\omega }\left(c_K\right):={\int }_{\omega c_K}^{+\mathrm{\infty }}\tilde{g}\left(\mathrm{\omega x}\right)f_{\xi }\left(x\right)\mathrm{d}x,$and using the definition of $\tilde{g}$ as a piecewise polynomial, we get: (A4) $I_{\omega }\left(c_K\right)={\int }_{\omega c_K}^{-\overline{\xi }\vee c_K}{g}_{-\omega }\left(\mathrm{\omega x}\right)f_{\xi }\left(x\right)\mathrm{d}x+{\int }_{-\overline{\xi }\vee c_K}^{\overline{\xi }\vee c_K}{g}_M\left(\mathrm{\omega x}\right)f_{\xi }\left(x\right)\mathrm{d}x+{\int }_{\overline{\xi }\vee c_K}^{+\mathrm{\infty }}{g}_{\omega }\left(\mathrm{\omega x}\right)f_{\xi }\left(x\right)\mathrm{d}x.$(A4) The thesis follows by applying Result Appendix A.2 at each term in (A4(A4) $I_{\omega }\left(c_K\right)={\int }_{\omega c_K}^{-\overline{\xi }\vee c_K}{g}_{-\omega }\left(\mathrm{\omega x}\right)f_{\xi }\left(x\right)\mathrm{d}x+{\int }_{-\overline{\xi }\vee c_K}^{\overline{\xi }\vee c_K}{g}_M\left(\mathrm{\omega x}\right)f_{\xi }\left(x\right)\mathrm{d}x+{\int }_{\overline{\xi }\vee c_K}^{+\mathrm{\infty }}{g}_{\omega }\left(\mathrm{\omega x}\right)f_{\xi }\left(x\right)\mathrm{d}x.$(A4) ) and exploiting the definition of $F_{\xi }$.

A.3. Almost exact simulation from the Heston model

In a MC framework, the most common scheme employed in the industry is the Euler-Maruyama discretization of the system of SDEs which describes the underlying process dynamics. For the stochastic volatility model of Heston, such a scheme can be improved, allowing for an exact simulation of the variance process $v\left(t\right)$ (see (10(10) $\begin{array}{rl}\mathrm{d}v\left(t\right)& =\kappa (\overline{v}-v(t\left)\right)\mathrm{d}t+\gamma \sqrt{v\left(t\right)}\mathrm{d}{W}_v\left(t\right),v\left(t_0\right)=v_0,\end{array}$(10) )), as shown in [5]. This results in increased accuracy, and avoids numerical issues due to the theoretical non-negativity of the process $v\left(t\right)$, leading to the so-called almost exact simulation of the Heston model [24].

Result Appendix A.3

Almost exact simulation from the Heston model

Given $X\left(t\right):=\log S\left(t\right)$, its dynamics¹¹ between the consequent times $t_i$ and $t_{i+1}$ is discretized with the following scheme: (A5) $\begin{array}{rl}{x}_{i+1}& \approx x_i+k_0+k_1{v}_i+k_2{v}_{i+1}+\sqrt{k_3{v}_i}\xi ,\\ v_{i+1}& =\overline{c}{\chi }^2(\delta ,\overline{\kappa }{v}_i)\end{array}$(A5) with the quantities: $\mathrm{\Delta t}:=t_{i+1}-t_1,\delta :=\frac{4\kappa \overline{v}}{{\gamma }^2},\overline{c}:=\frac{\overline{v}}{\delta }(1-{\mathrm{e}}^{-\kappa \mathrm{\Delta t}}),\overline{\kappa }:={\overline{c}}^{-1}{\mathrm{e}}^{-\kappa \mathrm{\Delta t}},$the noncentral chi-squared random variable ${\chi }^2(\delta ,\eta )$ with δ degrees of freedom and non-centrality parameter η, and $\xi \sim \mathcal{N}(0,1)$. The remaining constants are defined as: $k_0:=(r-\frac{\rho }{\gamma }\kappa \overline{v})\mathrm{\Delta t},k_1:=(\frac{\mathrm{\rho \kappa }}{\gamma }-\frac{1}{2})\mathrm{\Delta t}-\frac{\rho }{\gamma },k_2:=\frac{\rho }{\gamma },k_3:=(1-{\rho }^2)\mathrm{\Delta t}.$

Derivation.

Given $X\left(t\right)=\log S\left(t\right)$, by applying Itô's Lemma and Cholesky decomposition on the dynamics in (9(9) $\begin{array}{rl}\mathrm{d}S\left(t\right)& =rS\left(t\right)\mathrm{d}t+\sqrt{v\left(t\right)}S\left(t\right)\mathrm{d}{W}_x\left(t\right),S\left(t_0\right)=S_0,\end{array}$(9) ) and (10(10) $\begin{array}{rl}\mathrm{d}v\left(t\right)& =\kappa (\overline{v}-v(t\left)\right)\mathrm{d}t+\gamma \sqrt{v\left(t\right)}\mathrm{d}{W}_v\left(t\right),v\left(t_0\right)=v_0,\end{array}$(10) ), we get: (A6) $\begin{array}{rl}\mathrm{d}X\left(t\right)& =(r-\frac{1}{2}v(t\left)\right)\mathrm{d}t+\sqrt{v\left(t\right)}\left[\rho \mathrm{d}{\tilde{W}}_v\right(t)+\sqrt{1-{\rho }^2}\mathrm{d}{\tilde{W}}_x(t\left)\right],\end{array}$(A6) (A7) $\begin{array}{rl}\mathrm{d}v\left(t\right)& =\kappa (\overline{v}-v(t\left)\right)\mathrm{d}t+\gamma \sqrt{v\left(t\right)}\mathrm{d}{\tilde{W}}_v\left(t\right),\end{array}$(A7) where ${\tilde{W}}_x\left(t\right)$ and ${\tilde{W}}_v\left(t\right)$ are independent BMs.

By integrating (A6(A6) $\begin{array}{rl}\mathrm{d}X\left(t\right)& =(r-\frac{1}{2}v(t\left)\right)\mathrm{d}t+\sqrt{v\left(t\right)}\left[\rho \mathrm{d}{\tilde{W}}_v\right(t)+\sqrt{1-{\rho }^2}\mathrm{d}{\tilde{W}}_x(t\left)\right],\end{array}$(A6) ) and (A7(A7) $\begin{array}{rl}\mathrm{d}v\left(t\right)& =\kappa (\overline{v}-v(t\left)\right)\mathrm{d}t+\gamma \sqrt{v\left(t\right)}\mathrm{d}{\tilde{W}}_v\left(t\right),\end{array}$(A7) ) in a the time interval $[t_i,t_{i+1}]$, the following discretization scheme is obtained: (A8) $\begin{array}{rl}{x}_{i+1}& =x_i+{\int }_{t_i}^{t_{i+1}}(r-\frac{1}{2}v(t\left)\right)\mathrm{d}t+\rho {\int }_{t_i}^{t_{i+1}}\sqrt{v\left(t\right)}\mathrm{d}{\tilde{W}}_v\left(t\right)+\sqrt{1-{\rho }^2}{\int }_{t_i}^{t_{i+1}}\sqrt{v\left(t\right)}\mathrm{d}{\tilde{W}}_x\left(t\right),\end{array}$(A8) (A9) $\begin{array}{rl}{v}_{i+1}& =v_i+\kappa {\int }_{t_i}^{t_{i+1}}(\hat{v}-v(t\left)\right)\mathrm{d}t+\gamma {\int }_{t_i}^{t_{i+1}}\sqrt{v\left(t\right)}\mathrm{d}{\tilde{W}}_v\left(t\right),\end{array}$(A9) where $x_i:=X\left(t_i\right)$, $x_{i+1}:=X\left(t_{i+1}\right)$, $v_i:=v\left(t_i\right)$, $v_{i+1}:=v\left(t_{i+1}\right)$.

Given $v_i$, the variance $v_{i+1}$ is distributed as a suitable scaled noncentral chi-squared distribution [24]. Therefore, we substitute ${\int }_{t_i}^{t_{i+1}}\sqrt{v\left(t\right)}\mathrm{d}{\tilde{W}}_v^{\lambda }\left(t\right)$ in (A8(A8) $\begin{array}{rl}{x}_{i+1}& =x_i+{\int }_{t_i}^{t_{i+1}}(r-\frac{1}{2}v(t\left)\right)\mathrm{d}t+\rho {\int }_{t_i}^{t_{i+1}}\sqrt{v\left(t\right)}\mathrm{d}{\tilde{W}}_v\left(t\right)+\sqrt{1-{\rho }^2}{\int }_{t_i}^{t_{i+1}}\sqrt{v\left(t\right)}\mathrm{d}{\tilde{W}}_x\left(t\right),\end{array}$(A8) ) using (A9(A9) $\begin{array}{rl}{v}_{i+1}& =v_i+\kappa {\int }_{t_i}^{t_{i+1}}(\hat{v}-v(t\left)\right)\mathrm{d}t+\gamma {\int }_{t_i}^{t_{i+1}}\sqrt{v\left(t\right)}\mathrm{d}{\tilde{W}}_v\left(t\right),\end{array}$(A9) ), ending up with: $\begin{array}{rl}{x}_{i+1}& =x_i+{\int }_{t_i}^{t_{i+1}}(r-\frac{1}{2}v(t\left)\right)\mathrm{d}t+\frac{\rho }{\gamma }(v_{i+1}-v_i-\kappa {\int }_{t_i}^{t_{i+1}}(\overline{v}-v(t\left)\right)\mathrm{d}t)+\sqrt{1-{\rho }^2}{\int }_{t_i}^{t_{i+1}}\sqrt{v\left(t\right)}\mathrm{d}{\tilde{W}}_x\left(t\right).\end{array}$We approximate the integrals in the expression above employing the left integration boundary values of the integrand, as in the Euler-Maruyama discretization scheme. The scheme (A5(A5) $\begin{array}{rl}{x}_{i+1}& \approx x_i+k_0+k_1{v}_i+k_2{v}_{i+1}+\sqrt{k_3{v}_i}\xi ,\\ v_{i+1}& =\overline{c}{\chi }^2(\delta ,\overline{\kappa }{v}_i)\end{array}$(A5) ) follows collecting the terms and employing the property ${\tilde{W}}_x\left(t_{i+1}\right)-{\tilde{W}}_x\left(t_i\right)\stackrel{\mathrm{d}}{=}\sqrt{\mathrm{\Delta t}}\xi$, with $\xi \sim \mathcal{N}(0,1)$ and $\mathrm{\Delta t}:=t_{i+1}-t_i$.

A.4. SC error analysis for Chebyshev interpolation

The two following lemmas are useful to show that the conditional complex ChF ${\varphi }_{A|\mathbf{V}}\left(z\right)=\mathbb{E}\left[{\mathrm{e}}^{\mathrm{izA}}|\mathbf{V}\right]$ is an analytic function of $z\in \mathbb{C}$. The first one provides the law of the conditional stock-price distribution, whereas the second one is meant to give algebraic bounds for the target function $A\left(S\right)$.

Lemma Appendix A.2

Conditional distribution under Heston

Let $S\left(t\right)$ be the solution at time t of Equation (9(9) $\begin{array}{rl}\mathrm{d}S\left(t\right)& =rS\left(t\right)\mathrm{d}t+\sqrt{v\left(t\right)}S\left(t\right)\mathrm{d}{W}_x\left(t\right),S\left(t_0\right)=S_0,\end{array}$(9) ) and $I_v\left(t\right):={\int }_{t_0}^{t}v\left(\tau \right)\mathrm{d}\tau$, with v driven by the dynamics in Equation (10(10) $\begin{array}{rl}\mathrm{d}v\left(t\right)& =\kappa (\overline{v}-v(t\left)\right)\mathrm{d}t+\gamma \sqrt{v\left(t\right)}\mathrm{d}{W}_v\left(t\right),v\left(t_0\right)=v_0,\end{array}$(10) ). Then, the following equality in distribution holds: $S\left(t\right)|I_v(t_0,t)\stackrel{\mathrm{d}}{=}\exp \left(\mu \right(I_v(t_0,t))+\sigma (I_v(t_0,t)\left)\xi \right),$with $\xi \sim \mathcal{N}(0,1)$, μ and σ defined as $\mu \left(y\right):=\log S\left(t_0\right)+r(t-t_0)-y/2$ and $\sigma \left(y\right):=\sqrt{y}$. Furthermore, for any $k=0,1,\dots$, the following holds: (A10) $\mathbb{E}\left[S\right(t{)}^k|I_v\left(t\right)]=\exp \left(\mathrm{k\mu }\right(I_v\left(t\right))+\frac{1}{2}{k}^2{\sigma }^2(I_v\left(t\right)\left)\right)$(A10) In other words, the stock price given the time-integral of the variance process $I_v$ is log-normally distributed, with parameters dependent on the time-integral $I_v$, and its moments up to any order are given as in Equation (A10(A10) $\mathbb{E}\left[S\right(t{)}^k|I_v\left(t\right)]=\exp \left(\mathrm{k\mu }\right(I_v\left(t\right))+\frac{1}{2}{k}^2{\sigma }^2(I_v\left(t\right)\left)\right)$(A10) ).

Proof.

Writing (9(9) $\begin{array}{rl}\mathrm{d}S\left(t\right)& =rS\left(t\right)\mathrm{d}t+\sqrt{v\left(t\right)}S\left(t\right)\mathrm{d}{W}_x\left(t\right),S\left(t_0\right)=S_0,\end{array}$(9) ) in integral form we get: $S\left(t\right)=S\left(t_0\right)\exp \left(r\right(t-t_0)-\frac{1}{2}{\int }_{t_0}^{t}v(\tau )\mathrm{d}\tau +{\int }_{t_0}^t\sqrt{v\left(\tau \right)}\mathrm{d}{W}_x(\tau \left)\right).$By considering the conditional distribution $S\left(t\right)|I_v\left(t\right)$ (instead of $S\left(t\right)$) the only source of randomness is given by the Itô's integral (and it is due to the presence of the Brownian motion $W_x\left(t\right)$). The thesis follows since the Itô's integral of a deterministic argument is normally distributed with zero mean and variance given by the time integral of the argument squared (in the same interval). Therefore, $S\left(t\right)|I_v\left(t\right)$ is log-normally distributed, with moments given as in (A10(A10) $\mathbb{E}\left[S\right(t{)}^k|I_v\left(t\right)]=\exp \left(\mathrm{k\mu }\right(I_v\left(t\right))+\frac{1}{2}{k}^2{\sigma }^2(I_v\left(t\right)\left)\right)$(A10) ).

Lemma Appendix A.3

Algebraic bounds

Let us consider $\{s_1,\dots ,s_N\}$, with $s_n>0$ for each $n=1,\dots ,N$. Then, for any $k=1,2,\dots$, we have:

(1)	$(\sum _n{s}_n{)}^k\le 2^{(N-1)(k-1)}\sum _n{s}_n^k$.
(2)	$(\underset{n}{min}{s}_n{)}^k\le s_{n^{\ast }}^k$ for any $n^{\ast }=1,\dots ,N$.

Proof.

The second thesis is obvious. We prove here the first one. We recall that in general, given a, b>0 and any $k=1,2,\dots$, the following inequality holds: (A11) $(a+b{)}^k\le 2^{k-1}(a^k+b^k).$(A11) Then, applying (A11(A11) $(a+b{)}^k\le 2^{k-1}(a^k+b^k).$(A11) ) N−1 times we get: ${\left(\sum _{n=1}^N{s}_n\right)}^k\le 2^{k-1}(s_1^k+{\left(\sum _{n=2}^N{s}_n\right)}^k)\le \cdots \le 2^{(N-1)(k-1)}{s}_N^k+\sum _{n=1}^{N-1}{2}^{n(k-1)}{s}_n^k,$which can be further bounded by: ${\left(\sum _{n=1}^N{s}_n\right)}^k\le 2^{(N-1)(k-1)}{s}_N^k+\sum _{n=1}^{N-1}{2}^{(N-1)(k-1)}{s}_n^k=\sum _{n=1}^N{2}^{(N-1)(k-1)}{s}_n^k.$

We have all the ingredients to prove Proposition 4.4.

Proof of Proposition 4.4

To exploit the characterization for entire ChFs in Theorem 4.3 we need to show the finiteness of each absolute moment as well as that Equation (25(25) $\underset{k\to +\mathrm{\infty }}{lim}{\left(\frac{|\mathbb{E}\left[Y^k\right]|}{k!}\right)}^{\frac{1}{k}}=0.$(25) ) is satisfied. Both the conditions can be proved using Lemmas Appendix A.2 and A.3. For $k=0,1,\dots$, we consider the two cases:

1. If $A=\frac{1}{N}\sum _{n}S\left(t_n\right)$, then thanks to Lemma Appendix A.3 we have: $\begin{array}{rl}\mathbb{E}\left[|A{|}^k|\mathbf{V}\right]=\frac{1}{N^k}\mathbb{E}\left[{\left(\sum _{n}S\right(t_n\left)\right)}^k|\mathbf{V}\right]& \le \frac{2^{(N-1)(k-1)}}{N^k}\mathbb{E}\left[\sum _{n}S\left(t_n{)}^k\right|\mathbf{V}\right]\\ & =\frac{2^{(N-1)(k-1)}}{N^k}\sum _n\mathbb{E}\left[S\right(t_n{)}^k|\mathbf{V}],\end{array}$whereas from Lemma Appendix A.2 follows: (A12) $\begin{array}{rl}\mathbb{E}\left[|A{|}^k|\mathbf{V}\right]& \le \frac{2^{(N-1)(k-1)}}{N^k}\sum _n\mathbb{E}\left[S\right(t_n{)}^k|\mathbf{V}]\end{array}$(A12) (A13) $\begin{array}{rl}& =\frac{2^{(N-1)(k-1)}}{N^k}\sum _n\exp \left(k{\mu }_n\right(\mathbf{V})+\frac{1}{2}{k}^2{\sigma }_n^2(\mathbf{V}\left)\right),\end{array}$(A13) where ${\mu }_n\left(\mathbf{V}\right):=\mu \left(I_v\right(t_n\left)\right)$ and ${\sigma }_n\left(\mathbf{V}\right):=\sigma \left(I_v\right(t_n\left)\right)$, $n=1,\dots ,N$.

2. If $A=\underset{n}{min}S\left(t_n\right)$, then we immediately have: (A14) $\begin{array}{rl}\mathbb{E}\left[|A{|}^k|\mathbf{V}\right]& \le \mathbb{E}\left[S\right(t_{n^{\ast }}{)}^k|\mathbf{V}]\end{array}$(A14) (A15) $\begin{array}{rl}& =\exp \left(k{\mu }_{n^{\ast }}\right(\mathbf{V})+\frac{1}{2}{k}^2{\sigma }_{n^{\ast }}^2(\mathbf{V}\left)\right),\end{array}$(A15) for an arbitrary $n^{\ast }=1,\dots ,N$.

The finiteness of the absolute moments up to any order follows directly from Equations (A13(A13) $\begin{array}{rl}& =\frac{2^{(N-1)(k-1)}}{N^k}\sum _n\exp \left(k{\mu }_n\right(\mathbf{V})+\frac{1}{2}{k}^2{\sigma }_n^2(\mathbf{V}\left)\right),\end{array}$(A13) ) and (A15(A15) $\begin{array}{rl}& =\exp \left(k{\mu }_{n^{\ast }}\right(\mathbf{V})+\frac{1}{2}{k}^2{\sigma }_{n^{\ast }}^2(\mathbf{V}\left)\right),\end{array}$(A15) ) respectively, since $I_v\left(t_n\right)$ are finite (indeed, they are time-integrals on compact intervals of continuous paths).

Eventually, thanks to Jensen's inequality we have $|\mathbb{E}\left[A|\mathbf{V}\right]{|}^k\le \mathbb{E}\left[|A{|}^k|\mathbf{V}\right]$. This, together with the at-most exponential growth (in k) of the absolute moments of $A|\mathbf{V}$, ensures that the limit in Equation (25(25) $\underset{k\to +\mathrm{\infty }}{lim}{\left(\frac{|\mathbb{E}\left[Y^k\right]|}{k!}\right)}^{\frac{1}{k}}=0.$(25) ) holds. Then, by Theorem 4.3, ${\varphi }_{A|\mathbf{V}}\left(z\right)$ is an entire function of the complex variable $z\in \mathbb{C}$.

Proof of Proposition 4.5

The goal here is to apply Morera's theorem. Hence, let $\gamma \in {\mathcal{S}}_{y^{\ast }}$ be any piecewise $C^1$ closed curve in the strip ${\mathcal{S}}_{y^{\ast }}$. Then: $\begin{array}{rl}{\int }_{\gamma }{\varphi }_A\left(z\right)\mathrm{d}z& \stackrel{\left(26\right)}{=}{\int }_{\gamma }{\int }_{{\mathrm{\Omega }}_{\mathbf{V}}}{\varphi }_{A|\mathbf{V}=\mathbf{v}}\left(z\right)\mathrm{d}{F}_{\mathbf{V}}\left(\mathbf{v}\right)\mathrm{d}z\\ & \stackrel{\mathrm{Fubini}}{=}{\int }_{{\mathrm{\Omega }}_{\mathbf{V}}}{\int }_{\gamma }{\varphi }_{A|\mathbf{V}=\mathbf{v}}\left(z\right)\mathrm{d}z\mathrm{d}{F}_{\mathbf{V}}\left(\mathbf{v}\right)\stackrel{\mathrm{Cauchy}}{=}0,\end{array}$where in the first equality we exploited the representation of the unconditional ChF ${\varphi }_A$ in terms of conditional ChFs ${\varphi }_{A|\mathbf{V}}$, in the second equality we use Fubini's theorem to exchange the order of integration, and eventually in the last equation we employ the Cauchy's integral theorem on ${\int }_{\gamma }{\varphi }_{A|\mathbf{V}=\mathbf{v}}\left(z\right)\mathrm{d}z$.