#  A new approach for American option pricing: The Dynamic Chebyshev method

Kathrin Glau$\vphantom{l}{}^{1,2}$, Mirco Mahlstedt$\vphantom{l}{}^{2,*}$,
Christian Pötz$\vphantom{l}{}^{1,2,}$111The authors thank the KPMG Center of
Excellence in Risk Management for their support.

$\vphantom{l}{}^{\text{1}}$Queen Mary University of London, UK
$\vphantom{l}{}^{\text{2}}$Technical University of Munich, Germany

###### Abstract

We introduce a new method to price American options based on Chebyshev
interpolation. In each step of a dynamic programming time-stepping we
approximate the value function with Chebyshev polynomials. The key advantage
of this approach is that it allows to shift the model-dependent computations
into an offline phase prior to the time-stepping. In the offline part a family
of generalised (conditional) moments is computed by an appropriate numerical
technique such as a Monte Carlo, PDE or Fourier transform based method. Thanks
to this methodological flexibility the approach applies to a large variety of
models. Online, the backward induction is solved on a discrete Chebyshev grid,
and no (conditional) expectations need to be computed. For each time step the
method delivers a closed form approximation of the price function along with
the options’ delta and gamma. Moreover, the same family of (conditional)
moments yield multiple outputs including the option prices for different
strikes, maturities and different payoff profiles. We provide a theoretical
error analysis and find conditions that imply explicit error bounds for a
variety of stock price models. Numerical experiments confirm the fast
convergence of prices and sensitivities. An empirical investigation of
accuracy and runtime also shows an efficiency gain compared with the least-
square Monte-Carlo method introduced by Longstaff and Schwartz (2001).

Keywords American Option Pricing, Complexity Reduction, Dynamic Programming,
Polynomial Interpolation

2010 MSC 91G60, 41A10

##  1 Introduction

A challenging task for financial institutions is the computation of prices and
sensitivities for large portfolios of derivatives such as equity options.
Typically, equity options have an early exercise feature and can either be
exercised at any time until maturity (American type) or at a set of pre-
defined exercise dates (Bermudan type). In lack of explicit solutions,
different numerical methods haven been developed to tackle this problem. One
of the first algorithms to compute American put option prices in the Black-
Scholes model has been proposed by Brennan and Schwartz (1977). In this
approach, the related partial differential inequality is solved by a finite
difference scheme. A rich literature further developing the PDE approach has
accrued since, including methods for jump models (Levendorskiĭ (2004), Hilber
et al. (2013)), extensions to two dimensions (Haentjens and in’t Hout (2015))
and combinations with complexity reduction techniques (Haasdonk et al.
(2013)). Besides PDE based methods a variety of other approaches has been
introduced, many of which trace back to the solution of the optimal stopping
problem by the dynamic programming principle, see e.g. Peskir and Shiryaev
(2006). For Fourier based solution schemes we refer to Lord et al. (2008),
Fang and Oosterlee (2009). Simulation based approaches are of fundamental
importance, the most prominent representative of this group is the Least-
squares Monte-Carlo (LSM) approach of Longstaff and Schwartz (2001), we refer
to Glasserman (2003) for an overview of different Monte-Carlo methods. Fourier
and PDE methods typically are highly efficient, compared to simulation,
however, they are less flexible towards changes in the model and particularly
in the dimensionality.

In order to reconcile the advantages of the PDE and Fourier approach with the
flexibility of Monte Carlo simulation, we propose a new approach. We consider
a dynamic programming time-stepping. Let ${X_{t}}$ be the underlying Markov
process and the value function $V_{t}$ is given by,

$\displaystyle V_{T}(x)$ $\displaystyle=g(x)$ $\displaystyle V_{t}(x)$
$\displaystyle=f\left(g(t,x),\mathbb{E}[V_{t+1}(X_{t+1})|X_{t}=x]\right).$

with time steps $t<t+1<\ldots<T$ and payoff function $g$. The computational
challenge is to compute $\mathbb{E}[V_{t+1}(X_{t+1})|X_{t}=x]$ for for all
time steps $t$ and all states $x$, where $V_{t+1}$ depends on all previous
time steps.

In order to tackle this problem, we approximate the value function in each
time step by Chebyshev polynomial interpolation. We thus express the value
function as a finite sum of Chebyshev polynomials

(1.1) $\displaystyle\mathbb{E}[V_{t+1}(X_{t+1})|X_{t}=x]\,\approx\sum
c_{j}^{t+1}\mathbb{E}[T_{j}(X_{t+1})|X_{t}=x]$

The choice of Chebyshev polynomials is motivated by the promising properties
of Chebyshev interpolation such as

  * •

The vector of coefficients $(c_{j}^{t+1})_{j=0,\ldots,N}$ is explicitly given
as a linear combination of the values $V_{t}(x_{k})$ at the Chebyshev grid
points $x_{k}$. For this, equation (1.1) needs to be solved at the Chebyshev
grid points $x=x_{k}$ only.

  * •

Exponential convergence of the interpolation for analytic functions and
polynomial convergence of differential functions depending on the order.

  * •

The interpolation can be implemented in a numerically stable way.

The computation of the continuation value at a single time step coincides with
the pricing of a European option. Its interpolation with Chebyshev polynomials
is proposed in Gaß et al. (2018), where the method shows to be highly
promising and exponential convergence is established for a large set of models
and option types. Moreover, the approximation of the value function with
Chebyshev polynomials has already proven to be beneficial for optimal control
problems in economics, see Judd (1998) and Cai and Judd (2013).

The key advantage of our approach for American option pricing is that it
collects all model-dependent computations in the generalized conditional
moments $\Gamma_{j,k}=\mathbb{E}[T_{j}(X_{t+1})|X_{t}=x_{k}]$. If there is no
closed-form solution their calculation can be shifted into an offline phase
prior to the time-stepping. Depending on the underlying model a suitable
numerical technique such as Monte Carlo, PDE and Fourier transform methods can
be chosen, which reveals the high flexibility of the approach. Once the
generalized conditional moments $\Gamma_{j,k}$ are computed, the backward
induction is solved on a discrete Chebyshev grid. Which avoids any
computations of conditional expectations during the time-stepping. For each
time step the method delivers a closed form approximation of the price
function $x\mapsto\sum c_{j}^{t}T_{j}(x)$ along with the options’ delta and
gamma. Since the family of generalized conditional moments $\Gamma_{j,k}$ are
independent of the value function, they can be used to generate multiple
outputs including the option prices for different strikes, maturities and
different payoff profiles. The structure of the method is also beneficial for
the calculation of expected future exposure which is the computational
bottleneck in the computation of CVA, as investigated in Glau et al. (2018).

The offline-online decomposition separates model and payoff yielding a modular
design. We exploit this structure for a thorough error analysis and find
conditions that imply explicit error bounds. They reflect the modularity by
decomposing into a part stemming from the Chebyshev interpolation, from the
time-stepping and from the offline computation. Under smoothness conditions
the asymptotic convergence behaviour is deduced.

We perform numerical experiments using the Black-Scholes model, Merton’s jump
diffusion model and the Constant Elasticity of Variance (CEV) model as a
represenative of a local volatility model. For the computation of the
generalized conditional moments we thus use different techniques, namely
numerical integration based on Fourier transforms and Monte Carlo simulation.
Numerical experiments confirm the fast convergence of option prices along with
its delta and gamma. A comprehensive comparison with the LSM reveals the
potential efficiency gain of the new approach, particularly when several
options on the same underlying are priced.

The rest of the article is organized as follows. We introduce the problem
setting and the new method in Section 2 and provide the error analysis in
Section 3. Section 4 discusses general traits of the implementation and
Section 5 presents the numerical experiments. Section 6 concludes the article,
followed by an appendix with the proof of the main result.

##  2 The Chebyshev method for Dynamic programming problems

First, we present the Bellman-Wald equation as a specific form of dynamic
programming. Second, we provide the necessary notation for the Chebyshev
interpolation. Then we are in a position to introduce the new approach and its
application to American option pricing.

###  2.1 Optimal stopping and Dynamic Programming

Let $X=(X_{t})_{0\leq t\leq T}$ be a Markov process with state space
$\mathbb{R}^{d}$ defined on the filtered probability space
$(\Omega,\mathcal{F},(\mathcal{F}_{t})_{t\geq 0},\mathbb{P})$. Let
$g:[0,T]\times\mathbb{R}^{d}\longrightarrow\mathbb{R}$ be a continuous
function with $\mathbb{E}\left[\sup_{0\leq t\leq
T}\left|g(t,X_{t})\right|\right]<\infty.$ Then

$\displaystyle V(t,x):=\sup_{t\leq\tau\leq
T}\mathbb{E}\left[g(\tau,X_{\tau})|X_{t}=x\right]\qquad\text{for
all}\quad(t,x)\in[0,T]\times\mathbb{R}^{d}$

over all stopping times $\tau$, see $(2.2.2^{\prime})$ in Peskir and Shiryaev
(2006). In discrete time the optimal stopping problems can be solved with
dynamic programming.

Namely, with time stepping $t=t_{0}<\ldots<t_{n}=T$ the solution of the
optimal stopping problem can be calculated via backward induction

$\displaystyle V_{T}(x)$ $\displaystyle=g(T,x)$ $\displaystyle V_{t_{u}}(x)$
$\displaystyle=\max\left(g(t_{u},x),\mathbb{E}[V_{t_{u+1}}(X_{t_{u+1}})|X_{t_{u}}=x]\right).$

Note that $n$ refers to the number of time steps between $t$ and $T$. For
notational convenience, we indicate the value function at each time step with
subscript $t_{u}$ to directly refer to the time step $t_{u}$. For a detailed
overview of optimal control problems in discrete time we refer to Peskir and
Shiryaev (2006).

###  2.2 Chebyshev polynomial interpolation

The univariate Chebyshev polynomial interpolation as described in detail in
Trefethen (2013) has a tensor based extension to the multivariate case, see
e.g. Sauter and Schwab (2010). Usually, the Chebyshev interpolation is defined
for a function on a $[-1,1]^{D}$ domain. For an arbitrary hyperrectangular
$\mathcal{X}=[\underline{x}_{1},\overline{x}_{1}]\times\ldots\times[\underline{x}_{D},\overline{x}_{D}]$,
we introduce a linear transformation
$\tau_{\mathcal{X}}:[-1,1]^{D}\rightarrow\mathcal{X}$ componentwise defined by

(2.1)
$\displaystyle\tau_{\mathcal{X}}(z_{i})=\overline{x}_{i}+0.5(\underline{x}_{i}-\overline{x}_{i})(1-z_{i}).$

Let $\overline{N}:=(N_{1},\ldots,N_{D})$ with $N_{i}\in\mathds{N}_{0}$ for
$i=1,\ldots,D$. We define the index set

$\displaystyle\mathcal{J}:=\\{j\in\mathbb{N}^{D}:1\leq j_{i}\leq N_{D}\
\text{for}\ i=1,\ldots,d\\}.$

The Chebyshev polynomials are defined for $z\in[-1,1]^{D}$ and
$j\in\mathcal{J}$ by

