remarkRemark hypothesisHypothesis claimClaim assumptionAssumption Gradient-
enhanced LSM for American optionsJ. Yang and G. Li

# Gradient-enhanced sparse Hermite polynomial expansions for pricing and
hedging high-dimensional American options††thanks: Submitted to the editors
May XX, 2024. JY acknowledges support from the University of Hong Kong via the
HKU Presidential PhD Scholar Programme (HKU-PS). GL acknowledges the support
from GRF (project number: 17317122) and Early Career Scheme (Project number:
27301921), RGC, Hong Kong.

Jiefei Yang Department of Mathematics, University of Hong Kong, Pokfulam, Hong
Kong<EMAIL_ADDRESS>Guanglian Li Department of Mathematics,
University of Hong Kong, Pokfulam, Hong Kong<EMAIL_ADDRESS>
###### Abstract

We propose an efficient and easy-to-implement gradient-enhanced least squares
Monte Carlo method for computing price and Greeks (i.e., derivatives of the
price function) of high-dimensional American options. It employs the sparse
Hermite polynomial expansion as a surrogate model for the continuation value
function, and essentially exploits the fast evaluation of gradients. The
expansion coefficients are computed by solving a linear least squares problem
that is enhanced by gradient information of simulated paths. We analyze the
convergence of the proposed method, and establish an error estimate in terms
of the best approximation error in the weighted $H^{1}$ space, the statistical
error of solving discrete least squares problems, and the time step size. We
present comprehensive numerical experiments to illustrate the performance of
the proposed method. The results show that it outperforms the state-of-the-art
least squares Monte Carlo method with more accurate price, Greeks, and optimal
exercise strategies in high dimensions but with nearly identical computational
cost, and it can deliver comparable results with recent neural network-based
methods up to dimension 100.

###### keywords:

sparse Hermite polynomial expansion, least squares Monte Carlo, backward
stochastic differential equation, high dimensions, American option

60G40, 91G60, 91G20, 62J05, 65C30

##  1 Introduction

The early exercise feature of American or Bermudan options gives holders the
right to buy (call) or sell (put) underlying assets before the expiration
date, and their accurate numerical calculation is of great practical
importance. Meanwhile, the efficient estimation of Greeks (i.e., derivatives
of the price function, e.g., delta and gamma) is vital for hedging and risk
management, since the theory of option pricing builds on the assumption of the
absence of arbitrage. For example, when the asset price is on the rise, the
gain in the long position of a call writer’s asset may offset the potential
loss of the call option.

Nonetheless, the early exercise feature of American options poses significant
challenges for computing price and Greeks, especially in high dimensions. One
of the most popular methods for high-dimensional American option pricing is
the least squares Monte Carlo (LSM) method [15, 22]. Computing Greeks in high
dimensions is more involved, and further developments with LSM have been
proposed [24, 5]. In this work, building on LSM, we shall develop a simple,
fast, and accurate algorithm, termed as gradient-enhanced least squares Monte
Carlo (G-LSM) method, c.f. Algorithm 1, for computing price and Greeks
simultaneously at all time steps for dimensions up to 100. The key
methodological innovations includes using sparse Hermite polynomial space with
a hyperbolic cross index set as the ansatz space for approximating the
continuation value functions (CVFs), and incorporating the gradient
information for computing the expansion coefficients.

We elaborate on the two methodological innovations. First, the main obstacle
of using gradient information lies in the computational expense: a $d$-variate
price has $d$ partial derivatives, which grows quickly with $d$, especially
when the derivative evaluation is costly. With the proposed sparse Hermite
polynomial ansatz space, the derivatives of polynomial bases can be obtained
at almost no extra cost, cf. Eq. 9 below. This allows greatly reducing the
computational cost. Second, although using a polynomial ansatz space for the
CVFs as LSM, G-LSM constructs the expansion coefficients via solving a linear
least squares problem enhanced by the gradient (and hence the name G-LSM), by
minimizing the mean squared error between the approximate and exact value
functions at $t_{k+1}$, c.f. Eq. 15 below. This differs markedly from LSM,
which approximates the conditional expectations by projection and minimizes
the mean squared error of approximating CVF at $t_{k}$. Numerical experiments
show that this choice can achieve better accuracy in price, Greeks, and
optimal exercise strategies than LSM.

In G-LSM, the approximation of the terminal condition at $t_{k+1}$ is obtained
by discretizing the linear backward stochastic differential equation (BSDE)
for the CVF, c.f. Theorem 4.1, which was recently innovated in a deep neural
network-based method for American option pricing [6]. The idea of matching the
terminal condition has been widely applied in solving high-dimensional BSDEs
with deep neural networks (DNNs) [7, 12]. In practice, it involves computing
the gradient of the CVF, and in turn that of the basis functions in the ansatz
space (in addition to function evaluation). When $d\gg 1$, for a complicated
ansatz space, evaluating the derivatives at all time steps can be prohibitive.
We shall show that the extra cost is nearly negligible for the sparse Hermite
polynomial ansatz, and that the overall complexity of G-LSM with $N$ time
steps, $M$ sample paths, and $N_{b}$ basis functions is
$\mathcal{O}(NMN_{b})$, nearly identical with that for LSM. Numerical results
show that the accuracy of G-LSM is competitive with DNN-based methods for
dimensions up to $d=100$.

In theory, the CVF can be formulated as a smooth, high-dimensional function in
${L}^{2}_{\omega}(\mathbb{R}^{d})$ with a Gaussian weight function
$\omega(\mathbf{y})$ [25]. This regularity enables the use of normalized and
generalized Hermite polynomials, which form an orthonormal basis of
$L^{2}_{\omega}(\mathbb{R}^{d})$. Furthermore, drawing on the geometric
convergence rate of the hyperbolic cross approximation with Hermite
polynomials [17], we shall prove the global convergence of G-LSM using BSDE
technique, stochastic and Malliavin calculus, and establish an error bound in
terms of time step size, statistical error of the Monte Carlo approximation,
and the best approximation error in weighted Sobolev space
$H^{1}_{\omega}(\mathbb{R}^{d})$, c.f. Theorem 5.8. In sum, the algorithmic
development of G-LSM, its error analysis and extensive numerical evaluation
represent the main contributions of the present work.

Now we situate the present study in existing works. Currently, there are two
popular classes of methods to price American options in high dimensions: (i)
least-squares Monte Carlo-based (LSM) methods and (ii) DNN-based methods. The
LSM method has shown tremendous success for pricing American or Bermudan
options with more than one stochastic factors. The original LSM [15] uses
polynomials to approximate the CVF, and other choices have also been explored,
e.g., Gaussian process [16] and DNNs [14, 4]. Recently, LSM with the
hierarchical tensor train technique has been studied in [2], which
demonstrates the success of polynomial approximation for CVFs in very high
dimensions. The proposed G-LSM is a variant of LSM that incorporates gradient
information that comes nearly for free. Due to the excellent capability for
high-dimensional approximation of DNNs, several methods based on DNNs have
been proposed for pricing American or Bermudan options, based on optimal
stopping problem (parameterizing the stopping time by DNNs and then maximizing
the expected reward [3]), free boundary PDEs (parameterizing PDE solutions
with DNNs [21]), or BSDEs (parameterizing the solution pair of the associated
reflected BSDE [8] by DNNs [12]). Within the framework of BSDEs, Chen and Wan
[6] suggest approximating the difference of the CVF between adjacent time
steps by averaging several trained neural networks, which has a quadratic
complexity in the number of time steps. Wang et al. [23] extend the deep BSDE
method [7] from European option pricing to Bermudan one, with the loss
function being the variance of the initial value, and Gao et al. [10] analyze
its convergence. In comparison with DNN-based methods, G-LSM enjoys high
efficiency and robustness, which involves only least-squares problems and is
easy to implement.

The structure of this article is organized as follows. In Section 2, we
describe the mathematical framework of pricing and hedging high-dimensional
American or Bermudan options, and in Section 3, we recall several useful
properties of generalized Hermite polynomials and approximation with sparse
hyperbolic cross index set. Then in Section 4, we derive the main algorithm,
i.e., gradient-enhanced least squares Monte Carlo (G-LSM) method, and
establish its local and global error estimates in Section 5. In Section 6, we
present extensive numerical results including prices, Greeks, optimal stopping
time, and computing time. We also present a comparative study with existing
methods. Finally, we conclude in Section 7 with further discussions.

Throughout, bold and plain letters represent multi-variable and scalars,
respectively, and the capital and bold letters $\mathbf{S}$,
$\tilde{\mathbf{X}}$ and $\mathbf{W}$ denote random vectors. The notation
$\mathbf{a}\cdot\mathbf{b}$ denotes the dot product of two vectors
$\mathbf{a}$ and $\mathbf{b}$, $\operatorname{Tr}(A)$ the trace of a square
matrix $A$, and ⊤ vector transpose. The notation
$\nabla_{\mathbf{x}}f(t,\mathbf{x})$ and
$\operatorname{Hess}_{\mathbf{x}}f(t,\mathbf{x})$ denote respectively the
gradient and Hessian of $f$ with respect to $\mathbf{x}$. For multi-index
$\bm{\alpha}=(\alpha_{1},\ldots,\alpha_{d})^{\top}\in\mathbb{N}_{0}^{d}$,
$|\bm{\alpha}|=\alpha_{1}+\dots+\alpha_{d}$. For a positive-valued integrable
function $\omega:\mathbb{R}^{d}\to\mathbb{R}^{+}$, the weighted space
$L^{2}_{\omega}(\mathbb{R}^{d})$ is defined by

$L^{2}_{\omega}(\mathbb{R}^{d}):=\\{f:\mathbb{R}^{d}\to\mathbb{R}:\|f\|_{L^{2}_{\omega}(\mathbb{R}^{d})}<\infty\\},\quad\text{
with
}\|f\|_{L^{2}_{\omega}(\mathbb{R}^{d})}^{2}:=\int_{\mathbb{R}^{d}}f(\mathbf{x})^{2}\omega(\mathbf{x})\mathrm{d}\mathbf{x}.$

The weighted Sobolev space $H_{\omega}^{m}(\mathbb{R}^{d})$, $m\in\mathbb{N}$,
is defined by

$H_{\omega}^{m}(\mathbb{R}^{d}):=\\{f:\mathbb{R}^{d}\to\mathbb{R}:\|f\|_{H_{\omega}^{m}(\mathbb{R}^{d})}<\infty\\},\quad\text{with
}\|f\|_{H_{\omega}^{m}(\mathbb{R}^{d})}^{2}=\sum_{0\leq|\bm{\alpha}|\leq
m}\left\|\frac{\partial^{\bm{\alpha}}f}{\partial\mathbf{x}^{\bm{\alpha}}}\right\|_{L^{2}_{\omega}(\mathbb{R}^{d})}^{2}.$

##  2 Bermudan option pricing and hedging

Now we describe the valuation framework for American or Bermudan option
pricing and hedging.

###  2.1 Option pricing and Greeks

The fair price of American option $v^{A}(t)$ at time $t\in[0,T]$ is expressed
as the solution to the optimal stopping problem in a risk neutral probability
space $(\Omega,\mathcal{F},(\mathcal{F}_{t})_{0\leq t\leq T},\mathbb{Q})$,

$v^{A}(t)=\sup_{\tau_{t}\in[t,T]}\mathbb{E}[e^{-r(\tau_{t}-t)}g(\mathbf{S}_{\tau_{t}})|\mathcal{F}_{t}],$

where $\tau_{t}$ is an $\mathcal{F}_{t}$-stopping time, $T>0$ is the
expiration date, $(\mathbf{S}_{t})_{0\leq t\leq T}$ is a collection of
$d$-dimensional price processes, and $g(\mathbf{S}_{t})\in
L^{2}(\Omega,\mathcal{F}_{t},\mathbb{Q})$ is the payoff depending on the type
of the option.

Numerically, the price of Bermudan option is used to approximate the American
one. The Bermudan option can be exercised at finite discrete times
$0=t_{0}<t_{1}<\dots<t_{N}=T$ with $\Delta t:=t_{k+1}-t_{k}$ for all
$k=0,1,\dots,N-1$. Using dynamic programming principle or Snell envelope
theory [20, Section 1.8.4], the Bermudan price function ${v}_{t_{k}}$ at time
$t_{k}$ is given by the following backward induction:

