Machine Learning for Pricing American Options in High-Dimensional Markovian and non-Markovian models

Introduction

Pricing American options is clearly a crucial question of finance but also a challenging one since computing the optimal exercise strategy is not an evident task. This issue is even more exacting when the underling of the option is a multi-dimensional process, such as a baskets of $d$ assets, since in this case the direct application of standard numerical schemes, such as finite difference or tree methods, is not possible because of the exponential growth of the calculation time and the required working memory.

Common approaches in this field can be divided in four groups: techniques which rely on recombinant trees to discretize the underlyings (see , and ), techniques which employ regression on a truncated basis of $L^{2}$ in order to compute the conditional expectations (see and ), techniques which exploit Malliavin calculus to obtain representation formulas for the conditional expectation (see , , , and ) and techniques which make use of duality-based approaches for Bermudan option pricing (see , and ).

Recently, Machine Learning algorithms (Rasmussen and Williams ) and Deep Learning techniques (Nielsen ) have found great application in this sector of option pricing.

Neural networks are used by Kohler et al. to price American options based on several underlyings. Deep Learning techniques are nowadays widely used in solving large differential equations, which is intimately related to option pricing. In particular, Han et al. introduce a Deep Learning-based approach that can handle general high-dimensional parabolic PDEs. E et al. propose an algorithm for solving parabolic partial differential equations and backward stochastic differential equations in high dimension. Beck et al. introduce a method for solving high-dimensional fully nonlinear second-order PDEs. As far as American options in high dimension are concerned, Becker et al. develop a Deep Learning method for optimal stopping problems which directly learns the optimal stopping rule from Monte Carlo samples.

Also Machine Learning techniques have made their contribution. For example, Dixon and Crépey present a multi-Gaussian process regression for estimating portfolio risk, and in particular the associated CVA. De Spiegeleer et al. propose to apply Gaussian Process Regression (GPR) to predict the price of the derivatives from a training set made of observed prices for particular combinations of model parameters. Ludkovski proposes to use GPR meta-models for fitting the continuation values of Bermudan options. Similarly, GoudenÚge et al. propose the GPR-MC, which is a backward induction algorithm that employs Monte Carlo simulations and GPR to compute the price of American options in very high dimension (up to 100). In the insurance context, Gan studies the pricing of a large portfolio of Variable Annuities in the Black-Scholes model by using clustering and GPR. Moreover, Gan and Lin propose a novel approach that combines clustering technique and GPR to efficiently evaluate policies considering nested simulations.

In this paper we present two numerical techniques which upgrade the GPR-MC approach by replacing the Monte Carlo based computation of the continuation value respectively with a tree step and with an exact integration step. In particular, the algorithms we propose proceed backward over time and compute the price function only on a set of predetermined points. At each time step, a binomial tree step or a closed formula for integration are used together with GPR to approximate the continuation value at these points. The option price is then obtained as the maximum between the continuation value and the intrinsic value of the option and the algorithms proceed backward. For the sake of simplicity, we name these new approaches Gaussian Process Regression - Tree (GPR-Tree) and Gaussian Process Regression - Exact Integration (GPR-EI). We observe that the use of the GPR method to extrapolate the option value is particularly efficient in terms of computing time with respect to other techniques such as Neural Networks, especially because a small dataset is considered here. Moreover, Le Gratiet et Garnier developed recent convergence results about GPR, extending the outcomes of Rasmussen and Williams , and founding the convergence rate when different kernels are employed.

In order to demonstrate the wide applicability of the GPR methods, we also consider the rough Bergomi model, which is a non-Markovian model with stochastic volatility. Such a model, introduced by Bayer et al. stood out for explaining implied volatility smiles and other phenomena in the pricing of European options. The non-Markovian property of the model makes it difficult to implement a methodologically correct approach to address the valuation of American options. The literature in this framework is really poor. Horvat et al. propose an approach based on Donsker’s approximation for fractional Brownian motion and on a tree with exponential complexity. More recently, Bayer et al. introduce a method based on Monte Carlo simulation and exercise Rate Optimization.

Numerical results show that both the GPR-Tree and the GPR-EI methods are accurate and reliable in the multi-dimensional Black-Scholes model. Moreover the computational times with respect to the GPR-MC method are improved. The GPR-Tree and the GPR-EI methods prove its accuracy also when applied to the rough Bergomi model.