$\displaystyle T_{j}(z)=\prod_{i=1}^{D}T_{j_{i}}(z_{i}),\quad
T_{j_{i}}(z_{i})=\cos(j_{i}\cdot\text{acos}(z_{i})),$

and the $j$-th Chebyshev polynomial on $\mathcal{X}$ as
$p_{j}(x)=T_{j}(\tau^{-1}_{\mathcal{X}}(x))1_{\mathcal{X}}(x)$. The Chebyshev
points are given by

$\displaystyle z^{k}=(z_{k_{1}},\dots,z_{k_{D}}),\
z_{k_{i}}=\cos\left(\pi\frac{k_{i}}{N_{i}}\right)\text{ for
}k_{i}=0,\ldots,N_{i}\text{ and }i=1,\ldots,D.$

and the transformed Chebyshev points by $x^{k}=\tau_{\mathcal{X}}(z^{k}).$ The
Chebyshev interpolation of a function $f:\mathcal{X}\rightarrow\mathbb{R}$
with $\prod_{i=1}^{D}(N_{i}+1)$ summands can be written as a sum of Chebyshev
polynomials

(2.2) $\displaystyle
I_{\overline{N}}(f)(x)=\sum_{j\in\mathcal{J}}c_{j}T_{j}(\tau_{\mathcal{X}}^{-1}(x))=\sum_{j\in\mathcal{J}}c_{j}p_{j}(x)\qquad\text{for}\quad
x\in\mathcal{X}$

with coefficients $c_{j}$ for $j\in\mathcal{J}$

(2.3) $\displaystyle c_{j}$
$\displaystyle=\Big{(}\prod_{i=1}^{D}\frac{2^{\mathds{1}_{\\{0<j_{i}<N_{i}\\}}}}{N_{i}}\Big{)}\sum_{k\in\mathcal{J}}{}^{{}^{\prime\prime}}f(x^{k})T_{j}(z^{k})$

where $\sum{}^{{}^{\prime\prime}}$ indicates the summand is multiplied by
$1/2$ if $k_{i}=0$ or $k_{i}=N_{i}$.

###  2.3 The Dynamic Chebyshev method

In this section, we present the new approach to solve a dynamic programming
problem via backward induction using Chebyshev polynomial interpolation.

######  Definition 2.1.

We consider a Dynamic Programming Problem (DPP) with value function

(2.4) $\displaystyle V_{T}(x)$ $\displaystyle=g(T,x)$ (2.5) $\displaystyle
V_{t_{u}}(x)$
$\displaystyle=f\left(g(t_{u},x),\mathbb{E}[V_{t_{u+1}}(X_{t_{u+1}})|X_{t_{u}}=x]\right),$

where $t=t_{0}<\ldots<t_{n}=T$ and
$f:\mathbb{R}\times\mathbb{R}\rightarrow\mathbb{R}$ is Lipschitz continuous
with constant $L_{f}$.

At the initial time $T=t_{n}$, we apply Chebyshev interpolation to the
function $g(T,x)$, i.e. for $x\in\mathcal{X}$,

$\displaystyle
V_{T}(x)=g(T,x)\approx\sum_{j\in\mathcal{J}}c_{j}(T)p_{j}(x)=:\widehat{V}_{T}(x)$

At the first time step $t_{n-1}$, the derivation of
$\mathbb{E}[g(t_{n},X_{t_{n}})|X_{t_{n-1}}=x]$ is replaced by
$\mathbb{E}[\sum_{j}c_{j}(t_{n})p_{j}(X_{t_{n}})|X_{t_{n-1}}=x]$ yielding

$\displaystyle V_{t_{n-1}}(x)$
$\displaystyle=f\left(g(t_{n-1},x),\mathbb{E}[V_{t_{n}}(X_{t_{n}})|X_{t_{n-1}}=x]\right)$
$\displaystyle\approx
f\Big{(}g(t_{n-1},x),\mathbb{E}\Big{[}\sum_{j\in\mathcal{J}}c_{j}(t_{n})p_{j}(X_{t_{n}})\Big{|}X_{t_{n-1}}=x\Big{]}\Big{)}$
$\displaystyle=f\Big{(}g(t_{n-1},x),\sum_{j\in\mathcal{J}}c_{j}(t_{n})\mathbb{E}\Big{[}p_{j}(X_{t_{n}})\Big{|}X_{t_{n-1}}=x\Big{]}\Big{)}.$

At time step $t_{n-1}$ the value function $V_{t_{n-1}}$ needs only to be
evaluated at the specific Chebyshev nodes. Hence, denoting with
$x^{k}=(x_{k_{1}},\ldots,x_{k_{D}})$ the Chebyshev nodes, it suffices to
evaluate

(2.6) $\displaystyle V_{t_{n-1}}(x^{k})\approx
f\Big{(}g(t_{n-1},x^{k}),\sum_{j\in\mathcal{J}}c_{j}(t_{n})\mathbb{E}\Big{[}p_{j}(X_{t_{n}})\Big{|}X_{t_{n-1}}=x^{k}\Big{]}\Big{)}=:\widehat{V}_{t_{n-1}}(x^{k}).$

A linear transformation of $(\widehat{V}_{t_{n-1}}(x^{k}))_{k\in\mathcal{J}}$
yields the Chebyshev coefficients according to (2.3) which determines the
Chebyshev interpolation $\widehat{V}_{t_{n-1}}=\sum_{j}c_{j}(t_{n-1})p_{j}$.
We apply this procedure iteratively as described in detail in Algorithm 1.

The stochastic part is gathered in the expectations of the Chebyshev
polynomials conditioned on the Chebyshev nodes, i.e.
$\Gamma_{j,k}(t_{u})=\mathbb{E}[p_{j}(X_{t_{u+1}})|X_{t_{u}}=x^{k}]$.
Moreover, if an equidistant time stepping is applied the computation can be
further simplified. If for the underlying stochastic process

(2.7)
$\displaystyle\Gamma_{j,k}(t_{u})=\mathbb{E}[p_{j}(X_{t_{u+1}})|X_{t_{u}}=x^{k}]=\mathbb{E}[p_{j}(X_{t_{1}})|X_{t_{0}}=x^{k}]=:\Gamma_{j,k}$

for $u=0,\ldots,n-1$, then the conditional expectations need to be computed
only for one time step, see Algorithm 2. One can pre-compute these conditional
expectations and thus, the method allows for an offline/online decomposition.

Algorithm 1 Dynamic Chebyshev algorithm

1:$\overline{N}\in\mathbb{N}^{D}$,
$\mathcal{X}=[\underline{x}_{1},\overline{x}_{1}]\times\ldots\times[\underline{x}_{D},\overline{x}_{D}]$,
$0=t_{0},\ldots,t_{n}=T$

2:Determine index set $\mathcal{J}$ and nodal points
$x^{k}=(x_{k_{1}},\ldots,x_{k_{D}})$

3:Pre-computation step:

4: For all $j,k\in\mathcal{J}$ and all $t_{u},\ u=0,\ldots,n-1$

5: Compute
$\Gamma_{j,k}(t_{u})=\mathbb{E}[p_{j}(X_{t_{u+1}})|X_{t_{u}}=x^{k}]$

6:Time $T$

7: $\widehat{V}_{T}(x^{k})=g(T,x^{k}),\ k\in\mathcal{J}$, derive

8:
$c_{j}(T)=D_{\overline{N}}(j)\sum_{k\in\mathcal{J}}\\!{}^{{}^{\prime\prime}}\widehat{V}_{T}(x^{k})T_{j}(z^{k})$

9: Obtain Chebyshev interpolation
$\widehat{V}_{T}(x)=\sum_{j\in\mathcal{J}}c_{j}(T)p_{j}(x)$ of $V_{T}(x)$

10:Iterative time stepping from $t_{u+1}\rightarrow t_{u},\ u=n-1,\ldots,1$

11: Given Chebyshev interpolation of
$\widehat{V}_{t_{u+1}}(x)=\sum_{j\in\mathcal{J}}c_{j}(t_{u+1})p_{j}(x)$

12: Derivation of $\widehat{V}_{t_{u}}(x^{k}),\ k\in\mathcal{J}$ at the nodal
points

13:
$\widehat{V}_{t_{u}}(x^{k})=f(g(t_{u},x^{k}),\sum_{j\in\mathcal{J}}c_{j}(t_{u+1})\Gamma_{j,k}(t_{u}))$

14: Derive
$c_{j}(t_{u})=D_{\overline{N}}(j)\sum_{k\in\mathcal{J}}\\!{}^{{}^{\prime\prime}}\widehat{V}_{t_{u}}(x^{k})T_{j}(z^{k})$

15: Obtain Chebyshev interpolation
$\widehat{V}_{t_{u}}(x)=\sum_{j\in\mathcal{J}}c_{j}(t_{u})p_{j}(x)$ of
$V_{t_{t}}(x)$

16:Deriving the solution at $t=0$

17: $\widehat{V}_{0}(x)=\sum_{j\in\mathcal{J}}c_{j}(0)p_{j}(x)$

Algorithm 2 Simplified Dynamic Chebyshev algrithm

1:Time steps $0=t_{1},\ldots,t_{n}=T$ with $\Delta t:=t_{u}-t_{u-1}$

2:Replace in Algorithm 1 Lines 2-4 with:

3:Pre-computation step:

4: Compute $\Gamma_{j,k}=\mathbb{E}[p_{j}(X_{\Delta t})|X_{0}=x^{k}]$ for all
$j,k\in\mathcal{J}$

##  3 Error Analysis

In this section we analyse the error of Algorithm 1, i.e.

(3.1)
$\displaystyle\varepsilon_{t_{u}}:=\max_{x\in\mathcal{X}}|V_{t_{u}}(x)-\widehat{V}_{t_{u}}(x)|.$

Two different error sources occur at $t_{u}$, the classical interpolation
error of the Chebyshev interpolation and a distortion error at the nodal
points. The latter comes from the fact that the values
$\widehat{V}_{t_{u}}(x^{k})$ are approximations of $V_{t_{u}}(x^{k})$. The
behaviour of the interpolation error depends on the regularity of the value
function. Here, we assume analyticity of the value function. The concept can
be extended to further cases such as assuming differentiability or piecewise
analyticity. The latter is discussed in preliminary form in (Mahlstedt, 2017,
Section 5.3) and is further investigated in a follow-up paper. Hence, we need
a convergence result for the Chebyshev interpolation which incorporates a
distortion error at the nodal points.

First, we introduce the required notation. A Bernstein ellipse
$\mathcal{B}([-1,1],\varrho)$ with $\varrho>1$ is defined as the open region
in the complex plane bounded by an ellipse with foci $\pm 1$ and semiminor and
semimajor axis lengths summing to $\varrho$. We define a generalized Bernstein
ellipse $\mathcal{B}(\mathcal{X},\varrho)$ around the hyperrectangle
$\mathcal{X}$ with parameter vector $\varrho\in(1,\infty)^{D}$ as

$\displaystyle\mathcal{B}(\mathcal{X},\varrho):=\mathcal{B}([\underline{x}_{1},\overline{x}_{1}],\varrho_{1})\times\ldots\times\mathcal{B}([\underline{x}_{D},\overline{x}_{D}],\varrho_{D})$