(1) $\displaystyle v_{t_{N}}(\mathbf{s})$
$\displaystyle=g_{t_{N}}(\mathbf{s}),$ $\displaystyle v_{t_{k}}(\mathbf{s})$
$\displaystyle=\begin{cases}g_{t_{k}}(\mathbf{s}),\quad\text{ if
}g_{t_{k}}(\mathbf{s})\geq\tilde{c}_{t_{k}}(\mathbf{s}),\\\
\tilde{c}_{t_{k}}(\mathbf{s}),\quad\text{ if
}g_{t_{k}}(\mathbf{s})<\tilde{c}_{t_{k}}(\mathbf{s}),\end{cases}\quad\text{
for }k=N-1:-1:0,$

where $g_{t_{k}}$ is the discounted payoff (exercise value) function and
$\tilde{c}_{t_{k}}$ is the CVF, defined by

(2)
$\tilde{c}_{t_{k}}(\mathbf{s})=\mathbb{E}[v_{t_{k+1}}(\mathbf{S}_{t_{k+1}})|\mathbf{S}_{t_{k}}=\mathbf{s}].$

The options delta and gamma are defined to be the first and second order
derivatives of the price function $v_{t_{k}}$ with respect to the price of
underlying assets, i.e.,

(3) $\displaystyle\Delta_{t_{k}}$ $\displaystyle:=\nabla
v_{t_{k}}(\mathbf{s})=\left(\frac{\partial v_{t_{k}}(\mathbf{s})}{\partial
s_{1}},\dots,\frac{\partial v_{t_{k}}(\mathbf{s})}{\partial
s_{d}}\right)^{\top}\quad\mbox{and}\quad\Gamma^{ij}_{t_{k}}$
$\displaystyle:=\frac{\partial^{2}v_{t_{k}}(\mathbf{s})}{\partial
s_{j}\partial s_{i}}.$

We consider all exercise and continuation values of Bermudan option discounted
to the time $t=0$.

###  2.2 Multi-asset model and transformation

One of the most classical models for high-dimensional American option pricing
is the multi-asset Black-Scholes model. Under the risk-neutral probability
$\mathbb{Q}$, the prices of $d$ underlying assets,
$\mathbf{S}_{t}=\left(S^{1}_{t},\dots,S^{d}_{t}\right)^{\top}$, follow the
correlated geometric Brownian motions

(4)
$\mathrm{d}S^{i}_{t}=(r-\delta_{i})S^{i}_{t}\,\mathrm{d}t+\sigma_{i}S^{i}_{t}\,\mathrm{d}\tilde{W}^{i}_{t}\quad\text{
with }S^{i}_{0}=s^{i}_{0},\quad i=1,2,\dots,d,$

where $\tilde{W}^{i}_{t}$ are correlated Brownian motions with correlation
$\mathbb{E}[\mathrm{d}\tilde{W}^{i}_{t}\mathrm{d}\tilde{W}^{j}_{t}]=\rho_{ij}\,\mathrm{d}t$,
and $r$, $\delta_{i}$ and $\sigma_{i}$ are the riskless interest rate,
dividend yields, and volatility parameters, respectively. We denote the
correlation matrix by $P=(\rho_{ij})_{d\times d}$, the volatility matrix by
$\Sigma$ (which is a diagonal matrix with volatility $\sigma_{i}$ on the
diagonal), and write the dividend yields as a vector
$\bm{\delta}=[\delta_{1},\dots,\delta_{d}]^{\top}$. Using the spectral
decomposition $\Sigma P\Sigma^{\top}=Q\Lambda Q^{\top}$, the rotated log-price
$\tilde{\mathbf{X}}_{t}:=Q^{\top}\ln(\mathbf{S}_{t}./\mathbf{s}_{0})$
satisfies an independent Gaussian distribution

$\displaystyle\tilde{\mathbf{X}}_{t}\sim\mathcal{N}\left(Q^{\top}\left(r-\bm{\delta}-\frac{1}{2}\Sigma^{2}\mathbf{1}\right)t,\Lambda
t\right).$

Let
$\bm{\mu}=Q^{\top}\left(r-\bm{\delta}-\frac{1}{2}\Sigma^{2}\mathbf{1}\right)$
and $\lambda_{i}$ be the $i$-th diagonal element of $\Lambda$. This gives a
transformation between underlying asset prices $\mathbf{S}_{t}$ and
independent Brownian motions
$\mathbf{W}_{t}=[W^{1}_{t},\dots,W^{d}_{t}]^{\top}$, i.e.,

(5)
$\mathbf{S}_{t}=\mathbf{s}_{0}\odot\exp\left(Q\left(\bm{\mu}t+\sqrt{\Lambda}\mathbf{W}_{t}\right)\right),$

where $\odot$ denotes componentwise product.

From Eq. 5 and Eq. 2, the CVF $c_{t_{k}}$ with respect to the independent
Brownian motions is

(6)
$c_{t_{k}}(\mathbf{w})=\mathbb{E}[u_{t_{k+1}}(\mathbf{W}_{t_{k+1}})|\mathbf{W}_{t_{k}}=\mathbf{w}],$

where $u_{t_{k+1}}$ represents the value function at time $t_{k+1}$ with
respect to the independent Brownian motions. Our aim is to develop a gradient-
enhanced least squares Monte Carlo method that can efficiently approximate
$c_{t_{k}}$, and thus provide accurate prices and their derivatives.

##  3 Sparse Hermite polynomial expansion and gradient

Sparse polynomial chaos expansion can serve as a surrogate model of unknown
stochastic variables with finite second-order moments. The motivations of
using sparse Hermite polynomial expansion for pricing and hedging American
options are twofold:

  1. 1.

Let $\omega_{t}$, $t>0$, be the Gaussian density function defined by
$\omega_{t}(\mathbf{y}):=\prod_{j=1}^{d}\rho_{t}(y_{j})$ with
$\rho_{t}(y):=\frac{1}{\sqrt{2\pi t}}\exp(-\frac{y^{2}}{2t})$. The CVF
$c_{t_{k}}$ satisfies $c_{t_{k}}\in{L}^{2}_{\omega_{t_{k}}}(\mathbb{R}^{d})$,
since the payoff $g(\mathbf{S}_{t_{k}})\in
L^{2}(\Omega,\mathcal{F}_{t_{k}},\mathbb{Q})$ has a finite second-order
moment. Thus, as an orthonormal basis for
$L^{2}_{\omega_{t_{k}}}(\mathbb{R}^{d})$, the set of normalized and
generalized Hermite polynomials emerges as a natural choice.

  2. 2.

The CVF $c_{t_{k}}\in{H}^{m}_{\omega_{t_{k}}}(\mathbb{R}^{d})$ for any
positive integer $m$ by [25, Lemmas 4.2 and 4.3]. The smoothness of the
function implies efficient polynomial approximation.

Next, we recall one-dimensional normalized and generalized Hermite polynomials
$H_{n}^{(t)}(y)$ and the tensorized $d$-dimensional polynomials
$H_{\bm{\alpha}}^{(t)}(\mathbf{y})$, $\bm{\alpha}\in\mathbb{N}_{0}^{d}$. For
$d=1$, using the standard Hermite polynomials $H_{n}$, we define the $n$th-
order normalized and generalized Hermite polynomial $H_{n}^{(t)}$ by

$H_{n}^{(t)}(y):=\frac{H_{n}(\frac{y}{\sqrt{t}})}{\sqrt{n!}},\quad\text{ for
}n\in\mathbb{N}_{0},t>0,y\in\mathbb{R}.$

Then $\\{H_{n}^{(t)}\\}_{n\in\mathbb{N}_{0}}$ forms a complete orthonormal
basis for ${L}^{2}_{\rho_{t}}(\mathbb{R})$:

$\mathbb{E}\left[H_{n}^{(t)}(W_{t})H_{m}^{(t)}(W_{t})\right]=\int_{\mathbb{R}}H_{n}^{(t)}(y)H_{m}^{(t)}(y)\rho_{t}(y)\,\mathrm{d}y=\delta_{nm},$

with $\delta_{nm}$ being the Kronecker delta. For $d>1$, the tensorized
Hermite polynomial with the multi-index
$\bm{\alpha}=(\alpha_{1},\ldots,\alpha_{d})^{\top}\in\mathbb{N}_{0}^{d}$
defined by

$H_{\bm{\alpha}}^{(t)}(\mathbf{y}):=\prod_{j=1}^{d}H_{\alpha_{j}}^{(t)}(y_{j})$

forms a complete orthonormal basis for ${L}^{2}_{\omega_{t}}(\mathbb{R}^{d})$,
satisfying

$\mathbb{E}\left[H_{\bm{\alpha}}^{(t)}(\mathbf{W}_{t})H_{\bm{\gamma}}^{(t)}(\mathbf{W}_{t})\right]=\int_{\mathbb{R}^{d}}H_{\bm{\alpha}}^{(t)}(\mathbf{y})H_{\bm{\gamma}}^{(t)}(\mathbf{y})\omega_{t}(\bm{y})\,\mathrm{d}\mathbf{y}=\delta_{\bm{\alpha},\bm{\gamma}},\quad\bm{\alpha},\bm{\gamma}\in\mathbb{N}_{0}^{d}.$

Here,
$\delta_{\bm{\alpha},\bm{\gamma}}=\prod_{j=1}^{d}\delta_{\alpha_{j},\gamma_{j}}$
and $\omega_{t}(\bm{y})=\prod_{j=1}^{d}\rho_{t}(y_{j})$.

For any fixed multi-index set $I\subset\mathbb{N}_{0}^{d}$, the Hermite
polynomial ansatz space $P_{I,k}$ is defined by

(7) $\displaystyle
P_{I,k}:=\text{span}\left\\{H_{\bm{\alpha}}^{(t_{k})}:\bm{\alpha}\in
I\right\\}.$

Then we aim to approximate the CVF $c_{t_{k}}$ in $P_{I,k}$ by

$c_{t_{k}}(\mathbf{W}_{t_{k}})\approx\sum_{\bm{\alpha}\in
I}\beta_{\bm{\alpha}}H_{\bm{\alpha}}^{(t_{k})}(\mathbf{W}_{t_{k}}).$

It is well-known that computing polynomial approximations in high dimensions
suffers from the notorious curse of dimensionality with tensor product-type
multi-index sets. Fortunately, the smoothness of $c_{t_{k}}$ implies a fast
decay of the coefficients in the polynomial expansion. The large coefficients
usually occur in a lower multi-index set $I$, that is, if $\bm{\alpha}\in I$
and $\bm{\gamma}\leq\bm{\alpha}$, then $\bm{\gamma}\in I$ [1, Section 1.5]. To
circumvent the curse of dimensionality, the decay property of the expansion
coefficients can be exploited and a hyperbolic cross sparse index set can be
used to construct an approximation. The hyperbolic cross multi-index set $I$
with maximum order $p\in\mathbb{N}$ is defined by

(8)
$I:=\Big{\\{}\bm{\alpha}=(\alpha_{j})_{j=1}^{d}\in\mathbb{N}_{0}^{d}:\prod_{j=1}^{d}\max(\alpha_{j},1)\leq
p\Big{\\}},$

which has a cardinality $\mathcal{O}(p(\ln p)^{d-1})$. The best approximation
error of the sparse Hermite polynomial approximation with hyperbolic cross
index set was analyzed in [17]; see Section 5 for details.

Now, we introduce a property of derivatives of Hermite polynomials, which
plays an important role in reducing the computational cost. The first-order
derivative of the one-dimensional normalized and generalized Hermite
polynomial $H_{n}^{(t)}(y)$ satisfies

$\frac{\mathrm{d}}{\mathrm{d}y}\left(H_{n}^{(t)}(y)\right)=\sqrt{\frac{n}{t}}H_{n-1}^{(t)}(y).$

For the $d$-dimensional Hermite polynomial, we have

(9) $\frac{\partial}{\partial
y_{j}}\left(H_{\bm{\alpha}}^{(t)}(\mathbf{y})\right)=\sqrt{\frac{\alpha_{j}}{t}}H_{\bm{\alpha}-\mathbf{e}_{j}}^{(t)}(\mathbf{y}),$