The reminder of the paper is organized as follows. Section 2 presents American options in the multi-dimensional Black-Scholes model. Section 3 and Section 4 introduce the GPR-Tree and the GPR-EI methods for the multi-dimensional Black-Scholes model respectively. Section 5 presents the American options in the rough Bergomi model. Section 6 and Section 7 introduce the GPR-Tree and the GPR-EI methods for the rough Bergomi model. Section 8 reports some numerical results. Section 9 draws some conclusions.

American options in the multi-dimensional Black-Scholes model

An American option with maturity $T$ is a derivative instrument whose holder can exercise the intrinsic optionality at any moment before maturity. Let $\mathbf{S}=(\mathbf{S}_{t})_{t\in[0,T]}$ denote the $d$-dimensional underlying process, which is supposed to randomly evolve according to the multi-dimensional Black-Scholes model: under the risk neutral probability, such a model is given by the following equation

For simulation purposes, the $d-$dimensional Black-Scholes model can be written alternatively using the Cholesky decomposition. Specifically, for $i=1,\dots,d$ we can write

where $\mathbf{B}$ is a $d$-dimensional uncorrelated Brownian motion and $\Sigma_{i}$ is the $i$-th row of the matrix $\Sigma$ defined as a square root of the correlation matrix $\Gamma$, given by

The GPR-Tree method in the multi-dimensional Black-Scholes model

The GPR-Tree method is similar to the GPR-MC method but the diffusion of the underlyings is performed through a step of a binomial tree. In particular, the algorithm proceeds backward over time, approximating the price of the American option with the price of a Bermudan option on the same basket. At each time step, the price function is evaluated only on a set of predetermined points, through a binomial tree step together with GPR to approximate the continuation value. Finally, the optionality is exploited by computing the option value as the maximum between the continuation value and the exercise value.

Let $N$ denote the number of time steps, $\Delta t=T/N$ be the time increment and $t_{n}=n\,\Delta t$ represent the discrete exercise dates for $n=0,1,\ldots,N$. At any exercise date $t_{n}$, the value of the option is determined by the vector of the underlying prices $\mathbf{S}_{t_{n}}$ as follows:

where $C$ denotes the continuation value of the option and it is given by the following relation:

We observe that, if the function $v\left(t_{n+1},\cdot\right)$ is known, then it is possible to compute $v\left(t_{n},\cdot\right)$ by approximating the expectation in (3.2). In order to obtain such an approximation, we consider a set $X$ of $P$ points whose elements represent certain possible values for the underlyings $\mathbf{S}$:

The GPR-Tree approximation $v_{N-2}^{GPR-Tree}\left(\cdot\right)$ of the value function $v\left(t_{N-2},\cdot\right)$ at time $t_{N-2}$ can be computed as follows:

Then, the function $v_{n}^{GPR-Tree}$ can be obtained as

The GPR-EI method in the multi-dimensional Black-Scholes model

The GPR-EI method differs from both the GPR-MC and GPR-Tree methods for two reasons. First of all, the predictors employed in the GPR step are related to the logarithms of the predictors used in the GPR-Tree method. Secondly, the continuation value at these points is computed through a closed formula which comes from an exact integration.

where $\omega_{1},\dots,\omega_{P}$ are weights that are computed by solving a linear system. The continuation value can be computed by integrating the function $u^{GPR}$ against a $d$-dimensional probability density. This calculation can be done easily by means of a closed formula.

Specifically, the GPR-EI method relies on the following Proposition.

Let $n\in\left\{0,\dots,N-1\right\}$ and suppose the function $u\left(t_{n+1},\cdot\right)$ at time $t_{n+1}$ to be known at $Z$. The GPR-EI approximation of the option value $u\left(t_{n},\cdot\right)$ at time $t_{n}$ at $\mathbf{z}^{p}$ is given by

$\sigma_{f}$, $\sigma_{l}$, and $\omega_{1},\dots,\omega_{P}$ are certain constants determined by the GPR approximation of the function $\mathbf{z}\mapsto u\left(t_{n+1},\mathbf{z}\right)$ for $k=1,\dots,P$, considering $Z$ as the predictor set, and $\Pi=\left(\Pi_{i,j}\right)$ is the $d\times d$ covariance matrix of the log-increments defined by $\Pi_{i,j}=\rho_{i,j}\sigma_{i}\sigma_{j}\Delta t$.