with
$\mathcal{B}([\underline{x},\overline{x}],\varrho):=\tau_{[\underline{x},\overline{x}]}\circ\mathcal{B}([-1,1],\varrho)$,
where for $x\in\mathds{C}$ we have the transform
$\tau_{[\underline{x},\overline{x}]}\big{(}\Re(x)\big{)}:=\overline{x}+\frac{\underline{x}-\overline{x}}{2}\big{(}1-\Re(x)\big{)}$
and
$\tau_{[\underline{x},\overline{x}]}\big{(}\Im(x)\big{)}:=\frac{\overline{x}-\underline{x}}{2}\Im(x)$
where the sets $\mathcal{B}([-1,1],\varrho_{i})$ are Bernstein ellipses for
$i=1,\ldots,D$.

######  Proposition 3.1.

Let $\mathcal{X}\ni x\mapsto f(x)$ be a real-valued function with an analytic
extension to some generalized Bernstein ellipse
$\mathcal{B}(\mathcal{X},\varrho)$ for $\varrho\in(1,\infty)^{D}$ with
$\sup_{x\in\mathcal{B}(\mathcal{X},\varrho)}|f(x)|\leq b$. Assume distorted
values $f^{\varepsilon}(x^{k})=f(x^{k})+\varepsilon(x^{k})$ with
$|\varepsilon(x^{k})|\leq\overline{\varepsilon}$ at all nodes $x^{k}$. Then

$\displaystyle\max_{x\in\mathcal{X}}\big{|}f(x)-I_{\overline{N}}(f^{\varepsilon})(x)\big{|}\leq\varepsilon_{int}(\varrho,N,D,B)+\overline{\varepsilon}\Lambda_{\overline{N}}.$

with

(3.2) $\displaystyle\varepsilon_{int}(\varrho,N,D,B):=2^{\frac{D}{2}+1}\cdot
B\cdot\bigg{(}\sum_{i=1}^{D}\varrho_{i}^{-2N_{i}}\prod_{j=1}^{D}\frac{1}{1-\varrho_{j}^{-2}}\bigg{)}^{\frac{1}{2}}$

and Lebesgue constant
$\Lambda_{\overline{N}}\leq\prod_{i=1}^{D}\big{(}\frac{2}{\pi}\log(N_{i}+1)+1\big{)}$.

###### Proof.

Using the linearity of the interpolation operator we obtain for the Chebyshev
interpolation of $f^{\varepsilon}$ with
$f^{\varepsilon}(x^{k})=f(x^{k})+\varepsilon(x^{k})$ that

$\displaystyle
I_{\overline{N}}(f^{\varepsilon})(x)=I_{\overline{N}}(f)(x)+I_{\overline{N}}(\varepsilon)(x).$

The tensor-based multivariate Chebyshev interpolation
$I_{\overline{N}}(\varepsilon)$ can be written in Lagrange form

$\displaystyle
I_{\overline{N}}(\varepsilon)(x)=\sum_{j\in\mathcal{J}}\varepsilon(x^{j})\lambda^{j}(x)\qquad\text{with}\quad\lambda^{j}(x)=\prod_{i=1}^{D}\ell_{j_{i}}(\tau^{-1}_{[\underline{x}_{i},\overline{x}_{i}]}(x_{i}))$

where $\ell_{j_{i}}(z)=\prod_{k\neq j_{i}}\frac{z-z_{k}}{z_{j_{i}}-z_{k}}$ is
the $j_{i}-$th Lagrange polynomial. This yields

$\displaystyle\max_{x\in\mathcal{X}}\left|I_{N}(\varepsilon)(x)\right|=\max_{x\in\mathcal{X}}\Big{|}\sum_{j\in\mathcal{J}}\varepsilon(x^{j})\lambda^{j}(x)\Big{|}\leq\overline{\varepsilon}\max_{x\in\mathcal{X}}\sum_{j\in\mathcal{J}}|\lambda^{j}(x)|=:\overline{\varepsilon}\Lambda_{\overline{N}}.$

The term $\Lambda_{\overline{N}}$ is the Lebesgue constant of the
(multivariate) Chebyshev nodes which is given by

$\displaystyle\Lambda_{\overline{N}}=\max_{x\in\mathcal{X}}\sum_{j\in\mathcal{J}}\big{|}\lambda^{j}(x)\big{|}=\max_{x\in\mathcal{X}}\sum_{j\in\mathcal{J}}\prod_{i=1}^{D}\big{|}\ell_{j_{i}}(x_{i})\big{|}=\prod_{i=1}^{D}\max_{x_{i}\in[\underline{x}_{i},\overline{x}_{i}]}\sum_{j_{i}=0}^{N_{i}}\Big{|}\ell_{j_{i}}\big{(}\tau^{-1}_{[\underline{x}_{i},\overline{x}_{i}]}(x_{i})\big{)}\Big{|}.$

Since
$\max_{x_{i}\in[\underline{x}_{i},\overline{x}_{i}]}\sum_{j_{i}=0}^{N_{i}}|\ell_{j_{i}}(\tau^{-1}_{[\underline{x}_{i},\overline{x}_{i}]}(x_{i}))|=\max_{z\in[-1,1]}\sum_{j_{i}=0}^{N_{i}}|\ell_{j_{i}}(z)|=\Lambda_{N_{i}}$,
which is the Lebesgue constant of the univariate Chebyshev interpolation, we
have $\Lambda_{\overline{N}}=\prod_{i=1}^{D}\Lambda_{N_{i}}$. From (Trefethen,
2013, Theorem 15.2) we obtain for the univariate Chebyshev interpolation
$\Lambda_{N}\leq\frac{2}{\pi}\log(N+1)+1$ and hence

(3.3)
$\displaystyle\Lambda_{\overline{N}}\leq\prod_{i=1}^{D}\Big{(}\frac{2}{\pi}\log(N_{i}+1)+1\Big{)}.$

For the distorted Chebyshev interpolation holds

$\displaystyle\big{|}f(x)-I_{\overline{N}}(f^{\varepsilon})(x)\big{|}\leq|f(x)-I_{\overline{N}}(f)(x)\big{|}+|I_{\overline{N}}(\varepsilon)(x)\big{|}.$

Therefore, the proposition follows directly from (3.3) and Sauter and Schwab
(2010). ∎

We use this result to investigate the error of the Dynamic Chebyshev method.
First, we introduce the following assumption.

######  Assumptions 3.2.

We assume $\mathcal{X}\ni x\mapsto V_{t_{u}}(x)$ is a real valued function
that has an analytic extension to a generalized Bernstein ellipse
$\mathcal{B}(\mathcal{X},\varrho_{t_{u}})$ with
$\varrho_{t_{u}}\in(1,\infty)^{D}$ and
$\sup_{x\in\mathcal{B}(\mathcal{X},\varrho_{t_{u}})}|V_{t_{u}}(x)|\leq
B_{t_{u}}$ for $u=1,\ldots,n$.

Proposition 3.5 provides conditions on the process $X$ and the functions $f$
and $g$ that guaranty Assumptions 3.2. Under this assumptions, we can apply
Proposition 3.1 to obtain an error bound for the Dynamic Chebyshev method at
each time step. This error bound has a recursive structure, since the values
of $V_{t_{u}}$ depend on the conditional expectation of $V_{t_{u+1}}$. The
interpolation error of the final time step is of form (3.2). At any other time
step $t_{u}$ an additional distortion error by approximating the function
values at the nodal points by

$\displaystyle V_{t_{u}}(x^{k})\approx
f\bigg{(}g(t_{u},x^{k}),\sum_{j\in\mathcal{J}}c_{j}(t_{u+1})\mathbb{E}[p_{j}(X_{t_{u+1}})|X_{t_{u}}=x^{k}]\bigg{)}=\widehat{V}_{t_{u}}(x^{k})$

comes into play. Proposition 3.1 yields

$\displaystyle\varepsilon_{t_{u}}:=\max_{x\in\mathcal{X}}|V_{t_{u}}(x)-\widehat{V}_{t_{u}}(x)|\leq\varepsilon_{int}(\varrho_{t_{u}},N,D,B_{t_{u}})+\Lambda_{\overline{N}}F_{t_{u}},$

where
$F_{t_{u}}:=\max_{j\in\mathcal{J}}|V_{t_{u}}(x_{j})-\widehat{V}_{t_{u}}(x_{j})|$.
The term $F_{t_{u}}$ depends on the function $f$ and the interpolation error
at the previous time step $t_{u+1}$.

Moreover, two additional error sources can influence the error bound. If there
is no closed-form solution for the generalized moments
$\mathbb{E}[p_{j}(X_{t_{u+1}})|X_{t_{u}}=x^{k}]$ a numerical technique, e.g.
numerical quadrature or Monte Carlo methods, introduces an additional error.
The former is typically deterministic and bounded whereas the latter is
stochastic. In order to incorporate this error in the following error analysis
we introduce some additional notation. The conditional expectation can be seen
as a linear operator which operates on the vector space of all continuous
functions $\mathcal{C}(\mathbb{R}^{D})$ with finite $L^{\infty}$-norm

$\displaystyle\Gamma_{t_{u}}^{k}:\mathcal{C}(\mathbb{R}^{D})\rightarrow\mathbb{R}\quad\text{with}\quad\Gamma_{t_{u}}^{k}(f):=\mathbb{E}[f(X_{t_{u+1}})|X_{t_{u}}=x^{k}].$

Define the subspace of all $D$ variate polynomials
$\mathcal{P}_{\overline{N}}(\mathcal{X}):=\text{span}\\{p_{j},\
j\in\mathcal{J}\\}$ equipped with the $L^{\infty}$-norm. We assume the
operator $\Gamma_{t_{u}}^{k}$ is approximated by a linear operator
$\widehat{\Gamma}_{t_{u}}^{k}:\mathcal{P}_{\overline{N}}(\mathcal{X})\rightarrow\mathbb{R}$
on $\mathcal{P}_{\overline{N}}(\mathcal{X})$ which fullfills one of the two
following conditions. For all $u=0,\ldots,n$ the approximation is either
deterministic and the error is bounded by a constant $\overline{\delta}$,

(GM)
$\displaystyle||\Gamma_{t_{u}}^{k}-\widehat{\Gamma}_{t_{u}}^{k}||_{op}:=\sup_{\begin{subarray}{c}p\in\mathcal{P}_{\overline{N}}\\\
||p||=1\end{subarray}}\left|\Gamma_{t_{u}}^{k}(p)-\widehat{\Gamma}_{t_{u}}^{k}(p)\right|\leq\overline{\delta}\quad\forall
k\in\mathcal{J}$

or the approximation is stochastic and uses $M$ samples of the underlying
process and the polynomials $p$ may have stochastic coefficients. In this case
we assume the error bound

($\text{GM}^{*}$)
$\displaystyle||\Gamma_{t_{u}}^{k}-\widehat{\Gamma}_{t_{u}}^{k}||_{op}:=\sup_{\begin{subarray}{c}p\in\mathcal{P}_{\overline{N}}\\\
||p||_{\infty}^{\star}=1\end{subarray}}\mathbb{E}\left[\left|\Gamma_{t_{u}}^{k}(p)-\widehat{\Gamma}_{t_{u}}^{k}(p)\right|\right]\leq\delta^{\star}(M)\quad\forall
k\in\mathcal{J}$

