Variance Reduction Applied to Machine Learning for Pricing Bermudan/American Options in High Dimension
Ludovic Goudenège Andrea Molent Antonino Zanette

	
Quantum enhancement of sensitivity achieved by photon-number-resolving detection in the dark port of a two-path interferometer operating at high intensities

    Holger F. Hofmann
    Submitted 6 July 2018 / Accepted 19 March 2019
==============================================================================================================================================================



	

		11pt
		

	
	

		Abstract
		
 	
In this paper we propose an efficient method to compute the price of multi-asset American  options, based on Machine Learning, Monte Carlo simulations and variance reduction technique. Specifically, the options we consider are written on a basket of assets, each of them following a Black-Scholes dynamics.
 In the wake of Ludkovski's approach <cit.>, we implement here a backward dynamic programming  algorithm which considers a finite number of uniformly distributed exercise dates.  On these dates, the option value is computed as the maximum between the exercise value and the continuation value, which is obtained by means of Gaussian process regression technique and Monte Carlo simulations. 
Such a method performs well for low dimension baskets but it is not accurate for very high dimension baskets.  In order to improve the dimension range, we employ the European option price as a control variate, which allows us to treat very large baskets and moreover to reduce the variance of  price estimators.
Numerical tests show that the proposed algorithm is fast and reliable, and it can handle also American options on very large baskets of assets, overcoming the problem of the curse of dimensionality. 
	

	
	Keywords: 
 Finance; Gaussian process regression; Control  variate; American options; Monte Carlo methods.

	11pt
	

§ INTRODUCTION

	
In this paper we consider one of the most compelling problems among the still open issues in the field of  computational finance: pricing and hedging American options in high dimension. 
From a practical point of view, the efficient numerical evaluation of American options which consider as underlying a baskets of d assets is very challenging because of the so-called “curse of dimensionality”, which avoids the direct application of standard numerical schemes such as finite difference or tree methods. Specifically, this curse of dimensionality means that the computational cost and the memory requirement  increase exponentially with the dimension of the problem.

Several new ideas have appeared in this research area, which can be divided into five groups. The first type of approach consists in employing a recombinant tree in order to obtain a discretization of the underlying diffusion. An example of this mode is given by
the stochastic mesh method of Broadie and Glasserman <cit.>, the quantization
algorithms of Bally, et al.
<cit.>, the stochastic grid method of Jain and Oosterlee <cit.>. 
The second idea makes use of regression on a truncated basis of L^2 in order to compute the
conditional expectations. This is done in Longstaff and Schwartz
<cit.> and in
Tsisiklis and Van Roy <cit.>. 
The third concept consists in exploiting the representation formulas for the conditional expectation using Malliavin calculus.
This has been done by Lions and Reigner <cit.>, Bouchard and Touzi <cit.>, Bally et al. <cit.> Caramellino and Zanette <cit.> and Abbas-Turki and Lapeyre<cit.>. 
Another group of ideas relies on duality-based approaches for Bermudan option pricing, which are proposed by Rogers <cit.>, Haugh and Kogan <cit.>,  Andersen and Broadie <cit.>, Schoenmakers et al. <cit.> and Lelong <cit.>, which can be used to construct bounds on the option value.
Finally, the last group consists of methods that employ Machine Learning techniques to learn the continuation value or the stopping rules. This has been proposed by Becker et al. <cit.>, Kohler et al. <cit.> and Ludkowski <cit.>.


European prices can be used as control variate while pricing American options, as done, for example, by Bally et al. <cit.> and by Caramellino and Zanette <cit.>.  Since multi-asset products are considered, efficiently computing European prices is not trivial and many authors developed valid methods in this field.  
	Some of them focused on computing lower and upper bounds, such as Deelstra et al. <cit.>,   Carmona and Durrleman <cit.>, Caldana et al. <cit.>.
	Other approaches for basket options are based on the approximation of the sum of the log-normal distributions with a simple  distribution by matching some moments, as done by Levy <cit.>, Milevsky and Posner <cit.>, Zhou and Wang <cit.>, Korn and Zeytun <cit.>. Moreover, an approximation approach is also proposed  by Li and Wu <cit.> for options on several mean-reverting assets.
	Recently, Glau et al. <cit.> and Glau et al. <cit.> consider Chebyshev based methods for pricing. Deep Learning techniques are nowadays widely used in solving large differential equations, which is intimately related to option pricing: recent progresses in this field have been achieved by  Han et al. <cit.>, E et al. <cit.> and Beck et al.  <cit.>.	Finally, efficient Monte Carlo approaches are developed by Jourdain and Lelong <cit.> and more recently by Bayer et al. <cit.>. 


In this paper, we propose a new method that combines  Machine Learning, Monte Carlo simulations and variance reduction control variate technique. In particular, the use of a control variate makes the method more stable and extends its applicability range to high very large baskets. Moreover, the variance of price estimator is significantly reduced.