The proof of Proposition 1 is reported in the Appendix A. Equation (4.5) allows one to compute the option price at time $t=0$ by proceeding backward. In fact, the function $u\left(t_{N},\cdot\right)$ is known at time $t_{N}=T$ throught (4.2) since the price function $v\left(t_{N},\cdot\right)$ is equal to the payoff function $\Psi\left(\cdot\right)$. Moreover, if an approximation of $u\left(t_{n+1},\cdot\right)$ is available, then one can approximate $u\left(t_{n},\cdot\right)$ at $Z$ by means of relation (4.5). Finally, the option price at time $t=0$ is approximated by $u_{0}^{GPR-EI}\left(\log\left(\mathbf{\mathbf{S}_{0}}\right)\right)$.

American options in the rough Bergomi model

The rough Bergomi model, introduced by Bayer et al. , shapes the underlying process $S_{t}$ and its volatility $V_{t}$ through the following relations:

with $r$ the (constant) interest rate, $\eta$ a positive parameter and $H\in\left(0,1\right)$ the Hurst parameter. The deterministic function $\xi_{0}\left(t\right)$ represents the forward variance curve and following Bayer et al. we consider it as constant. The process $W_{t}^{1}$ is a Brownian motion, whereas $\widetilde{W}_{t}^{H}$ is a Riemann-Liouville fractional Brownian motion that can be expressed as a stochastic integral:

with $W_{t}^{2}$ a Brownian motion and $\rho$ the instantaneous correlation coefficient between $W_{t}^{1}$ and $W_{t}^{2}$.

The rough Bergomi model stood out for its ability to explain implied volatility and other phenomena related to European options. Moreover, it is particularly interesting from a computational point of view as it is a non-Markovian model and therefore it is not possible to apply standard techniques for American options.

where $\mathcal{F}_{t}$ is the natural filtration generated by the couple $\left(W_{s}^{1},\widetilde{W}_{s}^{H}\right)$ for $s\in\left[0,t\right]$. We point out that, as opposed to the multi-dimensional Brownian motion, in this case, the stopping time $\tau$ does not only depend from the actual values of $S$ and $V$ but, since these are non-Markovian processes, it depends on the whole filtration, that is from the whole observed history of the processes.

The GPR-Tree method in the rough Bergomi model

The GPR-Tree method can be adapted to price American options in the rough Bergomi model. Despite the dimension of the model is only two, it is a non-Markovian model which obliges one to take into account the past history when evaluating the price of an option. So, the price of an option at a certain moment depends on all the filtration at that moment. Clearly, evaluating an option by considering the whole history of the process (a continuous process) is not possible. To overcome such an issue, we simulate the process on a finite number of dates and we consider the sub-filtration induced by these observations. First of all, we consider a finite number $N$ of time steps that determines the time increment $\Delta t=\frac{T}{N}$, and we employ the scheme presented in Bayer et. al to generate a set of $P$ simulations of the couple $\left(S_{t},V_{t}\right)$ at $t_{n}=n\,\Delta t$ for $n=1,\ldots,N$. In particular, if we set $\Delta W_{n}^{1}=W_{t_{n}}^{1}-W_{t_{n-1}}^{1}$, then the $2N$-dimensional random vector $\mathbf{R}$, given by

follows a zero-mean Gaussian distribution. Moreover, using the relations stated in Appendix B, one can calculate the covariance matrix $\Upsilon$ of $\mathbf{R}$ and its lower triangular square root $\Lambda$ by using the Cholesky factorization. The vector $\mathbf{R}$ can be simulated by computing $\Lambda\mathbf{G}$, where $\mathbf{G}=\left(G_{1},\dots,G_{2N}\right)^{\top}$ is a vector of independent standard Gaussian random variables. Finally, a simulation for $\left(S_{t_{n}},V_{t_{n}}\right)_{n=0,\dots,N}$ can be obtained from $\mathbf{R}$ by considering the initial values

First of all, the GPR-Tree method simulates $P$ different samples for the vector $\mathbf{G}$, namely $\mathbf{G}^{p}$ for $p=1,\dots,P$, and it computes the corresponding paths $\left(S_{t_{1}}^{p},V_{t_{1}}^{p},\dots,S_{t_{N}}^{p},V_{t_{N}}^{p}\right)$ according to (6.2), (6.3) and (6.4). To summarize the values assumed by $S$ and $V$, let us define the vector

for $i,j\in\left\{0,\dots,N\right\}$ and $i<j$. Moreover, we also define

where $\log$ stands for the natural logarithm.