with norm $||p||_{\infty}^{\star}=\max_{x\in\mathcal{X}}\mathbb{E}[|p(x)|]$.
In order to incorporate stochasticity of $\widehat{V}_{t_{u}}(x)$, we replace
(3.1) by

(3.4)
$\displaystyle\varepsilon_{t_{u+1}}:=\max_{x\in\mathcal{X}}\mathbb{E}\left[\left|V_{t_{u}}(x)-\widehat{V}_{t_{u}}(x)\right|\right].$

Note that in the deterministic case (3.1) and (3.4) coincide. Additionally, a
truncation error is introduced by restricting to a compact interpolation
domain $\mathcal{X}$. We assume that the conditional expectation of the value
function outside this set is bounded by a constant

(TR)
$\displaystyle\mathbb{E}[V_{t_{u+1}}(X_{t_{u+1}})\mathds{1}_{\mathbb{R}^{D}\setminus\mathcal{X}}|X_{t_{u}}=x^{k}]\leq\varepsilon_{tr}.$

The following theorem provides an error bound for the Dynamic Chebyshev
method.

######  Theorem 3.3.

Let the DPP be given as in Definition 2.1. Assume the regularity Assumptions
3.2 hold and the boundedness of the truncation error (TR). Then we have

(3.5)
$\displaystyle\varepsilon_{t_{u}}\leq\sum_{j=u}^{n}C^{j-u}\varepsilon_{int}^{j}\
+\
\Lambda_{\overline{N}}L_{f}\sum_{j=u+1}^{n}C^{j-(u+1)}(\varepsilon_{tr}+\varepsilon_{gm}\overline{V}_{j})$

with with $\varepsilon_{gm}=\overline{\delta}$ if assumption (GM) holds and
$\varepsilon_{gm}=\delta^{\star}(M)$ if assumption ($\text{GM}^{*}$) holds and
$C=\Lambda_{\overline{N}}L_{f}(1+\varepsilon_{gm})$,
$\overline{V}_{j}=\max_{x\in\mathcal{X}}|V_{t_{j}}(x)|$ and
$\varepsilon_{int}^{j}=\varepsilon_{int}(\varrho_{t_{j}},N,D,B_{t_{j}})$.

###### Proof.

The proof of the theorem can be found in the appendix. ∎

The following corollary provides a simplified version of the error bound (3.5)
presenting its decomposition into three different error sources (interpolation
error $\varepsilon_{int}$, truncation error $\varepsilon_{tr}$ and the error
from the numerical computation of the generalized moments $\varepsilon_{gm}$).

######  Corollary 3.4.

Let the setting be as in Theorem 3.3. Then the error is bounded by

(3.6) $\displaystyle\varepsilon_{t_{u}}$
$\displaystyle\leq\left(\varepsilon_{int}(\underline{\varrho},N,D,\overline{B})+\varepsilon_{tr}+\varepsilon_{gm}\overline{V}\right)\tilde{C}^{n+1-u}$

with $\tilde{C}=\max\\{2,C\\}$, $\underline{\varrho}=\min_{1\leq u\leq
n}\varrho_{t_{u}}$, $\overline{B}=\max_{1\leq u\leq n}B_{t_{u}}$,
$\overline{V}=\max_{u\leq j\leq n}\overline{V}_{j}$.

Moreover, if $\varepsilon_{tr}=0$, $L_{f}=1$ and $N=N_{i}$, $i=1,\ldots,D$ the
error bound can be simplified further. Under ($\text{GM}^{*}$)
$\delta^{\star}(M)\leq c/\sqrt{M}$, $c>0$ yields

$\displaystyle\varepsilon_{t_{u}}\leq\tilde{c}_{1}\varrho^{-N}\log(N)^{D\,n}\
+\ \tilde{c}_{2}\log(N)^{D\,n}M^{-0.5}.$

for some constants $\tilde{c}_{1},\tilde{c}_{2}>0$. Under (GM) the term
$M^{-0.5}$ is replaced by $\overline{\delta}$.

###### Proof.

Assuming $C>2$ and using the geometric series, the first term in the error
bound (3.5) can be rewritten as

$\displaystyle\sum_{j=u}^{n}C^{j-u}\varepsilon_{int}^{j}\leq\overline{\varepsilon}_{int}\sum_{j=u}^{n}C^{j-u}=\overline{\varepsilon}_{int}\sum_{k=0}^{n-u}C^{k}=\overline{\varepsilon}_{int}\left(\frac{1-C^{n+1-u}}{1-C}\right)\leq\overline{\varepsilon}_{int}\,C^{n+1-u},$

where
$\overline{\varepsilon}_{int}=\max_{j}\varepsilon_{int}^{j}=\max_{j}\varepsilon_{int}(\varrho_{t_{j}},N,D,B_{t_{j}})\leq\varepsilon_{int}(\underline{\varrho},N,D,\overline{B})$
for $\underline{\varrho}=\min_{1\leq u\leq n}\varrho_{u}$ and
$\overline{B}=\max_{1\leq u\leq n}B_{t_{u}}$. For $C\leq 2$ the sum is bounded
by $\overline{\varepsilon}_{int}\,2^{n+1-u}$. Similar, we obtain for the
second term in the error bound (3.5) with
$\beta=(\varepsilon_{tr}+\varepsilon_{gm}\overline{V}_{j})$

$\displaystyle\Lambda_{\overline{N}}L_{f}\sum_{j=u+1}^{n}C^{j-(u+1)}\beta_{j}\leq\Lambda_{\overline{N}}L_{f}\,\overline{\beta}\sum_{k=0}^{n-(u+1)}C^{k}\leq\Lambda_{\overline{N}}L_{f}\,\overline{\beta}C^{n-u}\leq\overline{\beta}C^{n+1-u}$

where $\overline{\beta}=\max_{j}\beta_{j}$. Moreover, we used
$\Lambda_{\overline{N}}L_{f}\leq\Lambda_{\overline{N}}L_{f}(1+\varepsilon_{gm})=C$
in the last step. Thus, we obtain the following error bound (3.5)

$\displaystyle\varepsilon_{t_{u}}$
$\displaystyle\leq(\overline{\varepsilon}_{int}\ +\
\overline{\beta})\,\tilde{C}^{n+1-u}=\left(\varepsilon_{int}(\underline{\varrho},N,D,\overline{B})+\varepsilon_{tr}+\varepsilon_{gm}\overline{V}\right)\tilde{C}^{n+1-u},$

where $\tilde{C}=\max\\{2,C\\}$ and $\overline{V}=\max_{j}\overline{V}_{j}$,
which shows (3.6).

Furthermore, using the definition of the error bound (3.2) and $N=N_{i}$,
$i=1,\ldots,D$ we conclude that
$\varepsilon_{int}(\underline{\varrho},N,D,\overline{B})\leq
c_{1}\varrho^{-N}$ for a constant $c_{1}>0$. For the Lebesgue constant of the
Chebyshev interpolation exists a constant $c_{2}>0$ such that

$\displaystyle\Lambda_{\overline{N}}\leq\prod_{i=1}^{D}\big{(}\frac{2}{\pi}\log(N+1)+1\big{)}\leq\prod_{i=1}^{D}\big{(}\frac{4}{\pi}+1\big{)}\log(N)\leq
c_{2}\log(N)^{D}.$

Under ($\text{GM}^{*}$), $\delta^{\star}(M)\leq c/\sqrt{M}$, $c>0$ yields with
$\varepsilon_{tr}=0$, $L_{f}=1$

$\displaystyle\varepsilon_{t_{u}}$
$\displaystyle\leq\left(\varepsilon_{int}(\underline{\varrho},N,D,\overline{B})+\varepsilon_{tr}+\varepsilon_{gm}\overline{V}\right)\left(\Lambda_{\overline{N}}L_{f}(1+\varepsilon_{gm})\right)^{n+1-u}$
$\displaystyle\leq\left(c_{1}\varrho^{-N}+c\overline{V}M^{-0.5}\right)\left(c_{2}\log(N)^{D}(1+cM^{-0.5})\right)^{n}$
$\displaystyle\leq\tilde{c}_{1}\varrho^{-N}\log(N)^{D\,n}\ +\
\tilde{c}_{2}\log(N)^{D\,n}M^{-0.5}$

and this converges towards zero for $N\rightarrow\infty$ if
$\sqrt{M}>\log(N)^{D\,n}$. If (GM) holds we have
$\varepsilon_{gm}=\overline{\delta}$ and the term $M^{-0.5}$ is replaced by
$\overline{\delta}$. ∎

The following proposition provides conditions under which the value function
has an analytic extension to some generalized Bernstein ellipse and
Assumptions 3.2 hold.

######  Proposition 3.5.

Consider a DPP as defined in (2.4) and (2.5) with equidistant time-stepping
and $g(t,x)=g(x)$. Let $X=(X_{t})_{0\leq t\leq T}$ be a Markov process with
stationary increments. Assume $e^{\langle\eta,\cdot\rangle}g(\cdot)\in
L^{1}(\mathbb{R}^{D})$ for some $\eta\in\mathbb{R}^{D}$ and $g$ has an
analytic extension to the generalized Bernstein ellipse
$\mathcal{B}(\mathcal{X},\varrho_{g})$. Furthermore, assume
$f:\mathbb{R}\times\mathbb{R}\rightarrow\mathbb{R}$ has an analytic extension
to $\mathbb{C}^{2}$. If

  * (i)

the characteristic function $\varphi^{x}$ of $X_{\Delta t}$ with $X_{0}=x$ is
in $L^{1}(\mathbb{R}^{D})$ for every $x\in\mathcal{X}$,

  * (ii)

for every $z\in\mathbb{R}^{D}$ the mapping $x\mapsto\varphi^{x}(z-i\eta)$ has
an analytic extension to $B(\mathcal{X},\varrho_{\varphi})$ and there are
constants $\alpha\in(1,2]$ and $c_{1},c_{2}>0$ such that $\sup_{x\in
B(\mathcal{X},\varrho_{\varphi})}|\varphi^{x}(z)|\leq
c_{1}e^{-c_{2}|z|^{\alpha}}$ for all $z\in\mathbb{R}^{D}$,

then the value function $x\mapsto V_{t_{u}}(x)$ of the DPP has an analytic
extension to $B(\mathcal{X},\varrho)$ with $\varrho=\varrho_{g}$.

###### Proof.

At $T$ the value function $x\mapsto V_{T}(x)$ is analytic since
$V_{T}(x)=g(x)$ and $g$ has an analytic extension by assumption. Moreover,
$e^{\langle\eta,\cdot\rangle}g(\cdot)\in L^{1}(\mathbb{R}^{D})$ for some
$\eta\in\mathbb{R}^{D}$. We assume
$e^{\langle\eta,\cdot\rangle}V_{t_{u+1}}(\cdot)\in L^{1}(\mathbb{R}^{D})$ and
$V_{t_{u+1}}$ has an analytic extension to $B(\mathcal{X},\varrho)$. Then the
function

