text:
\section{Introduction}

In recent years, mixed Monte Carlo-partial differential equation (MC-PDE) methods for European options have seen an increase in research activity. In the context of stochastic volatility (SV) models with one-way coupling, these methods revolve around simulating the SV process, computing an expectation by solving a PDE conditional on the volatility path, and averaging over paths. The approach has been around for some time as in \cite{HW} and \cite{Lewis} but has seen renewed interest in \cite{Ang:2013zl}, \cite{LLP}, \cite{LP}, \cite{DJM}, \cite{DJS} and \cite{cozmareis}.

For high-dimensional European option pricing problems, under stochastic volatility, finite difference methods cannot be readily applied, and the correlations between the underlying processes often make the system non-affine which rules out Fourier-based quadrature methods. Also, full Monte-Carlo (MC) methods applied to such systems suffer from high variance and computational costs. An alternative is to find a middle ground between the two approaches where one simulates the underlying volatility processes and solves the resulting lower dimensional conditional PDEs, which may often be handled efficiently. This mixed method results in dimension reduction from the PDE perspective and variance reduction from the MC perspective.

As previous research which utilizes this strategy is focused on pricing European style options and addresses the dimension and variance reduction in computing relevant expected values, we analyse the mixed MC-PDE framework for Bermudan style options. The Bermudan context requires dealing with a high dimensional PDE between exercise dates, along with a high dimensional grid for accurately locating the exercise region; the latter being an issue that arises when moving from European to Bermudan options. When pricing Bermudan options the primary object of interest is the optimal stopping policy and the exercise boundaries that it defines. We note that one can always approximate the price of an American style option by considering a Bermudan option with high number of exercise dates as discussed in \cite{Bouchard}.

To deal with these issues, we develop a hybrid method which mixes the least squares Monte Carlo (LSMC) approaches of \cite{LS} and \cite{TV} with PDE techniques. The basic idea of our version of a mixed LSMC-PDE algorithm is to simulate paths of the underlying SV process, solve the PDE along each path, and regress these conditional expectations onto a family of basis functions over the volatility state-space. The algorithm can be viewed as an extension of \cite{TV} and reduces the monitoring of a high dimensional grid by replacing the volatility dimensions with a few regression coefficients. The approach provides variance reduction from the Monte-Carlo perspective, dimension reduction from the PDE perspective and alters the regression problems that are solved at each time step such that they are simpler than in standard LSMC. It is interesting to note that LSMC, when applied to SV problems, tends to be fairly inaccurate in determining the optimal exercise boundaries. In the literature, there have been a few direct modifications to LSMC, applied to SV problems, such as in \cite{Hilpisch}, \cite{Ludkovski2}, and \cite{ludkov} that address this issue. There have also been other types of probabilistic approaches such as in \cite{Ait} and \cite{Sircar}. The approaches in \cite{Hilpisch} and \cite{Ait} are specific to the Heston model and appear to be non-applicable to models outside the affine class. \cite{Sircar} develops a highly efficient method for multi-scale SV models. The work of \cite{Ludkovski2}, \cite{ludkov} cast LSMC as a classification problem using various experimental designs, more sophisticated regression methods, and apply their method to a one dimensional mean-reverting SV model as an example. The approach of \cite{SGBM} also looks promising, however, it typically requires a choice of basis functions for which one can compute (or approximate) expectations in closed form, which need to be developed on a case-by-case basis. It's also worth noting the work of \cite{particle}, which deals with the related problem of Bermudan option pricing under unobservable SV.