Then, the GPR-Tree method computes the option value for each of these $P$ trajectories, proceeding backward in time and considering the past history coded into the filtration. Since we consider only a finite number of steps, we approximate the filtration $\mathcal{F}_{t_{n}}$ with the natural filtration $\hat{\mathcal{F}_{t_{n}}}$ generated by the $2n$ variables $W_{t_{1}}^{1},\widetilde{W}_{t_{1}}^{H},\dots,W_{t_{n}}^{1},\widetilde{W}_{t_{n}}^{H}$. Moreover, $\hat{\mathcal{F}}_{t_{n}}$ is equal to the filtration generated by $S_{t_{1}},V_{t_{1}},\dots,S_{t_{n}},V_{t_{n}}$ because there exists a deterministic bijective function that allows one to obtain $W_{t_{1}}^{1},\widetilde{W}_{t_{1}}^{H},\dots,W_{t_{n}}^{1},\widetilde{W}_{t_{n}}^{H}$ from $S_{t_{1}},V_{t_{1}},\dots,S_{t_{n}},V_{t_{n}}$ and vice versa. Therefore, when we calculate the option value conditioned by filtration $\hat{\mathcal{F}}_{t_{n}}$, it is enough to conditioning with respect to the knowledge of the variables $S_{t_{1}},V_{t_{1}},\dots,S_{t_{n}},V_{t_{n}}$.

The GPR-Tree method proceeds backward in time, using a tree method and the GPR to calculate the option price with respect to the initially simulated trajectories. As opposed to the multi-dimensional Black-Scholes model, here we perform more than one single tree step, so as to reduce the number of GPR regressions and thus increasing the computational efficiency. In particular, we consider $N=N^{Tree}\cdot m$ with $N^{Tree}$ and $m$ natural numbers that represent how many times the tree method is used and the number of time steps employed, respectively.

After simulating the $P$ random paths $\left\{\mathbf{SV}_{1:N}^{p},\ p=1,\dots,P\right\}$, we compute the tree approximation of the option value $v\left(t_{N-m},\mathbf{SV}_{1:\left(N-m\right)}^{p}\right)$ at time $t_{N-m}$ for each path as follows:

with $C_{N-m}^{Tree}$ stands for the the approximation of the continuation value function at time $t_{N-m}$ obtained by means of a tree approach, which discretizes each component of the Gaussian vector $\mathbf{G}_{\left[2\left(N-m\right)+1\right]:2N}$ that generates the process. As opposed to the multi-dimensional Black-Scholes model, the approximation of the independent Gaussian components of $\mathbf{G}$ through the equiprobable couple $\left\{-1,+1\right\}$ is not suitable since the convergence to the right price is too slow. So, we propose to use the same discrete approximation employed by Alfonsi in , which is stated in the following Lemma.

So, for each path $p$, we consider a quadrinomial tree with $m$ time steps, and we use it to compute the continuation value. In particular, we consider the discrete time process $\left(\hat{S}_{k}^{p},\hat{V}_{k}^{p}\right)_{k\in\left\{N-m,\dots,N\right\}}$defined through

where $\Lambda_{2k+1}$ is the $2k+1$-th rows of the matrix $\Lambda$ and $\Lambda_{2k+2}$ the $2k+2$-th row. Moreover, $\hat{G}_{j}^{p}=G_{j}^{p}$ for $j=1,\dots,2\left(N-m\right)$ and the other components, that is $\hat{G}_{j}^{p}$ for $j=2\left(N-m\right)+1,\dots,2N$, are sampled by using the random variable $A$ of Lemma 2.

An option value is assigned to each node of the tree: at maturity, that is for $k=N,$ it is equal to the payoff $\Psi\left(\hat{S}_{N}^{p}\right)$, and for $k=N-m,\dots,N-1$ it can be obtained as the maximum between the exercise value and the discounted mean value at the future nodes, weighted according to the transition probabilities determined by the probability distribution of $A$.

This approach allows us to compute the function $v_{N-m}^{GPR-Tree}\left(\mathbf{SV}_{1:\left(N-m\right)}^{p}\right)$ for $p=1,\dots,P$. We point out that, since the quadrinomial tree is not recombinant, the number of nodes grows exponentially with the number of time steps $m$. Therefore, $m$ must be small. A similar problem arises with the tree approach proposed by Horvat et al. . In order to overcome such an issue, we apply the GPR method to approximate the function $u_{N-m}^{GPR-Tree}\left(\log\left(\mathbf{SV}_{1:\left(N-m\right)}^{p}\right)\right)=v_{N-m}^{GPR-Tree}\left(\mathbf{SV}_{1:\left(N-m\right)}^{p}\right)$. Specifically, consider a natural number $J$ and define $d_{n}=2\min\left(n,J+1\right)$. We train the GPR method considering the predictor set given by