$\displaystyle x\mapsto
V_{t_{u}}(x)=f\left(g(x),\mathbb{E}[V_{t_{u+1}}(X_{t_{u+1}})|X_{t_{u}}=x]\right)$

is analytic if
$x\mapsto\mathbb{E}[V_{t_{u+1}}(X_{t_{u+1}})|X_{t_{u}}=x]=\mathbb{E}[V_{t_{u+1}}(X_{\Delta
t}^{x})]$ has an analytic extension. From (Gaß et al., 2018, Conditions 3.1)
we obtain conditions (A1)-(A4) under which a function of the form
$(p^{1},p^{2})\mapsto\mathbb{E}[f^{p^{1}}(X^{p^{2}})]$ is analytic. In our
case we only have the parameter $p^{2}=x$ and so $X^{p^{2}}=X_{\Delta t}^{x}$.
Condition (A1) is satisfied since
$e^{\langle\eta,\cdot\rangle}V_{t_{u+1}}(\cdot)\in L^{1}(\mathbb{R}^{D})$ and
for (A2) we have to verify that $|\widehat{V}_{t_{u+1}}(-z-i\eta)|\leq
c_{1}e^{c_{2}|z|}$ for constants $c_{1},c_{2}>0$.

$\displaystyle|\widehat{V}_{t_{u+1}}(-z-i\eta)|$
$\displaystyle=\Big{|}\int\limits_{\mathbb{R}^{D}}e^{i\langle
y,-z-i\eta\rangle}V_{t_{u+1}}\text{d}y\Big{|}$
$\displaystyle\leq\int\limits_{\mathbb{R}^{D}}|e^{-i\langle
y,z\rangle}|\left|e^{\langle y,\eta\rangle}V_{t_{u+1}}(y)\right|\text{d}y$
$\displaystyle\leq||e^{\langle\eta,\cdot\rangle}V_{t_{u+1}}(\cdot)||_{L^{1}}$

and thus (A2) holds. The remaining conditions (A3)-(A4) are equivalent to our
conditions (i)-(ii) and (Gaß et al., 2018, Theorem 3.2) yields the analyticity
of $x\mapsto\mathbb{E}[V_{t_{u+1}}(X_{\Delta t}^{x})]$ on the Bernstein
ellipse $\mathcal{B}(\mathcal{X},\varrho_{\varphi})$. Hence, $x\mapsto
V_{t_{u}}(x)$ is a composition of analytic functions and therefore analytic on
the intersection of the domains of analyticity
$\mathcal{B}(\mathcal{X},\varrho_{\varphi})\cap\mathcal{B}(\mathcal{X},\varrho_{g})=\mathcal{B}(\mathcal{X},\varrho)$
with $\varrho=\min\\{\varrho_{g},\varrho_{\varphi}\\}$.

It remains to prove that $e^{\langle\eta,\cdot\rangle}V_{t_{u}}(\cdot)\in
L^{1}(\mathbb{R}^{D})$. Here the Lipschitz continuity of $f$ yields

$\displaystyle||e^{\langle\eta,\cdot\rangle}V_{t_{u}}(\cdot)||_{L^{1}}\leq
L_{f}\left(||e^{\langle\eta,\cdot\rangle}g(\cdot)||_{L^{1}}+||e^{\langle\eta,\cdot\rangle}V_{t_{u+1}}(\cdot)||_{L^{1}}\right)<\infty.$

∎

Often, the discrete time problem (2.4) and (2.5) is an approximation of a
continuous time problem and thus, we are interested in the error behaviour for
$n\rightarrow\infty$.

######  Remark 3.6.

Assume the setup of Corollary 3.4. Moreover, assume that
$\varepsilon_{tr}=\varepsilon_{gm}=0$. If we let $N$ and $n$ go to infinity,
we have to ensure that the error bound tends to zero. We use that
$\varepsilon_{int}(\underline{\varrho},N,D,\overline{B})\leq
C_{1}\varrho^{-\underline{N}}$ for a constant $C_{1}>0$ and
$\underline{N}=\min_{i}N_{i}$. The following condition on the relation between
$n$ and $N$ ensures convergence

$\displaystyle
n<\frac{\log(\varrho)}{C_{1}D}\cdot\frac{\underline{N}}{\log(\Lambda_{\overline{N}})+\log(L_{f})}\,+\,1.$

##  4 Implementational aspects of the Dynamic Chebyshev method

In this section we discuss several approaches to compute the generalized
moments (2.7) which contain the model dependent part. Moreover, preparing the
numerical experiments we tailor the Dynamic Chebyshev method to the pricing of
American put options.

###  4.1 Derivation of generalized moments

Naturally, the question arises how the generalized moments (2.7) can be
derived. Here, we present four different ways and illustrate all approaches in
the one-dimensional case $\mathcal{X}=[\underline{x},\overline{x}]$. Similar
formulas can be obtained for a multidimensional domain.

Probability density function

For the derivation of $\mathbb{E}[p_{j}(X_{t_{u+1}})|X_{t_{u}}=x_{k}]$, let
the density function of the random variable $X_{t_{u+1}}|X_{t_{u}}=x_{k}$ be
given as $f^{u,k}(x)$. Then, the conditional expectation can be derived by
solving an integral,

$\displaystyle\mathbb{E}[p_{j}(X_{t_{u+1}})|X_{t_{u}}=x_{k}]=\int_{\underline{x}}^{\overline{x}}T_{j}(\tau^{-1}_{[\underline{x},\overline{x}]}(y))\,f^{u,k}(y)\text{d}y$

using $p_{j}(y)=T_{j}(\tau^{-1}_{\mathcal{X}}(y))1_{\mathcal{X}}(y)$. This
approach is rather intuitive and easy to implement.

Fourier Transformation

Assume the process $X$ has stationary increments and the characteristic
function $\varphi$ of $X_{\Delta t}$ is explicitly available. We apply
Parseval’s identity, see Rudin (1973), and use Fourier transforms

$\displaystyle\mathbb{E}[p_{j}(X_{t_{u+1}})|X_{t_{u}}=x_{k}]=\int_{-\infty}^{\infty}p_{j}(x+x_{k})F(\text{d}x)=\frac{1}{2\pi}\int_{-\infty}^{\infty}\widehat{p_{j}^{x_{k}}}(\xi)\varphi(-\xi)\text{d}\xi,$

where $p_{j}^{x_{k}}(x)=p_{j}(x+x_{k})$. Using the definition of
$\tau_{[\underline{x},\overline{x}]}$, we can express the Fourier transform of
$p_{j}^{x_{k}}(x)$ with the help of the Chebyshev polynomial $T_{j}(y)$. This
yield

(4.1)
$\displaystyle\mathbb{E}[p_{j}(X_{t_{u+1}})|X_{t_{u}}=x_{k}]=\frac{1}{2\pi}e^{-i\xi
x^{k}}e^{i\xi(\overline{x}-\frac{\overline{x}-\underline{x}}{2})}\frac{\overline{x}-\underline{x}}{2}\int_{-\infty}^{\infty}\widehat{T}_{j}\Big{(}\frac{\overline{x}-\underline{x}}{2}\xi\Big{)}\varphi(-\xi)\text{d}\xi.$

The Fourier transform of the Chebyshev polynomials $\widehat{T_{j}}$ are
presented in Dominguez et al. (2011) and the authors also provide a Matlab
implementation.

Truncated moments

In this approach, we use that each one-dimensional Chebyshev polynomial can be
represented as a sum of monomials, i.e.

$\displaystyle T_{j}(x)=\sum_{l=0}^{j}a_{l}x^{l},\ j\in\mathbb{N}.$

The coefficients $a_{l},\ l=0,\ldots,j$, can easily be derived using the
chebfun function poly(), see Driscoll et al. (2014). Then,

$\displaystyle\mathbb{E}[p_{j}(X_{t_{u+1}})|X_{t_{u}}=x_{k}]$
$\displaystyle=\mathbb{E}[T_{j}(\tau^{-1}_{\mathcal{X}}(X_{t_{u+1}}))1_{\mathcal{X}}(X_{t_{u+1}})|X_{t_{u}}=x_{k}]$
$\displaystyle=\sum_{l=0}^{j}a_{l}\mathbb{E}[(\tau^{-1}_{\mathcal{X}}(X_{t_{u+1}}))^{l}1_{\mathcal{X}}(X_{t_{u+1}})|X_{t_{u}}=x_{k}]$

As $\tau_{\mathcal{X}}$ is linear the computation of the generalized moments
has thus been reduced to deriving truncated moments.

Monte-Carlo simulation

Lastly, especially in cases for which neither a probability density function,
nor a characteristic function of the underlying process is given, Monte-Carlo
simulation is a suitable choice. For every nodal point $x_{k}$ one simulates
$N_{MC}$ paths $X_{t_{u+1}}^{i}$ of $X_{t_{u+1}}$ with starting value
$X_{t_{u}}=x_{k}$. These simulations can then be used to approximate

$\displaystyle\Gamma_{t_{u},t_{u+1}}(p_{j})(x^{k})=\mathbb{E}[p_{j}(X_{t_{u+1}})|X_{t_{u}}=x_{k}]\approx\frac{1}{N_{MC}}\sum_{i=1}^{N_{MC}}p_{j}(X_{t_{u+1}}^{i})$

for every $j\in\mathcal{J}$. For an overview of Monte-Carlo simulation from
SDEs and variance reduction techniques we refer to Glasserman (2003).

###  4.2 American Put Option

In the numerical section we use the Dynamic Chebyshev method to price American
put options. Assuming an asset model of the form $S_{t}=e^{X_{t}}$, the DPP
becomes

$\displaystyle V_{T}(x)$ $\displaystyle=(K-e^{x})^{+}$ $\displaystyle
V_{t_{u}}(x)$
$\displaystyle=\max\left\\{(K-e^{x})^{+},e^{-r(t_{u+1}-t_{u})}\mathbb{E}[V_{t_{u+1}}(X_{t_{u+1}})|X_{t_{u}}=x]\right\\}.$

Typically, the support of the underlying process $X_{t}$ is $\mathbb{R}$ and
restricting the domain to $\mathcal{X}$ introduces a truncation error. We
reduce this error by exploiting the asymptotic behaviour of the payoff. If
$X_{t_{u}}$ is below the exercise boundary the option is exercised at the
value $K-e^{X_{t_{u}}}$ which we exploit for $X_{t_{u}}<\underline{x}$. The
function $x\mapsto V_{t_{u}}$ tends to zero from above for
$x\rightarrow\infty$ and thus for $\overline{x}$ large enough the truncation
to zero for $x>\overline{x}$ is justified. Hence, we introduce the following
modification of the Dynamic Chebyshev method:

$\displaystyle V_{t_{u+1}}(x)$
$\displaystyle=V_{t_{u+1}}(x)1_{\\{x<\underline{x}\\}}+V_{t_{u+1}}(x)1_{\\{x\in\mathcal{X}\\}}+V_{t_{u+1}}(x)1_{\\{x>\overline{x}\\}}$
$\displaystyle\approx(K-e^{x})1_{\\{x<\underline{x}\\}}+\widehat{V}_{t_{u+1}}(x)1_{\\{x\in\mathcal{X}\\}}$