where $\mathbf{e}_{j}$ is the $j$-th canonical basis vector. This implies that
for a lower multi-index set $I\subset\mathbb{N}_{0}^{d}$, we have
$\bm{\alpha}-\mathbf{e}_{j}\in I$ if $\bm{\alpha}\in I$. Thus, once the
evaluations of $H_{\bm{\alpha}}^{(t)}(\mathbf{y})$ for $\bm{\alpha}\in I$ are
available, the gradients of $H_{\bm{\alpha}}^{(t)}(\mathbf{y})$ for
$\bm{\alpha}\in I$ can be evaluated cheaply.

##  4 Algorithm and complexity

Now we derive the main methodology to approximate the CVF $c_{t_{k}}$ by
matching values of $u_{t_{k+1}}$, and analyze its computational complexity.
Below we abbreviate the notations $c_{t_{k}}$, $u_{t_{k}}$ and
$\mathbf{W}_{t_{k}}$ to $c_{k}$, $u_{k}$ and $\mathbf{W}_{k}$, etc, for
$k=0,1,\dots,N$.

###  4.1 Gradient-enhanced Least Squares (G-LS)

First we derive a linear backward stochastic differential equation (BSDE) for
the CVF $c_{k}$.

######  Theorem 4.1.

The CVF $c_{k}(\mathbf{W}_{k})$ satisfies the linear BSDE

$c_{k}(\mathbf{W}_{k})=u_{k+1}(\mathbf{W}_{k+1})-\int_{t_{k}}^{t_{k+1}}\nabla_{\mathbf{w}}f(t,\mathbf{W}_{t})\cdot\mathrm{d}\mathbf{W}_{t},$

where $f(t,\mathbf{w}),t\in[t_{k},t_{k+1}]$, is defined by
$f(t,\mathbf{w}):=\mathbb{E}[u_{k+1}(\mathbf{W}_{k+1})|\mathbf{W}_{t}=\mathbf{w}].$

######  Proof 4.2.

This is a direct consequence of martingale representation theorem [26, Theorem
2.5.2]. For the sake of completeness, we provide a brief proof. By the
Feynman-Kac formula, $f(t,\mathbf{w})$ satisfies a $d$-dimensional parabolic
PDE subject to a terminal condition