We term $u_{N-m}^{GPR}$ the function obtained by the aforementioned regression, which depends on $\log\left(\mathbf{SV}_{\max\left\{1,N-m-J\right\}:N-m}^{p}\right)$. We stress out that if we consider $J=N-m-1$ (or greater), then the function $u_{N-m}^{GPR}$ would consider all the observed values of $S$ and $V$ as predictors. Anyway, numerical tests show that it is enough to consider smaller values of $J$, which reduces the dimension $d_{N-m}$ of the regression and thus improves the numerical efficiency. A similar approach is taken by Bayer et al. .

Once we have obtained $u_{N-m}^{GPR}$, we can approximate the option value $v\left(t_{N-2m},\mathbf{SV}_{1:\left(N-2m\right)}^{p}\right)$ at time $t_{N-2m}$ by means of the tree approach again. The only difference in this case is that the value attributed to the terminal nodes is not determined by the payoff function, but through the function $u_{N-m}^{GPR}$. We term $v_{N-2m}^{GPR-Tree}$ the function obtained after this backward tree step. If we train the GPR method considering the predictor set given by

then we obtain the function $u_{N-2m}^{GPR}$, which can be employed to repeat the tree step and the GPR step, proceeding backward up to obtaining the initial option price by backward induction.

The GPR-EI method in the rough Bergomi model

The GPR-EI method can be adapted to price American options in the rough Bergomi model. Just like the GPR-Tree approach, the GPR-EI method starts by simulating $P$ different paths $\left(S_{t_{1}}^{p},V_{t_{1}}^{p},\dots,S_{t_{N}}^{p},V_{t_{N}}^{p}\right)$ for the processes $S$ and $V$, and it goes on by solving a backward induction problem, through the use of the GPR method and a closed formula for integration.

As opposed to the multi-dimensional Black-Scholes model, in the rough Bergomi case the use of the squared exponential kernel is not suitable because it is a isotropic kernel and the predictors employed have different nature (prices and volatilities at different times) and thus changes in each predictor impact differently on the price. So, we employ the Automatic Relevance Determination (ARD) Squared Exponential Kernel, that has separate length scale for each predictor and it is given by

with $d$ the number of the considered predictors. Specifically, the GPR-EI method relies on the following Propositions.

The GPR-EI approximation of the option value at time $t_{N-1}$ at $\mathbf{SV}_{\max\left\{1,N-1-J\right\}:\left(N-1\right)}^{p}$ is given by:

where $\sigma_{f}$, $\sigma_{j}$, and $\omega_{1},\dots,\omega_{P}$ are certain constants determined by the GPR approximation of the function $\log\left(S_{T}\right)\mapsto\Psi\left(S_{T}\right)$. Moreover,

The proof of Proposition 3 is reported in the Appendix C. Therefore, we can compute the value of the option at time $t_{N-1}$ for each simulated path by using (7.2).

Let $n\in\left\{0,\dots,N-2\right\}$ and suppose the option price function $v\left(t_{n+1},\cdot\right)$ at time $t_{n+1}$ to be known for all the simulated paths $\left\{\mathbf{SV}_{1:N}^{p},p=1,\dots,P\right\}$. Define

where $\Lambda_{2n+2}$ is the $2n+2$-th row of the matrix $\Lambda$ and $\underline{\mathbf{G}}^{p}=\left(G_{1}^{p},\dots,G_{2n}^{p},0\dots,0\right)^{\top}$, and

where $\sigma_{d_{n+1}-1}$, $\sigma_{d_{n+1}}$, $\sigma_{f}$ and $\omega_{1},\dots,\omega_{P}$ are certain constants determined by the GPR approximation of the function $\log\left(\mathbf{SV}_{1:n+1}\right)\mapsto v\left(t_{n+1},\mathbf{SV}_{1:n+1}\right)$ considering $\left\{\mathbf{SV}_{\max\left\{1,n+1-J\right\}:n+1}^{p},p=1,\dots,P\right\}$ as the predictor set. Moreover, $h_{q}^{p}$ and $f_{q}^{p}$ are two factors given by

where $z_{i}^{p}=\log\left(S_{n+1-\left(i-1\right)/2}^{p}\right)$ if $i$ is even and $z_{i}^{p}=\log\left(V_{n+1-i/2}^{p}\right)$ if $i$ is odd, for $i=1,\dots,d_{n+1}$.

The proof of Proposition 4 is reported in the Appendix D. Relations (7.2) and (7.7) can be used to compute the option price at time $t=0$ by backward induction.

Numerical results