and thus

$\displaystyle\mathbb{E}[V_{t_{u+1}}(X_{t_{u+1}})|X_{t_{u}}=x_{k}]$
$\displaystyle\approx\mathbb{E}[(K-e^{X_{t_{u+1}}})1_{\\{X_{t_{u+1}}<\underline{x}\\}}|X_{t_{u}}=x_{k}]$
$\displaystyle\quad+\sum_{j=0}^{N}c_{j}(t_{u+1})\Gamma_{j,k}(t_{u})$

for $\underline{x}$ small and $\overline{x}$ large enough. One can precompute
$\mathbb{E}[(K-e^{X_{t_{u+1}}})1_{\\{X_{t_{u+1}}<\underline{x}\\}}|X_{t_{u}}=x_{k}]$.
We emphasize that similar modifications to reduce the truncation error can be
found for other payoff profiles, e.g. for digitals, butterfly options or any
other combination of different put options.

Moreover, we also modify the first time step. Instead of approximating the
payoff with Chebyshev polynomials at $t_{n}=T$ we just use the price of a
European option to compute the continuation value at $t_{n-1}$. The kink of
the payoff is in this case ”smoothed” and convergence is improved.

The option’s sensitivities Delta and Gamma can be computed by tanking the
first or second derivative of

$\displaystyle
S\mapsto\widehat{V}_{0}(\log(S))=\sum_{j\in\mathcal{J}}c_{j}(t_{0})p_{j}(\log(S)).$

Thus Delta and Gamma are expressed as the sum of derivatives of Chebyshev
polynomials. In particular, their derivation comes without any additional
computational costs.

##  5 Numerical experiments

In this section, we use the Dynamic Chebyshev method to price American put
options and we numerically investigate the convergence of the method.
Moreover, we compare the method with the Least-squares Monte-Carlo method of
Longstaff and Schwartz (2001).

###  5.1 Stock price models

For the convergence analysis we use three different stock price models.

The Black-Scholes model:
In the classical model of Black and Scholes (1973) the stock price process is
modelled by the SDE

$\displaystyle\text{d}S_{t}=rS_{t}\text{d}t+\sigma S_{t}\text{d}W_{t}$

where $r$ is the risk-free interest rate and $\sigma>0$ is the volatility. In
this model the log-returns $X_{t}=\log(S_{t})$ are normally distributed and
for the double truncated moments

$\displaystyle\mathbb{E}[X^{m}1_{[a,b]}(X)]\qquad\text{for}\qquad
X\sim\mathcal{N}(\mu,\sigma^{2})$

explicit formulas are available.Kan and Robotti (2017) present results for the
(multivariate) truncated moments and provide a Matlab implementation.

The Merton jump diffusion model:
The jump diffusion model introduced by Merton (1976) adds jumps to the
classical Black-Scholes model. For $S_{t}=S_{0}e^{X_{t}}$ the log-returns
$X_{t}$ follow a jump diffusion with volatility $\sigma$ and added jumps
arriving at rate $\lambda>0$ with normal distributed jump sizes according to
$\mathcal{N}(\alpha,\beta^{2})$. The characteristic function of $X_{t}$ is
given by

$\displaystyle\varphi(z)=exp\left(t\left(ibz-\frac{\sigma^{2}}{2}z^{2}+\lambda\left(e^{iz\alpha-\frac{\beta^{2}}{2}z^{2}}-1\right)\right)\right)$

with risk-neutral drift
$b=r-\frac{\sigma^{2}}{2}-\lambda\Big{(}e^{\alpha+\frac{\beta^{2}}{2}}-1\Big{)}$.

The Constant Elasticity of Variance model:
The Constant Elasticity of Variance model (CEV) as stated in Schroder (1989)
is a local volatility model based on the stochastic process

(5.1) $\displaystyle\text{d}S_{t}=rS_{t}\text{d}t+\sigma
S_{t}^{\beta/2}\text{d}W_{t}\quad\text{for}\quad\beta>0.$

Hence the stock volatility $\sigma S_{t}^{(\beta-2)/2}$ depends on the current
level of the stock price. For the special case $\beta=2$ the model coincides
with the Black-Scholes model. However, from market data one typically observes
a $\beta<2$. The CEV-model is one example of a model which has neither a
probability density, nor a characteristic function in closed-form.

###  5.2 Convergence analysis

In this section we investigate the convergence of the Dynamic Chebyshev
method. We price American put options along with the options’ Delta and Gamma
in the Black-Scholes and the Merton jump diffusion model, where we can use the
COS method of Fang and Oosterlee (2009) as benchmark. The COS method is based
on the Fourier-cosine expansion of the density function and provides fast and
accurate results for the class of Lévy models.

For the experiments, we use the following parameter sets in the Black-Scholes
model

$\displaystyle K=100,\quad r=0.03,\quad\sigma=0.25,\quad T=1,$

and for the Merton jump diffusion model

$\displaystyle K=100,\quad
r=0.03,\quad\alpha=-0.5,\quad\beta=0.4,\quad\sigma=0.25,\quad\lambda=0.4$

and we use $32$ time steps.

For both models the generalized moments are computed by the Fourier approach
as stated in (4.1). We truncate the integral at $|\xi|\leq 250$ and use
Clenshaw-Curtis with $500$ nodes for the numerical integration. For the
Fourier transform of the Chebyshev polynomials the implementation of Dominguez
et al. (2011) is used. We run the Dynamic Chebyshev method for an increasing
number of Chebyshev nodes $N=50,100,\ldots,750$. Then, option prices and their
sensitivities delta and gamma are calculated on a grid of different values of
$S_{0}$ equally distributed between $60\%$ and $140\%$ of the strike $K$. The
resulting prices and Greeks are compared using the COS method as benchmark and
the maximum error over the grid is calculated. Here we use the implementation
provided in von Sydow et al. (2015).

Figure 5.1 shows the error decay for the Black-Scholes model (left hand side)
and the Merton model (right hand side). We observe that the method converges
and an error below $10^{-3}$ is reached for $N=300$ Chebyshev nodes. The
experiments confirm that the method can be used for an American put option.

Figure 5.1: Error decay prices Dynamic Chebyshev in the BS model (left) and
the Merton model (right). The conditional expectation of the Chebyshev
polynomials are calculated with the Fourier transformation.

###  5.3 Dynamic Chebyshev with Monte-Carlo

So far we empirically investigated the error decay of the method for option
price and their derivatives. In this section we compare the Dynamic Chebyshev
method with the Least Square Monte-Carlo approach of Longstaff and Schwartz
(2001) in terms of accuracy and runtime.

####  5.3.1 The Black-Scholes model

As a first benchmark, we use the Black-Scholes model with an interest rate of
$r=0.03$ and volatility $\sigma=0.25$. Here, we look at a whole option price
surface with varying maturities and strikes. We choose $12$ different
maturities between one month and four years given by

$\displaystyle T\in\\{1/12,2/12,3/12,6/12,9/12,1,15/12,18/12,2,30/12,3,4\\}$

and strikes equally distributed between $80\%$ and $120\%$ of the current
stock price $S_{0}=100$ in steps of $5\%$. We fix $n=504$ time steps (i.e.
exercise rights) per year.

We compare the Dynamic Chebyshev method to the Least Squares Monte-Carlo
approach. We run both methods for an increasing number of Monte-Carlo paths
according to

(5.2) $\displaystyle M\in\left\\{2500,5000,10000,20000,40000,80000\right\\}.$

The convergence of the DC method depends on both, the number of nodes $N$ and
the number of Monte-Carlo paths $M$. For an optimal convergence behaviour one
needs to find a reasonable relationship between these factors. For the
following experiments, we fix the number of Chebyshev nodes as
$N=\sqrt{2}\sqrt{M}$.

Figure 5.2 shows the price surface and the error surface for $N=400$ and
$M=80000$. The error was estimated by using the COS method as benchmark. We
reach a maximal error below $4*10^{-2}$ on the whole option surface.

In Figure 5.3 the $log_{10}$-error is shown as a function of the
$log_{10}$-runtime for both methods. The left figure compares the total
runtimes and the right figure compares the offline runtime. For the Dynamic
Chebyshev method the total runtime includes the offline-phase and the online
phase. The offline-phase consists of the simulation of one time step of the
underlying asset process $X_{\Delta t}$ for $N+1$ starting values
$X_{0}=x_{k}$ and of the computation of the conditional expectations
$E[p_{j}(X_{\Delta t})|X_{0}=x_{k}]$ for $j,k=0,\ldots,N$. The online phase is
the actual pricing of the American option for all strikes and maturities.
Similar, the total runtime of the Least-Square Monte-Carlo method includes the
simulation of the Monte-Carlo paths (offline-phase) and the pricing of the
option via backward induction (online-phase).

We observe that the Dynamic Chebyshev method reaches the same accuracy in a
much lower runtime. For example, a maximum error of $0.1$ is reached in a
total runtime of $0.5s$ with the Dynamic Chebyshev method whereas the LSM
approach needs $98s$. This means the Dynamic Chebyshev method is nearly $200$
times faster for the same accuracy. For the actual pricing in the online
phase, the gain in efficiency is even higher. We observe that the Dynamic
Chebyshev method outperforms the Least-Square Monte-Carlo method in terms of
the total runtime and the pure online runtime. Moreover, we observe that the
performance gain from splitting the computation into an offline and an online
phase is much higher for the Dynamic Chebyshev method. For instance, in the
example above the online runtime of the Dynamic Chebyshev method is $0.05s$
whereas the LSM takes $95s$, a factor of $1900$ times more.

Figure 5.2: Price surface and corresponding error of the Dynamic Chebyshev
method in the Black-Scholes model. The conditional expectations are calculated
with Monte-Carlo.

Figure 5.3: Log-Log plot of the total/online runtime vs. accuracy. Comparison
of the Dynamic Chebyshev method with the Least-Square Monte-Carlo algorithm.

The main advantage of the Dynamic Chebyshev method is that once the
conditional expectations are calculated, they can be used to price the whole
option surface. The pure pricing, i.e. the online phase, is highly efficient.
Furthermore, one only needs to simulate one time step $\Delta t$ of the
underlying stochastic process instead of the complete path. We investigate
this efficiency gain by varying the number of options and the number of time
steps (exercise rights). Figure 5.4 compares the total runtime of the DC
method with the total runtime of the LSM method for an increasing number of
options and for an increasing number time steps. As expected, we can
empirically confirm that the efficiency gain by the Dynamic Chebyshev methods
increases with number of options and the number of exercise rights. In both
cases, the runtime of the DC method stays nearly constant whereas the runtime
of the LSM method increases approximately linearly.

Figure 5.4: Total runtime of the DC and the LSM method for an increasing
number of options (left) and an increasing number of timesteps (right).

####  5.3.2 The CEV model

Next, we use the constant elasticity of variance (CEV) model for the
underlying stock price process. We perform the same experiments as in the last
section. The parameters in the CEV model, as in (5.1) are the following