First of all, we implement a version  of the Ludkovski's algorithm <cit.>. Such an algorithm  proceeds backward over time by computing the price function on a set of prearranged points which represents possible values of the underlying. In particular, at each time step, it uses a set of Monte Carlo simulations together with Gaussian Process Regression (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. We term such an algorithm GPR Monte Carlo (GPR-MC).
The GPR-MC algorithm works very well for small baskets (in his paper, Ludkovski considers up to 5 dimensional basket), but it does not for large ones. In this paper, we show that, if one considers the European price as a control variate, the algorithm improves significantly and the variance of the price estimator is reduced. We term GPR Monte Carlo Control Variate (GPR-MC-CV) this new algorithm. Moreover, in order to compute the European prices, we suggest to use a semi-analytical formula, named  GPR-EI formula, introduced by Goudenège et al. in <cit.>, which proves to be efficient when many repeated computations of European prices have to be performed, or alternatively, Quasi-Monte Carlo simulations. 
 Finally, we investigate   the benefits brought by control variate technique to the GPR-Tree and GPR-EI approaches introduced by Goudenège et al. <cit.>.
	The paper is organized as follows. In Section 2 we
	present American options for the Black-Scholes d-dimensional model. In Section
	3 we briefly review Gaussian Process Regression, we present the GPR-EI formula, the GPR-MC method and the GPR-MC-CV method. Furthermore, we also investigate the use of control variate technique for the GPR-Tree and GPR-EI methods. In Section 4 we report
    some numerical results about pricing and variance reduction. Finally, Section 5 draws the 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, from inception up to maturity.
Let 𝐒=(𝐒_t)_t∈ [0,T] denote the d-dimensional underlying  process.
Such a stochastic process is assumed to randomly evolve according to the multidimensional Black-Scholes model: under the risk neutral measure, such a model is given by the following equation

    dS^i_t=(r-η_i)   S^i_t  dt +σ_i  S^i_t  dW^i_t,    i=1,…, d,

with 𝐒_0=( s_0,1,…,s_0,d)^⊤∈_+^d the spot price, r  the (constant) spot interest rate, η=(η_1, …,
η_d)^⊤ the vector of (constant) dividend rates, 
 σ = ( σ_1, …,σ_d)^⊤  the vector of (constant) volatilities, 𝐖 a d-dimensional correlated Brownian motion and ρ_ij the
instantaneous correlation coefficient between W^i_t and W^j_t.
Moreover, let Ψ(𝐒_T) denote the cash-flow associated with the option.
Thus, the price at time t of an American option having maturity T
and payoff function Ψ  : _+^d→ is then


    v^AM(t,𝐱)=sup_τ∈ T_t,T_t,𝐱[ e^-r (τ-t) Ψ(𝐒_τ)],

	 where T_t,T stands for the set of all the stopping times
taking values on [t,T] and _t,𝐱[ ·] is the expectation given all the information at time t and assuming 𝐒_t=𝐱.

For simulation purposes, the d-dimensional Black-Scholes model can be written alternatively using the Cholesky decomposition. 
Specifically, for i ∈{1, …, d} we can write

    dS^i_t = S^i_t ((r-η_i)  dt + σ_i Σ_i d𝐁_t),

where 𝐁 is a d-dimensional Brownian motion and Σ_i is the i-th row of the matrix Σ defined as a square root of the
correlation matrix Γ, given by

    Γ = [    1 ρ_12      ρ_1d; ρ_21    1    ⋱    ⋮;    ⋮    ⋱    ⋱    ⋮; ρ_d1              1 ].

	