In this Section we present some numerical results about the effectiveness of the proposed algorithms. The first section is devoted to the numerical tests about the multi-dimensional Black-Scholes model, while the second is devoted to the rough Bergomi model. The algorithms have been implemented in MATLAB and computations have been preformed on a server which employs a $2.40$ GHz Intel® Xenon® processor (Gold 6148, Skylake) and 20 GB of RAM.

Following GoudenÚge et al. , we consider an Arithmetic basket Put, a Geometric basket Put and a Call on the Maximum of $d$-assets.

In particular, we use the following parameters $T=1$, $S_{0}^{i}=100$, $K=100$, $r=0.05$, constant volatilities $\sigma_{i}=0.2$, constant correlations $\rho_{ij}=0.2$ and $N=10$ exercise dates. Moreover, we consider $P=250,\ 500$ or $1000$ points. As opposed to the other input parameters, we vary the dimension $d$, considering $d=2,\,5,\,10,\,20,\,40$ and $100$.

We present now the numerical results obtained with the GPR-Tree and the GPR-EI methods for the three payoff examples.

Geometric basket Put is a particularly interesting option since it is possible to reduce the problem of pricing it in the $d$-dimensional model to a one dimensional American Put option in the Black-Scholes model which can be priced straightforwardly, for example using the CRR algorithm with $1000$ steps (see Cox et al. ). Therefore, in this case, we have a reliable benchmark to test the proposed methods. Moreover, when $d$ is smaller than $10$ we can also compute the price by means of a multi-dimensional binomial tree (see Ekvall ). In particular, the number of steps employed for the multi-dimensional binomial tree is equal to $200$ when $d=2$ and to $50$ when $d=5$. For values of $d$ larger than $5$, prices cannot be approximated via such a tree, because the memory required for the calculations would be too large. Furthermore, we also report the prices obtained with the GPR-MC method, employing $P=1000$ points and $M=10^{5}$ Monte Carlo simulations, for comparison purposes. As far as the GPR-Tree is concerned, we compute the prices only for the values of $d$ smaller than $40$ since for higher values of $d$ the tree step becomes over time demanding. In fact, the computation of the continuation value with the tree step grows exponentially with the dimension $d$ and for $d=40$ it requires the evaluation of the GPR approximation at $2^{40}\approx 10^{12}$ points for every times step and for every point of $X$.

Results are reported in Table 1. We observe that the two proposed methods provide accurate and stable results and the computational time is generally very small, except for the GPR-Tree method at $d=20$. Moreover, the computer processing time of the GRP-EI method increases little with the size of the problem and this makes the method particularly effective when the dimension of the problem is high. This is because the computation of the expected value and the training of the GPR model are minimally affected by the dimension of the problem.

Figure 8.1 investigates the convergence of the GPR methods changing the dimension $d$. As we can see, the relative error is small with all the considered methods, but the computational time required by the GPR-Tree method and the GPR-EI method is generally smaller with respect to the GPR-MC method.

1.2 Arithmetic basket Put option

As opposed to the Geometric basket Put option, in this case we have no method to obtain a fully reliable benchmark. Therefore we only consider the prices obtained by means of the GPR-MC method, employed with $P=1000$ points and $M=10^{5}$ Monte Carlo simulations. Moreover, for small values of $d$, a benchmark can be obtained by means of a multi-dimensional tree method (see Boyle et al. ), just as shown for the Geometric case. Results are reported in Table 2. Similarly to the Geometric basket Put, the prices obtained are reliable and they do not change much with respect to the number $P$ of points. As opposed to the GPR-Tree method, which can not be applied for high values of $d$, the GPR-EI method requires a small computational time for all the values concerned of $d$.

1.3 Call on the Maximum

As for the Arithmetic basket Put, in this case we have no numerical methods to obtain a fully reliable benchmark. However, for small values of $d$, we can approximate the price obtained by means of a multi-dimensional tree method. Moreover, we also consider the price obtained with the GPR-MC method. Results, which are shown in Table 3, have an accuracy comparable to the one obtained for the Arithmetic basket Put option.

2 Rough Bergomi model