(10) $\left\\{\begin{aligned} \frac{\partial f}{\partial
t}(t,\mathbf{w})+\frac{1}{2}\operatorname{Tr}\left(\operatorname{Hess}_{\mathbf{w}}f(t,\mathbf{w})\right)&=0,\quad
t\in[t_{k},t_{k+1}),\\\
f(t_{k+1},\mathbf{w})&=u_{t_{k+1}}(\mathbf{w}).\end{aligned}\right.$

Let $Y_{t}:=f(t,\mathbf{W}_{t})$. By Itô’s formula, we have

$\mathrm{d}Y_{t}=\left(\frac{\partial f}{\partial
t}(t,\mathbf{W}_{t})+\frac{1}{2}\operatorname{Tr}\left(\operatorname{Hess}_{\mathbf{w}}f(t,\mathbf{W}_{t})\right)\right)\,\mathrm{d}t+\nabla_{\mathbf{w}}f(t,\mathbf{W}_{t})\cdot\mathrm{d}\mathbf{W}_{t}.$

In view of Eq. 10, the drift term vanishes. After taking the stochastic
integral and using the terminal condition in Eq. 10, we obtain the desired
assertion.

Next, we construct an approximation of the CVF ${c}_{k}$ by matching the
terminal condition over the interval $[t_{k},t_{k+1}]$. By Theorem 4.1, the
terminal value over $[t_{k},t_{k+1}]$ can be approximated by the Euler
discretization, i.e.,

(11) $\bar{u}_{k+1}(\mathbf{W}_{k+1})={c}^{{\rm
CLS}}_{k}(\mathbf{W}_{k})+\nabla{c}^{\operatorname{CLS}}_{k}(\mathbf{W}_{k})\cdot\Delta\mathbf{W}_{k},$

where $\Delta\mathbf{W}_{k}:=\mathbf{W}_{k+1}-\mathbf{W}_{k}$ is the Brownian
increment and ${c}^{\operatorname{CLS}}_{k}$ denotes an approximation to the
CVF ${c}_{k}$. Let $\hat{u}_{k+1}$ be the value function computed in the last
time step. Then using $P_{I,k}$ defined in Eq. 7 and Eq. 8 as the ansatz space
for ${c}^{\operatorname{CLS}}_{k}$, we solve for
${c}^{\operatorname{CLS}}_{k}$ using the least squares regression:

(12)
${c}_{k}^{\operatorname{CLS}}(\mathbf{W}_{k})=\operatorname*{argmin}_{\psi\in
P_{I,k}}E_{k}(\psi),$

with $E_{k}(\cdot):P_{I,k}\to\mathbb{R}^{+}$ being the quadratic loss (i.e.,
mean squared error) defined by

$\displaystyle
E_{k}(\psi):=\mathbb{E}\left[\left(\hat{u}_{k+1}(\mathbf{W}_{k+1})-\psi(\mathbf{W}_{k})-\nabla\psi(\mathbf{W}_{k})\cdot\Delta\mathbf{W}_{k}\right)^{2}\right].$

Finally, the numerical value $\hat{u}_{k}$ at time $t_{k}$ is updated to be
the discounted exercise or continuation value at $t_{k}$:

(13)
$\displaystyle\hat{u}_{k}(\mathbf{W}_{k})=\begin{cases}g_{k}(\mathbf{S}_{k}),\quad&\text{
if exercise},\\\ c_{k}^{\operatorname{CLS}}(\mathbf{W}_{k}),\quad&\text{ if
continue}.\end{cases}$

Since the option is only profitable when exercised in the in-the-money region
$\Omega_{\operatorname{ITM}}=\\{\mathbf{s}\in\mathbb{R}^{d}:g_{k}(\mathbf{s})>0\\}$,
we make the decision of exercising the option when
$c_{k}^{\operatorname{CLS}}(\mathbf{W}_{k})<g_{k}(\mathbf{S}_{k})$ and
$\mathbf{S}_{k}\in\Omega_{\operatorname{ITM}}$; otherwise, the option will be
continued.

###  4.2 Gradient-enhanced Least Squares Monte Carlo (G-LSM)

Now we derive the methodology based on the Monte Carlo method to solve Eq. 12
numerically and analyze its computational complexity. Let $N_{b}=|I|$ and let
$\\{\phi_{n}^{k}(\mathbf{W}_{k})\\}_{n=1}^{N_{b}}$ be the set of Hermite
polynomials in a scalar-indexed form. Then any function $\psi\in P_{I,k}$ can
be expressed as

$\displaystyle\psi(\mathbf{W}_{k})=\sum_{n=1}^{N_{b}}\beta_{n}\phi_{n}^{k}(\mathbf{W}_{k}),$

with its gradient given by

(14)
$\nabla\psi(\mathbf{W}_{k})=\sum_{n=1}^{N_{b}}\beta_{n}\nabla\phi_{n}^{k}(\mathbf{W}_{k}).$

In practice, the continuous least squares problem Eq. 12 is solved by
minimizing its Monte Carlo approximation:

(15) $\hat{c}_{k}:=\operatorname*{argmin}_{\psi\in
P_{I,k}}\frac{1}{M}\sum_{m=1}^{M}\big{(}\hat{u}_{k+1}(\mathbf{W}^{m}_{k+1})-\psi(\mathbf{W}^{m}_{k})-\nabla\psi(\mathbf{W}^{m}_{k})\cdot\Delta\mathbf{W}^{m}_{k}\big{)}^{2},$

where $\\{\mathbf{W}^{m}_{k}\\}_{m=1}^{M}$ are $M$ independent paths of the
Gaussian random process $\mathbf{W}_{k}$ and
$\Delta\mathbf{W}^{m}_{k}:=\mathbf{W}^{m}_{k+1}-\mathbf{W}^{m}_{k}$ is the
$m$th path of increment $\Delta\mathbf{W}_{k}$. Let
$A_{k}\in\mathbb{R}^{M\times N_{b}}$, $\bm{\beta}_{k}\in\mathbb{R}^{N_{b}}$
and $\mathbf{\hat{u}}_{k+1}\in\mathbb{R}^{M}$ with their components defined by

$\displaystyle(A_{k})_{mn}$
$\displaystyle=\phi_{n}^{k}(\mathbf{W}^{m}_{k})+\nabla\phi_{n}^{k}(\mathbf{W}^{m}_{k})\cdot\Delta\mathbf{W}^{m}_{k},$
$\displaystyle(\bm{\beta}_{k})_{n}$
$\displaystyle=\beta_{n},\quad(\mathbf{\hat{u}}_{k+1})_{m}=\hat{u}_{k+1}(\mathbf{W}^{m}_{k+1})$

for $m=1,\dots,M$, $n=1,\dots,N_{b}$.

Then finding the optimal polynomial in Eq. 15 amounts to solving the classical
least squares problem

$\bm{\beta}_{k}=\operatorname*{argmin}_{\bm{\beta}}\|A_{k}\bm{\beta}-\mathbf{\hat{u}}_{k+1}\|_{2}^{2}=(A_{k}^{\top}A_{k})^{-1}A_{k}^{\top}\mathbf{\hat{u}}_{k+1}.$

The proposed algorithm is summarized in Algorithm 1.

######  Remark 4.3.

Note that the well-known least squares Monte Carlo (LSM) method [15, 22] also
approximates the CVF with a finite number of basis functions. The orthonormal
Hermite polynomials is one possible choice. However, LSM computes the
coefficients by projecting the value, $u_{t_{k+1}}(\mathbf{W}_{t_{k+1}})$ in
Eq. 6, onto a finite-dimensional space spanned by the basis functions due to
the projection nature of conditional expectation:

$\hat{c}_{k}^{LSM}:=\operatorname*{argmin}_{\psi\in
P_{I,k}}\frac{1}{M}\sum_{m=1}^{M}\big{(}\hat{u}_{k+1}(\mathbf{W}^{m}_{k+1})-\psi(\mathbf{W}^{m}_{k})\big{)}^{2}.$

In contrast, G-LSM computes the coefficients by matching the approximate and
exact value of $u_{t_{k+1}}(\mathbf{W}_{t_{k+1}})$, see Eq. 15. Their
numerical performance will be compared in Section 6.

Algorithm 1 G-LSM

0: Market parameters: $\mathbf{S}_{0},r,\delta_{i},\sigma_{i},P$ Option
parameters: payoff function $g(\mathbf{s})$ Algorithm parameters: $N,M,p$

0: Option price $v_{0}$

1: Compute the hyperbolic cross multi-index set $I$ with maximum order $p$

2: Generate $M$ sample paths

3: Initialize values
$\mathbf{\hat{u}}_{N}=(e^{-rT}g(\mathbf{S}_{N}^{m}))_{m=1}^{M}$

4: Initialize stopping times $\tau_{\star}=T$

5: for $k=N-1:-1:1$ do

6: $\Phi_{k}=$ basis matrix$(\mathbf{W}_{k},I)$ with
$(\Phi_{k})_{m,n}=\phi_{n}^{k}(\mathbf{W}_{k}^{m})$

7: Compute matrix $A_{k}$ with Algorithm 2

8: Solve system of linear equations: $A_{k}\bm{\beta}_{k}=\mathbf{u}_{k+1}$

9: Update $(\mathbf{\hat{u}}_{k})_{m}=\begin{cases}e^{-rk\Delta
t}g(\mathbf{S}^{m}_{k}),&\text{if exercise}\\\
(\Phi_{k}\bm{\beta}_{k})_{m},&\text{if continue}\end{cases}\text{ and
}\tau_{\star}=\begin{cases}k\Delta t,&\text{if exercise}\\\
\tau_{\star},&\text{if continue}\end{cases}$

10: end for

11:
$v_{0}=\max\left\\{\frac{1}{M}\sum_{m=1}^{M}e^{-r\tau_{\star}^{m}}g(\mathbf{S}^{m}_{\tau_{\star}^{m}}),g(\mathbf{S}_{0})\right\\}$

For the efficient computation of the matrix $A_{k}$, using Eq. 9, i.e.,

$\frac{\partial\phi_{n}^{k}}{\partial
w_{j}}(\mathbf{W}^{m}_{k})=\sqrt{\frac{\alpha_{j}}{t_{k}}}\phi_{n^{\prime}}^{k}(\mathbf{W}^{m}_{k})\quad\text{
for
}\phi_{n}^{k}=H_{\bm{\alpha}}^{(t_{k})},\phi_{n^{\prime}}^{k}=H_{\bm{\alpha}-\mathbf{e}_{j}}^{(t_{k})},$

we can compute the matrix $A_{k}$ from the basis matrix $\Phi_{k}$ and the
increment $\Delta\mathbf{W}_{k}$. This routine is summarized in Algorithm 2,
where $(A_{k})_{n}$ represents the $n$-th column of the matrix $A_{k}$.

Finally, we analyze the computational complexity of Algorithm 1. We employ
backward induction Eq. 1 with $N$ time steps to price American options, which
has a complexity $\mathcal{O}(N)$. For each fixed $t_{k}$, the computation
consists of steps 6-9 in Algorithm 1. Step 6 has a complexity
$\mathcal{O}(MN_{b})$, when using $M$ samples and $N_{b}$ basis functions.
Step 7 is detailed in Algorithm 2, which has a linear complexity in the number
of nonzero $\alpha_{j}$ for all $\bm{\alpha}\in I$ and $j=1,\dots,d$,

$\displaystyle\|I\|_{0}:=\sum_{\bm{\alpha}\in
I}\\#\left\\{j=1,\cdots,d:\alpha_{j}\neq 0\right\\}.$

Note that $\|I\|_{0}$ can be calculated by $(N_{b}-N_{b,p-1})d$ for maximum
order $p$, where $N_{b,p-1}$ represents the number of basis functions with
maximum order $p-1$. Figure 1 shows that $\|I\|_{0}$ scales nearly linearly
with respect to $N_{b}$. Thus, the complexity of Algorithm 2 is linear in
$N_{b}$. Step 8 involves solving a linear system of size $M\times N_{b}$.
Using an iterative solver, the cost is $\mathcal{O}(MN_{b})$ if the matrix $A$
is well-conditioned. Step 9 involves matrix-vector multiplication, which has a
complexity $\mathcal{O}(MN_{b})$. Hence, with a fixed sample size $M$ and
number $N$ of time steps, the total cost of the proposed G-LSM is
$\mathcal{O}(NMN_{b})$. We will numerically verify the complexity analysis in
Section 6.3.

Algorithm 2 Compute the matrix $A_{k}$ using $I,\Phi_{k}$ and
$\Delta\mathbf{W}_{k}$

1: Initialize $A_{k}=\Phi_{k}$

2: for $j=1:d$ do

3: for $\bm{\alpha}\in I$ do

4: if $\alpha_{j}\geq 1$ then

5: $(A_{k})_{n}=(A_{k})_{n}+(\Delta
W_{k})_{j}\odot(\Phi_{k})_{n^{\prime}}\sqrt{\frac{\alpha_{j}}{t_{k}}}$.

6: end if

7: end for

8: end for

|  |
---|---|---
(a) $p=4$ | (b) $p=6$ | (c) $p=10$
Figure 1: $\|I\|_{0}$ scales linearly with respect to $N_{b}$.

###  4.3 Computing deltas

Now we present the approximation of deltas for constructing hedging strategies
based on Algorithm 1. In the backward induction loop, Eq. 1 and Eq. 3 imply
the deltas at time $t_{k}$ is given by

$\frac{\partial v_{k}(\mathbf{s})}{\partial s_{j}}=\left\\{\begin{aligned}
\frac{\partial g_{k}(\mathbf{s})}{\partial s_{j}},\text{ if exercise},\\\
\frac{\partial{\tilde{c}}_{k}(\mathbf{s})}{\partial s_{j}},\text{ if
continue},\end{aligned}\right.\quad j=1,\dots,d.$

In the continuation region, $\tilde{c}_{k}(\mathbf{s})=c_{k}(\mathbf{w})$ with
$\mathbf{s}=\mathbf{s}_{0}\odot\exp(Q(\mu t_{k}+\sqrt{\Lambda}\mathbf{w}))$.
Once we have obtained the coefficients $\bm{\beta}_{k}$ of the expansion, the
gradient $\nabla c_{k}(\mathbf{w})$ of the continuation value
$c_{k}(\mathbf{w})$ is approximated by $\nabla\hat{c}_{k}(\mathbf{w})$ in Eq.
14. Hence, we calculate the deltas in the continuation region by

$\frac{\partial{\tilde{c}}_{k}(\mathbf{s})}{\partial
s_{j}}\approx\sum_{i=1}^{d}\frac{\partial\hat{c}_{k}(\mathbf{w})}{\partial
w_{i}}\frac{\partial w_{i}}{\partial s_{j}},\quad k=1,\dots,{N-1}.$

The gammas $\Gamma_{t_{k}}^{ij}$ can be approximated in a similar way.

Particularly, the Greeks at time $t=0$ involves solving an additional linear
system, which arises from minimizing the
$L_{\omega_{t_{1}}}^{2}(\mathbb{R}^{d})$ error of approximating
$u_{1}(\mathbf{W}_{1})$. We employ a slightly different ansatz from the case
of $t_{k}>0$, since the expansions in $H_{\bm{\alpha}}^{(0)}$ is not
applicable. Instead, we take

$\hat{c}_{0}(\mathbf{w})=\sum_{\bm{\alpha}\in
I}\beta_{\bm{\alpha}}H_{\bm{\alpha}}^{(\Delta t)}(\mathbf{w}).$

Since $c_{0}(\mathbf{W}_{0})$ is a deterministic value given
$\mathbf{W}_{0}=\mathbf{0}$, we consider the Euler approximation

$u_{1}(\mathbf{W}_{1})\approx\hat{c}_{0}(\mathbf{0})+\nabla\hat{c}_{0}(\mathbf{0})\cdot\mathbf{W}_{1}.$

Similar to Eq. 15, the coefficients $\beta_{\bm{\alpha}}$ can be approximated
by solving an $M\times N_{b}$ linear system. After that, the delta at $t=0$ is
approximated by

$\frac{\partial{\tilde{c}}_{0}(\mathbf{s}_{0})}{\partial
s_{j}}\approx\sum_{i=1}^{d}\frac{\partial\hat{c}_{0}(\mathbf{0})}{\partial
w_{i}}\frac{\partial w_{i}}{\partial s_{j}}(\mathbf{s}_{0}),\quad\mbox{with
}\mathbf{s}=\mathbf{s}_{0}\odot\exp(Q\sqrt{\Lambda}\mathbf{w}).$

##  5 Convergence analysis

Now we analyze the convergence of Algorithm 1. We assume the Lipschitz
continuity of the discounted payoff function $g_{k}(\cdot)$. The discounted
payoff function $g_{k}(\cdot)$ is $L$-Lipschitz continuous: for
$k=0,1,\dots,N$

(16) $|g_{k}(\mathbf{s})-g_{k}(\mathbf{s}^{\prime})|\leq
L\|\mathbf{s}-\mathbf{s}^{\prime}\|,\quad\forall\mathbf{s},\mathbf{s}^{\prime}\in\mathbb{R}^{d}.$

Recall that $\hat{u}_{k}(\mathbf{w})$ is the numerical value function at time
$t_{k}$ computed by Algorithm 1, which approximates the exact value
${u}_{k}(\mathbf{w})$ for $k=0,1,\dots,N-1$. We aim to establish an upper
bound for $\max_{0\leq k\leq
N-1}\|u_{k}-\hat{u}_{k}\|_{{L}^{2}_{\omega_{k}}(\mathbb{R}^{d})}^{2}$.

###  5.1 One-step error estimation

First, we analyze the one step error over the interval $[t_{k},t_{k+1}]$ for
$k=0,1,\dots,N-1$, i.e., the numerical value function
$\hat{u}_{k}(\mathbf{w})$ as an approximation to the exact one
${u}_{k}(\mathbf{w})$. Given the previous value function $\hat{u}_{k+1}$, we
define the current CVF $\check{c}_{k}(\mathbf{W}_{k})$ by

(17)
$\check{c}_{k}(\mathbf{W}_{k})=\mathbb{E}[\hat{u}_{k+1}(\mathbf{W}_{k+1})|\mathbf{W}_{k}].$

The classic backward Euler scheme for BSDE approximate $\nabla\check{c}_{k}$
by

(18) $\check{Z}_{k}:=\frac{1}{\Delta
t}\mathbb{E}[\hat{u}_{k+1}(\mathbf{W}_{k+1})\Delta\mathbf{W}_{k}|\mathbf{W}_{k}].$

With the sparse Hermite polynomial ansatz space $P_{I,k}$, let $c_{k}^{*}\in
P_{I,k}$ be the best approximation to $\check{c}_{k}$ in
$H^{1}_{\omega_{k}}(\mathbb{R}^{d})$ defined by

(19) $c^{*}_{k}:=\operatorname*{argmin}_{\psi\in
P_{I,k}}\|\check{c}_{k}-\psi\|_{H^{1}_{\omega_{k}}(\mathbb{R}^{d})}.$

Also we define the best approximation error $\mathcal{E}_{k}^{\rm best}$ by

(20) $\mathcal{E}_{k}^{\rm
best}:=\|\check{c}_{k}-c_{k}^{*}\|_{H^{1}_{\omega_{k}}(\mathbb{R}^{d})}^{2}.$

Next, we denote the statistical error of solving the discrete least squares
problem Eq. 15 by

(21) $\displaystyle\mathcal{E}^{\rm
stat}_{k}:=\|\hat{c}_{k}-c_{k}^{\operatorname{CLS}}\|_{{L}^{2}_{\omega_{k}}(\mathbb{R}^{d})}^{2},$

where $c_{k}^{\rm CLS}$ solves the continuous least squares problem Eq. 12 for
$k=0,1,\dots,N-1$. In Theorem 5.4 below, we derive the convergence of
$c_{k}^{\operatorname{CLS}}$ to the exact $c_{k}$ provided that the previous
approximation error
$\|u_{k+1}-\hat{u}_{k+1}\|_{L_{\omega_{k+1}}^{2}(\mathbb{R}^{d})}^{2}$,
$\mathcal{E}_{k}^{\rm best}$ in Eq. 20, and time step size $\Delta t$ are
small. We first provide in Lemma 5.1 an estimate of $\mathcal{E}_{k}^{\rm
best}$ in Eq. 20, which is a direct application of [17, Theorem 4.2].

######  Lemma 5.1.

For the hyperbolic cross index set $I$ with maximum order $p$ defined in Eq.
8, given $\check{c}_{k}\in\mathcal{K}^{m}(\mathbb{R}^{d},\omega_{k})$, we have

$\mathcal{E}_{k}^{\rm best}\leq
C(m,d)\frac{1}{p^{m-1}}|\check{c}_{k}|_{\mathcal{K}^{m}_{\omega_{k}}(\mathbb{R}^{d})}^{2},$

where $\mathcal{K}^{m}_{\omega_{k}}(\mathbb{R}^{d})$ is the weighted Koborov-
type space defined by

$\mathcal{K}^{m}_{\omega_{k}}(\mathbb{R}^{d})=\\{u:\partial^{\bm{\alpha}}u\in{L}^{2}_{\omega_{k}}(\mathbb{R}^{d}),0\leq|\bm{\alpha}|_{\infty}\leq
m\\},$

and $|\cdot|_{\mathcal{K}^{m}_{\omega_{k}}(\mathbb{R}^{d})}$ is the seminorm
defined by

$|u|_{\mathcal{K}^{m}_{\omega_{k}}(\mathbb{R}^{d})}=\left(\sum_{|\bm{\alpha}|_{\infty}=m}\|\partial^{\bm{\alpha}}u\|_{{L}^{2}_{\omega_{k}}(\mathbb{R}^{d})}^{2}\right)^{1/2}.$

The next result provides a one-step error estimate of the $Z$ component in the
BSDE by the backward Euler scheme Eq. 18. The proof is inspired by the idea of
[11, pages 816-817]. By martingale representation theorem, the definition of
$\check{c}_{k}$ Eq. 17 implies the existence of a square integrable process
$\tilde{Z}_{s}$, $t_{k}\leq s\leq t_{k+1}$, such that

(22)
$\hat{u}_{k+1}(\mathbf{W}_{k+1})=\check{c}_{k}(\mathbf{W}_{k})+\int_{t_{k}}^{t_{k+1}}\tilde{Z}_{t}\cdot\mathrm{d}\mathbf{W}_{t}.$

######  Lemma 5.2.

Let $\check{Z}_{k}$ and $\tilde{Z}_{t}$ be defined in Eq. 18 and Eq. 22. Then
for some constant $C_{1}>0$,

$\mathbb{E}\left[(\check{Z}_{k}-\tilde{Z}_{t_{k}})^{2}\right]\leq C_{1}\Delta
t^{2}.$

######  Proof 5.3.

By the integration by parts formula of Malliavin calculus, we obtain

$\check{Z}_{k}=\frac{1}{\Delta
t}\mathbb{E}\left[\hat{u}_{k+1}(\mathbf{W}_{k+1})\Delta\mathbf{W}_{k}|\mathbf{W}_{k}\right]=\frac{1}{\Delta
t}\mathbb{E}\left[\int_{t_{k}}^{t_{k+1}}D_{t}\mathbf{W}_{k+1}\nabla\hat{u}_{k+1}(\mathbf{W}_{k+1})\,\mathrm{d}t\bigg{|}\mathbf{W}_{k}\right],$

where $D_{t}\mathbf{W}_{k+1}$ is an identity matrix of size $d\times d$.
Together with the identity

$\tilde{Z}_{t_{k}}=\nabla\check{c}_{k}(\mathbf{W}_{k})=\frac{1}{\Delta
t}\int_{t_{k}}^{t_{k+1}}\nabla\check{c}_{k}(\mathbf{W}_{k})\,\mathrm{d}t,$

it further leads to

$\displaystyle\check{Z}_{k}-\tilde{Z}_{t_{k}}$ $\displaystyle=\frac{1}{\Delta
t}\mathbb{E}\left[\int_{t_{k}}^{t_{k+1}}\nabla\hat{u}_{k+1}(\mathbf{W}_{k+1})-\nabla\check{c}_{k}(\mathbf{W}_{k})\,\mathrm{d}s\bigg{|}\mathbf{W}_{k}\right].$

Then, by [11, Equation (26)], we obtain

$\check{Z}_{k}-\tilde{Z}_{t_{k}}=\int_{t_{k}}^{t_{k+1}}\mathbb{E}\left[G(t,\mathbf{W}_{t})|\mathbf{W}_{k}\right]\,\mathrm{d}t,$

for a bounded function $G$. Hence, there holds
$|\check{Z}_{k}-\tilde{Z}_{t_{k}}|=\mathcal{O}(\Delta t)$, and there exists a
constant $C_{1}>0$ such that
$\mathbb{E}\left[(\check{Z}_{k}-\tilde{Z}_{t_{k}})^{2})\right]\leq C_{1}\Delta
t^{2}$.

Next, we provide an error estimate of solving the continuous least squares
problem Eq. 12 in terms of the error of the previous value function, the best
approximation error in the sparse Hermite polynomial ansatz space Eq. 20 and
the time step size $\Delta t$. The proof is inspired by the foundational work
[12].

######  Theorem 5.4.

For $\Delta t\leq 1/4$, with the constant $C_{1}>0$ in Lemma 5.2, there holds

$\|c_{k}-c_{k}^{\operatorname{CLS}}\|_{{L}^{2}_{\omega_{k}}(\mathbb{R}^{d})}^{2}\leq(1+4\Delta
t)\|u_{k+1}-\hat{u}_{k+1}\|_{{L}^{2}_{\omega_{k+1}}(\mathbb{R}^{d})}^{2}+\tfrac{1}{2\Delta
t}\mathcal{E}_{k}^{\rm best}+C_{1}\Delta t^{2}.$

######  Proof 5.5.

By the inequality $(a+b)^{2}\leq(1+4\Delta t)a^{2}+(1+\frac{1}{4\Delta
t})b^{2}$, we have

$\displaystyle\|c_{k}-c_{k}^{\operatorname{CLS}}\|_{{L}^{2}_{\omega_{k}}(\mathbb{R}^{d})}^{2}$
$\displaystyle\leq(1+4\Delta
t)\|c_{k}-\check{c}_{k}\|_{{L}^{2}_{\omega_{k}}(\mathbb{R}^{d})}^{2}+(1+\tfrac{1}{4\Delta
t})\|\check{c}_{k}-c_{k}^{\rm
CLS}\|_{{L}^{2}_{\omega_{k}}(\mathbb{R}^{d})}^{2}$
$\displaystyle\leq(1+4\Delta
t)\|c_{k}-\check{c}_{k}\|_{{L}^{2}_{\omega_{k}}(\mathbb{R}^{d})}^{2}+\tfrac{1}{2\Delta
t}\|\check{c}_{k}-c_{k}^{\operatorname{CLS}}\|_{{L}^{2}_{\omega_{k}}(\mathbb{R}^{d})}^{2},$

since $\Delta t\leq 1/4$. For the first term, the definitions Eq. 6 and Eq.
17, the tower property, and Jensen’s inequality lead to

$\displaystyle\mathbb{E}\left[\left(c_{k}(\mathbf{W}_{k})-\check{c}_{k}(\mathbf{W}_{k})\right)^{2}\right]$
$\displaystyle=\mathbb{E}\left[\left(\mathbb{E}\left[u_{k+1}(\mathbf{W}_{k+1})-\hat{u}_{k+1}(\mathbf{W}_{k+1})|\mathbf{W}_{k}\right]\right)^{2}\right]$
$\displaystyle\leq\mathbb{E}\left[\left(u_{k+1}(\mathbf{W}_{k+1})-\hat{u}_{k+1}(\mathbf{W}_{k+1})\right)^{2}\right]=\|u_{k+1}-\hat{u}_{k+1}\|_{{L}^{2}_{\omega_{k+1}}(\mathbb{R}^{d})}^{2}.$

Hence, it remains to show

(23)
$\|\check{c}_{k}-c_{k}^{\operatorname{CLS}}\|_{{L}^{2}_{\omega_{k}}(\mathbb{R}^{d})}^{2}\leq\mathcal{E}_{k}^{\rm
best}+2C_{1}\Delta t^{3}.$

By plugging Eq. 22 into the quadratic loss function $E_{k}(\cdot)$ in Eq. 12,
we obtain

$\displaystyle E_{k}(\psi)$
$\displaystyle=\mathbb{E}\left[\left(\check{c}_{k}(\mathbf{W}_{k})-\psi(\mathbf{W}_{k})+\int_{t_{k}}^{t_{k+1}}\tilde{Z}_{s}\cdot\mathrm{d}\mathbf{W}(s)-\nabla\psi(\mathbf{W}_{k})\cdot\Delta\mathbf{W}_{k}\right)^{2}\right]$
$\displaystyle=\mathbb{E}\left[\left(\check{c}_{k}(\mathbf{W}_{k})-\psi(\mathbf{W}_{k})\right)^{2}\right]+\mathbb{E}\left[\left(\int_{t_{k}}^{t_{k+1}}\tilde{Z}_{s}\cdot\mathrm{d}\mathbf{W}(s)-\nabla\psi(\mathbf{W}_{k})\cdot\Delta\mathbf{W}_{k}\right)^{2}\right].$

Now Itô isometry implies the identity

$\displaystyle\mathbb{E}\left[\left(\int_{t_{k}}^{t_{k+1}}\tilde{Z}_{s}\cdot\mathrm{d}\mathbf{W}(s)-\nabla\psi(\mathbf{W}_{k})\cdot\Delta\mathbf{W}_{k}\right)^{2}\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\left(\int_{t_{k}}^{t_{k+1}}\tilde{Z}_{s}\cdot\mathrm{d}\mathbf{W}(s)-\int_{t_{k}}^{t_{k+1}}\check{Z}_{k}\cdot\mathrm{d}\mathbf{W}(s)+\check{Z}_{k}\cdot\Delta\mathbf{W}_{k}-\nabla\psi(\mathbf{W}_{k})\cdot\Delta\mathbf{W}_{k}\right)^{2}\right]$
$\displaystyle=$
$\displaystyle\mathbb{E}\left[\int_{t_{k}}^{t_{k+1}}|\tilde{Z}_{s}-\check{Z}_{k}|^{2}\,\mathrm{d}s\right]+\Delta
t\mathbb{E}\left[\left|\check{Z}_{k}-\nabla\psi(\mathbf{W}_{k})\right|^{2}\right]$
$\displaystyle+2\mathbb{E}\left[\int_{t_{k}}^{t_{k+1}}(\tilde{Z}_{s}-\check{Z}_{k})\,\mathrm{d}s\right]\cdot\mathbb{E}[\check{Z}_{k}-\nabla\psi(\mathbf{W}_{k})].$

The definitions of $\check{Z}_{k}$ and $\tilde{Z}_{s}$ in Eq. 18 and Eq. 22,
together with Itô isometry, yield

$\displaystyle\check{Z}_{k}=\frac{1}{\Delta
t}\mathbb{E}\left[\int_{t_{k}}^{t_{k+1}}\tilde{Z}_{s}\,\mathrm{d}s\bigg{|}\mathbf{W}_{k}\right],$

which implies

$\displaystyle\mathbb{E}\left[\int_{t_{k}}^{t_{k+1}}(\tilde{Z}_{s}-\check{Z}_{k})\,\mathrm{d}s\right]=0.$

Hence, we can express the loss function $E_{k}(\psi)$ as

$E_{k}(\psi)=\mathbb{E}\left[\left(\check{c}_{k}(\mathbf{W}_{k})-\psi(\mathbf{W}_{k})\right)^{2}\right]+\mathbb{E}\left[\int_{t_{k}}^{t_{k+1}}|\tilde{Z}_{s}-\check{Z}_{k}|^{2}\,\mathrm{d}s\right]+\Delta
t\mathbb{E}\left[\left|\check{Z}_{k}-\nabla\psi(\mathbf{W}_{k})\right|^{2}\right].$

The equality Eq. 12 implies that $E_{k}(c_{k}^{\rm CLS})\leq E_{k}(\psi)$ for
all $\psi\in P_{I,k}$. Taking $\psi:=c_{k}^{*}$ as defined in Eq. 19, we
obtain

$\displaystyle\quad\mathbb{E}\left[\left(\check{c}_{k}(\mathbf{W}_{k})-c_{k}^{\rm
CLS}(\mathbf{W}_{k})\right)^{2}\right]+\Delta
t\mathbb{E}\left[\left(\check{Z}_{k}-\nabla c_{k}^{\rm
CLS}(\mathbf{W}_{k})\right)^{2}\right]$
$\displaystyle\leq\mathbb{E}\left[\left(\check{c}_{k}(\mathbf{W}_{k})-c_{k}^{*}(\mathbf{W}_{k})\right)^{2}\right]+\Delta
t\mathbb{E}\left[\left(\check{Z}_{k}-\nabla
c_{k}^{*}(\mathbf{W}_{k})\right)^{2}\right].$

Note that $\tilde{Z}_{t_{k}}=\nabla\check{c}_{k}(\mathbf{W}_{k})$. By
subtracting and adding $\tilde{Z}_{t_{k}}$, we have

$\displaystyle\quad\|\check{c}_{k}-c_{k}^{\operatorname{CLS}}\|_{{L}^{2}_{\omega_{k}}(\mathbb{R}^{d})}^{2}$
$\displaystyle\leq\|\check{c}_{k}-c_{k}^{*}\|_{{L}^{2}_{\omega_{k}}(\mathbb{R}^{d})}^{2}+2\Delta
t\mathbb{E}\left[(\check{Z}_{k}-\tilde{Z}_{t_{k}})^{2}\right]+2\Delta
t\mathbb{E}\left[|\nabla\check{c}_{k}(\mathbf{W}_{k})-\nabla
c_{k}^{*}(\mathbf{W}_{k})|^{2}\right]$
$\displaystyle\leq\|\check{c}_{k}-c_{k}^{*}\|_{H^{1}_{\omega_{k}}(\mathbb{R}^{d})}^{2}+2C_{1}\Delta
t^{3},$

where the last inequality follows from the condition $2\Delta t\leq 1$ and
Lemma 5.2. This shows the estimate Eq. 23, and completes the proof.

Now, we can give the one-step error propagation of
$\|u_{k}-\hat{u}_{k}\|_{{L}^{2}_{\omega_{k}}(\mathbb{R}^{d})}^{2}$ given
$\|u_{k+1}-\hat{u}_{k+1}\|_{{L}^{2}_{\omega_{k+1}}(\mathbb{R}^{d})}^{2}$.

######  Corollary 5.6.

For $\Delta t\leq 1/6$, there holds

$\|u_{k}-\hat{u}_{k}\|_{{L}^{2}_{\omega_{k}}(\mathbb{R}^{d})}^{2}\leq(1+6\Delta
t)\|u_{k+1}-\hat{u}_{k+1}\|_{{L}^{2}_{\omega_{k+1}}(\mathbb{R}^{d})}^{2}+\left(1+\frac{1}{2\Delta
t}\right)\mathcal{E}_{k}^{\rm best}+\frac{6}{5}C_{1}\Delta
t^{2}+\frac{6}{5\Delta t}\mathcal{E}_{k}^{\rm stat},$

where $\mathcal{E}_{k}^{\rm best}$ and $\mathcal{E}_{k}^{\rm stat}$ are
defined in Eq. 20 and Eq. 21, and $C_{1}>0$ is the constant in Lemma 5.2.

######  Proof 5.7.

First, using inequality $(a+b)^{2}\geq(1-\Delta t)a^{2}-\frac{1}{\Delta
t}b^{2}$ for any $a,b\in\mathbb{R}$, leads to

$\|c_{k}-c_{k}^{\operatorname{CLS}}\|_{{L}^{2}_{\omega_{k}}(\mathbb{R}^{d})}^{2}\geq(1-\Delta
t)\|c_{k}-\hat{c}_{k}\|_{{L}^{2}_{\omega_{k}}(\mathbb{R}^{d})}^{2}-\frac{1}{\Delta
t}\|\hat{c}_{k}-c_{k}^{CLS}\|_{{L}^{2}_{\omega_{k}}(\mathbb{R}^{d})}^{2}.$

This and Theorem 5.4 yield

$\displaystyle\quad\|c_{k}-\hat{c}_{k}\|_{{L}^{2}_{\omega_{k}}(\mathbb{R}^{d})}^{2}$
$\displaystyle\leq\frac{1+4\Delta t}{1-\Delta
t}\|u_{k+1}-\hat{u}_{k+1}\|_{{L}^{2}_{\omega_{k+1}}(\mathbb{R}^{d})}^{2}+\frac{1}{2\Delta
t(1-\Delta t)}\mathcal{E}_{k}^{\rm best}+\frac{C_{1}\Delta t^{2}}{1-\Delta
t}+\frac{1}{\Delta t(1-\Delta t)}\mathcal{E}_{k}^{\rm stat}$
$\displaystyle\leq(1+6\Delta
t)\|u_{k+1}-\hat{u}_{k+1}\|_{{L}^{2}_{\omega_{k+1}}(\mathbb{R}^{d})}^{2}+\left(1+\frac{1}{2\Delta
t}\right)\mathcal{E}_{k}^{\rm best}+\frac{6}{5}C_{1}\Delta
t^{2}+\frac{6}{5\Delta t}\mathcal{E}_{k}^{\rm stat},$

provided that $\Delta t\leq 1/6$. Finally, using inequality
$|\max(a,b)-\max(a,c)|\leq|b-c|$, we derive

$\|u_{k}-\hat{u}_{k}\|_{{L}^{2}_{\omega_{k}}(\mathbb{R}^{d})}^{2}=\|\max(G_{k},c_{k})-\max(G_{k},\hat{c}_{k})\|_{{L}^{2}_{\omega_{k}}(\mathbb{R}^{d})}^{2}\leq\|c_{k}-\hat{c}_{k}\|_{{L}^{2}_{\omega_{k}}(\mathbb{R}^{d})}^{2}.$

This completes the proof.

###  5.2 Global error estimation

Finally, we prove a global error estimate.

######  Theorem 5.8.

For $\Delta t\leq 1/6$, there holds

$\max_{0\leq k\leq
N-1}\|u_{k}-\hat{u}_{k}\|_{{L}^{2}_{\omega_{k}}(\mathbb{R}^{d})}^{2}\leq
C\left(\Delta t+N\sum_{k=0}^{N-1}\left(\mathcal{E}_{k}^{\rm
best}+\mathcal{E}_{k}^{\rm stat}\right)\right),$

where the constant
$C=\max\left(\frac{6}{5T}e^{6T},\frac{6T}{5}C_{1}e^{6T}\right)$ only depends
on the finite time horizon $T$ and the constant $C_{1}>0$ defined in Lemma
5.2.

######  Proof 5.9.

For $\Delta t\leq 1/6$, Corollary 5.6 implies

$\|u_{k}-\hat{u}_{k}\|_{{L}^{2}_{\omega_{k}}(\mathbb{R}^{d})}^{2}\leq(1+6\Delta
t)\|u_{k+1}-\hat{u}_{k+1}\|_{{L}^{2}_{\omega_{k+1}}(\mathbb{R}^{d})}^{2}+\frac{6}{5\Delta
t}\left(\mathcal{E}_{k}^{\rm best}+\mathcal{E}_{k}^{\rm
stat}\right)+\frac{6}{5}C_{1}\Delta t^{2}.$

Using the discrete Gronwall’s inequality and $\Delta t=T/N$, we obtain

$\max_{0\leq k\leq
N-1}\|u_{k}-\hat{u}_{k}\|_{{L}^{2}_{\omega_{k}}(\mathbb{R}^{d})}^{2}\leq
e^{6T}\|u_{N}-\hat{u}_{N}\|_{{L}^{2}_{\omega_{N}}(\mathbb{R}^{d})}^{2}+e^{6T}\frac{6N}{5T}\sum_{k=0}^{N-1}\left(\mathcal{E}_{k}^{\rm
best}+\mathcal{E}_{k}^{\rm stat}\right)+e^{6T}\frac{6T}{5}C_{1}\Delta t.$

Since the terminal condition $u_{N}$ is known, the result follows by
$\|u_{N}-\hat{u}_{N}\|_{{L}_{\omega_{N}}^{2}(\mathbb{R}^{d})}=0$.

Theorem 5.8 implies that the numerical value function $\hat{u}_{k}$ computed
by Algorithm 1 at each time $t_{k}$ approximates the exact one ${u}_{k}$ well
if $\Delta t$, $\mathcal{E}_{k}^{\rm best}$, and $\mathcal{E}_{k}^{\rm stat}$
are small. Using Lemma 5.1, the best approximation error $\mathcal{E}_{k}^{\rm
best}$ decays geometrically with the increase of the order $p$ of Hermite
polynomials. By the law of large numbers, the statistical error
$\mathcal{E}_{k}^{\rm stat}$ between the discrete and the continuous least
squares will be small when using a large number of sample paths.

##  6 Numerical examples

In this section, we present several examples of high-dimensional Bermudan
option pricing to show the efficiency and accuracy of Algorithm 1. The codes
for the numerical experiments can be founded in the GitHub repository
https://github.com/jiefeiy/glsm-american. The accuracy of the computed prices
$\hat{v}_{0}(\mathbf{s}_{0})$ and deltas $\nabla\hat{v}_{0}(\mathbf{s}_{0})$
are measured by the relative errors defined by

$\frac{|\hat{v}_{0}(\mathbf{s}_{0})-v_{0}^{\dagger}|}{|v_{0}^{\dagger}|}\times
100\%\text{ and
}\frac{\|\nabla\hat{v}_{0}(\mathbf{s}_{0})-\Delta_{0}^{\dagger}\|}{\|\Delta_{0}^{\dagger}\|}\times
100\%,$

respectively, where $v_{0}^{\dagger}$ and $\Delta_{0}^{\dagger}$ are exact
price and delta at time $t=0$. Unless otherwise stated, the results are the
average of 10 independent runs. The parameter settings of the examples are
summarized in Table 1, which has been considered previously [25, 21, 6, 3, 9].

Table 1: The parameters used in Examples 1-5. Examples 3 and 4 share the parameters except the volatility $\sigma_{i}$. Example | Parameters
---|---
1\. Geometric basket put | $K=100,T=0.25,r=0.03,\delta_{i}=0,\sigma_{i}=0.2,\rho_{ij}=0.5,N=50$
2\. Geometric basket call | $K=100,T=2,r=0,\delta_{i}=0.02,\sigma_{i}=0.25,\rho_{ij}=0.75,N=50$
3\. Max-call with $d$ symmetric assets | $K=100,T=3,r=0.05,\delta_{i}=0.1,\sigma_{i}=0.2,\rho_{ij}=0,N=9$
4\. Max-call with $d$ asymmetric assets | $\sigma_{i}=\begin{cases}0.08+0.32\times(i-1)/(d-1),\text{ if }d\leq 5\\\ 0.1+i/(2d),\text{ if }d>5\end{cases}$
5\. Put option under Heston model |  $K=10,T=0.25,r=0.1,v_{0}=0.0625,\rho=0.1,$ $\kappa=5,\theta=0.16,\nu=0.9,N=50$

###  6.1 Example 1: Bermudan geometric basket put

Geometric basket options are benchmark tests for high-dimensional option
pricing problems, since they can be reduced to one-dimensional problems, and
thus highly accurate prices are available. Indeed, the price of the
$d$-dimensional problem equals that of the one-dimensional American option
with initial price $\hat{s}_{0}$, volatility $\hat{\sigma}$, and dividend
yield $\hat{\delta}$ given respectively by

$\hat{s}_{0}=\Big{(}\prod_{i=1}^{d}s_{0}^{i}\Big{)}^{1/d},\quad\hat{\sigma}=\frac{1}{d}\sqrt{\sum_{i,j}\sigma_{i}\sigma_{j}\rho_{ij}},\quad\hat{\delta}=\frac{1}{d}\sum_{i=1}^{d}\Big{(}\delta_{i}+\frac{\sigma_{i}^{2}}{2}\Big{)}-\frac{\hat{\sigma}^{2}}{2}.$

We consider the example of Bermudan geometric basket put from [13, 25]. The
exact prices are computed by solving the reduced one-dimensional problem via a
quadrature and interpolation-based method [25] for Bermudan options. We
present in Table 2 the computed option prices and their relative errors using
G-LSM and LSM, with the same ansatz space for the CVF. The results show that
G-LSM achieves higher accuracy than LSM in high dimensions: G-LSM has a
relative error $0.55\%$ for $d=15$, which is almost ten times smaller than
that by LSM. That is, by incorporating the gradient information, the accuracy
of LSM can be substantially improved.

Table 2: The price and relative errors of Bermudan geometric basket put computed by LSM and G-LSM using $M=100,000$ samples, with $s_{0}^{i}=100$ and $p=10$. $d$ | $N_{b}$ | LSM | error | G-LSM | error | $v_{0}^{\dagger}$
---|---|---|---|---|---|---
$1$ | 11 | 3.6715 | 0.16% | 3.6532 | 0.34% | 3.6658
$2$ | 29 | 3.1886 | 0.18% | 3.1791 | 0.12% | 3.1831
$3$ | 56 | 3.0113 | 0.28% | 2.9896 | 0.44% | 3.0029
$5$ | 141 | 2.8700 | 0.70% | 2.8381 | 0.41% | 2.8499
$10$ | 581 | 2.8030 | 2.71% | 2.7160 | 0.48% | 2.7290
$15$ | 1446 | 2.8258 | 5.15% | 2.6726 | 0.55% | 2.6874

Table 3 gives the computed option prices and their relative errors by G-LSM
with different maximum polynomial orders $p$ for the $20$-dimensional
geometric basket put. The relative error decays steadily as the order $p$ of
Hermite polynomials increases, which agrees with Theorem 5.8.

Table 3: The computed prices and their relative errors by G-LSM for $d=20$ with polynomial order $p$ and $M=100,000$ samples. $v_{0}^{\dagger}=2.6664$. $p$ | 2 | 5 | 10 | 12
---|---|---|---|---
$\hat{v}_{0}(\mathbf{s}_{0})$ | 2.6389 | 2.6414 | 2.6529 | 2.6634
error | 1.03% | 0.94% | 0.50% | 0.11%

###  6.2 Example 2: American geometric basket call

Now we consider the example of American geometric basket call option from [21,
6] to demonstrate that the proposed G-LSM can achieve the same level of
accuracy as the DNN-based method [6], and take $M=720,000$ samples as in [6].
The prices and deltas are given in Table 4 and Table 5, respectively, where
the results of [6] are the average of 9 independent runs. From Table 4, both
methods have similar accuracy for the price. From Table 5, the relative error
of delta using G-LSM and DNN varies slightly with the dimension $d$. This is
probably because the DNN-based method computes the delta via a sample average,
while G-LSM uses the derivatives of the value function directly. The relative
error of the delta computed by G-LSM increases slightly for larger dimensions,
possibly due to the small magnitude of the exact delta values.

Table 4: Prices of American geometric basket call at $t=0$ using $M=720,000$ samples. |  |  |  | G-LSM |  |  DNN [6]
---|---|---|---|---|---|---
$d$ | $s_{0}^{i}$ | $p$ | exact | price | error |  | price | error
7 | 100 | 10 | 10.2591 | 10.2475 | 0.11% |  | 10.2468 | 0.12%
13 | 100 | 10 | 10.0984 | 10.0781 | 0.20% |  | 10.0822 | 0.16%
20 | 100 | 10 | 10.0326 | 10.0141 | 0.18% |  | 10.0116 | 0.21%
100 | 100 | 6 | 9.9345 | 9.8980 | 0.37% |  | 9.9163 | 0.18%
Table 5: Deltas of American geometric basket call at $t=0$ using $M=720,000$ samples. $\mathbf{1}=(1,1,\dots,1)^{\top}$. |  |  |  | G-LSM |  |  DNN [6]
---|---|---|---|---|---|---
$d$ | $s_{0}^{i}$ | $p$ | exact | delta | error |  | delta | error
7 | 100 | 10 | $0.0722*\mathbf{1}$ | $(0.0724,0.0725,\dots,0.0724)^{\top}$ | 0.32% |  | $0.0717*\mathbf{1}$ | 0.69%
13 | 100 | 10 | $0.0387*\mathbf{1}$ | $(0.0389,0.0387,\dots,0.0388)^{\top}$ | 0.39% |  | $0.0384*\mathbf{1}$ | 0.78%
20 | 100 | 10 | $0.0251*\mathbf{1}$ | $(0.0253,0.0253,\dots,0.0254)^{\top}$ | 0.59% |  | $0.0249*\mathbf{1}$ | 0.80%
100 | 100 | 6 | $0.00502*\mathbf{1}$ | $(0.00501,0.00508,\dots,0.00498)^{\top}$ | 1.45% |  | $0.00498*\mathbf{1}$ | 0.80%

With a priori knowledge of the regularity of the CVF, we can approximate
functions with polynomials in high dimensions. We briefly compare the
complexity of the two approaches. Compared with the DNN approximator, G-LSM
can involve fewer unknown parameters, and involves a simpler numerical task
(solving least-squares problems versus minimizing nonconvex losses). Indeed,
at each fixed time step, [6] suggests training the neural network with $L=7$
hidden layers and width $d+5$ in each layer, leading to more than $L(d+5)^{2}$
parameters. In contrast, the number of undetermined parameters in G-LSM is the
number $N_{b}$ of basis functions, which has a cardinality $\mathcal{O}(p(\ln
p)^{d-1})$. For example, for $d=100$, the neural network approach involves
more than $77175$ parameters, whereas G-LSM with $p=6$ involves only
$N_{b}=15451$ basis functions.

Next, Fig. 2 shows the classification results of continued and exercised data
using G-LSM and LSM with the number of simulated paths $M=100,000$ in $d=7$ or
20. Compared with the exact exercise boundary, G-LSM achieves better accuracy
in determining the exercise boundary than LSM despite of using the same ansatz
and number of paths. Thus, even with the right ansatz space, LSM might fail
the task of finding exercise boundary in high dimensions using only a limited
number of samples. Compared with [6, Figure 5] and [18, Figure 6], Fig. 2
domonstrates that G-LSM can detect accurate exercise boundary with fewer
number of paths than the DNN-based method.

(a) G-LSM  (b) LSM  (c) G-LSM  (d) LSM

Figure 2: Classification of the simulated continued and exercised data using
$M=100,000$ samples, with $p=10$. Black star dots represent the exact exercise
boundary.

###  6.3 Example 3: Bermudan max-call with symmetric assets

To benchmark G-LSM on high-dimensional problems without exact solutions and to
validate the complexity analysis in Section 4, we test Bermudan max-call
option and report the computing time. The computing time is calculated as
follows. For a fixed time step, $T_{\rm bas}$ is the time for generating basis
matrix $\Phi$, $T_{\rm mat}$ is the time for assembling matrix $A$, $T_{\rm
lin}$ is the time for solving linear system, and $T_{\rm up}$ is the time for
updating values. The overall computing time is $T_{\rm tot}\approx(N-1)(T_{\rm
bas}+T_{\rm mat}+T_{\rm lin}+T_{\rm up})$.

Table 6 presents the prices and computing time (in seconds) for Bermudan max-
call options with $d$ symmetric assets. The reference 95% confidence interval
(CI) is taken from [3]. The reference CI is computed with more than $3000$
training steps and a batch of $8192$ paths in each step, which in total
utilizes more than $10^{7}$ paths. The last five columns of the table report
the computing time for the step 6, 7, 8, 9 in Algorithm 1, and the total time,
respectively. All the computation for this example was performed on an Intel
Core i9-10900 CPU 2.8 GHz desktop with 64GB DDR4 memory using MATLAB R2023b.
It is observed that the prices computed by G-LSM fall into or stay very close
to the reference 95% CI, confirming the high accuracy of G-LSM. Furthermore,
the time for generating basis matrix, $T_{\rm bas}$, dominates the overall
computing time. Hence, the cost mainly arises from evaluating Hermite
polynomials on sampling paths, which is also required by LSM. In comparison
with LSM , $T_{\rm mat}$ is the extra cost to incorporate the gradient
information and takes only a small fraction of the total time. Therefore,
G-LSM has nearly identical cost with LSM.

We present in Table 7 the computing time for each basis function to verify
that the complexity is almost linear in $N_{b}$ as analyzed in Section 4. The
ratio does not vary much with the dimension $d$. Since
$N_{b}=\mathcal{O}(p(\ln p)^{d-1})$ with polynomials up to order $p$, the
total computing cost of G-LSM exhibits a polynomial growth, which overcomes
the so-called curse of dimensionality. Figure 3 shows the growth of $N_{b}$
with respect to the dimension for $1\leq d\leq 200$ and $p=6$, indicating that
$N_{b}$ exhibits a nearly quadratic growth in dimensions for $d\leq 200$.

Table 6: The results for Bermudan max-call options with $d$ symmetric assets using $M=100,000$ samples. $s_{0}^{i}=100$ for $i=1,\dots,d$. $d$ | $p$ | $N_{b}$ | reference 95% CI | G-LSM | $T_{\rm bas}$ | $T_{\rm mat}$ | $T_{\rm lin}$ | $T_{\rm up}$ | $T_{\rm tot}$
---|---|---|---|---|---|---|---|---|---
2 | 10 | 29 | $[13.880,13.910]$ | 13.8970 | 0.1678 | 0.0128 | 0.0024 | 0.0039 | 1.5725
3 | 10 | 56 | $[18.673,18.699]$ | 18.6715 | 0.2890 | 0.0265 | 0.0048 | 0.0060 | 2.6347
5 | 10 | 141 | $[26.138,26.174]$ | 26.0553 | 0.5718 | 0.0758 | 0.0149 | 0.0132 | 5.4662
10 | 10 | 581 | $[38.300,38.367]$ | 38.1738 | 1.9932 | 0.3164 | 0.1397 | 0.0497 | 19.9308
20 | 10 | 2861 | $[51.549,51.803]$ | 51.6508 | 11.3278 | 1.7153 | 3.2000 | 0.2434 | 131.9362
30 | 5 | 1456 | $[59.476,59.872]$ | 59.5475 | 7.8154 | 0.7750 | 0.7044 | 0.1247 | 75.3066
50 | 5 | 3926 | $[69.560,69.945]$ | 69.7216 | 29.3088 | 2.0516 | 5.6871 | 0.3281 | 299.4528
100 | 4 | 5351 | $[83.357,83.862]$ | 83.6777 | 71.7678 | 2.7185 | 10.8018 | 0.4306 | 684.9927
Table 7: Ratio of the computing time and number $N_{b}$ of basis functions.  $d$ | $T_{\rm bas}/N_{b}$ | $T_{\rm mat}/N_{b}$ | $T_{\rm lin}/N_{b}$ | $T_{\rm up}/N_{b}$ | $T_{\rm tot}/N_{b}$
---|---|---|---|---|---
2 | 0.0058 | 0.0004 | 0.0001 | 0.0001 | 0.0542
3 | 0.0052 | 0.0005 | 0.0001 | 0.0001 | 0.0470
5 | 0.0041 | 0.0005 | 0.0001 | 0.0001 | 0.0388
10 | 0.0034 | 0.0005 | 0.0002 | 0.0001 | 0.0343
20 | 0.0040 | 0.0006 | 0.0011 | 0.0001 | 0.0461
30 | 0.0054 | 0.0005 | 0.0005 | 0.0001 | 0.0517
50 | 0.0075 | 0.0005 | 0.0014 | 0.0001 | 0.0763
100 | 0.0134 | 0.0005 | 0.0020 | 0.0001 | 0.1280
Figure 3: The number $N_{b}$ of basis functions in dimension $1\leq d\leq 200$
with hyperbolic cross index set versus linear and quadratic scale in $d$, with
$p=6$.

Table 6 reports the time $T_{\rm lin}$ of solving linear system using the
built-in MATLAB function cgs. The main memory requirement is storing the
matrix $A\in\mathbb{R}^{M\times N_{b}}$, which takes $M\times N_{b}\times 8$
bytes in double-precision floating-point format. For $d=100$, the memory for
storing $A$ is about 4GB. This limits the application of G-LSM with direct
solvers in higher dimensions. To remedy the issue, one can solve the linear
system on a large RAM server, or use single-precision floating-point, or with
stochastic gradient descent.

Figure 4 shows the classification of continued and exercised sample points
computed by G-LSM and LSM in the example of two-dimensional max-call. G-LSM
yields a smoother exercise boundary than LSM. Compared with the exercise
boundary computed in literature [19, Figure 3], G-LSM exhibits higher
accuracy, albeit that the same ansatz space for the CVF is employed.

(a) G-LSM, at $t_{1}$ (b) G-LSM, at $t_{4}$ (c) G-LSM, at $t_{7}$ (d) LSM, at
$t_{1}$ (e) LSM, at $t_{4}$ (f) LSM, at $t_{7}$

Figure 4: Classification of simulated continued and exercised data for 2-$d$
Bermudan max-call option at different times $t_{k}$, $k=1,4,7$. Here, we take
$s_{0}^{i}=100$ for $i=1,2$, $M=100,000$ and $p=20$.

###  6.4 Example 4: Bermudan max-call with asymmetric assets

Example 3 assumes that all underlying assets follow the same dynamic. To
further show the robustness of G-LSM, we consider the Bermudan max-call option
but each asset has different volatility. The reference 95% confidence interval
(CI) is taken from [3]. The pricing results for different initial price
$s_{0}^{i}$ and dimension $d$ are listed in Table 8. The standard error (s.e.)
is calculated by
$\sqrt{\frac{1}{10(10-1)}\sum_{i=1}^{10}(v_{0}^{(i)}-\bar{v})^{2}}$ with
$\bar{v}$ being the average of $10$ independent runs. Similar to Example 3,
the prices computed by G-LSM always fall into or stay very close to the
reference 95% CI.

Table 8: Results for Bermudan max-call options with $d$ asymmetric assets, with $M=100,000$.  $d$ | $p$ | $s_{0}^{i}$ | reference 95% CI | G-LSM | s.e.
---|---|---|---|---|---
2 | 10 | 90 | $[14.299,14.367]$ | 14.3472 | 0.0230
2 | 10 | 100 | $[19.772,19.829]$ | 19.8019 | 0.0371
2 | 10 | 110 | $[27.138,27.163]$ | 27.1041 | 0.0213
3 | 10 | 90 | $[19.065,19.104]$ | 19.0266 | 0.0241
3 | 10 | 100 | $[26.648,26.701]$ | 26.6931 | 0.0411
3 | 10 | 110 | $[35.806,35.835]$ | 35.8363 | 0.0472
5 | 10 | 90 | $[27.630,27.680]$ | 27.6032 | 0.0259
5 | 10 | 100 | $[37.940,38.014]$ | 37.9309 | 0.0405
5 | 10 | 110 | $[49.445,49.533]$ | 49.3711 | 0.0473
10 | 10 | 90 | $[85.857,86.087]$ | 85.8221 | 0.0376
10 | 10 | 100 | $[104.603,104.864]$ | 104.7052 | 0.1029
10 | 10 | 110 | $[123.570,123.904]$ | 123.4777 | 0.0686
20 | 10 | 90 | $[125.819,126.383]$ | 126.4276 | 0.0980
20 | 10 | 100 | $[149.480,150.053]$ | 150.4028 | 0.1194
20 | 10 | 110 | $[173.144,173.937]$ | 173.8584 | 0.1349
30 | 5 | 90 | $[154.378,155.039]$ | 154.6913 | 0.1128
30 | 5 | 100 | $[181.155,182.033]$ | 181.6733 | 0.1385
30 | 5 | 110 | $[208.091,209.086]$ | 208.1267 | 0.1160
50 | 5 | 90 | $[195.793,196.963]$ | 196.6921 | 0.0890
50 | 5 | 100 | $[227.247,228.605]$ | 227.7831 | 0.1385
50 | 5 | 110 | $[258.661,260.092]$ | 259.7261 | 0.1413
100 | 4 | 90 | $[263.043,264.425]$ | 263.0543 | 0.1515
100 | 4 | 100 | $[301.924,303.843]$ | 302.0623 | 0.2130
100 | 4 | 110 | $[340.580,342.781]$ | 340.8581 | 0.2233

###  6.5 Example 5: Bermudan put under Heston model

The previous experiments and theoretical analysis are concerned with the most
frequently used multi-asset Black-Scholes model, cf. Section 2.2. We now
generalize G-LSM to price Bermudan option under the Heston model. The Heston
model defines the dynamic of the log-price process,
$X_{t}^{1}:=\ln(S_{t}/S_{0})$, and the volatility process, $v_{t}$, by two-
dimensional SDEs:

$\displaystyle\mathrm{d}X_{t}^{1}$
$\displaystyle=(r-\tfrac{1}{2}v_{t})\,\mathrm{d}t+\sqrt{v_{t}}\left(\rho\,\mathrm{d}W_{t}^{1}+\sqrt{1-\rho^{2}}\,\mathrm{d}W_{t}^{2}\right),$
$\displaystyle\mathrm{d}v_{t}$ $\displaystyle=\kappa(\theta-
v_{t})\,\mathrm{d}t+\nu\sqrt{v_{t}}\,\mathrm{d}W_{t}^{1},$

where $W_{t}^{1}$ and $W_{t}^{2}$ are two independent Wiener processes, and
the model parameters $r$, $\kappa$, $\theta$, $\nu$ and $\rho$ represent the
interest rate, the speed of mean reversion, the mean level of variance, the
variance of volatility process, and the correlation coefficient, respectively.
Since the transition density of log-variance has better regularity [9], we
take the log-variance process $X_{t}^{2}:=\ln(v_{t})$ as the regression
variable. Let $\mathbf{X}_{t}=[X_{t}^{1},X_{t}^{2}]^{\top}$ be a two-
dimensional process. The discounted continuation value at time $t_{k}$ is
given by

$c_{k}(\mathbf{X}_{t_{k}})=\mathbb{E}[u_{k+1}(\mathbf{X}_{t_{k+1}})|\mathbf{X}_{t_{k}}],$

with $u_{k+1}$ being the discounted value function. In view of the martingale
representation theorem and backward Euler approximation, we obtain

$u_{k+1}(\mathbf{X}_{t_{k+1}})\approx
c_{k}(\mathbf{X}_{t_{k}})+\big{(}\sigma(\mathbf{X}_{t_{k}})^{\top}\nabla
c_{k}(\mathbf{X}_{t_{k}})\big{)}\cdot\Delta\mathbf{W}_{k},$

where $\sigma(\cdot)$ is the covariance matrix of $\mathbf{X}_{t}$ given by

$\sigma(\mathbf{X}_{t})=\begin{bmatrix}\rho\exp(X_{t}^{2}/2)&\sqrt{1-\rho^{2}}\exp(X_{t}^{2}/2)\\\
\nu\exp(-X_{t}^{2}/2)&0\end{bmatrix}.$

The reference prices are computed by the COS method [9], with $2^{7}$ cosine
basis functions to approximate the transition density and $2^{7}$ points for
the quadrature rule in log-variance dimension. The results of G-LSM are
calculated using 70 basis functions (polynomials up to order $20$) with
Hermite polynomials in log-price dimension and Chebyshev polynomials in log-
variance dimension. The results in Table 9 show that the computed prices by
G-LSM coincide with the reference prices up to the two places after the
decimal separator for most cases. In Fig. 5, we plot the classification
results of continued and exercised points using two approaches. Given that the
G-LSM method is simulation-based, it has great potential for stochastic
volatility models in higher dimensions.

Table 9: Prices of Bermudan put option under the Heston model. $M=100,000$. $p=20$. $s_{0}$ | 8 | 9 | 10 | 11 | 12
---|---|---|---|---|---
COS | 1.9958 | 1.1061 | 0.5186 | 0.2131 | 0.0818
G-LSM | 1.9949 | 1.0972 | 0.5146 | 0.2118 | 0.0807

(a) COS  (b) G-LSM

Figure 5: Classification of continued and exercised grid/simulated points
using COS and G-LSM under the Heston model at time $t_{25}$ with initial price
$s_{0}=8$.

##  7 Conclusions and outlook

In this work, we have proposed a novel gradient-enhanced least squares Monte
Carlo (G-LSM) method that employs sparse Hermite polynomials as the ansatz
space to price and hedge American options. The method enjoys low complexity
for the gradient evaluation, ease of implementation and high accuracy for
high-dimensional problems. We analyzed rigorously the convergence of G-LSM
based on the BSDE technique, stochastic and Malliavin calculus. Extensive
benchmark tests clearly show that it outperforms least squares Monte Carlo
(LSM) in high dimensions with almost the same cost and it can also achieve
competitive accuracy relative to the deep neural networks-based methods.

There are several avenues for further research. The superiority of G-LSM over
LSM in high dimensions indicates that matching option values at $t_{k+1}$
might be a better choice than at $t_{k}$ for approximating the continuation
value function. There are other variants of LSM with different ansatz spaces,
and it is natural to ask whether incorporating the gradient information will
also result in improved performance for these variants. Moreover, to solve
higher dimensional problems, e.g., $d=1000$, the hierarchical tensor train
technique can be applied, which has been combined with LSM in [2], and it is
natural to combine the technique with G-LSM. Numerical results also show the
potential of G-LSM for stochastic volatility models. It would be interesting
to investigate G-LSM for more challenging financial models, e.g., rough
volatility models.

## Acknowledgments

The authors acknowledge the support of research computing facilities offered
by Information Technology Services, the University of Hong Kong.

## References

  * [1] B. Adcock, S. Brugiapaglia, and C. G. Webster, Sparse Polynomial Approximation of High-Dimensional Functions, SIAM, Philadelphia, PA, 2022.
  * [2] C. Bayer, M. Eigel, L. Sallandt, and P. Trunschke, Pricing high-dimensional Bermudan options with hierarchical tensor formats, SIAM Journal on Financial Mathematics, 14 (2023), pp. 383–406.
  * [3] S. Becker, P. Cheridito, and A. Jentzen, Deep optimal stopping, The Journal of Machine Learning Research, 20 (2019), pp. 2712–2736.
  * [4] S. Becker, P. Cheridito, and A. Jentzen, Pricing and hedging American-style options with deep learning, Journal of Risk and Financial Management, 13 (2020), p. 158.
  * [5] B. Bouchard and X. Warin, Monte-Carlo valuation of American options: facts and new algorithms to improve existing methods, in Numerical Methods in Finance: Bordeaux, June 2010, Springer, Berlin, 2012, pp. 215–255.
  * [6] Y. Chen and J. W. Wan, Deep neural network framework based on backward stochastic differential equations for pricing and hedging American options in high dimensions, Quantitative Finance, 21 (2021), pp. 45–67.
  * [7] W. E, J. Han, and A. Jentzen, Deep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations, Communications in Mathematics and Statistics, 5 (2017), pp. 349–380.
  * [8] N. El Karoui, C. Kapoudjian, E. Pardoux, S. Peng, and M.-C. Quenez, Reflected solutions of backward SDE’s, and related obstacle problems for PDE’s, The Annals of Probability, 25 (1997), pp. 702–737.
  * [9] F. Fang and C. W. Oosterlee, A Fourier-based valuation method for Bermudan and barrier options under Heston’s model, SIAM Journal on Financial Mathematics, 2 (2011), pp. 439–463.
  * [10] C. Gao, S. Gao, R. Hu, and Z. Zhu, Convergence of the backward deep bsde method with applications to optimal stopping problems, SIAM Journal on Financial Mathematics, 14 (2023), pp. 1290–1303.
  * [11] E. Gobet and C. Labart, Error expansion for the discretization of backward stochastic differential equations, Stochastic Processes and their Applications, 117 (2007), pp. 803–829.
  * [12] C. Huré, H. Pham, and X. Warin, Deep backward schemes for high-dimensional nonlinear PDEs, Mathematics of Computation, 89 (2020), pp. 1547–1579.
  * [13] P. Kovalov, V. Linetsky, and M. Marcozzi, Pricing multi-asset American options: A finite element method-of-lines with smooth penalty, Journal of Scientific Computing, 33 (2007), pp. 209–237.
  * [14] B. Lapeyre and J. Lelong, Neural network regression for Bermudan option pricing, Monte Carlo Methods and Applications, 27 (2021), pp. 227–247.
  * [15] F. Longstaff and E. Schwartz, Valuing American options by simulation: a simple least-squares approach, The Review of Financial Studies, 14 (2001), pp. 113–147.
  * [16] M. Ludkovski, Kriging metamodels and experimental design for Bermudan option pricing, Journal of Computational Finance, 22 (2018), pp. 37–77.
  * [17] X. Luo, Error analysis of the Wiener–Askey polynomial chaos with hyperbolic cross approximation and its application to differential equations with random input, Journal of Computational and Applied Mathematics, 335 (2018), pp. 242–269.
  * [18] A. S. Na and J. W. L. Wan, Efficient pricing and hedging of high-dimensional American options using deep recurrent networks, Quantitative Finance, 23 (2023), pp. 631–651.
  * [19] A. M. Reppen, H. M. Soner, and V. Tissot-Daguette, Deep stochastic optimization in finance, Digital Finance, 5 (2023), pp. 91–111.
  * [20] R. Seydel and R. Seydel, Tools for computational finance, vol. 3, Springer, 2006.
  * [21] J. Sirignano and K. Spiliopoulos, DGM: A deep learning algorithm for solving partial differential equations, Journal of Computational Physics, 375 (2018), pp. 1339–1364.
  * [22] J. Tsitsiklis and B. Van Roy, Regression methods for pricing complex American-style options, IEEE Transactions on Neural Networks, 12 (2001), pp. 694–703.
  * [23] H. Wang, H. Chen, A. Sudjianto, R. Liu, and Q. Shen, Deep learning-based BSDE solver for LIBOR market model with application to Bermudan swaption pricing and hedging, arXiv preprint arXiv:1807.06622, (2018).
  * [24] Y. Wang and R. Caflisch, Pricing and hedging American-style options: a simple simulation-based approach, The Journal of Computational Finance, 13 (2009), pp. 95–125.
  * [25] J. Yang and G. Li, On sparse grid interpolation for American option pricing with multiple underlying assets, arXiv preprint arXiv:2309.08287, (2023).
  * [26] J. Zhang, Backward stochastic differential equations, Springer, 2017\.