Our approach has its roots in \cite{LLP} where it is very briefly mentioned, but not analysed. Our contribution is a precise development of the algorithm, proof of convergence, discussion of complexity and complexity reduction methods, followed by a series of numerical examples. We begin by formulating the algorithm and its basic mechanics. Next, we develop theoretical aspects of the algorithm such as its probabilistic convergence and how it fits within the LSMC paradigm. With mild boundedness assumptions on our regression architecture and certain `separability conditions' on our model, we show that the algorithm converges almost surely. On a side note, we also develop a theoretical framework for the connection between the sampling measure which governs the distribution of the paths that we simulate and the pricing measure which induces the PDEs we solve. Such a discussion is missing from the current literature on mixed MC-PDE methods. We discuss two methods to improve the algorithm's complexity. The first is a clustering method that allows us to reduce the number of paths for which we solve the conditional PDEs. The second is a multi-level MC (MLMC)/multi-grid approach which allows us to use the PDEs solved on a coarse grid as a control variate. Finally, we apply the algorithm to the single and multi-asset Heston model, and demonstrate its merit in comparison to a full LSMC approach. The purpose of this paper is to develop the connection between LSMC and PDE methods, and not to showcase a state-of-the-art approach to either.

The remainder of this paper is organized as follows. In Section \ref{section:setup} we describe our notation and set up our optimal stopping problem. In Sections \ref{section:LSMC-PDEalg} and \ref{section:algpf} we describe the hybrid algorithm, show it is well-defined and converges probabilistically to our approximation scheme. Section \ref{section:pathred} discusses the clustering and multigrid methods which improve the complexity. In Sections \ref{section:Heston1} and \ref{section:Heston2} we apply the algorithm to the $2d$ and $4d$ Heston models and compare our numerical results to a full LSMC simulation and, in the $2d$ case, a finite difference solution.

\section{Problem Set Up and Notation}
\label{section:setup}

\subsection{Probability Spaces, SDEs, and Conditional Expectations}

\subsubsection{Initial Probability Space}

We suppose the existence of a probability space $(\Omega,{\mathbb F},{\mathbb Q})$ which may accomodate a $d_S + d_v$ dimensional stochastic process $X = (S^{(1)}_{t},...,S^{(d_s)}_{t},v^{(1)}_{t},...v^{(d_v)}_{t})_{t\in[0,T]}$ satisfying a system of SDEs with a strong, unique solution. We further suppose this system of SDEs exhibits one-way coupling in a sense which we describe below. We begin by defining mappings
\begin{align*}
\mu_S & : [0,T] \times \dsR^{d_S + d_v}  \to \dsR^{d_S}\,, \\
\mu_v & : [0,T] \times \dsR^{d_v}        \to \dsR^{d_v}\,,  \\
\sigma_S & : [0,T] \times \dsR^{d_S + d_v} \to \mathbb M_{d_S \times d_W}(\dsR)\,,  \\
\sigma_v & : [0,T] \times \dsR^{d_v} \to \mathbb M_{d_v \times d_{W, v} }(\dsR)\,,  \
\end{align*}
and a $d_W$-dimensional Brownian motion, $W=(W_t)_{t\in[0,T]}$, with independent components where $d_{W,v} < d_W$. We also denote the final $d_{W,v}$ components of $W_t$ as $W^v_t$, i.e., $W^{v,i}_t=W_t^{i-d_W+d_{W,v}}$ $\forall\;i\in\{d_W-d_{W,v}+1,\dots,d_W\}$. The process $X=(X_t)_{t\in[0,T]}$ is assumed to satisfy the following system of SDEs
\begin{align*}
dS_t &= \mu_S(t,S_t,v_t) dt + \sigma_S(t,S_t,v_t) \cdot dW_t\,, \qquad\text{and} \\
dv_t &= \mu_v(t,v_t) dt + \sigma_v(t,v_t)  \cdot dW^v_t\,. \
\end{align*}

While our algorithm may be formulated for such a general system of SDEs, our theoretical results and numerical examples assume there exist mappings $\widetilde{\mu}_S, \widetilde{\sigma}_S$ such that
\begin{align}
\mu^{(i)}_S(t,S,v)     & = S^{(i)} \ \widetilde{\mu}^{(i)}_S(t,v) \nonumber \\
\sigma^{(i)}_S(t,S,v)  & = S^{(i)} \ \widetilde{\sigma}^{(i)}_S(t,v) \
\label{eqn:musigma}
\end{align}
where the superscript $i$ indicates the $i^{th}$ entry of the vectors $S, \mu$ and $\tilde{\mu}$ and $i^{th}$ row of the matrices $\sigma_S, \widetilde{\sigma}_S$. That is, $X_t$ is a pure SV model without a local-volatility (LV) component. Examples of such models may be found in \cite{Stein}, \cite{Heston}, \cite{dHeston}, \cite{fastMR} and, \cite{4o2Heston} among others.

Finally, let ${\mathcal F}^{v}_{s,t}=\sigma( v_u)_{u\in[s,t]}$, ${\mathcal F}^S_{s,t}=\sigma( S_u)_{u\in[s,t]}$ and ${\mathcal F}^{W^v}_{s,t}=\sigma(W^v_u)_{u \in[s,t]}$ i.e., the natural filtrations generated by $v,S$ and $W^v$, respectively. To extend this notation, we sometimes write ${\mathcal F}^{Z}_t := {\mathcal F}^{Z}_{0,t}$ for some process $Z$.

Given $ t_n \in [0,T]$, we define a new class of (conditional) probability measures ${\mathbb Q}_{t_n,S}$ via ${\mathbb Q}_{t_n,S}(B) = {\mathbb Q}( B \ | \ S_{t_n} = S)$ for $B \in {\mathcal F}^S_{t_n,T} \vee {\mathcal F}^v_{t_n,T}$.

For a realization of $v_t$ and $W^v_t$ on $[t_n,t_{n+1}]$, which we denote as $[v]_{t_n}^{t_{n+1}}$, we assume there exists a finite-dimensional statistic of the path,
\begin{equation*}
\Lambda_n : C_0([t_n,t_{n+1}])^{d_v + d_{W^v}} \rightarrow \dsR^{d_{\Lambda}},
\end{equation*}
such that the following Markovian-like relation holds
\begin{align}
&{\mathbb E}^{{\mathbb Q}}[ \ h(S_{t_{n+1}}, v_{t_{n+1}}  ) \ | \ S_{t_n} = S \ , \  {\mathcal F}^{W^v}_{t_{n} ,t_{n+1} } \   ]  \nonumber \\
&\quad= {\mathbb E}^{{\mathbb Q}}[ \  h(S_{t_{n+1}}, v_{t_{n+1}}   ) \ | \ S_{t_n} = S \ , \ \Lambda_n([v]^{t_{n+1}}_{t_n}) , \ v_{t_{n+1}} \ ] \
\label{eqn:markov}
\end{align}
where $h : \dsR^{ d_S + d_v   } \to \dsR$ is Borel measurable such that the expectations are well-defined. On a high level, the conditional expectation generates a PDE over $[t_n,t_{n+1}] \times \dsR^{d_S}$ which depends on the simulated path, $[v]_{t_n}^{t_{n+1}}$, as follows
\begin{equation*}
\left\{
\begin{aligned}
\frac{\partial g}{\partial t}(t,S) + \mathcal{L}^{[v]}g(t,S) & = 0\,, \\
											    g(t_{n+1},S) & = f_{n+1}(S,v_{t_{n+1}})\,
\end{aligned}
\right.
\end{equation*}
for some function $f_{n+1}$, possibly different from $h_{n+1}$. The vector $\Lambda_n([v]_{t_n}^{t_{n+1}})$ captures the dependency of $g$ on $[v]_{t_n}^{t_{n+1}}$ over $[t_n,t_{n+1})$, as induced by $\mathcal{L}^{[v]}$, and $v_{t_{n+1} }$ corresponds to the PDE's boundary conditions. For example, in the Heston model, it can be shown that $\Lambda_n$ takes the following form
\begin{equation*}
\Lambda_n([v]_n^{n+1}) = \left( \ \int_{t_n}^{t_{n+1}} \sqrt{v_s} \ dW^{v}_s, \;\int_{t_n}^{t_{n+1}} v_s \ ds \  \right).
\end{equation*}
It is worth noting that generating the conditional PDE for (\ref{eqn:markov}) is, in general, non-trivial and to the best of our knowledge there are two approaches in the literature: the drift discretization method in \cite{LLP}, \cite{cozmareis}, \cite{DJM} and conditionally-affine decomposition of \cite{DJS}. Since our numerical examples assume that (\ref{eqn:musigma}) holds, conditional on a volatility path, our model is affine, and we employ the latter technique followed by Fourier Space Time-Stepping (FST) (see \cite{JJS}) in order to solve the induced affine conditional PDEs. Also, in Section \ref{section:pathred}, we introduce a multi-level Monte-Carlo scheme to further increase the algorithm's efficiency.

In our algorithm we will often encounter expectations of the form
\begin{equation*}
{\mathbb E}\left[ \, \left.\phi(v_{t_n} ) h(S_{t_{n+1}},v_{t_{n+1}}) \, \right| \ S_{t_n} = S \ , \ {\mathcal F}^{W^v}_{t_{n},t_{n+1} } \, \right]
\end{equation*}
which leads us to define
\begin{equation*}
\Theta_n([v]_n^{n+1}) = (v_{t_n}, \Lambda_n([v]_n^{n+1}), v_{t_{n+1}})
\end{equation*}
where $\phi: \dsR^{d_v} \to \dsR.$ In later sections we will replace $\Theta_n([v]_n^{n+1})$ with $[v]_{t_{n}}^{t_{n+1}}$ and simply write
\begin{equation*}
  {\mathbb E}^{{\mathbb Q}}\left[ \left. \  \phi(v_{t_n})  h( \ S_{t_{n+1}},   v_{t_{n+1}} \  ) \ \right| \ S_{t_n} = S \ , \  [v]^{{t_{n+1} }}_{t_n}   \ \right]:=
  {\mathbb E}^{{\mathbb Q}}\left[ \  \left. \phi(v_{t_n})  h(S_{t_{n+1}}, v_{t_{n+1}} ) \ \right| \ S_{t_n} = S \ , \Theta_n([v]_n^{n+1})  \ \right].
\end{equation*}
Equation (\ref{eqn:markov}) gives rise to the mappings $G_{f,t_n,S}$ defined by
\begin{equation}
\begin{split}
G_{f,t_n,S} & : \dsR^{d_{\Theta}} \to \dsR \,, \\
G_{f,t_n,S}( \theta ) & = {\mathbb E}^{{\mathbb Q}} \left.[ \ f(S_{t_{n+1} }, v_{t_{n+1}},v_{n} ) \ | \ S_{t_n} = S, \  \Theta_n = \theta \ \right.] \
\end{split}
\label{eqn:Gfn}
\end{equation}
and the conditional probability measures ${\mathbb Q}_{t_n,S,\theta}$ defined via
\begin{equation*}
{\mathbb Q}_{t_n,S,\theta}(B) : = {\mathbb Q}( \ B \mid S_{t_n} = S, \  \Theta_n = \theta  )  \
\end{equation*}
for $B \in \sigma( S_{t_{n+1}} ) \vee \sigma(v_{t_{n+1}}) \vee \sigma(v_{t_n})$. Letting $\widetilde{{\mathbb Q}}_{ \Theta_n}$ denote the distribution of $\Theta_n$ on $\dsR^{d_{\Theta}}$ we suppose our conditional expectation may be written as
\begin{equation*}
{\mathbb E}^{{\mathbb Q}}\left.[ f( S_{t_{n+1}}, v_{t_{n+1}}, v_{t_{n} }   ) \mid S_{t_n} = S, \Theta_n = \theta \right.] \  =   \int_{\Omega} f(S_{t_{n+1} }(\omega), v_{t_{n+1}}, v_{t_n}  ) \ d{\mathbb Q}_{t_n , S,\theta,n} (\omega)
\end{equation*}
so that we have the following relation
\begin{align}
& {\mathbb E}^{{\mathbb Q}}\left.[ \ f(S_{t_{n+1}},  v_{t_{n+1}} , v_{t_n}) \mid S_{t_n} = S \right.]  \nonumber \\
&\quad= \int_{\dsR^{d_{\Theta}} } \int_{\Omega} f(S_{t_{n+1} }(\omega), v_{t_{n+1}},v_{t_n}  ) \ d{\mathbb Q}_{t_n , S,\theta} (\omega ) \ d\widetilde{{\mathbb Q}}_{\Theta_n} (\theta) . \
\label{eqn:mrelation}
\end{align}

\subsubsection{Inherited Sampling Probability Space}
\label{sec:SamplingSpace}

A consequence of one-way coupling is the ability to simulate paths of $v_t$ independently of $S_t$. We now make sense of the notion of an iid collection of sample paths of $v_t$. Since we only realize $v_t$ through the statistics $\Theta_n$, we only describe how to generate iid copies of $\Theta_n$.

Denote the ordered subset $\{0 = t_0,...,t_n,t_{n+1},...,t_N = T \} \subset [0,T] $, and let $\{ [t_n,t_{n+1}] \}_{n=0}^{N-1}$ be the corresponding intervals. Given a path $v_t$ on $[0,T]$, define the $d_{\Theta} \times N$ dimensional matrix
\begin{equation*}
\Theta([v]) = \left.[ \ \Theta_0( [v]^{1}_0) \ , \ldots , \ \Theta_{N-1}( [v]^{N}_{N-1})  \right.].
\end{equation*}
This random matrix induces a measure $\widetilde{{\mathbb Q}}_{\Theta}$ on $\dsR^{Nd_{\Theta}}$. Given $\widetilde{{\mathbb Q}}_{\Theta}$, we introduce a new probability space $( \Omega', {\mathcal F}',{\mathbb Q}')$ equipped with a collection of independent random matrices
\begin{equation*}
\left.\{ \Theta([v^j] ) \right.\}_{j=1}^{\infty}\,,
\end{equation*}
such that each $\Theta([ v^j ]) $ has distribution $\widetilde{{\mathbb Q}}_{\Theta}$ on $\dsR^{Nd_{\Theta} }$. This construction follows from Kolmogorov's Extension Theorem applied to measures on $\dsR^{Nd_{\Theta} } $ (for a proof on $\dsR$ see \cite{Durrett}, for more general spaces see \cite{AB}). It then follows that each column $\Theta_n([v^j]_n^{n+1} )$ has distribution $\widetilde{{\mathbb Q}}_{\Theta_n}$ on $\dsR^{d_{\Theta} }$. Although the process $S_t$ is not defined on $\Omega'$, we can still compute relevant expectations involving this process using $G_{f,t_n,S}(\Theta)$ as defined in (\ref{eqn:Gfn}).

When taking limits in our algorithm, we consider expressions of the form
\begin{equation*}
\lim_{N\to\infty} \frac{1}{N}\sum_{j=1}^N {\mathbb E}^{{\mathbb Q}}\left.[ \ f(S_{n+1},v_{t_{n+1} },v_{t_n }) \ \mid \ S_{t_n} = S, [v^j]_n^{n+1}  \right]
\end{equation*}
which, by the SLLN, will converge to ${\mathbb E}^{{\mathbb Q}'}\left.[    G_{f,t_n,S}\left.( \Theta\left.([v]_{t_n}^{t_{n+1}} \right.)   \right.)    \right.]$ a.s. under ${\mathbb Q}'$. For our purposes, however, we require convergence to
\begin{equation*}
{\mathbb E}^{{\mathbb Q}}\left[ \ f(S_{t_{n+1}},v_{t_{n+1}},v_{t_n} ) \ \mid \ S_{t_n} = S \right].
\end{equation*}
To establish the equivalence between these expressions:
\begin{align}
{\mathbb E}^{{\mathbb Q}'}\left.[    G_{f,t_n,S}\left.( \Theta\left.([v]_{t_n}^{t_{n+1}} \right.)   \right.)    \right.] & = \int_{\Omega'} G_{f,t_n,S}(\Theta_n([v(\omega')]_n^{n+1}  )   ) \  d{\mathbb Q}'(\omega') \nonumber\\
&= \int_{\dsR^{d_{\Theta}}} G_{f,t_n,S}(\theta) \ d\widetilde{{\mathbb Q}}_{\Theta_n}(\theta) \nonumber \\
& = \int_{\dsR^{d_{\Theta}}} {\mathbb E}^{{\mathbb Q}}\left.[ f( S_{t_{n+1}}, v_{t_{n+1}}, v_{t_n}   ) \mid S_{t_n} = S, \Theta_n = \theta \right.] \ d \widetilde{{\mathbb Q}}_{\Theta_n}(\theta) \nonumber   \\
& = \int_{\dsR^{d_{\Theta}}}   \int_{\Omega} f(S_{t_{n+1} }(\omega), v_{t_{n+1}}, v_{t_n}  ) \ d{\mathbb Q}_{t_n , S,\theta,n} (\omega)  \ d \widetilde{{\mathbb Q}}_{\Theta_n}(\theta) \nonumber  \\
& = {\mathbb E}^{{\mathbb Q}}\left.[ f(S_{t_{n+1}}, v_{t_{n+1}} , v_{t_n}  ) \mid S_{t_n} = S \right.] \label{eqn:measequiv} \
\end{align}
where the second equality follows from $\Theta([v]_{t_n}^{t_{n+1}}) $ being $ \widetilde{{\mathbb Q}}_{\Theta_n}$ distributed.

\subsection{Notation, Bermudan Option Pricing, and the Dynamic Programming Principle}

Our goal is to price a Bermudan style option on $X$ whose payoff at each exercise date depends only on $S$. To proceed we need to introduce some additional notation.

Given a function $f:\mathds{R}_+\times \mathds{R}^n \mapsto \mathds{R}$  (i.e., $f(t,Z)$) we often denote the canonical stochastic processes corresponding to it as $(f_t)_{t\in[0,T]}$ with $f_t := f(t,Z_t)$, where $Z_t$ is some $n$-dimensional stochastic process. We also denote Euclidean norms of elements $x \in \dsR^n$ or $ x \in \mathds{M}_{n\times m}(\dsR)$ by $|x|$. Given any function $h: \dsR^n \to \dsR$, we write $||h||_{\infty} := \sup_{x \in \dsR^n} |h(x)|$ and $\text{supp} \ h  :=  \text{cl}\{ x \in \dsR^n \ | \   |h(x)| > 0  \}  $. Letting $X$ be an open subset of $\dsR^{n}$ we define $C_0(X)$ to be the set of continuous functions on $f: X \to \dsR$ that vanish at infinity. By vanish at infinity, we mean that for every ${\varepsilon} > 0$, the set $\{ x \in X \ | \ |f(x)| > {\varepsilon}\}$ is compact. We also let $C_c(X)$ be the set of compactly supported, continuous functions on $X$.

Let $\{t_0,...,t_M\} \subset [0,T]$ be an ordered set of exercise dates with $\Delta t_k = t_{k+1} - t_{k}$ and $h_{t_i}(S) : \dsR^{d_S} \to \dsR$ be our exercise function at each date. Later on, we place bounds on $h_{t_i}(S_{t_i})$ to ensure the algorithm converges. We often suppress the subscript $k$ in $\Delta t_k$ and $h_k$ to simplify notation when the context is clear. We also suppose the risk free rate is a constant, $r > 0$.

Valuing a Bermudan option requires developing an algorithm for evaluating
\begin{equation*}
V_t = \sup_{\tau \in T_t} {\mathbb E}^{{\mathbb Q}}\left[ \left.e^{-r  \tau}h_{\tau}(S_{\tau} )\right. \mid {\mathcal F}_t \right]\;,
\end{equation*}
where ${\mathcal F}_t = {\mathcal F}^{S}_t  \vee {\mathcal F}^{v}_t $ and $T_t$ is the set of $\mathcal F$-stopping times taking values in $\{t_k: k\in\{1,\dots,M\} \text{ and } t_k > t\}$. By the Markov property, $V_t$ depends only on $(t,S_t,v_t)$, and we can write $V_t=V(t,S_t,v_t)$ for some function $V:\mathds R_+\times \mathds R^{d_S}\times \dsR^{d_v}\mapsto \mathds R$.

At time $T = t_M$ we have $V(t_M,S_{t_M},v_{t_M}) = h_{t_M}(S_{t_M})$. We then define a new function, $C(t,S,v)$, denoted as the continuation value, at times $t_k$ by
\begin{equation*}
C(t_k,S,v) = e^{-r \Delta t_k} \;{\mathbb E}^{{\mathbb Q}}\left[ \left.V(t_{k+1},S_{t_{k+1}},v_{t_{k+1}}  )\right. \mid \left.S_{t_k} = S, v_{t_k} = v\right. \right].
\end{equation*}
By the discrete dynamic programming principle (DPP), for $k < M$ , we may express $V(t_k,S_{t_k},v_{t_k})$ as the maximum of the continuation value and the immediate exercise value at $t_k$:
\begin{equation*}
V_{t_{k}} = \max\left( h_{t_k}(S_{t_k}) \;, \;C(t_k,S_{t_k},v_{t_k}) \right)\;.
\end{equation*}
From this point on, for notational simplicity, we condense notation and replace our $t_k$ subscripts with simply $k$. We also replace arguments depending on time by subscripts.

\section{A Hybrid LSMC/PDE algorithm}
\label{section:LSMC-PDEalg}

We now describe a hybrid-method for computing $V_k(S,v)$ which is based on the \cite{TV} approach, but uses conditional PDEs to incorporate dimensional and variance reduction. We begin by giving an intuitive explanation and provide a formal, pseudo-code based, description in \ref{app:algstate}.

\subsection{Algorithm Overview}

We simulate $N$ paths of $v$ starting from the initial value $v_0$. Each path of $v_t$ over $[t_k,t_{k+1}]$ is represented as $[v]_{k}^{k+1}$. Given a product set $\mathcal{S} \subset \dsR^{d_S}$, we compute $V_k$ over the domain $\mathcal{S} \times \dsR^{d_v}$. The set ${\mathcal S}$ is the domain of the conditional expectations that we compute; in practice, it is the grid for our numerical PDE solver. We suppose the discretized form of ${\mathcal S}$ has $N_{s}$ points in each dimension so that there are $N^{d_S}_s$ points in total. Given the value $V_{k+1}(S_{k+1},v_{k+1})$ of the option at time $t_{k+1}$ we proceed to compute the continuation value at $t_k$. The algorithm begins at time $t_M = T$ where $V_M(S_M,v_M) = h_M(S_M)$.

\subsubsection{Solving along $S$ to obtain the pre-surface}\label{sec:step1}
For each simulation path $j\in\mathfrak{N}$ ($\mathfrak{N}=\{1,\dots,N\}$), we compute (beginning with $k=M-1$)
\begin{equation*}
C_k^j(s_i):=e^{-r\Delta t}\;{\mathbb E}^{{\mathbb Q}} \left[ \; \left.V_{k+1}(S_{k+1},v_{k+1}) \; \right.| \;S_{k} = s_i\,, \;[v^j]_{k}^{k+1}  \; \right]
\end{equation*}
for all $s_i \in \mathcal{S}$. These may be computed simultaneously over $\mathcal{S}$ for each path using a numerical PDE solver, or if the model, conditional on $v^j$, admits a semi-analytic  solution, then this may be computed for each $S \in \mathcal{S}$.

\begin{figure}
\centering
  \centering
  \includegraphics[width=\textwidth]{graphics/surface_precomp-eps-converted-to.pdf}
  \caption{The pre and completed continuation surface. The pre-surface is generated by solving a PDE along each variance path. It is smooth along the $S$ axis, and noisy across the $v$-axis. The completed surface is generated by regressing across the $v$-axis.}
  \label{fig:test1}
\end{figure}

\subsubsection{Regress across $v$ to obtain the completed surface}

For each $s_j \in \mathcal{S}$, from the previous step, we have $N$ realizations of the continuation value along each volatility path, i.e., $\{C_k^j(s_i)\}_{j\in\mathfrak{N}}$.

Next, apply least-squares regression to project this onto a family  $\{\phi_m(\cdot)\}_{m=1}^{d_B}$ of linearly independent basis functions over our volatility space. This results in a vector of coefficients $a(s_i)$ of length $d_B$, and provides the continuation value at $S_{t_{k}} = s_i$ for any point in the volatility space as follows:
\begin{equation*}
C_{k}(s_i,v) = \sum_{m=1}^{d_B} a_{m,k}(s_i)\phi_m(v) = a_{k}(s_i)\cdot \phi(v).
\end{equation*}

\subsubsection{Obtaining the Option Price}

The price of the option is then given by $V_k(s_i,v) = \max( h_{k}(s_i), C_k(s_i,v) )$. These steps are repeated from \ref{sec:step1} for all times $t_k$ where $k = M-1,...,1$.

\subsubsection{A Direct Estimate on the Time Zero Price}

Since at time zero, there is only a single value for $v_0$, we obtain an estimate for our time-zero prices by
\begin{equation*}
V_{d,0}(s_i,v_0) = \max\left(h_{t_0}(s_i) \;,\; \tfrac{1}{N}\sum_{j=1}^{N}e^{-r\Delta t}\;{\mathbb E}^{{\mathbb Q}}\left[ \left. V_1(S_{1},v_1) \; \right.| \;  S_{0} = s_i\;,\; [v^j]_{0}^{1} \;   \right] \;\right)\,
\end{equation*}
which is often biased high. Following \cite{SGBM}, we call this the direct estimator.

\subsubsection{A Lower Estimate on the Time Zero Price}

Given our estimated regression coefficients, we obtain a sub-optimal exercise policy $\tau(t,S,v)$ defined on $\{t_1,...,t_{M-1} \} \times {\mathcal S} \times \dsR^{d_v}$. Thus, we may define a lower estimate via the expectation
\begin{equation}
{\mathbb E}^{{\mathbb Q}} \left[ \; e^{-r \tau} h( S_{\tau} ) \mid  S_0, v_0 \; \right].
\label{eqn:lest}
\end{equation}
In traditional LSMC, one simulates a new independent set of paths $(S_t,v_t)$ to approximate (\ref{eqn:lest}). In the class of models we study, simulating both $S_t$ and $v_t$ undermines the variance reduction obtained by the algorithm and we instead use a hybrid approach.

To this end, we denote the $t_k$ holding and exercise regions by $\Gamma_k$ and $\Gamma^c_k$, respectively. We then simulate $N$ new independent paths of $v_t$ on $[0,T]$, compute
\begin{equation}
{\mathbb E}^{{\mathbb Q}}\left[ \; e^{-r \tau } h( S_{\tau})  \mid  S_0, [v^j]_0^T   \; \right]
\label{eqn:lestcond}
\end{equation}
via a PDE approach for $j\in\mathfrak{N}$, and take the average. To compute (\ref{eqn:lestcond}), for each $j\in \mathfrak{N}$, first set $V^j_M(S,v) = h(S)$. Next, compute
\begin{equation*}
U^j_{M-1}(S) = e^{-r\Delta t}\,  {\mathbb E}^{{\mathbb Q}} \left[   V^j_M(S_M,v_M) \mid S_{t_{M-1}} = S , [v^j]_{M-1}^M  \right]
\end{equation*}
via a PDE method for all $S \in {\mathcal S}$. The option price at time $t = t_{M-1}$ is then given by
\begin{equation*}
V^j_{M-1}(S) = U^{j}_{M-1}(S) \cdot I_{(S,v^j_{M-1}) \in \Gamma_{M-1}} +  h(S) \cdot I_{(S,v^j_{M-1}) \in \Gamma^{c}_{M-1}}
\end{equation*}
After repeating this procedure for times $k = M-2,...,1$ we obtain the lower estimate
\begin{equation}
V_{l,0}(S) = \frac{1}{N} \sum_{j=1}^{N} e^{-r \Delta t}\,  {\mathbb E}^{{\mathbb Q}} \left[ \;   V^j_1(S_1,v_1)  \mid S_{t_{1}} = S , [v^j]_{0}^1 \;   \right]
\label{eqn:lowestim}
\end{equation}
for all $S \in {\mathcal S}$.

\subsection{Discussion}

We refer the reader to \ref{app:algstate} for a pseudo-code based formal description.

Although we solve a PDE over thousands of paths of $v_t$ over each time interval $[t_n,t_{n+1}]$ and solve a linear-regression problem for each $s_i \in {\mathcal S}$ and $t_n$, the computational costs and run time are not as high as they may seem. First, the PDEs over each volatility path, and the regressions at $s_i \in {\mathcal S}$, are independent and can be parallelized. Also, based on (\ref{eqn:estcoef}), the regression problems at time $t_n$ require only one matrix inversion, and its result is applied to each of the $N_s^{d_S}$ regression sites. Finally, in Section \ref{section:pathred}, we discuss two  methods which allow us to reduce the algorithm's complexity.

We immediately see that our algorithm provides dimensional reduction from the PDE and regression perspective, and variance reduction from the MC perspective for a fixed simulation budget. If one employs a fully numerical scheme to solve the conditional PDEs, the algorithm is, in principle, capable of handling $3 + n$ dimensional problems where one solves a PDE over three dimensional asset space and simulates $n$ volatility variables. Although our set up is described in the context of an asset-volatility space setting, the algorithm can be applied to any situation where certain variables appear in the payoff and others appear in the background. For these settings, one should simulate the background variables and solve PDEs over the variables that appear in the payoff function.

As we shall see, the algorithm tends to be fairly accurate, for a given computational budget, in determining the time-zero value surface and optimal exercise regions. This may be attributed to the stability provided by our PDE grid, ${\mathcal S}$. When the conditional PDEs are solved along each path, we obtain our pre-surface as described in Subsection \ref{sec:step1}. At this point one has two choices: global or local regression. Our regression approach can be viewed as a special type of dimension reduced, local regression which is tailored to the presence of ${\mathcal S}$ and is equivalent to local regression onto $N_s^{d_S}$ carefully chosen regions. If $N_s = 128$ and $d_S = d_v = 2$, we are regressing onto $16,384$ families of basis functions at the cost of inverting a single matrix of size $d_B \times d_B$ where $d_B$ is about $10$, and a single matrix-multiplication for each $S \in {\mathcal S}$. Also, $C_n(s_j,v)$ is typically simple to fit as a function of $v$ and often does not require more than three or four monomials in each volatility dimension, resisting the well-known Basis Selection Problem.

Working with ${\mathcal S}$ has other advantages as well. In comparison to standard approaches to LSMC, there is a fundamental shift in how we compare the continuation value to the exercise value and locate the exercise boundary. At time $n$, when setting the value of $V_n$ for each $s_i \in {\mathcal S}$ we have
\begin{equation*}
V^N_n(s_i,v) = \max \left( h_n(s_i), C^{N}_n\left(s_i,v\right) \right)\,,
\label{eqn:compare}
\end{equation*}
and note that $h_n(s_i)$ is a deterministic constant as opposed to a function of a random variable. Thus, we have reduced the problem of locating the boundary from a global problem over the variables $S,v$ to a sequence of lower dimensional problems which are simpler in nature and exhibit less noise. Our approach essentially stores a function in $v$-space within each element of ${\mathcal S}$ and allows us to write our continuation surface using the separating representation
\begin{equation*}
C^N_n(s_i,v) = a^N_n(s_i) \cdot \phi(v).
\end{equation*}
To compute sensitivities with respect to $S$, one may simply differentiate $V^d_0(S,v_0)$ as in a standard PDE approach or the method of \cite{SGBM}. To get sensitivities with respect to $v$, one may employ the method of \cite{Caf}. Details and results may be found in \cite{Farahany}, where it is seen that sensitivities with respect to $S$ tend to exhibit much lower variance compared to the approach of  \cite{Caf} and comparable methods.

From the design of experiments perspective, our algorithm may be viewed as a batched design nested within a probabilistic design. The simulation of $v_t$ stemming from $v_0$ on $[0,T]$ and the repeated use of its paths, corresponds to a probabilistic design for the variable $v$. The solving of the conditional PDEs may be viewed as a sort of batched design for $S$ in the following sense. At time $t_n$, we solve a conditional PDE over $[v^j]_n^{n+1}$ and obtain the pre-regression continuation value of the option along this path for each $S \in {\mathcal S}$. This procedure is equivalent to selecting batches at each $S \in {\mathcal S}$, simulating an ``infinite'' number of paths of $S_t$ on $[t_n,t_{n+1}]$ conditional on $S_{t_n} = S$ and $[v^j]_n^{n+1}$ and equating the pre-regression continuation value for $[v^j]_n^{n+1}$ to the average payoff of these ``infinite" number of paths.

\section{Theoretical Aspects of the Algorithm}
\label{section:algpf}

In this section, we show that the LSMC-PDE algorithm, as described in Section \ref{section:LSMC-PDEalg}, is well defined and converges probabilistically. We remind the reader that for each $n$, the random variables $\{[v^j]_n^{n+1}\}_{j=1}^{\infty}$ are defined on the space $\Omega'$ and are iid (see the discussion in Section \ref{sec:SamplingSpace}). For notational convenience, we suppose the risk free rate is $0$. As our algorithm is based on the \cite{TV} approach to LSMC, many of our expressions are similar.

\subsection{Idealized Continuation Functions}

For each $n$, we consider a family of \textit{idealized continuation functions}, $C_{n}$, which are constructed by means of backwards induction. We begin by writing $C_M \equiv 0$ and $C_n(S,v) = a_n(S)\cdot \phi(v)$ for $n < M$ where $a_n(S)$ results from regressing the random variable
\begin{equation*}
 {\mathbb E}^{{\mathbb Q}}\left[ \left.\max(h_{n+1}(S_{n+1}),C_{n+1}(S_{n+1},v_{n+1}) )\right. \mid \left.S_{n} = S, \ v_n  \right.  \right] \ \ \text{onto the basis} \ \  \{\phi_m(\cdot)\}_{m=1}^{d_B}
\end{equation*}
for each $S \in {\mathcal S}$.
The coefficient vector $a_n(S)$ is the vector that minimizes the mapping $H_{n,S}: \dsR^{d_B} \to \dsR$ defined by
\begin{equation}
H_{n,S}(a) = {\mathbb E}^{{\mathbb Q}}\left[ \left( {\mathbb E}^{{\mathbb Q}}\left[ \left.f_{n+1}(S_{n+1},v_{n+1})\right. \mid \left. S_n = S, v_n\right. \right]  - a \cdot \phi(v_n)  \right)^2 \mid  \left.S_n = S\right. \right]
\label{eqn:functionH}
\end{equation}
where
\begin{equation}
\begin{split}
f_{n}(S,v) & = \max( h_n(S), C_{n}(S,v) ) \ \text{ for } \ n < M\;, \text{ and }\,  \\
f_{M}(S,v) & = h_M(S)\;. \
\end{split}
\label{eqn:functionf}
\end{equation}
To minimize $H$, we obtain the first order conditions and obtain the normal equations, resulting in the coefficients
\begin{equation*}
a_n(S) = A_n^{-1}  \cdot   {\mathbb E}^{{\mathbb Q}}\left[  \left.\phi(v_n) f_{n+1}(S_{n+1},v_{n+1})\right. \mid \left. S_{n} = S\right.   \right]
\end{equation*}
where $A_n = {\mathbb E}^{{\mathbb Q}}[ \phi(v_{t_n}) \phi(v_{t_n})^{\intercal} ] $.

\subsection{Almost-Idealized Continuation Functions}

Next, for a fixed $t_n$, we define a new type of continuation value, called the \textit{almost-idealized continuation functions}, $\widetilde{C}_n^{N}(S,v) = \widetilde{a}_n^N(S) \cdot \phi_n(v)$. These random variable are obtained by running the dynamic programming algorithm with the idealized continuation value at all times $k = M,...,n+1$. At time step $n$ we then estimate $\widetilde{a}^{N}_n(S)$ using our $N$ paths of $v_t$ and future idealized continuation values. This gives us the following regression coefficients for each $\omega \in \Omega'$
\begin{equation*}
\widetilde{a}^{N}_n(S,\omega') = \big[ A^{N}_n (\omega')\big]^{-1}   \frac{1}{N} \sum_{j=1}^{N}   \ \phi(v^j_n(\omega')) \cdot {\mathbb E}^{{\mathbb Q}}\left[\left.f_{n+1}(S_{n+1},v_{n+1}(\omega')) \right. \mid \left. S_{n} = S, [v^j(\omega')]_n^{n+1} \right. \right]
\end{equation*}
where $A^{N}_n = \frac{1}{N} \sum_{j=1}^{N} \phi(v_n(\omega'))\phi(v_n(\omega')  )^T $. Note that  $f_{n+1}$ involves the idealized continuation value at time $n+1$.

\subsection{Estimated Continuation Functions}

The \textit{estimated continuation functions} are the continuation functions produced from our algorithm: ${C}_n^{N}(S,v) = {a}_n^N(S) \cdot \phi_n(v)$. The regression coefficients are given by
\begin{equation}
a^{N}_n(S,\omega') =  \big[ A^{N}_n(\omega')\big]^{-1}   \cdot   \frac{1}{N} \sum_{j=1}^{N}   \  \phi(v^j_n(\omega')) \cdot  {\mathbb E}^{{\mathbb Q}}\left[\left. f_{n+1}^{N}(S_{n+1},v^j_{n+1}(\omega'),\omega' )\right.  \mid \left. S_{n} = S, [v^j(\omega')]_n^{n+1} \right. \right] \
\label{eqn:estcoef}
\end{equation}
where
\begin{align*}
f^{N}_{n}(S,v,\omega') & = \max( h_n(S), C^{N}_{n}(S,v,\omega') ) \ \text{ for } \ n < M\,, \text{ and } \\
f^{N}_{M}(S,v,\omega') & = h_M(S). \
\end{align*}

\subsection{Truncation Scheme}
We now state the following truncation scheme for our least-squares regression. It ensures that the coefficients produced by the algorithm are well defined and converge in a sense to be described later on.

\begin{assumption}[Truncation Conditions]
\
\label{assump:trunc}
\begin{enumerate}
\item The basis functions $\{\phi_m\}_{m=1}^{d_B}$ are bounded and supported on a compact rectangle.
\label{item:cond1}
\item The norm of the matrix
\label{item:cond2}
\begin{equation*}
\left[A^N_{n}(\omega')\right]^{-1} =  \left[  \frac{1}{N} \sum_{j=1}^{N} \phi(v^j_n(\omega')) \phi(v^j_n(\omega')  )^T    \right]^{-1}
\end{equation*}
is uniformly bounded for all $N$, $n$, provided the inverse is defined.
\item For each $i=1,...,M$, the exercise values $h_{t_i}(\cdot)$ are bounded with compact support in ${\mathcal S} \subset \dsR^{d_S}$.
\label{item:cond3}
\end{enumerate}
\end{assumption}

Condition (\ref{item:cond1}) may be imposed by limiting the support of $\{\phi_l\}$ on a bounded domain as they are typically smooth. By making $\text{supp} \ \phi \subset \dsR^{d_v}$ to be a very large rectangle, the value function is essentially unaffected.

Condition (\ref{item:cond2}) is imposed by replacing $[A^N_n]^{-1}$ with $[A^N_n]^{-1} I^{n,N}_{R}$ where $I^{n,N}_{R}$ is the indicator of the event that $[A^N_n]^{-1}$ is uniformly bounded by some constant $R$. If $R > |A^{-1}_n|$ then we have $[A^N_n]^{-1} I^{n,N}_{R} \to A^{-1}_n$ ${\mathbb Q}'$-a.s. Again, by making $R$ a very large constant, this has essentially no effect on the values obtained by the algorithm.

Condition (\ref{item:cond3}) on the functions $h$ are always satisfied in practice as numerically solving a PDE involves truncation of $h$'s domain.

\begin{lemma} Given the truncation conditions, the functions $H_{n,S}$ defined in (\ref{eqn:functionH}) are finite valued for all $n = 1,...,M-1$.
\label{Hlemma}
\end{lemma}
The proof is omitted due to its simplicity. The next lemma establishes a useful relationship between the idealized, almost idealized and estimated coefficients.
\begin{lemma}
Let $n \in \{1,...,M-2\}$, $S \in {\mathcal S}$. There exists a constant, $c$, which depends on our truncation conditions, such that
\begin{equation*}
|a_n^{N}(S) - a_n(S)| \leq c \cdot \beta^{N}_n(S) + \delta^{N}_n(S)
\end{equation*}
where
\begin{equation*}
\beta^{N}_n(S) = \frac{1}{N} \sum_{j=1}^{N}  {\mathbb E}^{{\mathbb Q}}\left[ \; \left. | a_{n+1}^{N}(S_{n+1}) - a_{n+1}(S_{n+1})    |    \right.  \mid  \left. S_{n} = S, [v^j]^{n+1}_n \right. \; \right] .
\end{equation*}
and
\begin{equation}
\delta^{N}_n(S) = | \widetilde{a}^{N}_n(S) - a_n(S)|
\label{eqn:deltadef}
\end{equation}
\label{lem:identity}
\end{lemma}

\begin{proof}
Given $n \in \{1,...,M-1\}$ and $S \in \mathcal{S}$ we have
\begin{equation*}
|a^{N}_n(S) - a_n(S)  | \leq | a^{N}_n(S) - \widetilde{a}^{N}_n(S)   |  + | \widetilde{a}_n^{N}(S)   - a_n(S) |.
\end{equation*}
After simplifying, we find
\begin{align}
|a^{N}_n(S) - \widetilde{a}^{N}_n(S)| &  \nonumber \\
\leq |[A^{N}_n]^{-1}_{R_A} |  \cdot & \frac{1}{N}  \sum_{j=1}^{N}   {\mathbb E}^{{\mathbb Q}}\left[ \left. | f_{n+1}^{N}(S_{n+1},v^j_{n+1}) - f_{n+1}(S_{n+1},v^j_{n+1}) |  \right. \mid \left. S_{n} = S, [v^j]_{n}^{n+1} \right. \right]   \cdot | \phi(v^j_n) | \nonumber \\
\leq c'  \cdot  \frac{1}{N}  &\sum_{j=1}^{N}   {\mathbb E}^{{\mathbb Q}}\left[ \left. | f_{n+1}^{N}(S_{n+1},v^j_{n+1}) - f_{n+1}(S_{n+1},v^j_{n+1}) |  \right. \mid \left. S_{n} = S, [v^j]_{n}^{n+1} \right. \right] \label{eqn:betasum} \
\end{align}
We then focus on the difference within the expectation. Using the inequality $|\max(a,b) - \max(a,c)| \leq |b - c|$ we find
\begin{align}
& {\mathbb E}^{{\mathbb Q}}\left[ \left. | f_{n+1}^{N}(S_{n+1},v^j_{n+1}) - f_{n+1}(S_{n+1},v^j_{n+1}) | \right. \mid \left. S_{n} = S, [v^j]^{n+1}_{n} \right.  \right]    \nonumber \\
& = {\mathbb E}^{{\mathbb Q}}\left[ \left.| \max( h_{n+1}(S_{n+1}), C_{n+1}^{N}(S_{n+1},v^j_{n+1}) ) - \max( h_{n+1}(S_{n+1}), C_{n+1}(S_{n+1},v^j_{n+1}) ) | \right.  \mid \left. S_{n} = S, [v^j]^{n+1}_{n} \right. \right] \nonumber    \\
&  \leq {\mathbb E}^{{\mathbb Q}}\left[ \left. | C^{N}_{n+1}(S_{n+1},v^j_{n+1})  - C_{n+1}(S_{n+1},v^j_{n+1})  |\right.  \mid \left. S_{n} = S, [v^j]_n^{n+1}\right. \right]    \nonumber \\
& \leq   {\mathbb E}^{{\mathbb Q}}\left[ \left. |  a^{N}_{n+1}(S_{n+1})- a_{n+1}(S_{n+1})| \right. \mid  \left. S_{n} = S, [v^j]_n^{n+1} \right. \right] \cdot  |\phi(v^j_{n+1}) | \nonumber \\
& \leq c'' \cdot  {\mathbb E}^{{\mathbb Q}}\left[ \left. |  a^{N}_{n+1}(S_{n+1})- a_{n+1}(S_{n+1})|    \right. \mid  \left. S_{n} = S, [v^j]^{n+1}_n\right. \right]  \label{eqn:finalterm}  \
\end{align}
where $c', c''$ depend on our truncation conditions. Substituting (\ref{eqn:finalterm}) into (\ref{eqn:betasum}), we obtain the result.
 \end{proof}
\subsection{Convergence Results for Models with Separability Properties}
In this section we prove that the coefficients also converge almost surely for a class of models with certain separability conditions. These separability properties are satisfied by pure SV models of the form (\ref{eqn:musigma}).

\begin{assumption}[Separability Conditions]
\
Let $n \in \{1,...,M-1\}$
	\begin{enumerate}
	\item The process $S_t$, for all $t\in[0,T]$, takes values in
	\[
	{\mathcal S} = \left\{ (S^{(1)},...,S^{(d_S)}) \in \dsR^{d_S} \ | \ S^{(i)} > 0, \forall i = 1,...,d_S     \right\}, \qquad \text{a.s.}
	\]
	\label{item:pdegrid}
	\item If $S_n = S$ almost surely, then $S_{n+1} = S \odot R_n  = (S_n^{(1)}R_n^{(1)},..., S_n^{(d_s)}R_n^{(d_S)} )$ where $R_n$ is adapted to $\vee_{t = t_n}^{t_M} {\mathcal F}_t$ and does not depend on the value of $S_n$. $R_n$ takes values in ${\mathcal S}$ a.s.
	\label{item:return}
	\item The exercise function, $h$, is continuous with compact support in ${\mathcal S}$. The basis functions $\phi$ are compactly supported and continuous on $\dsR^{d_v}$.
	\label{item:ctny}
\end{enumerate}
	\label{assump:separabilityconditions}
\end{assumption}

Condition (\ref{item:pdegrid}) limits our analysis to assets which take only positive values such as equities and foreign exchange rates.

Condition (\ref{item:return}) allows us to separate our future asset price as a product of its current price and return. The assumption that $R_n$ take values in ${\mathcal S}$ implies that they are finite valued a.s. As a result, letting ${\varepsilon}>0$, there exists tuples $\{(r^{(i)}_{l},r^{(i)}_{h})\}_{i=1}^{d_S}$, where $r^{(i)}_l > 0$ such that
\begin{equation}
\Omega_{r_l,r_h}^{{\varepsilon}} = \left\{ \  r^{(i)}_{l} \leq R^{(i)}_n \leq r^{(i)}_{h} \mid \forall \ n \in \{1,\ldots,M-1\} , \ \forall \ i \in \{1,\dots,d_S\} \  \right\} \subset \Omega
\label{eqn:posquad1}
\end{equation}
satisfies ${\mathbb Q}(\Omega_{r_l,r_h}^{{\varepsilon}}) > 1 - {\varepsilon}.$ We also write
\begin{equation*}
E^{{\varepsilon}}_{r_l,r_h} = \left \{  R \in \dsR^{d_S}  \mid  r^{(i)}_l \leq R^{(i)} \leq r^{(i)}_h, \forall \ i \in \left\{1,\ldots,d_S\right\} \   \right\}.
\label{eqn:posquad2}
\end{equation*}
Given $E_{r_l,r_h}^{{\varepsilon}}$, we may find an open set $U^{{\varepsilon}}$ such that $E_{r_l,r_h}^{{\varepsilon}} \subset U^{{\varepsilon}}$ and $ \overline{U^{{\varepsilon}}} \subsetneq {\mathcal S}$. By Urysohn's Lemma/Tietze's Extension theorem (see \cite{Munkres}), there exists a map $\eta_{E} : {\mathcal S} \to \dsR$ such that $\eta_{E} = 1$ on $E^{{\varepsilon}}_{r_l,r_h}$, $\eta_{E} = 0$ on $ {\mathcal S} \setminus U^{{\varepsilon}}$ and $||\eta_E||_{\infty} \leq 1$, i.e. a bump function supported on $E$. In most cases, our notation will suppress dependence on $r_l, r_h, {\varepsilon}$.

Condition (\ref{item:ctny}) allows us to apply the Stone-Weierstrass (SW) theorem which underlies the `separation technique' that will be demonstrated in upcoming Lemmas. We apply the version of SW for functions on unbounded domains that vanish at infinity (See \cite{Folland}). Suppose we are given the payoff function for a call option, i.e. $g : \dsR \to \dsR$ where $g(x) = (x-K)_+$. To modify $g$ such that it falls within our assumption, we first truncate its support to obtain a function $f(x) = (x-K)_+ I_{(0,R_1)}(x)$ where $R_1$ is large number. Finally, we continuously extend $f$ on $\dsR$ such that $f = 0$ on $(R_2,\infty)$ where $R_2 > R_1$. A similar construction may be done for a put option payoff near 0, and payoffs on higher dimensional domains.