Following Bayer et al. , we consider an American Put option and we use the same parameters: $T=1$, $H=0.07$, $\rho=-0.90$, $\xi_{0}=0.09$, $\eta=1.9$, $S_{0}=100$, $r=0.05$ and strike $K=70,80,\dots,120,130$ or $140$. As far as the GPR-Tree is concerned, we employ $N=50$ or $N=100$ time steps with $m=2$, $P=500,1000,2000$ or $4000$ random paths, and $J=0,1,3,7$ or $15$ past values. As far as the GPR-EI is concerned, we employ $N=50$ or $N=100$ time steps, $P=1000,2000,4000$ or $8000$ random paths, and $J=0,1,3,7$ or $15$ past values. Similar to what observed by Bayer et al. , the difference changing the value of $J$ does not impact significantly on the price, which indicates that considering the non-Markovian nature of the processes in the formulation of the exercise strategies is not particularly relevant. Conversely, using a large number of predictors significantly increases computational time. Numerical results are reported in Tables 4 and 5, together with the results reported by Bayer et al. in . Prices are very close to the benchmark, except for the case $K=120$: in this case with both the two GPR methods we obtain a price which is close to $20.20$ while Bayer et al. obtain $20.00$. Anyway, it is worth noticing that the relative gap between these two results is less than $1\%$ .

Conclusions

In this paper we have presented two numerical methods to compute the price of American options on a basket of underlyings following the Black-Scholes dynamics. These two methods are based on the GPR-Monte Carlo method and improve its results in terms of accuracy and computational time. The GPR-Tree method can be applied for dimensions up to $d=20$ and it proves to be very efficient when $d\leq 10$. The GPR-Exact Integration method proves to be particularly flexible and stands out for the small computational cost which allows one to obtain excellent estimates in a very short time. The two methods also turns out to be an effective tool to address non-Markovian problems such as the pricing of American options in the rough Bergomi model. These two methods are thus a step forward in overcoming the curse of dimensionality.

References

Appendix A Proof of Proposition 1

Let $n\in\left\{0,\dots,N-1\right\}$ and suppose the function $u\left(t_{n+1},\cdot\right)$ at time $t_{n+1}$ to be known at $Z$. Let us define the quantity

for $p=1,\dots,P$. The function $u\left(t_{n},\cdot\right)$ at time $t_{n}$ at $\mathbf{z}^{p}$ follows

where $\mathbf{Z}_{t_{n+1}}$ is the random variable defined as

Let us define $\Pi=\left(\Pi_{i,j}\right)$ as the $d\times d$ covariance matrix of the log-increments, that is $\Pi_{i,j}=\rho_{i,j}\sigma_{i}\sigma_{j}\Delta t$ . Moreover, let $\Lambda$ be a square root of $\Pi$ and $\mathbf{G}$ as a vector that follows a standard Gaussian law. Then, we observe that $\mathbf{Z}_{t_{n+1}}$ has the following conditional law

Moreover, relation (A.8) can also be stated as

Let $f_{\mathbf{z}^{p}}\left(\mathbf{z}\right)$ denote the density function of $\mathbf{Z}_{t_{n+1}}$ given $\log\left(\mathbf{S}_{t_{n}}\right)-\left(r-\frac{1}{2}\boldsymbol{\sigma}^{2}\right)t_{n}=\mathbf{z}^{p}$ . Specifically,

In particular, with reference to (A.13), the additional parameters $\sigma_{l}$ and $\sigma_{f}$ are called hyperparameters and are obtained by means of a maximum likelihood estimation. So let

be the GPR approximation of the function $u\left(t_{n+1},\mathbf{z}\right)$, where $\boldsymbol{\omega}=\left(\omega_{1},\dots,\omega_{q},\dots\omega_{P}\right)^{\top}$ in (A.14) is a vector of weights that can be computed by solving a linear system (see Rasmussen and Williams ). The GPR-EI approximation $C_{n}^{GPR-EI}$ of the continuation value is then given by

To compute each integral in (A.16), we observe that

where $\ast$ is the convolution product and $g_{\mathbf{-z}^{q}}$ is the density function of a Gaussian random vector which has law given by $\mathcal{N}\left(-\mathbf{z}^{q},\sigma_{l}^{2}I_{d}\right)$. Moreover, the convolution product of the densities of two independent random variables is equal to the density of their sum (see Hogg et al. ) and we can obtain the following relation which allows one to exactly compute the integrals in (A.16):

Therefore, the GPR-EI approximation $C_{n}^{GPR-EI}$ at $\hat{\mathbf{x}}^{p}$ reads

and the GPR-EI approximation $u_{n}^{GPR-EI}$ of the option value $u\left(t_{n},\cdot\right)$ at time $t_{n}$ and at $\mathbf{z}^{p}$ is given by

Appendix B Covariance of the vector R𝑅R in (6.1)

Let us report the formulas for the covariance of the components of the vector $R$ in (6.1). For all $n=1,\dots,N$, and $m=1,\dots,n-1$, the following relations hold:

Appendix C Proof of Proposition 3