$\displaystyle\sigma=0.25,\quad r=0.03,\quad\beta=1.5.$

Similarly, we compare the Dynamic Chebyshev and the LSM method by computing
the prices of an option price surface. We use the same parameter
specifications for $K$, $T$ and $n$. We run both methods for an increasing
number of Monte-Carlo simulations $M$ and fix $N=\sqrt{2}\sqrt{M}$. Figure 5.5
shows the price surface and the error surface for $N=400$ and $M=80000$. The
error is calculated using a binomial tree implementation of the CEV model
based on Nelson and Ramaswamy (1990).

Figure 5.5: Price surface and corresponding error of the Dynamic Chebyshev
method in the CEV model. The conditional expectations are calculated with
Monte-Carlo.

In Figure 5.6 the $log_{10}$-error is shown as a function of the
$log_{10}$-runtime for both methods. The left figure compares the total
runtimes and the right figure compares the offline runtimes. Again, we observe
that the Dynamic Chebyshev method is faster for the same accuracy and it
profits more from an offline-online decomposition. For example, the total
runtime of the DC method to reach an accuracy below $0.03$ is $3.5s$ whereas
LSM takes $136s$. For the online runtimes this out-performance is $1s$ to
$122s$.

Figure 5.6: Log-Log plot of the total/online runtime vs. accuracy. Comparison
of the Dynamic Chebyshev method with the Least-Square Monte-Carlo algorithm.

Investigating this efficiency gain further, we look at the performance for
different numbers of options and time steps (exercise rights). Similarly to
the last section, Figure 5.7 compares the total runtime of the DC method with
the total runtime of the LSM method for an increasing number of options and
time steps. In both cases, the runtime of the DC method stays nearly constant
whereas the runtime of the LSM method increases approximately linearly. This
observation is consistent with the findings for the Black-Scholes model in the
previous section.

Figure 5.7: Total runtime of the DC and the LSM method for an increasing
number of options (left) and an increasing number of timesteps (right).

##  6 Conclusion and Outlook

We have introduced a new approach to price American options via backward
induction by approximating the value function with Chebyshev polynomials.
Thereby, the computation of the conditional expectation of the value function
in each time step is reduced to the computation of conditional expectations of
polynomials. The proposed method separates the pricing of an option into a
part which is model-dependent (the computation of the conditional
expectations) and the pure pricing of a given payoff which becomes independent
of the underlying model. The first step, the computation of the conditional
expectation of the Chebyshev polynomials, is the so-called offline phase of
the method. The design of the method admits several qualitative advantageous:

  * •

If the conditional expectations are set-up once, we can use them for the
pricing of many different options. Thus, the actual pricing in the online step
becomes very simple and fast.

  * •

In the pre-computation step one can combine the method with different
techniques, such as Fourier approaches and Monte-Carlo simulation. Hence the
method can be applied in a variety of models.

  * •

The proposed approach is very general and flexible and thus not restricted to
the pricing of American options. It can be used to solve a large class of
optimal stopping problems.

  * •

We obtain a closed-form approximation of the option price as a function of the
stock price at every time step. This approximation can be used to compute the
option’s sensitivities Delta and Gamma at no additional costs. This holds for
all models and payoff profiles, even if Monte-Carlo is required in the offline
phase.

  * •

The method is easy to implement and to maintain. The pre-computation step is
well-suited for parallelization to speed-up the method.

We have investigated the theoretical error behaviour of the method and
introduced explicit error bounds. We put particular emphasis on the
combination of the method with Monte-Carlo simulation. Numerical experiments
confirm that the method performs well for the pricing of American options. A
detailed comparison of the method with the Least-Square Monte-Carlo approach
proposed by Longstaff and Schwartz (2001) confirmed a high efficiency gain.
Especially, when a high number of options is priced, for example a whole
option price surface. In this case, the Dynamic Chebyshev method highly
profits from the offline-online decomposition. Once the conditional
expectations are computed they can be used to price options with different
maturities and strikes. Besides the efficiency gain, the closed-form
approximation of the price function is a significant advantage because it
allows us to calculate the sensitivities. Since Longstaff and Schwartz (2001)
introduced their method different modifications have been introduced. Either
to increase efficiency or to tackle the sensitivities. For example the
simulation algorithm of Jain and Oosterlee (2015) is comparable to LSM in
terms of efficiency but is able to compute the Greeks at no additional costs.
Moreover dual approaches were developed to obtain upper bounds for the option
price, see Rogers (2002) and more recently Belomestny et al. (2018).

The presented error analysis of the method under an analyticity assumption is
the starting point for further theoretical investigations in the case of
piecewise analyticity and (piecewise) differentiability. The former allows to
cover rigorously the American option pricing problem and a preliminary version
is presented in Mahlstedt (2017). The qualitative merits of the method can be
exploited in a variety of applications. Glau et al. (2018) take advantage of
the closed-form approximation to efficiently compute the expected exposure of
early-exercise options as a step in CVA calculation. Moreover, the method can
be used to price different options such as Barrier options, Swing options or
multivariate American options.

##  Appendix A Proof of Theorem 3.3

###### Proof.

Consider a DPP as defined in Definition 2.1, i.e. we have a Lipschitz
continuous function

$\displaystyle|f(x_{1},y_{1})-f(x_{2},y_{2})|\leq
L_{f}(|x_{1}-x_{2}|+|y_{1}-y_{2}|).$

Assume that the regularity Assumption 3.2 and the assumption on the truncation
error (TR) hold. Then we have to distinguish between the deterministic case
(GM) and the stochastic case ($\text{GM}^{*}$). In the first case, the
expectation in the error bound can simply be ignored. First, we apply
Proposition 3.1. At time point $T$ there is no random part and no distortion
error. Thus,

$\displaystyle\max_{x\in\mathcal{X}}\mathbb{E}\left[|V_{T}(x)-\widehat{V}_{T}(x)|\right]=\max_{x\in\mathcal{X}}|V_{T}(x)-\widehat{V}_{T}(x)|\leq\varepsilon_{int}(\varrho_{t_{n}},N,D,M_{t_{n}}).$

For the ease of notation we will from know on write
$\varepsilon^{j}_{int}=\varepsilon_{int}(\varrho_{t_{j}},N,D,M_{t_{j}})$. We
obtain for the error at $t_{u}$

(A.1)
$\displaystyle\max_{x\in\mathcal{X}}\mathbb{E}\left[|V_{t_{u}}(x)-\widehat{V}_{t_{u}}(x)|\right]\leq\varepsilon^{u}_{int}+\Lambda_{\overline{N}}F(f,t_{u})$

with maximal distortion error
$F(f,t_{u})=\max_{k\in\mathcal{J}}\mathbb{E}\left[|V_{t_{u}}(x^{k})-\widehat{V}_{t_{u}}(x^{k})|\right]$.
Note that whether (GM) or ($\text{GM}^{*}$) hold an approximation error of the
conditional expectation of $\widehat{V}_{t_{u+1}}$ is made, i.e.
$\mathbb{E}[\widehat{V}_{t_{u+1}}(X_{t_{u+1}})|X_{t_{u}}=x^{k}]=\Gamma_{t_{u}}^{k}(\widehat{V}_{t_{u+1}})\approx\widehat{\Gamma}_{t_{u}}^{k}(\widehat{V}_{t_{u+1}})$.
The Lipschitz continuity of $f$ yields

$\displaystyle\left|V_{t_{u}}(x^{k})-\widehat{V}_{t_{u}}(x^{k})\right|$
$\displaystyle=\Big{|}f\left(g(t_{u},x^{k}),\Gamma_{t_{u}}^{k}(V_{t_{u+1}})\right)\
-f\Big{(}g(t_{u},x^{k}),\widehat{\Gamma}_{t_{u}}^{k}(\widehat{V}_{t_{u+1}})\Big{)}\Big{|}$
$\displaystyle\leq L_{f}\Big{(}\Big{|}g(t_{u},x^{k})-g(t_{u},x^{k})\Big{|}\
+\Big{|}\Gamma_{t_{u}}^{k}(V_{t_{u+1}})-\
\widehat{\Gamma}_{t_{u}}^{k}(\widehat{V}_{t_{u+1}})\Big{|}\Big{)}$
$\displaystyle=L_{f}\Big{(}\Big{|}\Gamma_{t_{u}}^{k}(V_{t_{u+1}})-\
\widehat{\Gamma}_{t_{u}}^{k}(\widehat{V}_{t_{u+1}})\Big{|}\Big{)}$
$\displaystyle\leq
L_{f}\Big{(}\Big{|}\Gamma_{t_{u}}^{k}(V_{t_{u+1}}1_{\mathcal{X}})-\Gamma_{t_{u}}^{k}(\widehat{V}_{t_{u+1}})\Big{|}+\Big{|}\Gamma_{t_{u}}^{k}(V_{t_{u+1}}1_{\mathbb{R}^{D}\setminus\mathcal{X}})\Big{|}$
$\displaystyle\qquad+\Big{|}\Gamma_{t_{u}}^{k}(\widehat{V}_{t_{u+1}})-\widehat{\Gamma}_{t_{u}}^{k}(\widehat{V}_{t_{u+1}})\Big{|}\Big{)}.$

Next, we consider the expectation for each of the three error terms. For the
first term we obtain

$\displaystyle\mathbb{E}\left[\Big{|}\Gamma_{t_{u}}^{k}(V_{t_{u+1}}1_{\mathcal{X}})-\Gamma_{t_{u}}^{k}(\widehat{V}_{t_{u+1}})\Big{|}\right]$
$\displaystyle=\mathbb{E}\left[\Big{|}\mathbb{E}[V_{t_{u+1}}(X_{t_{u+1}})\mathds{1}_{\mathcal{X}}-\widehat{V}_{t_{u+1}}(X_{t_{u+1}})|X_{t_{u}}=x^{k}]\Big{|}\right]$
$\displaystyle\leq\max_{x\in\mathcal{X}}\mathbb{E}\left[|V_{t_{u+1}}(x)-\widehat{V}_{t_{u+1}}(x)|\right]=\varepsilon_{t_{u+1}}$

and for the second term we have

$\displaystyle\mathbb{E}\left[\Big{|}\Gamma_{t_{u}}^{k}(V_{t_{u+1}}1_{\mathbb{R}^{D}\setminus\mathcal{X}})\Big{|}\right]\leq\mathbb{E}\left[\varepsilon_{tr}\right]=\varepsilon_{tr}.$

For the last term we have to distinguish two cases. If we assume (GM) holds,
the operator norm yields

$\displaystyle\Big{|}\Gamma_{t_{u}}^{k}(\widehat{V}_{t_{u+1}})-\widehat{\Gamma}_{t_{u}}^{k}(\widehat{V}_{t_{u+1}})\Big{|}$
$\displaystyle=\Big{|}\left(\Gamma_{t_{u}}^{k}-\widehat{\Gamma}_{t_{u}}^{k}\right)\Big{(}\widehat{V}_{t_{u+1}}\Big{)}\Big{|}$
$\displaystyle\leq\Big{|}\Big{|}\Gamma_{t_{u}}^{k}-\widehat{\Gamma}_{t_{u}}^{k}\Big{|}\Big{|}_{op}\
\Big{|}\Big{|}\widehat{V}_{t_{u+1}}\Big{|}\Big{|}_{\infty}$
$\displaystyle\leq\overline{\delta}\
\Big{|}\Big{|}\widehat{V}_{t_{u+1}}\Big{|}\Big{|}_{\infty}.$