§ MACHINE LEARNING FOR AMERICAN OPTIONS IN THE MULTI-DIMENSIONAL BLACK-SCHOLES MODEL

	
	

 §.§ Gaussian Process Regression

	
	In this Section, we present a brief review of Gaussian Process Regression
	and for a comprehensive treatment we refer to Rasmussen and Williams
	<cit.>.
	
	Gaussian Process Regression (GPR),  also known as Kriging (see Matheron <cit.>, Journel and Huijbregts <cit.>), is a class of non-parametric
	kernel-based probabilistic models which represents the input data
	as the random observations of a Gaussian stochastic process. The most
	important advantage of this approach in relation to other parametric
	regression techniques is that it is possible to effectively exploit
	a complex dataset which may consist of points sampled randomly in
	a multidimensional space.
	
	In general, a Gaussian process 𝒢 is a collection of random
	variables defined on a common probability space (Ω,ℱ,P),
	any finite number of which have consistent joint Gaussian distributions.
	We are interested in Gaussian processes for which the random variables
	in 𝒢 are indexed by a point 𝐱∈ℝ^d,
	d∈ℕ. Therefore, for all 𝐱∈ℝ^d,
	𝒢(𝐱):Ω→ℝ
	is a Gaussian random variable and if X={𝐱_p,p=1,…,P}⊂ℝ^d
	then (𝒢(𝐱_1),…,𝒢(𝐱_P))^⊤
	is a random Gaussian vector. Moreover, a Gaussian process is fully
	specified by its mean function μ(𝐱):ℝ^d→ℝ
	(which is usually assumed to be zero) and by its covariance function k(𝐱,𝐱'):ℝ^d×ℝ^d→ℝ. 
	
	Now, let us consider a training set 𝒟 of P observations
	(the input data), 𝒟={(𝐱_p,y_p), p=1,…,P}
	where X={𝐱_p,p=1,…,P}⊂ℝ^d
	denotes the set of input vectors and Y={ y_p,p=1,…,P}⊂ℝ
	denotes the set of scalar outputs. These observations are modeled
	as the realization of the sum of a Gaussian process and a noise source.
	Specifically,
	
    y_p=f_p+ε_p,

	where { f_p=𝒢(𝐱_p),p=1,…,P}
	is a Gaussian process and {ε_p,p=1,…,P}
	are i.i.d. random variables such that ε_p∼𝒩(0,σ_P^2).
	Moreover, the distribution of 𝐟=(f_1… f_P)^⊤
	is assumed to be given by
	
    𝐟∼𝒩(0,K(X,X)),

	where K(X,X) is a P× P matrix with K(X,X)_p_1,p_2=k(𝐱_p_1,𝐱_p_2) for p_1,p_2=1,…,P with k:ℝ^d×ℝ^d→ℝ   the so called   kernel function.
	Thus
	
    𝐲∼𝒩(0,K(X,X)+σ_P^2I_P),

	where I_P is the P× P identity matrix.
	
	Now, in addition, let us consider a test set X̃ of M
	points {𝐱̃_m,m=1,…,M}. The
	realizations f̃_m=𝒢(𝐱̃_m)
	are not known but rather we want to estimate them by exploiting the
	observed realizations of 𝒢 in 𝒟. The a
		priori joint distribution of 𝐲 and 𝐟̃=(f̃_1,…,f̃_M)^⊤
	is given by 
	
    [[  𝐲; 𝐟̃ ]]∼𝒩([[ 0_P; 0_M ]],[[ K(X,X)+σ_P^2I_P         K(X,X̃);         K(X̃,X)        K(X̃,X̃) ]])

	where K(X̃,X̃) is a M× M matrix
	given by K(X̃,X̃)_m_1,m_2=k(𝐱̃_m_1,𝐱̃_m_2) for m_1,m_2=1,…,M,
	K(X,X̃) is a P× M matrix given by K(X,X̃)_p,m=k(𝐱_p,𝐱̃_m) for p=1,…,P, m=1,…,M
	and K(X̃,X) is a M× P matrix given by
	K(X̃,X)_m,p=k(𝐱̃_m,𝐱_p) for m=1,…,M, p=1,…,P. 
	
	Since we know the values for the training set, we can consider the
	conditional distribution of 𝐟̃ given 𝐲.
	It is possible to prove that 𝐟̃|X̃,𝐲,X
	follows a Gaussian distribution given by 
	
    𝐟̃|X̃,𝐲,X∼𝒩(𝔼[𝐟̃|X̃,𝐲,X],Cov[𝐟̃|X̃,𝐲,X]),

	where 
	
    𝔼[𝐟̃|X̃,𝐲,X]=K(X̃,X)[K(X,X)+σ_P^2I_P]^-1𝐲

	and
	
    Cov[𝐟̃|X̃,𝐲,X]=K(X̃,X̃)-K(X̃,X)[K(X,X)+σ_P^2I_P]^-1K(X,X̃).

	Therefore, a natural choice consists in predicting the values 𝐟̃
	through 𝔼[𝐟̃|X̃,𝐲,X]. 
	Moreover, by using equation  (<ref>),  one can define a function f^GPR:ℝ^d→ℝ that approximates the function 𝐱_p↦ y_p by setting
	
    f^GPR( 𝐱̃)    =𝔼[𝐟̃|{𝐱̃},𝐲,X]
       =∑_p=1^Pk(𝐱̃,𝐱^p)ω_p,

	where ω=( ω_1,…,ω_1)^⊤ is a vector of weights determined by 
	
    ω=[K(X,X)+σ_P^2I_P]^-1𝐲.
 
	
	The computation in (<ref>) requires the knowledge of
	the covariance function K and of the noise variance σ_P^2.
	A commonly used covariance function is the Matern 3/2  kernel k_Ma:ℝ^d×ℝ^d→ℝ, which is given by
	
    k_Ma(𝐱,𝐱')=σ_f^2(1+√(3)‖𝐱-𝐱'‖ _2/σ_l)exp(-√(3)‖𝐱-𝐱'‖ _2/σ_l)  for 𝐱,𝐱'∈ℝ^d,

	where σ_f^2 is called the signal variance and σ_l
	is called the length-scale. Another possible choice is the Squared Exponential kernel k_SE:ℝ^d×ℝ^d→ℝ,
	which is given by
	
    k_SE(𝐱,𝐱')=σ_f^2exp(-‖𝐱-𝐱'‖ _2^2/ 2σ_l^2)    for 𝐱,𝐱'∈ℝ^d.

	In general, the choice of kernel function is performed by using a log-likelihood criterion.  	 The parameters
	σ_f^2, σ_l of the kernel function and
	σ_P^2 of the noise are called hyperparameters and need
	to be estimated. A common approach is to consider the maximum likelihood
	estimates which can be obtained by maximizing the log-likelihood function
	of the training data, that is by maximizing the following function:
	
    -1/2log((K(X,X)+σ_P^2I_P))-1/2𝐲^⊤[K(X,X)+σ_P^2I_P]^-1𝐲.



	
	The development of the GPR model can be divided in the training step
	and the evaluation step (also called testing step). The training step
	only requires the knowledge of the training set 𝒟 and
	it consists in estimating  the hyperparameters
	and computing the vector of weights ω.
	The evaluation step can be computed only after the training step has
	been accomplished and it consists in obtaining the predictions via
	the computation of K(X̃,X) ω.
	We stress out that the training step is independent of the test set
	X̃. Thus one can store the values computed during the training
	step and perform the evaluation step many times with a small computational
	cost, which is 𝒪(P· M).
	
 We observe that the computation time depends only marginally on the size d of the space where the points lie, as the value of d only impacts in the time taken to calculate distances between the points which appears in the covariance matrix K, that is ‖𝐱-𝐱'‖ _2.

	


 §.§ Machine Learning Exact Integration for European options

In order to improve the GPR-MC approach, we employ the European option price as a control variate. 
	Here, we propose to compute such a price by means of the semi-analytical formula introduced by Goudenège et al. <cit.>, that we term GPR-EI formula. This computation is based on two steps. First of all, the payoff function is approximated  by means of  GPR. Then, the European price is computed as the discounted expected value of the final cash flow, that is a multidimensional integral of the payoff function with respect to the log-underlying process density. Such an integral can be computed by means of a closed formula when replacing the true payoff function with its GPR approximation. 

Let us consider a set Z={𝐳^q,q=1,…,Q} consisting of Q points in ℝ^d
 quasi-randomly distributed
according to the law of the vector (σ_1W_T^1,…,σ_dW_T^d)^⊤.
In particular, we define 

    𝐳_i^q=√(T)σ_iΣ_i𝐡^q,

where Σ_i is i-th row of the matrix Σ and 𝐡^q
is the q-th point of the Halton sequence in ℝ^d (other
low-discrepancy sequence can be considered, such as Solob's or Faure's
ones).
Let u:Z→ℝ be the function defined by

    u(𝐳):=Ψ(𝐒_0exp((r-η-1/2σ^2)T+𝐳)).

In a nutshell, the main idea is to approximate the function u by
training the GPR method on the  set Z. In particular, we employ
the Squared Exponential kernel defined in (<ref>). Equation (<ref>) allows one to approximate the function u(·) by

    u^GPR(𝐳)=∑_q=1^Qk_SE(𝐳^q,𝐳)ω_q,

where ω_1,…,ω_P are weights. The continuation value can be computed
by integrating the function u^GPR against a d-dimensional
probability density. The use of the Squared Exponential kernel  allows one to easily perform such a calculation  by means of a closed formula. 
Specifically, the GPR-EI method relies on the following Proposition.

 Let us consider an European option with payoff function Ψ, inception t=0, maturity T, and multidimensional underlying following the dynamics in (<ref>) with spot price 𝐒_0.	The  price of such an  option at t=0 can be approximated	by
	
    v^EU=e^-rTE_0,𝐒_0[Ψ(𝐒_T)]≈ e^-rT∑_q=1^Qω_qσ_f^2σ_l^de^-1/2(𝐳^q)^⊤(T·Π+σ_l^2I_d)^-1(𝐳^q)/√((T·Π+σ_l^2I_d))

	where σ_f, σ_l, and ω_1,…,ω_Q
	are certain constants determined by the GPR approximation of the function
	𝐳↦ u(𝐳) considering Z as
	the predictor set, and Π=(Π_i,j) is the d× d
	covariance matrix of the vector (σ_1W_T^1,…,σ_dW_T^d)^⊤, that is Π_i,j=ρ_i,jσ_iσ_j.
 The proof of this Proposition is very similar to the one reported
in <cit.>.

Despite the GPR-EI formula (<ref>) is adapted to compute the option price supposing the spot price to be 𝐒_0 and the time to maturity to be T, it works quite well also for   spots close to 𝐒_0 and time to maturity smaller than T. The following Proposition states how to do that.


	Let us consider and European option with payoff function Ψ, inception 0<t̃<T, maturity T, and multidimensional underlying following the dynamics in (<ref>). Let  𝐒̃ be the vector of the spot prices at time t̃ and define 𝐳̃∈ℝ^d such that 
	
    𝐒̃=𝐒_0exp((r-η-1/2σ^2)t̃+𝐳̃) .
 
	The  price of such an  option at t̃ can be approximated	by
	
    v^EU=e^-r(T-t̃)E_t̃,𝐒̃[Ψ(𝐒_T)]≈ e^-r(T-t̃)∑_q=1^Qω_qσ_f^2σ_l^de^-1/2(𝐳^q-𝐳̃)^⊤((T-t̃)·Π+σ_l^2I_d)^-1(𝐳^q-𝐳̃)/√(((T-t̃)·Π+σ_l^2I_d))

	where σ_f, σ_l, and ω_1,…,ω_Q
	 and Π=(Π_i,j) are defined according to Proposition <ref>.

The proof of   Proposition <ref> derives directly from Proposition <ref> by considering 𝐒̃ in place of 𝐒_0. The hyperparameters   σ_f, σ_l, and the weights   ω_1,…,ω_Q
need to be computed only once and then we can use formulas (<ref>) and (<ref>) to compute the European prices. The resolution of the linear systems within the exponential factors and the computation of the matrix determinant in (<ref>) and (<ref>) can be done quite fast by computing the Cholesky decomposition of the matrices (T-t̃)·Π+σ_l^2I_d for each of the few possible values of t, that is t=0,t_1,…,t_N-1.
 For this reason, turns out to be faster than repeated Monte Carlo simulations to compute the many European prices to be used as control variate. 
 



 §.§ Machine Learning Control Variate algorithm for  American options



  §.§.§ The GPR Monte Carlo Method

Let us introduce the GPR Monte Carlo approach.
We approximate the price of an American option with the price of a Bermudan option on the same basket.  Specifically, let N be the number of time steps and Δ t=T/N the time increment. The discrete exercise dates are  t_n=n Δ t, as
n=1,…,N.
If 𝐱 represents the vector of the underlying prices at the exercise date t_n, then the price of the Bermudan option v^BM is given by

    v^BM(t_n,𝐱)=max(Ψ(𝐱),E_t_n,𝐱[e^-rΔ tv(t_n+1,𝐒_t_n+1)]).

First of all, by knowing the function v^BM(t_n+1,·), one can compute   v^BM(t_n,·) by approximation 
of the expectation in  (<ref>).
In order to do that, we consider a set X^n of P points whose coordinates represent certain possible values for the underlyings at time t_n: 

    X^n={𝐱^n,p=(x_1^n,p,…,x_d^n,p),p=1,…,P}⊂ℝ^d.

 Suppose now we want to compute  v^BM(t_n,·) but only for 𝐱^n,p∈ X^n.
This goal can be achieved by means of a one step Monte Carlo simulation. In particular, for each 𝐱^n,p∈ X^n, we simulate a set of points X̃^n_p   

    X̃^n_p={𝐱̃^n,p,m=(x̃_1^n,p,m,…,x̃_d^n,p,m),m=1,…,M}⊂ℝ^d

of M  possible values for 𝐒_t_n+1 according to the law of 𝐒_t_n+1|.. In particular, for i=1,…,d, n=1,…,N, p=1,…,P, m=1,…,M, we define

    x̃ ^n,p,m_i =x^n,p_i e^(r-η_i - 1/2σ^2_i) Δ t+√(Δ t)σ_iΣ_i  𝐆^n,p,m,

where 𝐆^n,p,m∼𝒩(0,I_d)  is a standard Gaussian random vector and Σ_i is the i-th row of the matrix Σ, just as in  (<ref>).
Then, the option price can be approximated for each 𝐱^n,p∈ X^n  by

    v̂^BM(t_n,𝐱^n,p)=max(Ψ(𝐱^n,p),e^-rΔ t/M∑_m=1^Mv^BM(t_n+1,𝐱̃^n,p,m)),

if the quantities v^BM(t_n+1,𝐱̃^n,p,m) are known for all of these simulated points 𝐱̃^n,p,m. If we proceed backward, the function v^BM(t,·) is known for t=T since it is equal to the payoff function Ψ(·) and thanks to (<ref>) it is known, through an approximation, also for t=t_N-1 and 𝐱^N-1,p∈ X^N-1. In order to assess v^BM(t_N-2,𝐱^N-2,p) for all 𝐱^N-2,p∈ X^N-2, and thus going on up to t = 0, it is necessary to evaluate the function v^BM(t_N-2,·) for all the points  in X̃^N-2=⋃_p=1^PX̃^N-2,p. This cannot be done directly since we know  v̂^BM(t_N-1,·) only for the points in X^N-1 and not for all those in X̃^N-2.
To overcome this issue, we compute the approximation of the function v̂^BM(t_N-1,·) by means of the GPR technique. In particular the set X^N-1 serves as the predictor set and {v̂^BM( t_N-1, 𝐱^N-1,p) ,p=1,…,P } as the response set.

 More generally, let v^BM,GPR_n(·) be the GPR approximation of v̂^BM(t_n,·) trained by considering X^n as the predictor set and {v̂^BM( t_n, 𝐱^n,p) ,p=1,…,P } as the response set, where v̂^BM is defined as in (<ref>). Then, we can proceed backward by computing 

    v̂^BM(t_n-1,𝐱^n-1,p)=max(Ψ(𝐱^n-1,p),e^-rΔ t/M∑_m=1^Mv^BM,GPR_n(𝐱̃^n-1,p,m)).

and by computing v^BM,GPR_n-1, that is the GPR approximation of v̂^BM( t_n-1,·). 
Finally, the option price at time t=0   is computed through

    v̂^BM(0,𝐒_0)=max(Ψ(𝐒_0),e^-rΔ t/M∑_m=1^Mṽ_1^BM,GPR(𝐱̃^0,m))

where the points 𝐱̃^0,1,…, 𝐱̃^0,M are random simulations of 𝐒_t_1 given by 


    x̃ ^0,m_i =S_0^i e^(r-η_i - 1/2σ^2_i) Δ t+√(Δ t)σ_iΣ_i  𝐆^0,m,

where 𝐆^0,m∼𝒩(0,I_d)  is a standard Gaussian random vector for any m∈{ 1,…,M}.

The choice of the sets X^n,n=1…,N-1 is a sensitive question. 
Similarly to what proposed by Ludkovski <cit.>, here we use a  deterministic space-filling sequence based on the Halton sequence.
Specifically, let 𝐡^p be the p-th point of the Halton quasi-random sequence in ℝ^d and Φ^-1 the  inverse cumulative distribution of a standard normal distribution.
We define the points 𝐱^n,p as follows: 

    𝐱^n,p_i =𝐒_0^i e^(r-η_i - 1/2σ^2_i) t_n+√(t_n)σ_iΣΦ^-1( H^p),

for i=1,…,d, n=1,…,N-1, and  p=1,…,P. 
This choice for the sets X^n proves to be the most effective, since the points used to train the GPR algorithm at time t_n are sampled according to the density function of the process 𝐒_t_n.



  §.§.§ The GPR Monte Carlo Control Variate Method

Let us present the GPR Monte Carlo Control Variate method (GPR-MC-CV), that is our proposed algorithm. 

The control variate technique  is commonly used to reduce the variance of Monte Carlo estimators, but it can also give its contribution in American pricing. Following Bally et al. <cit.> and  Caramellino and Zanette <cit.>, we employ the European price as a control variate for the American price. Let us consider an American and an European option  with the same payoff function Ψ and maturity T, and let     v^AM,v^EU denote their prices  respectively. For a fixed time t and underlying stocks 𝐱, we define the American-European price gap as:  

    v(t,𝐱)=v^AM(t,𝐱)-v^EU(t,𝐱).

Then 

    v(T,𝐱)=0,

and it is straightforward to see that

    v(t,𝐱)=sup_τ∈𝒯_t,TE_t,𝐱[e^-r(τ-t)Ψ̂(τ,𝐒_τ)].

where 𝒯_t,T stands for the set of all stopping times taking values in [ t,T] and Ψ̂ is defined by

    Ψ̂(t,𝐱)= Ψ(𝐱)-v^EU(t,𝐱).

We stress out that Ψ̂(T,𝐱)=0 and the function Ψ̂ depends on the time variable also.
Therefore, in order to numerically evaluate v^AM(0,𝐒_0), one can arrange a dynamic programming principle, based on Bermudan approximation, actually equal to the one in Section <ref> by replacing Ψ with Ψ̂. Once the initial price gap v (0,𝐒_0) has been calculated, one can retrieve the American price by computing

    v^AM(0,𝐒_0)=v(0,𝐒_0)+v^EU(0,𝐒_0).



		


 The sketch of the GPR-MC-CV algorithm is presented here.
	


[16.4cm][l]
[l]16.4cm


















 

	We remark that when using a quasi-Monte Carlo sequence it is important to consider leaping. This technique consists in considering only some uniformly subsampled points of the original sequence, which improves convergence. However,  the leap values,  must be chosen with care.  In fact, many values lead to sequences that do not touch on large sub-hyper-rectangles of the unit hypercube, failing to be a uniform quasi-random point set (see Kocis and Whiten <cit.>). A common rule for choosing the leap values for the Halton sequence consists in  setting the value to q-1, where q is a prime number that has not been used to generate the sequence. 


	We observe that the Monte Carlo evaluation of the continuation value    can be easily parallelized since the  summations in (<ref>), are independent of each other and can be calculated separately. Thus, this feature  allows one to significantly reduce the computational time.
	

As observed by Ludkovsi <cit.>, the main computational cost is	due to the training of the GPR model, which is proportional to the cube of the observation amount. In our case, this training has to be performed one time to compute the European prices with a cost 𝒪(Q^3 ) (with Q   the number of points employed in   European price computation), and N-2 times within the algorithm to approximate the American-European gap at a give time, thus 𝒪(N· P^3 ) (with N  the number of time steps and P the number of points used to train the GPR models at each time step). On the other hand, the cost of the Monte Carlo step depends on both the number of evaluations to be performed and to the number of points employed: the cost for such a step is 𝒪(N· P · M) (with M  the number of Monte Carlo simulations employed in estimating the continuation gap value). Finally, we observe that if we compute the European prices by using M'  Monte Carlo simulations instead of by using the GPR-EI formula, then the cost  𝒪(Q^3 ) is replaced by 𝒪(N· P · M').

	
  


  §.§.§ The Control Variate for GPR-Tree and GRP-EI

Although the control variable technique was initially conceived as a variance reduction techniques for Monte Carlo methods, it can also be a valid support in other contexts.
We investigate the benefits brought by this technique to the GPR-Tree and GPR-EI techniques introduced by Goudenège et al. <cit.> for pricing American options in high dimension. In particular, as proposed for the GPR-MC method, we use the European price as a control variate and we employ   GPR-Tree (or   GPR-EI) to compute the American-European price gap.
Let us give a brief introduction of these two numerical approaches. We  refer the interested reader to <cit.> for more details.
 
The GPR-Tree method   is similar to the GPR-MC method here proposed. The main difference consists in the use of a tree step  in place of random simulations to compute the continuation value. In particular, for each time step t_n and for each point 𝐱^p, 2^d future values are generated according to the tree method proposed by Ekvall <cit.>, in place of Monte Carlo simulations. Such a method is particularly efficient when the dimension d is low (that is, indicatively, it does not exceed 10).
 
The GPR-EI method differs from both the GPR-MC and GPR-Tree methods for three reasons. 
First of all, the predictors employed in the GPR step are related to the logarithms of the underlying value. Then, the continuation value at these points is computed through a  closed formula which comes from an  exact integration. Finally, the GPR-EI method employs the Squared Exponential kernel, which is given by
 	
    k(𝐱,𝐱')=σ_f^2exp(-1/2 σ_l^2𝐱-𝐱' ^2_2),

 	for 𝐱,𝐱'∈ℝ^d, where d is the dimension of the regression problem. 
 

	

§ NUMERICAL RESULTS


	
In this Section we report some numerical results in order to investigate
the effectiveness of the proposed Machine Learning algorithm for pricing
American options in the multi-dimensional Black-Scholes model. 

First of all, we compare the  GPR-MC and GPR-MC-CV methods considering Geometric and Arithmetic basket put options and then we focus on a Call on the Maximum option. Moreover, we study the benefits of using the control variable also for GPR-Tree and GPR-EI methods. Finally, we investigate the variance of the price estimators about the two methods.
We stress out that the GPR-Tree method is interesting only for low dimension options: when d exceeds 10, the method still works, but computational times grow exponentially.



 §.§ Geometric and Arithmetic Basket Put  Options

In this test we focus on two payoff that depend on the mean of the underlyings. Specifically, we consider the following payoff examples:


  * Geometric basket Put

    Ψ(𝐒_T)=(K - ( ∏_i=1^d S_T^i) ^1/d)_+,


  * Arithmetic basket Put 

    Ψ(𝐒_T)=(K - 1/d∑_i=1^d S_T^i)_+.






We consider both the GPR-MC  and the GPR-MC-CV method in order to investigates the benefits induced by the control variate technique.
We consider the same parameters as in <cit.>: T=1, S_i=100, K=100, r =0.05, equal (null) dividend rates η_i=0.0, equal volatilities σ_i=0.2, equal correlations
ρ_ij=0.2 and N=10 exercise dates. Moreover, we consider P=250,500 or 1000 points, M=10^3,10^4 or 10^5 Monte Carlo simulations and Q=10000 points for the computation of the European prices with the GPR-EI formula.
As opposed to the other input parameters, we vary the dimension d, considering d=2, 5, 10, 20, 40 and 100.
The algorithm has 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.


We present now the numerical results for the two payoff examples.

 First of all, let us present the European results, obtained by means of the GPR-EI formula.
Table <ref> reports the prices, changing the dimension d and the number of employed points Q. Moreover, we also report a Benchmark price computed by Monte Carlo simulation considering 10^6 samples (95% confidence intervals are ±0.01 for all the benchmark values.). As we can see with only 1000 points we can obtain accurate results in any considered dimension. 





Let us now focus on the American results. As far as the Geometric basket Put is considered, it is possible to reduce the problem of
pricing in the d-dimensional model to a one dimensional American Put option in the
Black-Scholes model with opportune parameters. The price of such a one dimensional American option   can be computed in a easy way, for example by using the CRR algorithm with 1000 steps (see Cox et al. <cit.>).
Therefore in this case we have a reliable  benchmark to test the algorithm. Moreover, when d is smaller than 10  we can also compute the price by means of a multi-dimensional binomial tree (see Ekvall et al. <cit.>). In particular, the number of steps employed for the 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.
Results are reported in Tables <ref> and <ref>. 
 We observe that both the two algorithms are very accurate in low dimension, despite we are approximating an American option with a Bermudan one. When larger baskets are considered, say d≥40, the prices obtained with the GPR-MC are less accurate and less stable while changing the number of points P and the number of Monte Carlo simulations M.  The computer processing time of the GPR-MC-CV method are a little higher than those of the GPR-MC  because European prices need to be computed. 

We also stress out that the computer processing time increase little with the size of the problem. This is due to the fact that the dimension affects significantly only the computational time of the Monte Carlo step while the GPR step is only minimally distressed (see Remark <ref>).

Table <ref> and <ref> report the  results for the GPR-Tree and GPR-EI methods employing or not the control variate technique.  By comparing the results of the two Tables, we observe that the option prices for d≤10 are very similar: in this case variate control technique is not crucial to improve convergence. As opposed to that, GPR-EI benefits sensitively from control variate technique when high values of d are considered.



 
 







 As opposed to the Geometric basket Put option,  we have no method to obtain a fully reliable benchmark when dealing with an Arithmetic basket Put  option. However, for small  values of d, a reference  price can be obtained by means of a multidimensional tree method (see Ekvall et al. <cit.>), just as shown for the Geometric case. 
Results are reported in Tables <ref>  and <ref>. The conclusions that we can draw in this case are similar to those for the Geometric case: both the two methods are accurate in low dimension,  while the control variate method is more effective in high dimension.
 
Table <ref> and <ref> report the  results for the GPR-Tree and GPR-EI methods employing or not the control variate technique.  Just as for the Geometric put option  we observe that the option prices for d≤10 are very similar: in this case variate control technique is not crucial to improve convergence. As opposed to that, control variate technique has an impact on GPR-EI results when high values of d are considered. Anyway, in this case, due to the lack of a benchmark price, it is difficult to draw clear cut conclusions.






 


 §.§ Call on the Maximum option

Let us consider a Call on the Maximum of d-assets American option, whose payoff is given by:

    Ψ(𝐒_T)=(max_i=1… d  S_T^i-K)_+.

 The Call on the Maximum setting is particularly interesting for investigating
scalability of our approaches in the dimension d of the problem. As observed by Ludkovski <cit.>, as opposed  to basket Put options, the stopping region of a  Call on the Maximum  consists of several disconnected pieces and this makes the pricing problem particularly challenging.
As done in the previous Section, we consider both the GPR-MC  and the GPR-MC-CV method in order to investigates the benefits induced by the use of this technique.
We consider the same parameters as those employed by Becker et al. <cit.>: T=3, S_i=100, K=100, r =
0.05, equal  dividend rates η_i=0.1, equal volatilities σ_i=0.2, equal (null) correlations
ρ_ij=0.0 and N=9 exercise dates. Moreover, we consider P=250,500 or 1000 points, M=10^3,10^4 or 10^5 Monte Carlo simulations. As opposed to the other input parameters, we vary the dimension d, considering d=2, 5, 10, 20, 30, 50 and 100.
In this particular case, because of the long maturity and unbounded payoff,  the GPR-EI formula is not very accurate when considering high dimension and initial points far from the spot 𝐒_0, and so we prefer computing the European price by means of Quasi-Monte Carlo simulation with 10^6 random simulations.


 First of all, let us present the European results, obtained by means of the GPR-EI formula.
Table <ref> reports the prices, changing the dimension d and the number of employed points Q. Moreover, we also report a Benchmark price computed by Monte Carlo simulation considering 10^6 samples (95% confidence intervals are ±0.01 for all the benchmark values). 



 The aforementioned testing set has also been considered by Becker et al. <cit.> and therefore  we report their results as reference prices.  Furthermore, for small  values of d, we can approximate the price obtained by means of a multidimensional tree method. 
Results, which are reported in Tables <ref> and <ref>, are quite meaningful. Both the two methods perform fine in low dimension, but when large baskets are considered outcomes are strongly different. As far as this particular dataset is considered, the GPR-MC approach gives several null results and others very high, which means that the GPR regression is not able to extrapolate the price surface correctly. In particular, this happens when d≥50. Increasing the number P of points   fixes things when d=50 and for P=1000 results are likely, although outside the confidence interval proposed by Becker et al. <cit.>. Anyway, when d=100 we always obtain null value, showing all the limits of the GPR-MC approach.
As opposed to the GPR-MC, the GPR-MC-CV method performs very well for all the considered dimensions and almost all the values obtained with P=1000 and M=10^5 are within the confidence intervals proposed by Becker et al. <cit.>.
 
Finally, Tables <ref> and <ref> report the results for the GPR-Tree and the GPR-EI method obtained by using or not the control variate technique. These two methods seems to be not very effective for the particular Bermudan option considered here. As far as the GPR-EI method is concerned, variates control technique improves the results. Such  a improvement is not evident with respect to the GPR-Tree method.  




 

 §.§ Variance Reduction


We conclude our numerical investigations by showing the effect of introducing a control variate on the variance of the estimated prices. 
In particular, we consider the same Geometric Put option as in Section <ref> and we price the same option 100 different times, changing the seed of the Monte Carlo generator. This allows us to estimate the variance of the price estimator and to make comparisons. 
Results are available in Tables <ref> and <ref>, that report the estimated standard deviations and their 95% confidence intervals, computed according to the method suggested by Sheskin <cit.>. It is evident that the the standard deviation (and thus the variance) of the prices obtained with the GPR-MC-CV method is several time lower than the one computed with the GPR-MC method. This is also confirmed for all the considered  combination of P and M, by the Hartley's F_max test (see Sheskin <cit.>) with a 99% confidence level.






§ CONCLUSIONS

In this paper we have proposed  a new approach   to price American options on baskets of assets, each of them following a Black-Scholes dynamics. The method employs Machine Learning technique,  Monte Carlo method and  variance reduction technique that exploits the European option price as a control variate. The European prices are computed by means of a semy-analitical formula or Quasi-Monte Carlo simulations. Numerical results show that the method is reliable and fast for baskets including up to 100 assets. The use of a control variate improves the algorithm accuracy and reduces the variance of the estimated prices. In certain cases, also the GPR-Tree and GRP-EI methods benefit from the use of a control variate. The computation time is small and shortly growing with respect to the dimension of the basket. Moreover,  the algorithm is partially parallelizable and therefore the computing time can be significantly reduced. Machine Learning seems to be a very promising tool for American option pricing in high dimension, overcoming the problem of the curse of dimensionality.
 	
	
	acm
	
	