Let us denote the random vector $\left(S_{t_{i}},V_{t_{i}},S_{t_{i+1}},V_{t_{i+1}},\dots,S_{t_{j}},V_{t_{j}}\right)^{\top}$ for $i,j\in\left\{0,\dots,N\right\}$ and $i<j$ with $\mathbf{SV}_{i:j}$ . We observe that the option value $v\left(t_{N},\cdot\right)$ at time $t_{N}$ is given by the payoff function $\Psi$, which only depends by the final value of the underlying. The option value $v\left(t_{N-1},\cdot\right)$ at time $t_{N-1}$ about the $p$-th path is given by

where $C$ stands for the continuation value and it is equal to

We approximate the continuation value in (C.2) by means of the GPR approximation of $\Psi$. In particular, let $\Psi^{GPR}\left(z\right)$ be the approximation of the function $z\mapsto\Psi\left(\exp\left(z\right)\right)$ by using the GPR method employing the Squared Exponential Kernel and considering the log-underlying values at maturity as predictors. Specifically, the predictor set is

where $k_{SE}$ is the Squared Exponential kernel, $\sigma_{l}$ is the characteristic length scale, $\sigma_{f}$ is the signal standard deviation and $\omega_{1},\dots,\omega_{P}$ are weights.

So we approximate the continuation value $C\left(t_{N-1},\mathbf{SV}_{1:\left(N-1\right)}^{p}\right)$ with the expression:

We observe that the law of $\log\left(S_{t_{N}}\right)$ given $S_{t_{1}}^{p},V_{t_{1}}^{p},\dots,S_{t_{N-1}}^{p},V_{t_{N-1}}^{p}$ is normal

Therefore, the GPR-EI approximation for the continuation value at time $t_{N-1}$ is as follows:

Taking advantage of the properties of the convolution between density functions, we obtain

Appendix D Proof of Proposition 4

In order to proceed backward, from $t_{N-2}$ up to $t_{1}$ we consider an integer positive value $J$ and train the GPR method considering the last $J+1$ observed values of the couple $\left(\log\left(S_{t_{n}}^{p}\right),\log\left(V_{t_{n}}^{p}\right)\right)$ as predictors, and the option price as response. Specifically, the predictor set is

approximates $v\left(t_{N-1},\mathbf{SV}_{1:\left(N-1\right)}^{p}\right)$.

Since the predictors have different nature (log-prices and log-volatilities at different times), we use the Automatic Relevance Determination (ARD) Squared Exponential Kernel $k_{ASE}$ to perform the GPR regression. In particular, if $d$ is the dimension of the space containing the predictors, it holds

As opposed to the Squared Exponential kernel, the ARD Squared Exponential kernel considers a different length scale $\sigma_{i}$ for each predictor that allows the regression to better learn the impact of each predictor on the response.

where $z_{i}^{q}=\log\left(S_{n+1-\left(i-1\right)/2}^{q}\right)$ if $i$ is even and $z_{i}^{q}=\log\left(V_{n+1-i/2}^{q}\right)$ if $i$ is odd, for $i=1,\dots,d_{n+1}$. This means that $z_{i}^{q}$ is the observed log-price at time $t_{n+1-\left(i-1\right)/2}$ of the $q$-th path if $i$ is even, and it is the observed log-volatility at time $t_{n+1-\left(i-1\right)/2}$ of the $q$-th path if $i$ is odd.

where $\Lambda_{2n+2}$ is the $2n+2$-th row of the matrix $\Lambda$ and $\underline{\mathbf{G}}^{p}=\left(G_{1}^{p},\dots,G_{2n}^{p},0\dots,0\right)^{\top}$. Moreover, the covariance matrix is given by

where $\Lambda_{i,j}$ stands for the element of $\Lambda$ in position $i,j$. Using a similar reasoning as done for the continuation value at time $t_{N-1}$, one can obtain the following GPR-EI approximation for the continuation value at time $t_{n-1}$:

where $h_{q}^{p}$ and $f_{q}^{p}$ are two factors given by

In particular, $h_{q}^{p}$ measures the impact of the past observed values on the price, whereas $f_{q}^{p}$ integrates the changes due to the diffusion of the underlying and its volatility.

Finally, we observe that, in order to compute $u_{n}^{GPR}$, we train the GPR method considering the predictor set given by

By induction we can compute the option price value for $n=N-2,\dots,0$ .

To conclude, we observe that the continuation value at time $t=0$ can be computed by using (D.9) and considering $h_{q}^{p}=1$ for $q=1,\dots,P$ since in this case, there are no past values to consider.