Next, we consider the second case and assume that ($\text{GM}^{*}$) holds.
Then we have

$\displaystyle\mathbb{E}\left[\Big{|}\Gamma_{t_{u}}^{k}(\widehat{V}_{t_{u+1}})-\widehat{\Gamma}_{t_{u}}^{k}(\widehat{V}_{t_{u+1}})\Big{|}\right]\leq\Big{|}\Big{|}\Gamma_{t_{u}}^{k}-\widehat{\Gamma}_{t_{u}}^{k}\Big{|}\Big{|}_{op}\
\Big{|}\Big{|}\widehat{V}_{t_{u+1}}\Big{|}\Big{|}_{\infty}^{\star}$
$\displaystyle\leq\delta^{\star}(M)\
\Big{|}\Big{|}\widehat{V}_{t_{u+1}}\Big{|}\Big{|}_{\infty}^{\star}.$

Hence in either case the following bound holds

$\displaystyle\mathbb{E}\left[\Big{|}\Gamma_{t_{u}}^{k}(\widehat{V}_{t_{u+1}})-\widehat{\Gamma}_{t_{u}}^{k}(\widehat{V}_{t_{u+1}})\Big{|}\right]\leq\varepsilon_{gm}\max_{x\in\mathcal{X}}\mathbb{E}\left[|\widehat{V}_{t_{u+1}}(x)|\right]$

with $\varepsilon_{gm}=\overline{\delta}$ if assumption (GM) holds and
$\varepsilon_{gm}=\delta^{\star}(M)$ if assumption ($\text{GM}^{*}$) holds. We
need an upper bound for the maximum of the Chebyshev approximation

$\displaystyle\max_{x\in\mathcal{X}}\mathbb{E}\left[|\widehat{V}_{t_{u+1}}(x)|\right]$
$\displaystyle\leq\max_{x\in\mathcal{X}}\mathbb{E}\left[|\widehat{V}_{t_{u+1}}(x)-V_{t_{u+1}}(x)|\right]+\max_{x\in\mathcal{X}}|V_{t_{u+1}}(x)|\leq\varepsilon_{t_{u+1}}+\overline{V}_{u+1}$

with $\overline{V}_{u+1}:=\max_{x\in\mathcal{X}}|V_{t_{u+1}}(x)|$. Hence, the
error bound (A.1) becomes

$\displaystyle\varepsilon_{t_{u}}\leq\varepsilon^{u}_{int}+\Lambda_{\overline{N}}L_{f}\big{(}(1+\varepsilon_{gm})\varepsilon_{t_{u+1}}+\varepsilon_{tr}+\varepsilon_{gm}\overline{V}_{u+1}\big{)}.$

By induction, we now show (3.5). For $u=n$ we have
$\varepsilon_{t_{n}}\leq\varepsilon^{n}_{int}$ and therefore (3.5) holds for
u=n. We assume that for $n,\ldots,u+1$ equation (3.5) holds. For the error
$\varepsilon_{t_{u}}$ we obtain

$\displaystyle\varepsilon_{t_{u}}$
$\displaystyle\leq\varepsilon^{u}_{int}+\Lambda_{\overline{N}}L_{f}\big{(}(1+\varepsilon_{gm})\varepsilon_{t_{u+1}}+\varepsilon_{tr}+\varepsilon_{gm}\overline{V}_{u+1}\big{)}$
$\displaystyle\leq\varepsilon^{u}_{int}+\Lambda_{\overline{N}}L_{f}\Big{(}(1+\varepsilon_{gm})\Big{(}\sum_{j=u+1}^{n}C^{j-(u+1)}\varepsilon^{j}_{int}\
+\
\Lambda_{\overline{N}}L_{f}\sum_{j=u+2}^{n}C^{j-(u+2)}(\varepsilon_{tr}+\varepsilon_{gm}\overline{V}_{j})\Big{)}$
$\displaystyle\quad+\varepsilon_{tr}+\varepsilon_{gm}\overline{V}_{u+1}\Big{)}$
$\displaystyle=\varepsilon^{u}_{int}+C\sum_{j=u+1}^{n}C^{j-(u+1)}\varepsilon^{j}_{int}\
+\Lambda_{\overline{N}}L_{f}\Big{(}C\sum_{j=u+2}^{n}C^{j-(u+2)}(\varepsilon_{tr}+\varepsilon_{gm}\overline{V}_{j})$
$\displaystyle\quad+\varepsilon_{tr}+\varepsilon_{gm}\overline{V}_{u+1}\Big{)}$
$\displaystyle=\varepsilon^{u}_{int}+\sum_{j=u+1}^{n}C^{j-u}\varepsilon^{j}_{int}\
+\Lambda_{\overline{N}}L_{f}\Big{(}\sum_{j=u+2}^{n}C^{j-(u+1)}(\varepsilon_{tr}+\varepsilon_{gm}\overline{V}_{j})+\varepsilon_{tr}+\varepsilon_{gm}\overline{V}_{u+1}\Big{)}$
$\displaystyle=\sum_{j=u}^{n}C^{j-u}\varepsilon^{j}_{int}\
+\Lambda_{\overline{N}}L_{f}\sum_{j=u+1}^{n}C^{j-(u+1)}(\varepsilon_{tr}+\varepsilon_{gm}\overline{V}_{j})$

which was our claim. ∎

## References

  * Belomestny et al. (2018) Belomestny, D., S. Häfner, and M. Urusov (2018).  Regression-based complexity reduction of the nested Monte-Carlo methods.  SIAM Journal on Financial Mathematics 9(2), 665–689.
  * Black and Scholes (1973) Black, F. and M. Scholes (1973).  The pricing of options and corporate liabilities.  Journal of Political Economy 81(3).
  * Brennan and Schwartz (1977) Brennan, M. J. and E. S. Schwartz (1977).  The valuation of American put options.  The Journal of Finance 2(32), 449–462.
  * Cai and Judd (2013) Cai, Y. and K. L. Judd (2013).  Shape-preserving dynamic programming.  Mathematical Methods of Operations Research 77(3), 407–421.
  * Dominguez et al. (2011) Dominguez, V., I. Graham, and V. Smyshlyaev (2011).  Stability and error estimates for Filon–Clenshaw–Curtis rules for highly oscillatory integrals.  IMA Journal of Numerical Analysis 31(4), 1253–1280.
  * Driscoll et al. (2014) Driscoll, T. A., N. Hale, and L. N. Trefethen (2014).  Chebfun guide.
  * Fang and Oosterlee (2009) Fang, F. and C. W. Oosterlee (2009).  Pricing early-exercise and discrete barrier options by Fourier-cosine series expansions.  Numerische Mathematik 114(1), 27.
  * Gaß et al. (2018) Gaß, M., K. Glau, M. Mahlstedt, and M. Mair (2018).  Chebyshev interpolation for parametric option pricing.  Finance and Stochastics.  https://doi.org/10.1007/s00780-018-0361-y.
  * Glasserman (2003) Glasserman, P. (2003).  Monte Carlo Methods in Financial Engineering, Volume 53.  Springer Science & Business Media.
  * Glau et al. (2018) Glau, K., R. Pachon, and C. Pötz (2018).  Chebyshev interpolation for the fast calculation of Credit Valuation Adjustments.  Working Paper.
  * Haasdonk et al. (2013) Haasdonk, B., J. Salomon, and B. Wohlmuth (2013).  A reduced basis method for the simulation of American options.  In Numerical Mathematics and Advanced Applications 2011, pp. 821–829. Springer.
  * Haentjens and in’t Hout (2015) Haentjens, T. and K. J. in’t Hout (2015).  ADI schemes for pricing American options under the Heston model.  Applied Mathematical Finance 22(3), 207–237.
  * Hilber et al. (2013) Hilber, N., N. Reich, C. Winter, and C. Schwab (2013).  Computational Methods for Qantitative Finance.  Springer.
  * Jain and Oosterlee (2015) Jain, S. and C. W. Oosterlee (2015).  The stochastic grid bundling method: Efficient pricing of Bermudan options and their Greeks.  Applied Mathematics and Computation 269, 412–431.
  * Judd (1998) Judd, K. L. (1998).  Numerical methods in economics.  MIT press.
  * Kan and Robotti (2017) Kan, R. and C. Robotti (2017).  On moments of folded and truncated multivariate normal distributions.  Journal of Computational and Graphical Statistics 26(4), 930–934.
  * Levendorskiĭ (2004) Levendorskiĭ, S. (2004).  Pricing of the American put under Lévy processes.  International Journal of Theoretical and Applied Finance 7(03), 303–335.
  * Longstaff and Schwartz (2001) Longstaff, F. A. and E. S. Schwartz (2001).  Valuing American Options by Simulation: A Simple Least-Squares Approach.  Review of Financial studies 14(1), 113–147.
  * Lord et al. (2008) Lord, R., F. Fang, F. Bervoets, and C. W. Oosterlee (2008).  A fast and accurate FFT-based method for pricing early-exercise options under Lévy processes.  SIAM Journal on Scientific Computing 30(4), 1678–1705.
  * Mahlstedt (2017) Mahlstedt, M. (2017).  Complexity Reduction for Option Pricing.  Ph. D. thesis, Technische Universität München.
  * Merton (1976) Merton, R. C. (1976).  Option pricing when underlying stock returns are discontinuous.  Journal of Financial Economics 3, 125–144.
  * Nelson and Ramaswamy (1990) Nelson, D. B. and K. Ramaswamy (1990).  Simple binomial processes as diffusion approximations in financial models.  The review of financial studies 3(3), 393–430.
  * Peskir and Shiryaev (2006) Peskir, G. and A. Shiryaev (2006).  Optimal Stopping and Free-Boundary Problems.  Springer.
  * Rogers (2002) Rogers, L. C. (2002).  Monte Carlo valuation of American options.  Mathematical Finance 12(3), 271–286.
  * Rudin (1973) Rudin, W. (1973).  Functional Analysis.  McGraw-Hill Book Co.
  * Sauter and Schwab (2010) Sauter, S. and C. Schwab (2010).  Boundary Element Methods, Translated and expanded from the 2004 German original, Volume 39.  Springer Series Computational Mathematics.
  * Schroder (1989) Schroder, M. (1989).  Computing the constant elasticity of variance option pricing formula.  The Journal of Finance 44(1), 211–219.
  * Trefethen (2013) Trefethen, L. N. (2013).  Approximation Theory and Approximation Practice.  SIAM books.
  * von Sydow et al. (2015) von Sydow, L., L. Josef Höök, E. Larsson, E. Lindström, S. Milovanović, J. Persson, V. Shcherbakov, Y. Shpolyanskiy, S. Sirén, J. Toivanen, et al. (2015).  BENCHOP–The BENCHmarking project in option pricing.  International Journal of Computer Mathematics 92(12), 2361–2379.