Under these assumptions, we have the following main result.

\begin{theorem}
Let $n \in \left\{1,\ldots,M-1\right\}$ and $S \in {\mathcal S}$ be fixed. Then we have
\begin{equation*}
\lim_{{N}\to \infty} | a_n^{N}(S) - a_n(S) | = 0
\end{equation*}
${\mathbb Q}'$-almost surely
\label{thm:as}
\end{theorem}
We recall that given $\omega' \in \Omega'$, we generate random, independent paths $\left\{ \left[v^i(\omega')\right]_0^T \right\}_{i=1}^{\infty}$. Theorem \ref{thm:as} tells us that for almost every choice of sequence of paths, our coefficients based on these paths converge to the true idealized coefficient.

The proof borrows ideas from \cite{TV} and \cite{Clement}, and takes the following steps
\begin{enumerate}
\item \underline{Lemma \ref{lemma:geom}}. Carry out a geometric construction that allows us to approximately separate functions $h$ of $S_{t_n}$ that are continuous and compactly supported in ${\mathcal S}$. The function that provides the approximate separation is denoted as $\psi$.
\item \underline{Lemma \ref{lemma:finalapprox}}. Use the geometric construction to show the explicit relationship between $h$ and the separating functions, $\psi$ and thus demonstrate what is referred to as a separating estimate.
\item \underline{Lemma \ref{lemma:mm1}}. Prove the theorem for $n = M-1$ and also obtain an almost-sure separating estimate for $|a^N_{M-1}(S) - a_{M-1}(S)|$.
\item \underline{Lemma \ref{lemma:mm2}}. Prove the theorem for $n = M-2$ and also obtain an almost-sure separating estimate for $|a^N_{M-2}(S) - a_{M-2}(S)|$. The separating estimate for $n = M-2$ involves the function $\delta^N_{M-2}(S)$ as defined in (\ref{eqn:deltadef}).
\item \underline{Lemma \ref{lemma:deltaest}}. Develop an almost-sure separating estimate for $\delta^N_n$ for all $n \in \{1,\ldots,M-2\}$.
\item \underline{Proposition \ref{prop:main} }. Prove the theorem for $n = \{1,\ldots,M-3\}$ using Lemma \ref{lemma:mm2} and lemma \ref{lemma:deltaest}. Also obtain an almost-sure separating estimate for $|a^N_{n}(S) - a_{n}(S)|$ which is used during the induction.
\end{enumerate}

\begin{lemma}
Let ${\varepsilon} > 0$ and $h: {\mathcal S} \to \dsR$ be continuous and compactly supported in ${\mathcal S}$. Let $\widetilde{h}: {\mathcal S} \times {\mathcal S} \to \dsR$ be defined via
\begin{equation*}
\widetilde{h}(S,R ) = h(S \odot R) \cdot \eta_{E}(R)
\end{equation*}
where $E$ results from (\ref{eqn:posquad1}). There exists a map $\psi : {\mathcal S} \times {\mathcal S} \to \dsR$ of the form
\begin{equation*}
\psi(S,R) = \sum_{i=1}^k \psi_{i,S}(S)\psi_{i,R}(R)
\end{equation*}
such that $\psi_{i,S}, \psi_{i,R} \in C_c({\mathcal S})$  and $|| \widetilde{h} - \psi ||_{\infty} < {\varepsilon}$
\label{lemma:geom}
\end{lemma}

\begin{proof}

It follows from the properties of the mapping $i(S,R) = S \odot R$ and compact support of $\eta_E$ and $h$, that $\widetilde{h}$ is compactly supported in ${\mathcal S} \times {\mathcal S}$. We also have that $\widetilde{h}$ is continuous on ${\mathcal S} \times {\mathcal S}$.

To construct $\psi$, we begin by defining the algebra of functions
\begin{equation*}
{\mathcal A } := \left\{ \sum_{i=1}^k \psi_{i,S}(S) \psi_{i,R}(R) \ | \ \psi_{i,S}, \psi_{i,R} \in C_c({\mathcal S}), \ k \in {\mathbb N}    \right\}
\end{equation*}
and show ${\mathcal A }$ is dense in $C_0({\mathcal S}\times{\mathcal S})$ under the uniform metric.

Given two distinct points $\{(S_i,R_i)\}_{i=1}^2$, without loss of $  $generality, we assume $S_1 \neq S_2$. We now find bump functions $\eta_{S,1}, \eta_{S,2} \in C_c({\mathcal S})$ that separate $S_1, S_2$, i.e., $\eta_{S,1}(S_1) = 1, \eta_{S,1}(S_2) = 0$, $\eta_{S,2}(S_1) = 0, \eta_{S,2}(S_2) = 1$, along with a bump function $\eta_{R} \in C_c({\mathcal S})$ that is supported on $R_1$ and $R_2$. Letting $\psi_{1}(S,R) = \eta_{S,1}(S) \cdot \eta_{R}(R)$ and $\psi_2(S,R) = \eta_{S,2}(S) \cdot \eta_{R}(R)$, we have that ${\mathcal A }$ separates points. Since ${\mathcal A }$ contains bump functions supported at each point $(S,R) \in {\mathcal S}\times{\mathcal S}$, it vanishes nowhere.
\end{proof}

\begin{lemma}
Let ${\varepsilon}>0$ and $h : {\mathcal S} \to \dsR $ be continuous and supported in ${\mathcal S}$. Let $\widetilde{h}$ be as in the statement of Lemma \ref{lemma:geom} and $\psi(S,R)$ be a separable ${\varepsilon}$-approximation of $\widetilde{h}$. There exists a random variable $F(S)$ on $\Omega'$ and function $G(S)$ such that
\begin{align*}
{\mathbb E}^{{\mathbb Q}}\left[ \phi(v_n) h(S_{n+1}) \mid S_n = S, [v(\omega')]_n^{n+1} \right] &= F(S,\omega') +  \phi(v_n(\omega'))\sum_{i=1}^k \psi_{S,i}(S)  {\mathbb E}^{{\mathbb Q}}\left[ \psi_{R,i}(R_n) \mid [v(\omega')]_n^{n+1}  \right] \\
 \ {\mathbb E}^{{\mathbb Q}}\left[ \phi(v_n) h(S_{n+1})  \mid S_n = S \right] \  &= G(S) + \sum_{i=1}^k \psi_{S,i}(S) {\mathbb E}^{{\mathbb Q}}\left[ \  \phi(v_n) \psi_{R,i}(R_n)  \right]
\end{align*}
where $|F(S,\omega')| \leq c \ {\varepsilon}, \  |G(S)| \leq c' \ {\varepsilon}$ for all $S \in {\mathcal S}$, $\omega' \in \Omega'$ where $c,c'$ depend on our truncation and separability conditions.
\label{lemma:finalapprox}
\end{lemma}

\begin{proof}
We will only show the first inequality as they are similar.
\begin{align*}
&  \hspace*{-1em} {\mathbb E}^{{\mathbb Q}} \left[ \phi(v_n) h(S_{n+1}) \mid S_n = S, [v(\omega')]_n^{n+1} \right]  \\
&  =  \ {\mathbb E}^{{\mathbb Q}}[ \phi(v_n)  h(S \odot R_n) (1-\eta_E(R_n) )   \ | \ S_n = S, [v(\omega')]_n^{n+1} ] + {\mathbb E}^{{\mathbb Q}}\left[ \phi(v_n)  h(S \odot R_n) \eta_{E}(R_n)    \mid S_n = S, [v(\omega')]_n^{n+1} \right]    \\
 & =  \   F_1(S,\omega') +  {\mathbb E}^{{\mathbb Q}}\left[ \phi(v_n)  \widetilde{h}(S \odot R_n) \mid S_n = S, [v(\omega')]_n^{n+1} \right]    \\
& =  \   F_1(S,\omega') +  {\mathbb E}^{{\mathbb Q}}\left[ \phi(v_n)  (  \widetilde{h}(S \odot R_n) - \psi(S,R_n) )  \mid S_n = S, [v(\omega')]_n^{n+1} \right]  + {\mathbb E}^{{\mathbb Q}}\left[ \phi(v_n)  \psi(S,R_n) \mid S_n = S, [v(\omega')]_n^{n+1} \right]       \\
& =  \ F_1(S,\omega')     +  F_2(S,\omega') + {\mathbb E}^{{\mathbb Q}}\left[ \phi(v_n)  \psi(S,R_n) \mid S_n = S, [v(\omega')]_n^{n+1} \right]   \\
& =  \ F_1(S,\omega')     +  F_2(S,\omega') +  {\mathbb E}^{{\mathbb Q}}\left[ \phi(v_n)  \psi(S,R_n) \mid S_n = S, [v(\omega')]_n^{n+1} \right]   \\
& =  \ F_1(S,\omega')     +  F_2(S,\omega') +  \phi(v_n(\omega'))\sum_{i=1}^k \psi_{S,i}(S)  {\mathbb E}^{{\mathbb Q}}\left[ \psi_{R,i}(R_n) \mid [v(\omega')]_n^{n+1}  \right] \
\end{align*}
It then follows that $|F_i(S,\omega')| < {\varepsilon}$ for all $S \in {\mathcal S}$,   $\omega' \in \Omega'$, $i = 1,2$. Finally, we set $F(S,\omega') = F_1(S,\omega') + F_2(S,\omega')$.
\end{proof}

\begin{lemma}
Let ${\varepsilon}>0$. There exists a set $\Omega'_{M-1} \subset \Omega'$ with ${\mathbb Q}'(\Omega'_{M-1}) = 1$ such that for $\omega' \in \Omega'_{M-1}$ we have the following
\label{lemma:sepest}
\begin{equation*}
|a^{N}_{M-1}(S,\omega') - a_{M-1}(S,\omega')| \leq c_{M-1} \  {\varepsilon} +  \sum_{i_{M-1} = 1}^{k_{M-1} } \alpha^{i_{M-1},N}_{M-1,M}(\omega') \cdot \psi_{S,i_{M-1}}(S) \\
\label{eqn:mclaim}
\end{equation*}
where $\alpha^{i_{M-1},N}_{M-1,M}$ are random variables that depend on $\{ [v^j]^{M}_{M-1} \}_{j=1}^{\infty}$ and
\begin{equation*}
\lim_{N\to \infty} \alpha^{i_{M-1},N}_{M-1,M}( \omega') = 0,
\end{equation*}
 $\psi_{S,i_{M-1}}$ satisfy the conditions of Lemma \ref{lemma:finalapprox} and $c_{M-1}$ depends on our truncation and separability conditions.
 \label{lemma:mm1}
 \end{lemma}

\begin{proof}

\

Given ${\varepsilon} > 0$ we have
\begin{align*}
|a^{N}_{M-1}(S) - a_{M-1}(S)| &=| \ [A^{N}_{M-1}]^{-1} \frac{1}{N} \sum_{j=1}^{N} \phi(v^j_{M-1}){\mathbb E}^{{\mathbb Q}}[ h(S_M) \ | \ [v^j]_{M-1}^M, S_{M-1}=S  ] \\
& \ \ \ \ - A_{M-1}^{-1}  {\mathbb E}^{{\mathbb Q}}[ h(S_M) \  \phi(v_{M-1}) \ | \ S_{M-1}=S  ] \  | .\
\end{align*}
By Lemma \ref{lemma:finalapprox} we can find a constant $c$ and a separable ${\varepsilon}$-approximation, $ \psi_{M-1}(S,R)$ such that
\begin{align*}
|a^{N}_{M-1}(S) - a_{M-1}(S)| & \leq | \  \left[A^{N}_{M-1}\right]^{-1} \frac{1}{N} \sum_{j=1}^{N} \phi(v^j_{M-1}) \sum_{i_{M-1} = 1}^{k_M} \psi_{S,{i_{M-1}  }}(S) {\mathbb E}^{{\mathbb Q}}\left[  \psi_{R,i_{M-1}  }(R_{M-1})  \mid [v^j]_{M-1}^M \ \right]  \\
& \ \ \ - A_{M-1}^{-1} \sum_{i_{M -1}=1}^{k_M} \psi_{S,{i_{M-1}}} (S){\mathbb E}^{{\mathbb Q}}\left[ \phi(v_{M-1}) \psi_{R,i_{M-1}  }(R_{M-1}) \ \right] \ | + c \ {\varepsilon} \\
& \leq  \sum_{i_{M -1 } =1}^{k_{M-1}  } |\psi_{S,i_{M -1} }(S)| \cdot \alpha^{i_{M-1},N}_{M-1,M} +  c \ {\varepsilon} \
\end{align*}
where the last line follows from interchanging the summations, $c$ depends on our truncation and separability conditions and
\begin{align*}
\alpha^{i_{M -1} ,N}_{M-1,M} &= | \  [A^{N}_{M-1}]^{-1} \frac{1}{N} \sum_{j=1}^{N} \phi(v^j_{M-1}){\mathbb E}^{{\mathbb Q}}\left[ \psi_{R,i_{M-1}  }(R_{M-1})  \mid [v^j]_{M-1}^M \ \right] \\
& \ \ \ \   - A_{M-1}^{-1} {\mathbb E}^{{\mathbb Q}}\left[  \phi(v_{M-1}) \psi_{R,i_{M-1}  }(R_{M-1}) \  \right] \ | .
\label{eqn:alphamm1}
\end{align*}
By the SLLN, there exists a set $\Omega'_{M-1}$ with ${\mathbb Q}'( \Omega'_{M-1} ) = 1$ such that $\lim_{N\to\infty} \alpha^{i_{M-1},N}_{M-1,M} = 0 $ on $\Omega'_{M-1}$ for all $i_{M-1} \in \{1,...,k_M\}$.
\end{proof}

\begin{lemma}
Let ${\varepsilon} > 0$. There exists a set $\Omega'_{M-2}$ such that ${\mathbb Q}'(\Omega'_{M-2}) = 1$ and that for $\omega' \in \Omega'_{M-2}$
\begin{align*}
|a^{N}_{M-2}(S,\omega') - a_{M-2}(S,\omega')| & \leq c^N_{M-2}(\omega') \ {\varepsilon} +  \sum_{i_{M-1} = 1}^{k_{M-1} } \alpha^{i_{M-1},N}_{M-1,M}(\omega')\sum_{i_{M-2}}^{k_{M-2}} \alpha^{i_{M-2},i_{M-1},N}_{M-2,M-1}(\omega')  \psi_{S,i_{M-2},i_{M-1}}(S) \\
& + \delta^N_{M-2}(S,\omega')\
\end{align*}
where $ \alpha^{i_{M-1},N}_{M-2,M-1}$ are random variables that depend on $\{ [v^j]_{M-2}^{M-1}\}_{j=1}^{\infty}$ and satisfy
\begin{equation*}
\limsup_{N\to\infty} \alpha^{i_{M-2},N}_{M-2,M-1}(\omega') < \infty.
\end{equation*}
Also, $\limsup_{N\to\infty} c^N_{M-2}(\omega') < \infty $ with a bound depending only on our truncation and separability conditions and $\psi_{S,i_{M-1},i_{M-2}}$ satisfy the conditions of Lemma \ref{lemma:finalapprox}. Lastly, $\delta^N_{M-2}$ is as defined in (\ref{eqn:deltadef}).
\label{lemma:mm2}
\end{lemma}

\begin{proof}
Given ${\varepsilon} > 0$, and using Lemma \ref{lem:identity} and \ref{lemma:mm1} we have that for $\omega' \in \Omega'_{M-1} $
\begin{align*}
|a^{N}_{M-2}(S,\omega') - a_{M-2}(S,\omega') | &\leq \widetilde{c} \frac{1}{N} \sum_{j=1}^{N}  {\mathbb E}^{{\mathbb Q}}\left[  \left. |a^{N}_{M-1}(S_{M-1},\omega') - a_{M-1}(S_{M-1},\omega') |   \right.  \mid  \left. S_{M-2} = S, [v^j(\omega')]^{M-1}_{M-2} \right. \right]   \\
&+  \delta^{N}_{M-2}(S,\omega') \\
& \leq c' {\varepsilon} + \sum_{i_{M-1}}^{k_{M-1}} \alpha^{i_{M-1},N}_{M-1,M}(\omega') \frac{1}{N}\sum_{j=1}^N {\mathbb E}^{{\mathbb Q}}\left[ \psi_{S,i_{M-1}}(S_{M-1}) \mid S_{M-2} = S  \ [v^j(\omega')]_{M-2}^{M-1}   \right] \\
& + \delta^{N}_{M-2}(S,\omega') \
\end{align*}
For each $\psi_{S,i_{M-1}}$, we apply Lemma \ref{lemma:finalapprox}, and obtain a separable ${\varepsilon}$-approximating function $\psi_{i_{M-1}}(S,R)$ of the form
\begin{equation*}
\psi_{i_{M-1}}(S,R) = \sum_{i_{M-2}=1}^{k_{M-2}} \psi_{S,i_{M-2},i_{M-1}}(S)\psi_{R,i_{M-2},i_{M-1}}(R)
\end{equation*}
and apply it to the expectation. This results in
\begin{align*}
|a^N_{M-2}(S,\omega') - a_{M-2}(S,\omega') | &\leq (c' + \widetilde{c}^{N}_{M-2}(\omega')) {\varepsilon} + \sum_{i_{M-1} = 1}^{k_{M-1} } \alpha^{i_{M-1},N}_{M-1,M}(\omega')\sum_{i_{M-2}=1}^{k_{M-2}} \alpha^{i_{M-2},i_{M-1},N}_{M-2,M-1}(\omega')  \psi_{i_{M-2},i_{M-1}}(S) \\
& + \delta^N_{M-2}(S,\omega') \
\end{align*}
where $\lim_{N\to \infty} \widetilde{c}^N_{M-2}(\omega') = 0$ and
\begin{equation}
\alpha^{i_{M-2},N}_{M-2,M-1} = \frac{1}{N} \sum_{j=1}^N {\mathbb E}^{{\mathbb Q}}\left[ \;  \psi_{R,i_{M-2},i_{M-1}}(R_{M-2}) \mid [v^j]_{M-2}^{M-1} \; \right] .
\label{eqn:alphamm2}
\end{equation}
By the SLLN, for each tuple $(i_{M-1},i_{M-2})$ there exists a set of full measure such that $\alpha^{i_{M-2},N}_{M-2,M-1}$ converges. Thus, we may find a set $\widetilde{\Omega'}_{M-2}$ with ${\mathbb Q}'(\widetilde{\Omega'}_{M-2}) = 1$ on which they all converge. Also, $\delta^N_{M-2}(S)$ converges to 0 on a set of full measure $\Omega_{M-2}^{'\delta}$. We then define $\Omega'_{M-2} = \widetilde{ \Omega'}_{M-2} \cap \Omega'_{M-1} \cap \Omega^{'\delta}_{M-2} $ and $c^N_{M-2} = c' + \widetilde{c}^N_{M-2}$ which implies the lemma.
\end{proof}

\begin{lemma}
Let $n \in \{1,...,M-2 \}$, $ \delta^{N}_n(S)$ as in (\ref{eqn:deltadef})  and ${\varepsilon} > 0$. There exists a set $\Omega^{'\delta}_n \subset \Omega'$ with ${\mathbb Q}'(\Omega_n^{'\delta}) = 1 $ such that for $\omega' \in \Omega_n^{'\delta}$ we have
\begin{equation*}
\delta^N_{n}(S,\omega') \leq c \cdot {\varepsilon} + \sum_{i_n=1}^{k_n} \beta^{i_n,N}_{n,n+1}(\omega') \psi_{i_n}(S)
\end{equation*}
where $\beta^{i_n,N}_{n,n+1}$ are random variables that depend on $\{ [v^j]_n^{n+1} \}_{j=1}^{\infty}$, $\lim_{N\to\infty} \beta^{i_n,N}_{n,n+1}(\omega') = 0$, $c$ depends on our truncation and separability conditions, and $\psi_{i_n}$ satisfy the conditions of Lemma \ref{lemma:finalapprox}.
\label{lemma:deltaest}
\end{lemma}

\begin{proof}
We break the proof up into multiple stages.

\textbf{Preliminary estimates on} $a_k, f_k$.

Let $f_k$ be as defined in (\ref{eqn:functionf}). We begin by showing that for each $j \in \{n,...,M-1\}$ that
\begin{align*}
a_j(S) &= F_j(S) + G_j(S) \\
f_j(S,v) &= H_j(S,v) + J_j(S,v)  \
\end{align*}
where $|F_k(S)| \leq c_{k} {\varepsilon}$, $|H_k(S,v)| \leq c'_{k} {\varepsilon}$ with $c_{k},c'_k$ depending on our truncation conditions and separability conditions. The functions $G_j$ admit the representation
\begin{equation*}
G_j(S) = \sum_{i_j =1}^{k_j} c_{i_j, j} \psi_{S,i_j}(S)
\end{equation*}
where $c_{i_j,j} \in \dsR^{d_B} $ and $\psi_{S,i_j} \in C_c({\mathcal S})$. Finally $J_j(S,v) = \max( h(S),G_{j}(S) \cdot \phi(v))$.
To this end, we let $j = M-1$ and have that
\begin{align}
a_{M-1}(S) =& A^{-1}_{M-1} {\mathbb E}^{{\mathbb Q}}\left[ \phi(v_{M-1}) h(S_M)  \mid S_{M-1} = S  \right] \nonumber \\
       =& A^{-1}_{M-1} {\mathbb E}^{{\mathbb Q}}[ \phi(v_{M-1}) h(S \odot R_{M-1})  (1-\eta_{E}(R_{M-1}))  \   ] \nonumber \\
       & + A^{-1}_{M-1} {\mathbb E}^{{\mathbb Q}}\left[ \phi(v_{M-1}) h(S \odot R_{M-1})  \eta_{E}(R_{M-1})  \right]  \
\label{eqn:am1}
\end{align}
where $E$ is as in (\ref{eqn:posquad2}) and focus on the second term. We let
\begin{align*}
\widetilde{\Xi}_{M-1} & : \dsR^{2d_S + d_v } \to \dsR^{d_B} \\
\widetilde{\Xi}_{M-1}(S,R,v_1) & =  \phi(v_1)  h(S \odot R) \eta_E(R) \\
                                   & =  \phi(v_1) \widetilde{h}(S ,R) \
\end{align*}
and apply Lemma \ref{lemma:geom} to find an ${\varepsilon}$-separating function of the form
\begin{equation*}
\psi_{M-1}(S,R) = \sum_{i_{M-1} = 1}^{k_{M-1}} \psi_{S,i_{M-1}}(S)\psi_{R,i_{M-1}}(R)
\end{equation*}
where $\psi_{S,i_{M-1}}, \psi_{R,i_{M-1}} \in C^0_c({\mathcal S})$ and $||\widetilde{h}-\psi_{M-1}||_{\infty} < {\varepsilon}$.
This allows us to write
\begin{equation*}
\widetilde{\Xi}_{M-1}(S,R,v_1) =  F_{2,M-1}(S,R,v_1) + \phi(v_1)\psi_{M-1}(S,R)
\end{equation*}
where
$ F_{2,M-1}(S,R,v_1) = \phi(v_1)(\widetilde{h}(S,R) - \psi_{M-1}(S,R) )$ and $
{\mathbb E}^{{\mathbb Q}}[ F_{2,M-1}(S,R_{M-1},v_{M-1}) ]  < c \ {\varepsilon}$ with $c$ depending on our truncation and separability conditions.

Returning to our expression in $(\ref{eqn:am1})$, we find
\begin{align*}
a_{M-1}(S) &= F_{M-1}(S) + \sum_{i_{M-1} =1}^{k_{M-1}} c_{i_{M-1},M-1} \psi_{S,i_{M-1} }(S) \\
           &=: F_{M-1}(S) + G_{M-1}(S) \
\end{align*}
where $|F_{M-1}(S)| \leq c_{M-1} \  {\varepsilon}$ for all $S \in {\mathcal S}$ , $c_{M-1}$ depends on our truncation and separability conditions, and
\begin{equation*}
c_{i_{M-1},M-1} = {\mathbb E}^{{\mathbb Q}}\left[  \phi(v_{M-1}) \psi_{R,i_{M-1}}(R_{M-1})   \right].
\end{equation*}
We now turn to $f_{M-1}(S,v)$
\begin{align*}
f_{M-1}(S,v) &= \max( h(S) , a_{M-1}(S)\cdot \phi(v) ) \\
&= \max(h(S)  , F_{M-1}(S)\cdot \phi(v) + G_{M-1}(S) \cdot \phi(v) )\\
&= \Big[\max(h(S), F_{M-1}(S)\cdot \phi(v) + G_{M-1}(S) \cdot \phi(v) ) - \max(h(S), G_{M-1}(S)\cdot \phi(v) ) \Big] \\
&+ \max(h(S), G_{M-1}(S) \cdot \phi(v) ). \
\end{align*}
Equating $H_{M-1}(S,v)$ to the term in brackets and
\begin{equation*}
J_{M-1}(S,v) := \max( h(S),G_{M-1}(S)\cdot \phi(v)),
\end{equation*}
gives us the required form and concludes the claim for $j = M-1$.

Now let $j \in \{1,...,M-2\}$, we have
\begin{align*}
a_j(S) & = A^{-1}_{j} {\mathbb E}^{{\mathbb Q}}\left[ \phi(v_j)f_{j+1}(S_{j+1},v_{j+1}) \right] \\
       & = A^{-1}_j {\mathbb E}^{{\mathbb Q}}\left[ \phi(v_j)H_{j+1}(S_{j+1},v_{j+1}) \right] \\
       & + A^{-1}_j {\mathbb E}^{{\mathbb Q}}\left[ \phi(v_j) \max(h(S\odot R_{j} ), G_{j+1}(S\odot R_j) \cdot \phi(v_{j+1})   ) (1 - \eta_{E_{j+1}}(R_n) )  )  \right]  \\
       & + A^{-1}_j {\mathbb E}^{{\mathbb Q}}\left[ \phi(v_j) \max(h(S\odot R_j), G_{j+1}(S\odot R_j) \cdot \phi(v_{j+1})   )  \eta_{E_{j+1} }(R_n)     \right],  \
\end{align*}
where $E_{j+1}$ is as in (\ref{eqn:posquad1}) so that ${\mathbb Q}(\Omega^{{\varepsilon}_{j+1},c}) < \frac{{\varepsilon}}{||J_{j+1}||_{\infty} }$, and focus on the final term. By assumption, $G_{j+1}(S) = \sum_{i_{j+1}=1}^{k_{j+1}} c_{i_{j+1} } \psi_{S,i_{j+1}  }(S)$ where $c_{i_{j+1} }  \in \dsR^{d_B}$. Letting $d_{i_{j+1}  }(v) = c_{i_{j+1} }\cdot\phi(v)$, we are led to consider the function
\begin{align*}
\widetilde{\Xi}_{j} &: \dsR^{2d_S + 2d_v} \to \dsR^{d_B} \\
\widetilde{\Xi}_{j}(S,R,v_1,v_2) &= \phi(v_1) \max( h(S\odot R)\eta_{E_{j+1}}(R), \sum_{i_{j+1}=1}^{k_{j+1}} d_{i_{j+1}}(v_2) \psi_{S,i_{j+1}}(S\odot R)\eta_{E_{j+1}}(R)  ) \\
&= \phi(v_1) \max( \widetilde{h}(S,R)  , \sum_{i_{j+1}=1}^{k_{j+1}} d_{i_{j+1}}(v_2) \widetilde{\psi}_{S,i_{j+1}}(S,R)  ) \
\end{align*}
where $\widetilde{\psi}_{S,i_{j+1}}(S,R) = \psi_{S,i_{j+1}}(S,R)\eta_{E_{j+1}}(R)$.
Before carrying out a construction very similar to the $j = M-1$ case, we let
\begin{equation*}
A_{j+1}  =   ||\widetilde{h}||_{\infty} + \sum_{i_{j+1} = 1}^{k_{j+1}} ||d_{i_{j+1}}||_{\infty} ||\widetilde{\psi}_{S,i_{j+1}} ||_{\infty} .
\end{equation*}

Next, we note that
\begin{equation*}
\max\left( \widetilde{h}   , \sum_{i_{j+1}=1}^{k_{j+1}} d_{i_{j+1}} \widetilde{\psi}_{S,i_{j+1}}  \right) \in C_c\left({\mathcal S} \times \left( {\mathcal S} \times \dsR^{d_v} \right) \right)
\end{equation*}
and find a function of the form
\begin{equation*}
\psi_j(S,R,v) = \sum_{i_{j}=1}^{k_{j}  } \psi_{S,i_j}(S)\psi_{R,v,i_j}(R,v)
\end{equation*}
such that $\psi_{S,i_j} \in C_c({\mathcal S})$, $\psi_{R,v,i_j} \in C_c({\mathcal S} \times \dsR^{d_v})$ and
\begin{equation*}
\left|\left|  \max \left( \widetilde{h} , \sum_{i_{j+1} =1}^{k_{j+1} } d_{i_{j+1}} \widetilde{\psi}_{S,i_{j+1} }  \right) - \psi_j \right|\right|_{\infty} < {\varepsilon}.
\end{equation*}
Finally, we write
\begin{equation*}
a_j(S) = F_{j,1}(S) + F_{j,2}(S) + F_{j,3}(S) + \sum_{i_j=1}^{k_j} c_{i_j,j} \psi_{S,i_j}(S)
\end{equation*}
where
\begin{align*}
F_{j,1}(S) &= A^{-1}_j {\mathbb E}^{{\mathbb Q}}\left[ \phi(v_j)H_{j+1}(S_{j+1},v_{j+1}) \ | \ S_j = S \right],\\
F_{j,2}(S) &= A^{-1}_j {\mathbb E}^{{\mathbb Q}}\left[ \phi(v_j) \max(h(S\odot R_{j} ), G_{j+1}(S\odot R_j) \cdot \phi(v_{j+1})   ) (1 - \eta_{E_{j+1}}(R_n) )  )  \right],  \\
F_{j,3}(S) &= {\mathbb E}^{{\mathbb Q}}\left[ \ \widetilde{\Xi}_j(S,R_j,v_j,v_{j+1}) - \phi(v_j) \psi_j(S,R_j,v_{j+1})   \ \right]\
\end{align*}
so that for  $i=1,2,3$, we have $||F_{j,i}||_{\infty} \leq c'_{j,i} {\varepsilon}  $ where $c'_{j,i}$ depend on our truncation and separability conditions and
\begin{equation*}
c_{i_j,j}  = {\mathbb E}\left[  \phi(v_j)\psi_{R,v,i_j}(R_j,v_{j+1} )  \right] .
\end{equation*}
Also we have that $\psi_{S,i_j }(S) \in C_c({\mathcal S}) $ . Writing $F_{j} = F_{j,1} + F_{j,2} + F_{j,3}$ and $G_j(S)$ to be the final term completes the claim for $a_j(S)$. Showing the result for $f_{j}(S,v)$ is analogous to the base case and so we omit the proof.

\textbf{Estimates of $\delta^N_n(S)$}

We write
\begin{align*}
\delta^N_n(S) &= | \ \left[A^N_n\right]^{-1} \frac{1}{N}\sum_{j=1}^N {\mathbb E}^{{\mathbb Q}}\left[ \phi(v_n) f_{n+1}(S_{n+1},v_{n+1})   \mid S_n = S, [v^j]_n^{n+1} \right] \\
& - A_n^{-1}{\mathbb E}^{{\mathbb Q}}\left[ \phi(v_n) f_{n+1}(S_{n+1},v_{n+1})   \mid S_n = S \right] \ | \\
& \leq | \ \left[A^N_n\right]^{-1} \frac{1}{N}\sum_{j=1}^N {\mathbb E}^{{\mathbb Q}}\left[ \phi(v_n) J_{n+1}(S_{n+1},v_{n+1})  \mid S_n = S, [v^j]_n^{n+1} \right] \\
&- A_n^{-1}{\mathbb E}\left[ \phi(v_n) J_{n+1}(S_{n+1},v_{n+1})  \mid S_n = S \right] \
+ c \ {\varepsilon} \\
& \leq | \ \left[A^N_n\right]^{-1} \frac{1}{N}\sum_{j=1}^N \phi(v^j_n){\mathbb E}^{{\mathbb Q}}\left[  \max( h(S_{n+1}) ,G_{n}(S_{n+1}) \cdot \phi(v_{n+1})) \eta_E(R_n) \mid S_n = S, [v^j]_n^{n+1} \right] \\
&- A_n^{-1}{\mathbb E}^{{\mathbb Q}} \left[ \phi(v_n)  \max( h(S_{n+1}),G_{n}(S_{n+1}) \cdot \phi(v_{n+1})) \eta_E(R_n)   \mid S_n = S \right] \ |  \
+ c' \ {\varepsilon} \
\end{align*}
where $c$ and $c'$ depend on our truncation conditions. We now focus on the expression within our expectations and define the function
\begin{align*}
\widetilde{\Xi}'_{n} &: \dsR^{2d_S  + 2d_v}  \to \dsR^{d_B} \\
\widetilde{\Xi}'_n(S,R,v_1,v_2) & = \phi(v_1) \max( h(S \odot R),G_{n}(S \odot R) \cdot \phi(v_2)) \ \eta_E(R) \\
                       & = \phi(v_1) \max( h(S \odot R)  \ \eta_E(R),G_{n}(S \odot R) \eta_E(R) \cdot \phi(v_2) ) \\
					   & = \phi(v_1) \max( \widetilde{h}(S,R)  ,\widetilde{G}_{n}(S,R) \cdot \phi(v_2) ). \
                       \end{align*}
Applying techniques that are exaclty analogous to previous steps, we obtain a separating estimate for $\widetilde{\Xi}'$ of the form
\begin{equation*}
\widetilde{\Xi}'_n(S,R,v_1,v_2) = F(S,R,v_1,v_2) + \phi(v_1)\sum_{i_n=1}^{k_n} \psi_{S,i_n}(S)\psi_{R,v,i_n}(R,v_2)
\end{equation*}
where $F$ is appropriately bounded and $\psi_{S,i_n} \in C^0_c({\mathcal S})$. This leads to
\begin{align*}
\delta^N_n(S) & \leq  \ c \ {\varepsilon}  + \sum_{i_n=1}^{k_n} |\psi_{S,i_n}(S)| \cdot \beta^{i_n,N}_{n,{n+1}}
\end{align*}
where $c$ again depends on our truncation conditions and
\begin{equation*}
\beta^{l,N}_{n,n+1} = | \  \left[A^N_n\right]^{-1} \ \frac{1}{N} \sum_{j=1}^N \phi(v^j_n) {\mathbb E}^{{\mathbb Q}}\left[ \psi_{R,v,i_n}(R_n,v_{n+1}) \mid [v^j]_n^{n+1} \ \right] -  A^{-1}_n \ {\mathbb E}^{{\mathbb Q}}\left[ \phi(v^j_n) \psi_{R,v,i_n}(R_n,v_{n+1}) \ \right] \ | .
\end{equation*}
By the SLLN, for each $i_n$, we have that $\lim_{N\to\infty} \beta^{i_n,N}_{n,n+1} = 0$ a.s. As a result, we may find a set $\Omega'^{\delta}_n$ with ${\mathbb Q}'(\Omega'^{\delta}_n) = 1$ such that they all converge which completes the proof.
\end{proof}

\begin{proposition}
Let $n \in \{1,...,M-3\}$ and ${\varepsilon}>0$. There exists a set $\Omega'_n$ such that ${\mathbb Q}'(\Omega'_n)=1$ and for  $\omega' \in \Omega'_n$ we have
\label{lemma:sepest}

\begin{equation}
|a^{N}_n(S,\omega') - a_n(S,\omega')|   \leq c_n^{N}(\omega') \ {\varepsilon} + \delta^N_n(S,\omega') + \alpha^N_n(S,\omega') + \sum_{l=n+1}^{M-2} \beta^N_{l,n}(S,\omega')
\label{eqn:main}
\end{equation}
where
\begin{equation*}
\alpha_n^N(S,\omega') = \sum_{i_{M-1} = 1}^{k_{M-1} } \alpha^{i_{M-1},N}_{M-1,M}(\omega') \sum_{i_{M-2}=1}^{k_{M-2}} \alpha^{i_{M-2},i_{M-1},N}_{M-2,M-1}(\omega') \ \dots \  \sum_{i_n=1}^{k_n} \alpha^{i_n,\dots,i_{M-1}, N}_{i_n,i_{n+1}}(\omega') \cdot \psi^{n,\alpha}_{i_n,\dots,i_{M-1}}(S),
\end{equation*}
\begin{equation*}
\beta_{l,n}^N(S,\omega') =    \sum_{i_{l} = 1}^{m_l } \beta^{l,i_{l},N}_{l-1,l}(\omega') \sum_{i_{l-1}=1}^{m_{l-1}} \beta^{l,i_l,i_{l-1},N}_{l-2,l-1}(\omega') \ \dots \  \sum_{i_{n}=1}^{m_{n}} \beta^{l,i_l,\dots,i_{n}, N}_{n,n+1}(\omega') \cdot \psi^{n,\beta}_{i_{l},\ldots,i_n}(S),
\end{equation*}
\begin{align*}
\lim_{N\to \infty} \alpha^{i_{M-1},N}_{M-1,M}( \omega') = 0, & \ \limsup_{N} \alpha^{ i_j,\ldots,i_l, N}_{k-1,k}(\omega')  < \infty, \ \forall \ k \in \{n,...,M-1\}  \
\end{align*}
and
\begin{align*}
 \forall \ l \in \{ n+1,\ldots,M-2\}, & \; \forall j \in \{ n,\ldots, l  \}  \\
 \lim_{N\to \infty} \beta^{l,i_l,N}_{l-1,l}( \omega') = 0 , & \ \limsup_{N} \beta^{ l,i_j,\ldots,i_l, N}_{k-1,k}(\omega') < \infty .\  \
\end{align*}
The bounds on $\limsup_{N} c^N_n$ depend on our truncation and separability conditions and $\psi^{k,\alpha}_{i_l,\ldots ,i_n },\psi^{k,\beta}_{i_l,\ldots ,i_n}$ satisfy the conditions of Proposition \ref{lemma:geom}.

\label{prop:main}
\end{proposition}

\begin{proof}
We begin by letting $n = M-3$ and by Lemma \ref{lem:identity} and Lemma \ref{lemma:mm2} we have on $\omega' \in \Omega'_{M-2} \cap \Omega^{'\delta}_{M-2}$
\begin{align}
&|a^{N}_{M-3}(S,\omega') - a_{M-3}(S,\omega') | \nonumber \\
&\leq \delta^{N}_{M-3}(S,\omega') + \widetilde{c} \frac{1}{N} \sum_{j=1}^{N}  {\mathbb E}^{{\mathbb Q}}\left[  \left. |a^{N}_{M-2}(S_{M-2},\omega') - a_{M-2}(S_{M-2},\omega') |   \right.  \mid  \left. S_{M-3} = S, [v^j(\omega')]^{M-2}_{M-3} \right. \right] \nonumber \\
&\leq  \delta^N_{M-3}(S,\omega') + c^{N}_{M-2}(\omega')  \nonumber \\
& + \widetilde{c}\sum_{i_{M-1} = 1}^{k_{M-1} } \alpha^{i_{M-1},N}_{M-1,M}(\omega')\sum_{i_{M-2}}^{k_{M-2}} \alpha^{i_{M-2},i_{M-1},N}_{M-2,M-1}(\omega') \frac{1}{N}\sum_{j=1}^N {\mathbb E}^{{\mathbb Q}}\left[ \psi_{S,i_{M-2},i_{M-1}}(S_{M-2}) \mid S_{M-3} = S, [v^j(\omega')]_{M-3}^{M-2}  \right] \label{eqn:C2} \\
& + \widetilde{c} \frac{1}{N} \sum_{j=1}^{N} {\mathbb E}^{{\mathbb Q}}[ \delta^N_{M-2}(S_{M-2},\omega')\ | \ S_{M-3} = S, [v^j(\omega')]_{M-3}^{M-2} \ ] \label{eqn:D2}   \
\end{align}
Line (\ref{eqn:C2}) may be handled just as in the proof of Lemma \ref{lemma:mm2}. As for line (\ref{eqn:D2}), we use Lemma \ref{lemma:deltaest} and write
\begin{align*}
& \frac{1}{N} \sum_{j=1}^{N} {\mathbb E}^{{\mathbb Q}}[ \delta^N_{M-2}(S_{M-2},\omega')\ | \ S_{M-3} = S, [v^j(\omega')]_{M-3}^{M-2} \ ] \\
& \leq c \ {\varepsilon} + \sum_{i_{M-2}=1}^{m_{M-2}} \beta_{M-2,M-1}^{M-2,i_{M-2},N}(\omega')\frac{1}{N}\sum_{j=1}^{N}{\mathbb E}^{{\mathbb Q}}[ \psi_{i_{M-2}}(S_{M-2}) \ | \ S_{M-3} = S, [v^j(\omega')]_{M-3}^{M-2}   ] \
\end{align*}
where $\lim_{N \to \infty} \beta_{M-3,M-2}^{M-2,i_{M-2},N}(\omega') = 0$ and we apply the usual separation technique, the details of which we omit. We then obtain a set $\Omega'_{M-3}$ such that ${\mathbb Q}'( \Omega'_{M-3} ) = 1$ and for $\omega' \in \Omega'_{M-3}$
\begin{align*}
&|a^{N}_{M-3}(S,\omega') - a_{M-3}(S,\omega') |  \\
&\leq c^{N}_{M-3}(\omega') +  \delta^N_{M-3}(S,\omega') \\
& + \sum_{i_{M-1} = 1}^{k_{M-1} } \alpha^{i_{M-1},N}_{M-1,M}(\omega')\sum_{i_{M-2}}^{k_{M-2}} \alpha^{i_{M-2},i_{M-1},N}_{M-2,M-1}(\omega') \sum_{i_{M-3}}^{k_{M-3}} \alpha^{i_{M-3},i_{M-2},i_{M-1},N}_{M-3,M-2}(\omega')\cdot \psi^{M-3,\alpha}_{S,i_{M-3},i_{M-2},i_{M-1} } (S)    \\
& + \sum_{i_{M-2}=1}^{k_{M-2}} \beta^{M-2,i_{M-2},N}_{M-2,M-1}(\omega') \sum_{i_{M-3}=1 }^{m_{M-3}} \beta^{M-2,i_{M-2},i_{M-3},N}_{M-3,M-2}(\omega') \psi^{M-2,\beta}_{i_{M-3},i_{M-2}} (S) \
\end{align*}
which corresponds to $l = M-2$ and the above random variables satisfy the necessary conditions. The remainder of the induction works exactly analogously to the base case and proofs of the previous Lemmas which establishes the Proposition.
\end{proof}

\section{Complexity Reduction and Discussion}
\label{section:pathred}

In the formluation of the algorithm, one simulates many paths of $v_t$ and solves a PDE over ${\mathcal S}$-space for each path. When $d_S \geq 2$ and the inner expectation is not available in closed form, the algorithm can be computationally demanding even with parallelization. For this purpose, we introduce two methods to greatly improve the algorithm's complexity by reducing the number of paths over which the PDE must be solved. The first is a MLMC/multi-grid method as outlined in \cite{Giles}, which is often used in the numerical analysis of stochastic PDEs (SPDEs). We may view our coefficients as arising from the expectations of solutions to SPDEs that are driven by the random variables $\Theta([v])$. MLMC approaches have also appeared in the mixed MC-PDE works of \cite{Ang:2013zl} and \cite{dmd3}. The second is a clustering algorithm where we group paths based on certain sufficient path statistics.

As noted in \cite{Giles}, due to the rapid growth in grid size with increases in level, MLMC methods when applied to SPDEs usually provide at least an order of magnitude in savings. Our clustering method, on the other hand, provides only minor savings. It provides mild improvements to the direct estimator, and its efficiency is sensitive to the number of paths. With that said, our back-tests show that one can be accurate to one or two decimal places with a number of paths $N_v \approx 10 \ 000$, due to the variance reduction from integrating over the variable $S_t$. If $d_S > 2$, for such a choice of $N_v$, it is often advantageous to apply clustering to reduce the number of paths by about a half.

In Sections \ref{section:Heston1} and \ref{section:Heston2}, we demostrate that the methods are effective when computing time-zero prices and exercise boundaries.

\subsection{MLMC/multi-grid Approximations }

We now discuss a MLMC/multi-grid approximation to (\ref{eqn:estcoef}) and (\ref{eqn:lowestim}). It is often the case that solving a PDE at a lower resolution adds a slight bias to the computed prices compared to solutions obtained at a higher resolution. This observation suggests that the PDE solved at a lower resolution can be used as a control variate. A pseudo-code description of the methodology based on this idea, and outlined below, is provided in Algorithm \ref{alg:MLMC} in \ref{app:algstate}.

\subsubsection{MLMC for Computing the Estimated Coefficients}

Assume we have $L$ independent simulations of $v_t$ on $[0,T]$, denoted $\{(v^{(l)}_t)\}_{l=0}^L$, with $N^{v}_l$ paths where $N^{v}_{l} > N^{v}_{l+1}$. Also, let $P_{n}^{(j,l)}(S)$ denote the numerical solution of the conditional PDE on the $j$th path with grid resolutions $N^{{\mathcal S}}_{l}$, in each dimension, over $[t_n,t_{n+1}]$ where $N^{{\mathcal S}}_l < N^{{\mathcal S}}_{l+1}$.

Using the usual multi-level MC appoach, we may write
\begin{equation}
a^{N}_n(S) \approx \left[A^{N_0}_n\right]^{-1} \Bigg[ \frac{1}{N^{v}_0} \sum_{j=1}^{N_0} \phi( v^{j,0}_n ) \cdot  P_n^{(j,0)}(S) + \sum_{k=0}^{L-1} \frac{1}{N_l}\sum_{j=1}^{N^{v}_l} \phi( v^{j,l}_n ) \cdot (P^{j,l+1}_n(S) - P^{j,l}_n(S)   ) \Bigg]\,,
\label{eqn:MLMCcoeff}
\end{equation}
where each of $P^{j,l}_n(S)$ is interpolated to have a resolution that matches the grid at level $L$.

This approach allows us to solve the PDE a low number of times at a fine resolution and many times on a very coarse resolution which can significantly reduces the algorithm's computational costs. A na\"ive implementation may end up carrying out many thousands of interpolations when computing (\ref{eqn:MLMCcoeff}), and also when generating terminal conditions for the PDEs. Our pseudo-code description (Algorithm \ref{alg:MLMC} in \ref{app:algstate}) shows how to keep the total number of interpolations to a minimum.

\subsubsection{MLMC for Low Estimates of the Time Zero Price}

We again carry out $L$ independent simulations of $v_t$ on $[0,T]$, denoted as $\{(v^{(l)}_t)\}_{l=0}^L$, with $N^{v}_l$ paths where $N^{v}_{l} > N^{v}_{l+1}$. Letting $P_0^{(j,l)}(S)$ denote the numerical solution of the conditional PDE on the $j$th path with grid resolution $N^{{\mathcal S}}_{l}$ over $[0,T]$ where $N^{{\mathcal S}}_l < N^{{\mathcal S}}_{l+1}$ we have
\begin{equation*}
V_{l,0}(S) \approx \frac{1}{N^{v}_0} \sum_{j=1}^{N^{v}_0} P_0^{(j,0)}(S) + \sum_{k=0}^{L-1} \frac{1}{N^{v}_l}\sum_{j=1}^{N^{v}_l}  (P_0^{j,l+1}(S) - P_0^{j,l}(S)   )\,,
 \end{equation*}
 where again we interpolate the lower resolution grids to match the highest resolution grid. We do not provide a pseudo-code description for this part as it is relatively simple.

\subsection{Clustering by Sufficient Path Statistics}

The idea is to approximate the sum within (\ref{eqn:estcoef}) using representatives of clusters that we now describe. In Section \ref{section:setup}, we supposed the existence of a finite dimensional functional of the path, $\Theta([v]_n^{n+1})$, such that (\ref{eqn:markov}) holds. We recall that this implies the conditional PDEs do not depend on the value of the path $[v]_n^{n+1}$ at each point in time, but rather on $\Theta_n([v_n]_{n}^{n+1})$ and, e.g., in the Heston model we have
\begin{equation}
\Theta([v]_n^{n+1}) =  \left( \ v_{t_n}, \Lambda([v]_n^{n+1}), v_{t_{n+1}} \ \right) =  \left( \ v_{t_n} , \int_{t_n}^{t_{n+1} } \sqrt{v_s} \ dW^v_s, \int_{t_n}^{t_{n+1}} v_s \ ds, v_{t_{n+1}}  \ \right).
\label{eqn:Hestonstat}
\end{equation}
Given a path statistic $\Theta([v]_n^{n+1} ) \in \dsR^{d_{\Theta}} $ we cluster our $N$ observations of $\Theta([v]_n^{n+1} )$ into $D_n$ clusters using hierarchical clustering. We use this clustering method as it can efficiently handle large values of $D_n$. Given the $k$th cluster we generate a representative $\Theta_k$ as follows
\begin{equation*}
\Theta^{*}_{n,k} = \frac{1}{D_{n,k}} \sum_{j=1}^{D_{n,k}} \Theta([v^{k,j}]^{n+1}_n)\,,
\end{equation*}
where $\{ [v^{k,j}]_{n}^{n+1} \}_{j=1}^{D_{n,k}}$ represents the paths in the $k$th cluster. We then obtain the following approximation for computing our coefficients:
\begin{equation}
a^{N}_n(S) \approx \big[ A^{N}_n \big]^{-1}   \cdot   \frac{1}{N} \sum_{k=1}^{D_n} D_{n,k}   \  \phi(v^{k*}_n) \cdot {\mathbb E}^{{\mathbb Q}}\left[\left.  f_{n+1}^{N}(S_{n+1},v^{k*}_{n+1})\right.  \mid \left. S_{n} = S, [v^{k*}]_n^{n+1} \right. \right]\,.
\label{eqn:pathred}
\end{equation}
Note that $[A_n^{N}]$ is a function of all $N$ paths and is not affected by the clustering algorithm. More formally, before solving a PDE over paths on $[t_n,t_{n+1}]$ we carry out Algorithm \ref{alg:clustering} listed in \ref{app:algstate}. In the approach outlined there, we decide on a fixed number of clusters at each time step (chosen arbitrarily). One may also choose the number of clusters adaptively at each time stage based on the dissimilarity matrix that is produced from their hierarchical tree.

\subsection{Combining MLMC and Clustering }

On a time increment $[t_n,t_{n+1}]$, when computing the coefficients in (\ref{eqn:estcoef}), we begin with the approximation (\ref{eqn:MLMCcoeff}). Next, since typically $N^v_0 \gg N^v_l$ for $l=1,...,L$, we cluster the paths of the $0$-level simulation of $v_t$ and apply approximation (\ref{eqn:pathred}) to compute $\frac{1}{N^v_0} \sum_{j=1}^{N^v_0} \phi( v^{j,0}_n ) \cdot  P_n^{j,0}(S)$.

\subsection{An Overview of Complexity}

While a complete complexity analysis lies outside the scope of this paper, we provide an overview of the sources of costs and error, along with a discussion of the complexity of our MLMC-FST scheme.

The direct estimator has three main stages: simulation of paths, solution of conditional PDEs, and the regression step. Assuming an Euler discretization for the simulation of paths, the FST method to solve the conditional PDEs, and Gauss-Jordan elimination for inverting $A^N_{n}$, we obtain the following expression for the algorithm's cost:
\begin{align*}
C &= C_{\text{Path. Sim.}} + C_{\text{PDE} } + C_{\text{Reg.} }  \\
  &\sim N_v \cdot N_t +  N_S^{d_S} \cdot (\log N_S )\cdot N_v + (N_v\cdot d_B + d_B^3 + d_B^2 \cdot N_S^{d_S} ).  \
\end{align*}
The terms in $C_{\text{Reg.}}$ correspond to the costs of constructing $A^N_n$, computing its inverse and multiplying $[A^N_n]^{-1}$ to a length $d_B$ vector at each of the $N_S^{d_S}$ regression sites.

We make the following rough estimate of the algorithm's error at a single point of our grid, $S \in {\mathcal S}$,
\begin{align*}
{\varepsilon}(S) &\sim {\varepsilon}_{\text{Path. Sim.}} + {\varepsilon}_{\text{PDE} } + {\varepsilon}_{\text{Reg.} }  \\
   &\sim \left(\frac{1}{\sqrt{N_v} } + \frac{1}{N_t} \right) + \left( \frac{1}{N_S^2} + \frac{1}{\sqrt{N_v}} \right) + E(d_B,N_v,N_t,N_S) \
\end{align*}
where determining $E$ is a topic for future papers. For other forms of LSMC, the structure of ${\varepsilon}_{\text{Reg.}}$ has been investigated in \cite{GlassermanYu}, \cite{Stentoft}, \cite{Egloff}. Another important question for future papers is how the quantities $N_v, N_t, N_S, d_B$ terms ought to be optimally balanced, especially in the context of MLMC that we discuss below.

In Sections \ref{section:Heston1} and \ref{section:Heston2}, alongside implementing the hybrid LSMC/PDE algorithm for Heston-type models, we also investigate how the bias, variance, and cost behave on the various levels of our MLMC-FST scheme in the context of single-period expected values. By numerically estimating the rates at which these quantities decay and grow, by the general multi-level theorem (GMLT) (see \cite{Giles2} and \cite{Giles}), we are able to suggest a bound for the complexity of the MLMC-FST component in isolation from the rest of the algorithm. As the number of points in the highest resolution level affects both the bias of the PDE solver and complexity of the regression stage, the choice of the highest level from the GMLT for the single period problem may not be optimal for the regression stage, and hence the entire hybrid algorithm. On the other hand, the GMLT sheds light on the potential improvements brought about by MLMC to the main bottle-neck of the algorithm.

The set up of our tests is a European option analogue of our Bermudan option with the same maturity and exercise function and we take the price of the option at a single initial value as our quantity of interest. We then proceed to plot $ \log_2\left(\left| {\mathbb E}\left[  P^{l}_0(S) - P_{0}(S) \right] \right|\right)$,   $\log_2\left(\mathbb{V}\left[Y_l\right] \right)$, and $\log_2\left( {\mathbb E}\left[C_l\right]\right) $ against level where \[ Y_l(S_0) =   P^{l}_0(S_0) - P^{l-1}_0(S_0)   \] and $C_l$ is the cost of computing $Y_l$, measured in CPU time. Based on the slopes of these plots, by the GMLT, we obtain suggested upper bounds for the complexity of the MLMC-FST scheme. Using the same notation as \cite{Giles}, we denote the slopes of the bias, variance, and cost plots as $\alpha, \beta,$ and $\gamma$, respectively.

Finally, the work of \cite{Giles} provides an optimal allocation scheme for the number of paths on each level. We note, however, that the traditional scheme is for a single period problem as opposed to a multi-period problem, as in our direct estimator. Heuristically, one could use the allocation from the single-period analogue for the multi-period problem, although we leave this for future investigation. Finally, the single-period allocation scheme may be applied to the low estimator.

\section{A $1d + 1d$ Example: Heston Model }
\label{section:Heston1}

\subsection{Model Set Up}
We work in the following form of the Heston model
\begin{align*}
dS_t = S_t(r \ dt + \sqrt{v_t} \ dW^S_t)\,, \qquad \text{and} \qquad
dv_t = \kappa (\theta - v_t ) \ dt + \gamma \sqrt{v_t} \ dW^{v}_t \
\end{align*}
where $d[W^S,W^{v}]_t = \rho \ dt $. Writing $W^{S}_t = \rho\ W^{v}_t + \sqrt{1 - \rho^2} \ W^{\perp,v}_t$, where $W^{\perp,v}$ is a Brownian motion independent of $W^{v}$, and setting $X_t=\log S_t$, we have
\begin{align*}
dX_t & = \big( r - \tfrac{1}{2}v_t \big) \ dt + \sqrt{1 - \rho^2}  \sqrt{v_t} \ dW^{\perp,v}_t + \rho \sqrt{v_t} \ dW^{v}_t \,.
\end{align*}
\subsection{Deriving the Conditional PDE}
To compute ${\mathbb E}\left[\left. e^{-r\Delta t_n} f(S_{n+1},v_{n+1} )\right. \mid \left.S_n = s_i , v_n= v\right. \right] $ we begin by iterating the total expectation and writing our expression in terms of $X_t = \log S_t$. Thus we consider a payoff of the form
\begin{equation*}
g(X,v) = f(\exp(X),v)\,,
\end{equation*}
so that
\begin{equation}
{\mathbb E}\left[\left. e^{-r\Delta t_n} f(S_{n+1},v_{n+1})\right. \mid \left.S_n , v_n\right. \right] = e^{-r\Delta t_n} {\mathbb E}\left[ \left. {\mathbb E}\left[  \left.g(X_{n+1},v_{n+1})\right. \mid \left.[v]^{n+1}_n, X_n\right. \right] \right. \mid \left.X_n , v_n\right. \right]\,.
\label{eqn:expectation-f-to-g}
\end{equation}
Focusing on the inner expectation, treating $[v]_n^{n+1}$ as a deterministic path on $[t_n,t_{n+1}]$, and representing it as a function of time $v_t$, we then write
\begin{align*}
X_t &= Y_t + Z_t\,,\ \quad \text{where}\\
Y_t &= X_{t_n} +  \int_{t_n}^t \big( r - \tfrac{1}{2}v_s  \big)\ ds + \int_{t_n}^t \sqrt{1 - \rho^2} \sqrt{v_s} dW^{v,\perp}_s\,, \quad \text{and} \quad
Z_t = \int_{t_{n}}^{t} \rho \sqrt{v_s} dW^v_s\,. \
\end{align*}
Re-writing the inner expectation in the rhs of \eqref{eqn:expectation-f-to-g} in terms of the process $Y_t$ we have
\begin{equation*}
{\mathbb E}\left[ \left. g(X_{n+1},v_{n+1})\right. \mid \left.[v]^{n+1}_n, X_n\right. \right] = {\mathbb E}\left[ \left.g(Y_{n+1} + Z_{n+1},v_{n+1})\right. \mid \left.[v]^{n+1}_n, Y_n\right. \right]:=u(t,Y_n).
\end{equation*}
By Feynman-Kac theorem, the function $u(t,y)$ can be written as the solution to the following PDE
\begin{equation}
\label{eq:mainpde}
\left\{
\begin{array}{rl}
\partial_t u(t,y) +  a_t \ \partial_y u(t,y) +  b_t \ \partial^2_yu(t,y) \!\!\!&= 0\,,  \\[0.25em]
u(t_{n+1},y) \!\!\!&= g(y + Z_{t_{n+1}},v_{t_{n+1}})\,,
\end{array}
\right.
\end{equation}
where $a_t = r - \frac{1}{2}v_t$ and $b_t = (1-\rho^2)v_t$ are deterministic functions of time.

\subsection{FFT Based Solution}
Taking the Fourier transform of (\ref{eq:mainpde}), as in the FST method of \cite{JJS}, and letting $\widehat{u}(t,\omega)$ denote the Fourier transform of $u(t,y)$ in the second argument, we have
\begin{align*}
\left\{
\begin{array}{rl}
\partial_t \widehat{u}(t,\omega) + (a_ti \omega - \frac{1}{2}b_t\omega^2) \widehat{u}(t,\omega) \!\!\!&= 0\,,  \\
\widehat{u}(t_{n+1},\omega) \!\!\! &= e^{i\omega Z_{t_{n+1}} }  \  \widehat{g}(\omega,v_{n+1}) \,,
\end{array}
\right.
\end{align*}
which is an ODE wrt time for $\widehat{u}$. Thus,
\begin{equation*}
\widehat{u}(t_{n},\omega) =  \widehat{u}(t_{n+1},\omega) \exp\left\{i\omega \int_{t_n}^{t_{n+1}} a_sds - \frac{1}{2} \omega^2 \int_{t_n}^{t_{n+1}} b_s ds   \right\}.
\end{equation*}
Writing $\Psi(\omega,t_n,t_{n+1}) = i\omega Z_{t_{n+1}} +  i\omega \int_{t_n}^{t_{n+1}} a_s \ ds  - \frac{1}{2}\omega^2 \int_{t_n}^{t_{n+1}} b_s \ ds  $  or more compactly $\Psi_{n,n+1}:=\Psi(\cdot,t_n,t_{n+1})$, and using the FST's discretization methodology with fast Fourier transforms (FFT), we have the recursion
\begin{equation*}
u_n = {\text{FFT}}^{-1}\Big[ {\text{FFT}}\Big[ g_{n+1} \Big] \exp(\Psi_{n,n+1})     \Big]\,,
\end{equation*}
where $g_k,u_k$ are discretizations of $g(t_k,y), u(t_k,y)$.

\subsubsection{Statistics for Clustering and Path Reduction}

As discussed in Section \ref{section:pathred} we use (\ref{eqn:Hestonstat}) for our clustering statistics.

\subsection{Numerical Results for the $2d$ Heston Model}
\label{section:results}
In this section, we compare the performance of our algorithm to the full LSMC algorithm and a finite difference implementation for a Bermudan put option written on the Heston model. After outlining our simulation methodology we show results for time-zero prices, optimal exercise boundaries, and MLMC-FST tests. For the finite difference approach, we choose our domain truncation and model parameters similar to \cite{Ikonen}. We differ in that we consider a Bermudan instead of an American option and set the risk free rate to be $2 \%$ as opposed to $10 \%$. We use a simple explicit scheme which we view as a reference solution as opposed to a benchmark that we are trying to beat.

\subsubsection{Model and Option Details}

We price a Bermudan style put option with equally spaced exercise dates.

\begin{table}[ht]
\begin{center}
 \begin{tabular}{||c c c c c c c c c c ||}
 \hline
 $T$  & $K$ & \text{Exercise Frequency} & $S_0$ & $v_0$  & $r$ & $\kappa$ & $\theta$ & $\eta$ & $\rho$  \\ [0.5ex]
 \hline\hline
 1   & 10  & $T/12$    & $10$ & $0.15$  & 0.02 & 5 & 0.16 & 0.9 & $0.1$ \\
 \hline
\end{tabular}
\end{center}
\caption{Parameters used in our option pricing model}
\label{table:parmod}
\end{table}

\subsubsection{Procedure for the LSMC tests}
We run $N_{trial}$ of a full LSMC algorithm as follows.
\begin{enumerate}
\item Simulate $N_{sim}$ paths of $(S,v)$ on $[0,T]$ with Euler discretization $1/N_{step}$.
\item Carry out the direct estimate approach for LSMC using the $N^{S,v}_{LSMC}$ simulation of $(S,v)$ and a two state-variable basis with dimension $d_{B,full}$.
\begin{itemize}
\item Retain the estimated coefficients
\end{itemize}
\item Independently simulate $N^{S,v}_{LSMC}$ paths of $(S,v)$ on $[0,T]$ with Euler discretization $1/N_{step}$.
\item Carry out the low estimate using the simulation of $(S,v)$ and the coefficients from the previous stages.
\end{enumerate}

\subsubsection{Procedure for the LSMC-PDE tests}
We run $N_{trial}$ of the LSMC/PDE and multi-level LSMC/PDE with clustering as follows. For our examples, each calculation has two levels with $N^{v}_0 > N^{v}_{1,i}$, $N^{{\mathcal S}}_{0}< N^{{\mathcal S}}_{1,i}$ with $i = 1,2,3$. That is, we employ three types of two-level MLMC approximations.

\begin{enumerate}
\item Independently simulate $N^v_{0}$ and $N^v_{1}$ paths of $(v,\Delta W^v)$ on $[0,T]$ with Euler discretization $1/N_{step}$.
\item For the $k=0$ paths, carry out the high biased LSMC-PDE algorithm with $ N^{{\mathcal S}}_{1,1},N^{{\mathcal S}}_{1,2}, N^{{\mathcal S}}_{1,3} $ points in the log-space grid. Use a one state-variable basis of dimension $d_{B,mix}$.
\begin{itemize}
\item Retain the estimated coefficients for each PDE resolution
\end{itemize}
\item Carry out the direct estimator of the multi-level LSMC/PDE with clustering algorithm using paths from levels $0,1$ of $(v,\Delta W^v)$ with resolutions $N^{{\mathcal S}}_{0}$ and $N^{{\mathcal S}}_{1,1}, N^{{\mathcal S}}_{1,2}, N^{{\mathcal S}}_{1,3}$ in the log-space grid. Use a one dimensional basis of dimension $d_{B,mix}$ and $N_{clust}$ clusters.
\begin{itemize}
\item Retain the estimated coefficients for each PDE resolution
\end{itemize}
\item Independently re-simulate $N^v_{0}$ and $N^v_{1}$ paths of $(v,\Delta W^v)$ on $[0,T]$ with Euler discretization $1/N_{step}$.
\begin{itemize}
\item Carry out the low estimator for LSMC/PDE and multi-level LSMC/PDE with clustering algorithms in the same manner for each grid resolution using estimated coefficients from the previous stages.
\end{itemize}
\end{enumerate}

\subsubsection{Procedure for the MLMC-FST Level Tests}
\begin{enumerate}
\item For each level $l \in \{l_1,\ldots,l_n \}$, compute $\left|{\mathbb E}\left[ P_0^l(S_0) - P_0(S_0) \right]\right|$ over $N_{trials}$. As mentioned before, the quantity $P_0(S_0)$ corresponds to the solution of the conditional PDE solved at a relatively high resolution, $N_{S,t}$.
\item For each level $l \in \{l_1,\ldots,l_n \}$, compute \[ {\mathbb V} \left[ Y_l(S_0) \right] = {\mathbb E}\left[ \left(Y_l(S_0) - {\mathbb E}\left[Y_l(S_0) \right]  \right)^2  \right] \] over $N_{trials}$.
\item For each level $l \in \{l_1,\ldots,l_n \}$, measure the expected CPU time for computing $C_l$, ${\mathbb E}\left[C_l\right]$ over $N_{trials}$.
\end{enumerate}

We choose $S_0$ to be our ATM-point, although our back-tests show that the results are not sensitive to this choice.

\subsubsection{Parameters used in LSMC and LSMC-PDE trials}

For LSMC-PDE, we use a monomial basis $\{\phi_l\}_{l=0}^{\infty}$ where $\phi_l(x) = x^l$. For LSMC we use monomials and powers of the payoff function $h$. By a two state-variable basis, we mean the family of functions of the form $\phi_k(x)\cdot \phi_l(y) $. We refer to the degree of the basis as the highest power of functions in our basis. We show our full LSMC and LSMC-PDE parameters in Table \ref{table:parsim} which also includes the other FST discretization parameters.
\begin{table}[ht]
\begin{center}
 \begin{tabular}{|| c c c c c c  ||}
 \hline
    $N_{trial}$ & $N^{S,v}_{LSMC}$ & $N^{v}_{sim,0}$ & $N^{v}_{sim,1}$ & $N_{clust}$ & $N_{step}$ \\ [0.5ex]
 \hline\hline
  100 &  100 000  & 10 000 & 100 & 4 500 & 1 000 \\
 \hline
 \hline
\end{tabular}
\end{center}
\begin{center}
 \begin{tabular}{|| c c c c c c  ||}
 \hline
    $N^{\cal_S}_{0}$ & $N^{\cal_S}_{1,1} $ & $N^{\cal_S}_{1,2} $ & $N^{\cal_S}_{1,3}$ & $(\log S/S_0)_{min}$  & $(\log S/S_0)_{max}$ \\ [0.5ex]
 \hline\hline
  $2^6$ & $2^7$  & $2^8$ & $2^9$ &  $-3$  & $3 $ \\
 \hline
 \hline
\end{tabular}
\end{center}
\begin{center}
 \begin{tabular}{|| c c ||}
 \hline
    $\text{deg}_{full}$ & $\text{deg}_{mix}$  \\ [0.5ex]
 \hline\hline
  3 & 5   \\
 \hline
 \hline
\end{tabular}
\caption{Parameters used in the LSMC and LSMC-PDE numerical simulations}
\end{center}
\label{table:parsim}
\end{table}

\begin{table}[ht]
\begin{center}
 \begin{tabular}{||c c c||}
 \hline
 $N_{S,l}$  & $N_{S,t}$ &  $N_{trials} $    \\ [0.5ex]
 \hline
 $2^5, \dots , 2^9 $   & $2^{10}$    & $1 \ 000$ \\
 \hline
\end{tabular}
\end{center}
\caption{Parameters used in our MLMC-FST tests  }
\end{table}

\subsubsection{Set up for the Finite Difference Method}
We employ an explicit finite difference scheme with equally spaced grid points. Our boundary conditions are as outlined in \cite{Ikonen}. Let $N_s$ and $N_v$ be the number of points in our spatial and volatility grids, $N_t$ be the number of time steps taken, and $S_{min}, S_{max}, v_{min}, v_{max}$ be the end points of our grids. These parameters have been found to be stable in determining the time-zero prices and exercise boundary.
\begin{table}[ht]
\begin{center}
 \begin{tabular}{||c c c c c c c ||}
 \hline
 $N_{s}$ & $N_{v}$ & $N_{t}$ & $v_{min}$ & $v_{max}$ & $S_{min}$ & $S_{max}$     \\ [0.5ex]
 \hline
 $2^9$ & $2^7$ & $100 \ 000$ & 0 & 1& 0 & 53   \\
 \hline
\end{tabular}
\end{center}
\caption{Parameters used for the finite difference computation. }
\label{table:parfd}
\end{table}

\subsubsection{Results for the Time Zero Price }
Tables \ref{table:resultslsmcpde}, \ref{table:resultslsmcpdeprml} and \ref{table:resultslsmc} compare the results of all three simulation algorithms.

\begin{table}[ht]
\begin{center}
 \begin{tabular}{||c c c c c||}
 \hline
 $N_{s}$   & Est. Type & Mean Price     & Std Dev  & Run Time (s)    \\ [0.5ex]
 \hline\hline
 $2^{9}$  &  Direct       &  1.4530       & 0.0014    &  10     \\
 		  &  Low          &  1.4528       & 0.0014    &  8    \\
 $2^{8}$  &  Direct       &  1.4532       & 0.0014    &  7 \\
 		  &  Low          &  1.4530       & 0.0014   &   6  \\
 $2^{7}$  &  Direct       &  1.4539       & 0.0014   &   5   \\
 		  &  Low          &  1.4537       & 0.0014   &   4  \\
 \hline
\end{tabular}
\end{center}
\caption{Resulting ATM price statistics for the LSMC-PDE algorithm. The reference value obtained from finite difference is 1.4528. 10 000 paths with 4th degree basis functions.  }
\label{table:resultslsmcpde}
\end{table}

\begin{table}[ht]
\begin{center}
 \begin{tabular}{||c c c c c||}
 \hline
 $N^{S}_{1}$   & Est Type & Mean Price     & Std Dev & Run Time (s)     \\ [0.5ex]
 \hline\hline
 $2^{9}$  &  Direct       &  1.4531       & 0.0014  &  23        \\
 		  &  Low          &  1.4530       & 0.0015  &   6       \\
 $2^{8}$  &  Direct       &  1.4534       & 0.0014  &  23     \\
 		  &  Low          &  1.4532       & 0.0015  &   6    \\
 $2^{7}$  &  Direct       &  1.4540       & 0.0014  &  23    \\
 		  &  Low          &  1.4538       & 0.0015  &   6   \\
 \hline
\end{tabular}
\end{center}
\caption{Resulting ATM price statistics for the LSMC-PDE algorithm with clustering and MLMC. $N^v_0 = 10 \ 000$ clustered to $4 \ 500$, $N^v_1 = 100$ and $N^{S}_0 = 2^6$ }
\label{table:resultslsmcpdeprml}
\end{table}

\begin{table}[ht]
\begin{center}
 \begin{tabular}{||c c c c||}
 \hline
 Estimate Type & Mean Price     & Std Dev  & Run Time (s)  \\ [0.5ex]
 \hline\hline
 Direct         &  1.4494       & 0.0020   & 48      \\
 Low            &  1.4487       & 0.0023   & 38      \\
 \hline
\end{tabular}
\end{center}
\caption{Resulting ATM price statistics for the LSMC algorithm. The reference value obtained from finite difference is 1.4528. 500,000 paths used in total with 5th degree polynomials and powers of the payoff functions.  }
\label{table:resultslsmc}
\end{table}

\subsubsection{Optimal Exercise Boundaries}
\label{sss:oeb}
We display our optimal exercise boundaries in \ref{section:OEB}.  Figures \ref{fig:lsmcpdebdy}, \ref{fig:lsmcpdeprmlbdy} and \ref{fig:tvrbdy} are obtained by computing the averages of the exercise boundaries for each trial and taking the absolute value of the difference between the average and the finite difference boundary. For each pair $(S,v)$ we have a probability that the algorithm will correctly decide on exercising. The yellow regions indicate that with probability 1, the algorithm is correct, whereas the blue regions indicate that with probability 1 the algorithm is incorrect.

The volatility axes for each exercise date was obtained by taking the mean of the 5th and 95th percentile.

\subsubsection{Results for MLMC-FST tests }
\label{subsubsec:mlmcfst}
Table \ref{table:fst1} features the results of MLMC-FST computation. Figures \ref{fig:fst1} provide plots of the estimates.

\begin{table}[ht]
\begin{center}
 \begin{tabular}{||c c c||}
 \hline
 Slope Type & Estimate     & $R^2$   \\ [0.5ex]
 \hline\hline
 $\alpha$      &  2.12    & 0.9993        \\
 $\beta$       &  4.16    & 0.9989         \\
 $\gamma$      &  0.30    & 0.9409       \\
 \hline
\end{tabular}
\end{center}
\caption{Resulting values for MLMC-FST tests for the Heston model. The parameters $\alpha, \beta, \gamma$ correspond to slopes of logarithmic bias, variance, and cost.  }
\label{table:fst1}
\end{table}

\section{A $2d + 2d$ Example: Multi-Dimensional Heston Model}
\label{section:Heston2}
In this section, we price a Bermudan put option on the max of two stocks, each with stochastic volatility. Since the option pricing PDE for this problem is 4-dimensional it is typically handled via LSMC.

\subsection{Model Set Up}

We work with the following version of the two stock price Heston model
\begin{align*}
dS^{(1)}_t = S^{(1)}_t(r dt + \sqrt{v^{(1)}_t} dW^{(1)}_t)\,,
&&
dv^{(1)}_t = \kappa_1(\theta_1 - v^{(1)}_t) dt + \eta_1 \sqrt{v^{(1)}_t} dW^{(3)}_t\,, \\
dS^{(2)}_t = S^{(2)}_t(r dt + \sqrt{v^{(2)}_t} dW^{(2)}_t)\,,
&&
dv^{(2)}_t = \kappa_2(\theta_2 - v^{(2)}_t) dt + \eta_2 \sqrt{v^{(2)}_t} dW^{(4)}_t \
\end{align*}
where $(W^{(i)})_{i=1}^4$ is a 4-dimensional Brownian motion with full correlation structure $\rho = [\rho_{i,j}]$. Switching to log-space for our asset prices $X^{(i)}_t = \log S_t^{(i)}$, $i=1,2$, and applying the Cholesky decomposition to the matrix $\rho$, we obtain the following system.
\begin{align*}
dX^{(1)}_t &= ( r - \tfrac{1}{2}v^{(1)}_t ) dt + \sqrt{v^{(1)}_t} \sum_{j=1}^4 a_{1,j} dB^{(j)}_t\,, \\
dX^{(2)}_t &= ( r - \tfrac{1}{2}v^{(2)}_t ) dt + \sqrt{v^{(2)}_t} \sum_{j=2}^4 a_{2,j} dB^{(j)}_t\,, \\
dv^{(1)}_t   &= \kappa_1(\theta_1 - v^{(1)}_t) dt + \eta_1 \sqrt{v^{(1)}_t} \sum_{j=3}^4 a_{3,j} dB^{(j)}_t\,,  \\
dv^{(2)}_t   &= \kappa_2(\theta_2 - v^{(2)}_t) dt + \eta_2 \sqrt{v^{(2)}_t} a_{4,4} dB^{(4)}_t\,,
\end{align*}
where $(B^{(i)}_t)_{i=1}^4$ are independent Brownian motions and $[a_{i,j}]$ is an upper triangular matrix satisfying the usual conditions with $a_{4,4} = 1$, and $\rho = aa^\intercal$.

\subsection{Conditional PDE and FFT Based Solutions}

We follow the same approach as in Section \ref{section:Heston1}, to compute ${\mathbb E}[e^{-r\Delta t_n}f(S_{n+1},v_{n+1} )\ | \ S_n = S, v_n = v ]$ where $S_t = (S^{(1)}_t,S^{(2)}_t)$, $v_t = (v^{(1)}_t, v^{(2)}_t)$. Again, writing $g(X,v) =  f(\exp(X),v)$ we have
\begin{equation}
\label{eqn:expectation-f-to-g-2d}
{\mathbb E}\left[\left.f(S_{n+1},v_{n+1})\right. \mid \left.S_n  , v_n \right. \right] =   {\mathbb E}\left[ \left. e^{-r\Delta t} {\mathbb E}\left[  \left.  g(X_{n+1},v_{n+1})\right. \mid \left.[v]^{n+1}_n, X_n\right. \right] \right. \mid \left.X_n , v_n\right. \right]\,. \\
\end{equation}
Next focusing on the inner expectation, and treating $[v]_n^{n+1} = ([v^{(1)}]_n^{n+1}, [v^{(2)}]_n^{n+1} )$ as a deterministic path in $\dsR^2$, leads to the following decompositions
\begin{align*}
X_t^{(1)} &= Y_t^{(1)} + Z_t^{(1)}\,,
\\
Y_t^{(1)} &= X^{(1)}_{t_n} + \int_{t_n}^t ( r - \tfrac{1}{2}v^{(1)}_s ) ds + \int_{t_n}^t \sqrt{v^{(1)}_t} \sum_{j=1}^2 a_{1,j} dB^{(j)}_t\,,
\\
Z_t^{(1)} &= \int_{t_n}^t \sqrt{v_t^{(1)}} \sum_{j=3}^4 a_{1,j} dB^{(j)}_t\,,
\end{align*}
and similarly
\begin{align*}
X_t^{(2)} &= Y_t^{(2)} + Z_t^{(2)}\,,
\\
Y_t^{(2)} &= X^{(2)}_{t_n} + \int_{t_n}^t ( r - \tfrac{1}{2}v^{(2)}_s ) ds + \int_{t_n}^t \sqrt{v^{(2)}_t}  a_{2,2} dB^{(2)}_t \,,
\\
Z_t^{(2)} &= \int_{t_n}^t \sqrt{v_t^{(2)}} \sum_{j=3}^4 a_{2,j} dB^{(j)}_t .
\end{align*}
Writing the inner expectation on the rhs of \eqref{eqn:expectation-f-to-g-2d} in terms of $(Y^{(1)}_t,Y^{(2)}_t)$ leads to
\begin{align*}
&{\mathbb E}\left[ \left. g(X^{(1)}_{n+1},X^{(2)}_{n+1},v^{(1)}_{n+1},v^{(2)}_{n+1})\right. \mid \left.[v]^{n+1}_n, X_n\right. \right] \\
& \quad = {\mathbb E}\left[ \left.g(Y^{(1)}_{n+1} + Z^{(1)}_{n+1},Y^{(2)}_{n+1} + Z^{(2)}_{n+1},v^{(1)}_{n+1},v^{(2)}_{n+1} )\right. \mid \left.[v]^{n+1}_n, Y_n = y_n \right. \right]
 =: u(t,y_1,y_2) \
\end{align*}
Applying the Feynman-Kac theorem, $u$ satisfies a PDE, which we then take the Fourier transform of, solve the transformed PDE in Fourier space, and then transform back to real space. With the FST's discretization methodology, we have
\begin{equation*}
u_n = \text{FFT}_2^{-1}[ \text{FFT}_2[g_{n+1}] \exp(\Psi_{n,n+1})  ].
\end{equation*}
where  $\text{FFT}_2$ denotes the 2$d$ FFT, and
\begin{align}
\Psi(\omega_1,\omega_2,t_n,t_{n+1}) =&\; i\omega_1 Z^{(1)}_{t_{n+1}} + i \omega_2 Z^{(2)}_{t_{n+1}} \nonumber \\
& + i \omega_1 \int_{t_n}^{t_{n+1}} ( r - \tfrac{1}{2} v^{(1)}_t )dt  + i \omega_2  \int_{t_n}^{t_{n+1}} ( r - \tfrac{1}{2} v^{(2)}_t )dt  \nonumber \\
& -\frac{1}{2} ( a^2_{1,1} + a^2_{1,2} )  \omega_1^2 \int_{t_n}^{t_{n+1}} v^{(1)}_t dt  - \tfrac{1}{2} a^2_{2,2}\ \omega^2_2  \int_{t_n}^{t_{n+1}} v^{(2)}_t dt \nonumber \\
& - a_{1,2} \  a_{2,2} \ \omega_1 \ \omega_2  \int_{t_n}^{t_{n+1}} \sqrt{v^{(1)}_t\ v^{(2)}_t} \ dt  .\
\label{eqn:charexp2d}
\end{align}

\subsubsection{Statistics for Clustering and Path Reduction}

From (\ref{eqn:charexp2d}), we identify the following path statistic for $[v]_n^{n+1} = ([v^{(1)}]_n^{n+1}, [v^{(2)}]_n^{n+1} ) $:
\begin{align*}
\Lambda_n([v]_n^{n+1}) =&\; \left(\; \int_{t_n}^t \sqrt{v_t^{(1)}} \sum_{j=3}^4 a_{1,j} dB^{(j)}_t, \;\int_{t_n}^t \sqrt{v_t^{(2)}} \sum_{j=3}^4 a_{2,j} dB^{(j)}_t , \right.\\
&\quad \left. \int_{t_n}^{t_{n+1}} v_t^{(1)} dt ,  \int_{t_n}^{t_{n+1}} v_t^{(2)} dt, \;\int_{t_n}^{t_{n+1}} \sqrt{v_t^{(1)}} \sqrt{v_t^{(2)}} dt \right)
\end{align*}
so that
\begin{equation*}
\Theta([v]_n^{n+1}) = ( v^{(1)}_{t_n}, v^{(2)}_{t_n} , \Lambda([v]_n^{n+1}), v^{(1)}_{t_{n+1}},v^{(2)}_{t_{n+1}} ).
\end{equation*}

\subsection{Numerical Results for the $4d$ Heston Model}
\subsubsection{Model and Option Details}
We price a Bermudan style max-put option with equally spaced exercise dates. Specifically, our payoff is
\begin{equation*}
h(S^{(1)},S^{(2)} ) = (K - \max(S^{(1)},S^{(2)} ))_+
\end{equation*}
where $K$ is the option's strike. We use the following parameters
\begin{table}[ht]
\begin{center}
 \begin{tabular}{||c || c || c || c || c || c || c || c || c || c || c || c || c || c   ||}
 \hline
 $T$  & $K$ & \text{Exercise Frequency} & $r$ & $S^{(1)}_0$  & $v^{(1)}_0$ & $\kappa_1$ & $\theta_1$ & $\eta_1$  & $S^{(2)}_0$  & $v^{(2)}_0$ & $\kappa_2$ & $\theta_2$ & $\eta_2$ \\ [0.5ex]
 \hline\hline
 1   & 10  & $T/12$  & 0.025 & 10   & 0.45  & 1.52 & 0.45 & 0.4  & 10   & 0.3  & 1.3 & 0.30 & 0.43   \\
 \hline
\end{tabular}
\[
\begin{matrix}
[\rho_{i,j}]
\end{matrix}
=
\begin{bmatrix}
    1       &   \rho_{S^{1},S^{2}} & \rho_{S^1,v^1}  & \rho_{S^1,v^2} \\
    \rho_{S^2,S^1}  & 1        & \rho_{S^2,v^1}  & \rho_{S^2,v^2} \\
    \rho_{v^1,S^1}  &  \rho_{v^1,S^2} &     1  & \rho_{v^1,v^2} \\
    \rho_{v^2,S^1}  & \rho_{v^2,S^2}   & \rho_{v^2,v^1} &  1 \
\end{bmatrix}
=
\begin{bmatrix}
    1       &   0.2 &   -0.3   & -0.15 \\
    0.2     &     1 &   -0.11  & -0.35 \\
    -0.3    & -0.11 &    1     &  0.2 \\
    -0.15   & -0.35 &    0.2   &  1 \
\end{bmatrix}
\]

\end{center}

\caption{Parameters used in our option pricing model. }
\label{table:parmod2}
\end{table}

\subsubsection{Procedure for Comparing the Simulation Algorithms and MLMC-FST Tests}
For comparing our results, we run $N_{trial}$ of a full LSMC and the multi-level LSMC-PDE with clustering algorithms. Our procedure has identical structure as in Section \ref{section:Heston1} except that we do not implement LSMC-PDE without complexity reduction methods. Our procedure for the MLMC-FST tests are also identical to Section \ref{subsubsec:mlmcfst}. Hence we omit the exact details.

\begin{table}[ht]
\begin{center}
 \begin{tabular}{|| c c c c c c  ||}
 \hline
    $N_{trial}$ & $N^{S,v}_{LSMC}$ & $N^{v}_{sim,0}$ & $N^{v}_{sim,1}$ & $N_{clust}$ & $N_{step}$ \\ [0.5ex]
 \hline\hline
  100 &  500 000  & 10 000 & 100 & 4 500 & 1 000 \\
 \hline
 \hline
\end{tabular}
\end{center}
\begin{center}
 \begin{tabular}{|| c c c c c c  ||}
 \hline
    $N^{\cal_S}_{0}$ & $N^{\cal_S}_{1,1} $ & $N^{\cal_S}_{1,2} $ & $N^{\cal_S}_{1,3}$ & $(\log S^{(i)}/S^{(i)}_0)_{min}$  & $(\log S^{(i)}/S^{(i)}_0)_{max}$ \\ [0.5ex]
 \hline\hline
  $2^6$ & $2^7$  & $2^8$ & $2^9$ &  -3  & 3  \\
 \hline
 \hline
\end{tabular}
\end{center}
\begin{center}
 \begin{tabular}{|| c c||}
 \hline
    $\text{deg}_{full}$ & $\text{deg}_{mix}$  \\ [0.5ex]
 \hline\hline
  4 & 4  \\
 \hline
 \hline
\end{tabular}
\end{center}
\caption{Parameters used in the LSMC and LSMC-PDE numerical simulations. }
\label{table:parsim2}
\end{table}

\subsubsection{Results for the Time Zero Price}
Table \ref{table:resultslsmcpdeprml2d} and \ref{table:resultslsmc2d} compare the results of a full LSMC simulation and the LSMC-PDE algorithm with MLMC and clustering.

\begin{table}[ht]
\begin{center}
 \begin{tabular}{||c c c c c||}
 \hline
 $N_{s,high}$ &  Est Type & Mean Price     & Std Dev & Run Times (s)            \\ [0.5ex]
 \hline\hline
$2^{9}$  &  Direct       &  1.1852        &  0.0057    & 157     \\
		 &  Low          &  1.1818        &  0.0061      & 130   \\
$2^{8}$  &  Direct       &  1.1858        &  0.0056    & 84    \\
		 &  Low          &  1.1823        &  0.0061      & 83    \\
$2^{7}$  &  Direct       &  1.1878        &  0.0056    & 68      \\
		 &  Low          &  1.1842        &  0.0056      & 74   \\
 \hline
\end{tabular}
\end{center}
\caption{Resulting ATM price statistics for LSMC-PDE algorithm w/ clustering and MLMC. We fix $N_{s,low} = 2^6, N_{sim,high} = 10 000, N_{sim,low} = 100$,  and $N_{cl} = 4500$. 3rd degree polynomials were used for a total of 10 basis functions. }
\label{table:resultslsmcpdeprml2d}
\end{table}

\begin{table}[ht]
\begin{center}
 \begin{tabular}{||c c c c||}
 \hline
 Estimate Type & Mean Price     & Std Dev   &  Run Times (s)              \\ [0.5ex]
 \hline\hline
 Direct         &  1.2121       &  0.0016  & 89        \\
 Low            &  1.1765       &  0.0018  & 74      \\
 \hline
\end{tabular}
\end{center}
\caption{Resulting ATM price statistics for the LSMC algorithm. 500,000 paths were used with 4th degree polynomials and powers of the payoff function}
\label{table:resultslsmc2d}
\end{table}

\subsubsection{Optimal Exercise Boundaries}
As our OEBs are now four-dimensional objects, we are not able to fully visualize them. To get some idea of their quality, however, we can analyze their slices. A natural choice are $v$-slices, as we already have some intuition about how they should look based on the payoff function. We consider the ``central slice" for the OEB, i.e. where $v_1 = \theta_1$ and $v_2 = \theta_2$, which correspond to the ``best" slice since we have the most amount of data around that point. In \cite{Farahany} other slices that are near the central slice are shown with similar results.

In \ref{section:OEB}, we show the central slice of the optimal exercise boundaries averaged over all 100 trials. Yellow regions indicate holding with probability 1, blue regions indicate exercising with probability 1. Our LSMC-PDE boundary corresponds to a highest-level resolution of $N_S = 2^8$.

\subsubsection{Results for MLMC-FST tests }
\label{subsubsec:mlmcfst2}
Table \ref{table:fst2} features the results of MLMC-FST computation. Figures \ref{fig:fst1} provide plots of the estimates.

\begin{table}[ht]
\begin{center}
 \begin{tabular}{||c c c||}
 \hline
 Slope Type & Estimate     & $R^2$   \\ [0.5ex]
 \hline\hline
 $\alpha$      &  2.03    & 0.9995        \\
 $\beta$       &  2.63    & 0.9640         \\
 $\gamma$      &  1.65    & 0.9785       \\
 \hline
\end{tabular}
\end{center}
\caption{Resulting values for MLMC-FST tests for the multi-dimensional Heston model. The parameters $\alpha, \beta, \gamma$ correspond to slopes of logarithmic bias, variance, and cost.  }
\label{table:fst2}
\end{table}

\section{Discussion of Results}

We provide a discussion of the results from the two examples we investigated above. We first provide example-specific comments for each model, followed by comments which apply to both examples.

\subsection{Heston model}

From Figures \ref{fig:tvrbdy}, we see that the LSMC boundary is fairly inaccurate for all exercise dates, with greater accuracy at the later dates. This is due to the accumulation of errors as one moves recursively through time and the lack of paths that are in-the-money at earlier exercise dates. Based on Figures \ref{fig:lsmcpdebdy} and \ref{fig:lsmcpdeprmlbdy}, we see the errors in the mixed algorithm appear only along the interface of the holding and exercise regions and it is able to closely approximate the true boundary across all exercise dates. The consistency arises from the comparison of the exercise value and continuation value being broken into a collection of cross-sectional comparisons for each $S \in {\mathcal S}$ as discussed in Section \ref{section:LSMC-PDEalg}. At time $n$, for each $S' \in {\mathcal S}$, we have an estimate of the function $C^N_n(S',v)$ which is monotonically increasing in $v$. Our goal is then to locate the points $v^* \in \dsR$ such that $C^N(S',v^*) = h(S')$ which is simple to compute since $h(S')$ is a deterministic constant. We see some noise in determining $v^*$ across trials which is due to the randomness in $C^N(S',\cdot)$, however, it is minimal, especially when comparing with the LSMC boundary.

As indicated in Tables \ref{table:resultslsmcpde}, \ref{table:resultslsmcpdeprml}, \ref{table:resultslsmc} the mixed algorithm performs considerably better than the full LSMC approach in determining the time zero price. This accuracy in our prices follows from the quality of our optimal exercise boundaries (OEBs). We note that the direct estimator and low estimator are fairly close to each other and suggest that the low estimator is somewhat redundant in this case. The full LSMC algorithm also provides a direct estimate of the price that is also close to the low estimator, however, this price is biased too low as evidenced by the reference value. The main point to note here is that the low estimates of the LSMC-PDE algorithms are higher than the low estimates of the LSMC algorithm and are, by definition, better.

Lastly, there is some discrepancy in the prices obtained by finite differences and the LSMC-PDE algorithms, mostly in the fourth decimal place. This discrepancy is likely due to the choice of type/number of basis functions, path discretization methods, interpolation used to obtain time-zero prices, and slight errors in the OEBs.

We also see LSMC-PDE algorithm with the complexity reduction methods is able to closely approximate the true exercise boundaries. Our run times show that the high estimate for the LSMC-PDE algorithm is more efficient than the complexity reduced version of the algorithm. This is owed to the 1$d$ FFT being faster than the use of hierarchical clustering. This tendency goes away when $d_S > 1$.

It is worth noting that the PDE approach is best for this $2d$ example as it has lower complexity. This example is used to demonstrate the value of the LSMC-PDE algorithm in comparison to standard LSMC and to show the effectiveness of the complexity reduction methods.

\subsection{Multi-dimensional Heston model}

Our slices pictured in Figure \ref{fig:2p2}, \ref{fig:lsmc2p2} show that the LSMC-PDE algorithm tends to be more precise than standard LSMC, with LSMC being only $80\%$ accurate in certain regions. We also note that they define fairly different regions for many time stages. As the low estimator produced by the LSMC-PDE algorithm is higher than that of LSMC, we are more inclined the trust the slice provided by the hybrid algorithm. We note, however, than showing the results of the central slice may only suggest that the hybrid algorithm's boundary is better.

The tightness of the direct and low estimator for LSMC-PDE suggest that, again, the low estimator is somewhat redundant. The same cannot be said for LSMC, in this case, as the direct estimator is considerably higher than the low estimator. We also note that our LSMC-PDE prices have also converged at a resolution of $2^7$, suggesting that the higher levels for this example are not necessary as well. On the other hand, the higher levels provide clearer pictures of the OEBs.

From Table \ref{table:resultslsmcpdeprml2d} we have that clustering reduces the costs of the direct estimator quite a bit, to the point that the regression step's run times are diminished for resolutions $2^7$ and $2^8$. Interesting, as the resolution increases from $2^8$ to $2^9$ the runtimes of the low and direct estimator jump up considerably compared to the increase from $2^7$ to $2^8$.

\subsection{Overall Comments}

Our choice of number of clusters and two-level scheme were for simplicity and conservativeness, and hence not optimized, as alluded to before.

One issue that arises from our results is the question of whether the hybrid algorithm is more efficient. For the Heston model, we see that the hybrid algorithm provides a great deal of variance reduction in time zero prices and OEBs especially given that LSMC was allowed a larger computational budget. It is interesting to note, however, that while the standard deviations for prices are relatively close to each other, the OEBs for LSMC ``seem" much noisier than the hybrid algorithm's, although this is based on subjective examination. For the multi-dimensional Heston model, we note that the variance of the hybrid algorithm's prices, for a comparable computational budget, is higher, although the hybrid algorithm's central slice is less noisy. This again points to a disconnect in the variance of the two estimated objects. On the other hand, it is important to note that for the hybrid algorithm, we obtain prices for all $S \in {\mathcal S}$ and not simply the ATM point, as in standard LSMC. Having the prices for these $S \in {\mathcal S}$ also allows us to compute the $\Delta$ and $\Gamma$ of the option without any extra simulations. As a result, the hybrid algorithm provides considerably more information than standard LSMC. Also, the information provided by the hybrid algorithm tends to be more accurate in terms of quality of the boundary and bias, as noted before.

Next, we turn to our MLMC-FST tests. From Table \ref{table:fst1} and \ref{table:fst2} and Figures \ref{fig:fst1} and \ref{fig:fst2}, by the GMLT it is suggested, since $\alpha > \frac{1}{2} \min(\beta,\gamma)$, and $\beta > \gamma$, that optimal complexity is attainable. That is, there exists a collection of levels and a path allocation such that one can obtain a RMSE of $O({\varepsilon})$ with cost $O({\varepsilon}^{-2})$ (up to a certain time-step resolution). This, again, is for our single-period computations in isolation from the entire hybrid algorithm. As one moves from the $2d$ example to the $4d$ example, we see that the accuracy across levels, measured by $\alpha$, remains unchanged. On the other hand, the variances and costs increase considerably as we increased the number of volatility processes and dimension of $S$. These results suggest that we may not conclude optimal complexity, via the GMLT, in the case of a $3d + 3d$ analogue, however, we leave this investigation to future study.

All run times were measured on a single core of a 3.60 GHz machine with 128 GB of RAM using  \textsc{Matlab} 2016. For the results in Tables \ref{table:resultslsmcpde} to \ref{table:resultslsmc2d} all run-times were measured by repeating the following process three times: clearing the memory, running the algorithm, and measuring the run time. The reported run time was then taken to be the minimum of the three.

\section{Conclusions}

By combining conditional PDEs with the LSMC approach of \cite{TV}, we have developed an algorithm that greatly improves traditional LSMC methods in the context of SV models for medium dimensional problems. We find it improves the OEBs and prices obtained by a full LSMC simulation for a given computational budget.

From a theoretical perspective we provide a proof of almost-sure convergence which uses geometric arguments and requires a few modifications to the basis, regression matrix, and payoff function. The proof is also highly dependent on separability of the model. Future research into theoretical properties of the algorithm ought to search for a proof with more natural assumptions and farther reach in terms of types of SV models. Another important theoretical question is establishing a Central Limit Theorem for the coefficients as in \cite{Clement}, although it is not clear that our separation technique will carry over.

There is also the question of developing duality-based high biased estimators as in \cite{duality2} and \cite{duality1}. Such estimates would allow us to form proper confidence intervals for our time-zero prices.

Our future work will revolve around applying the algorithm to other SV models and sharpening of our theoretical results. There is also the issue we raised before regarding a complete complexity analysis and strategies for optimization. In terms of computing sensitivities, as mentioned before, the details may be found in \cite{Farahany}. There, one can also find applications to the double Heston model (\cite{dHeston}) and a mean-reverting commodity model with jumps.

\clearpage

\section{Bibliography}

\bibliographystyle{chicago}

redpajama_set_name: RedPajamaArXiv