Finite Difference Solution Ansatz approach in Least-Squares Monte Carlo

We assume that the American-style option is exercisable after (excluding) $t=0$ and up to (including) the final expiry $T$ for every time interval $\Delta T$. We define the set of exercise times on or after $t$ as $\pi(t):=\{t_{k}\geq t|t_{k}=k\Delta T,1\leq k\leq M=\frac{T}{\Delta T}\}$. The exercisability is either on the option holder side or on the issuer side. Within the complete probability space $\Omega$, there exists a discrete filtration $\{\operatorname{\mathcal{F}}(t_{k})|1\leq k\leq M\}$ on $\pi(0)$. For notational clarity, we use $XYZ(\omega,t)$ to denote a stochastic variable under $\operatorname{\mathcal{F}}(t)$ with $\omega\in\Omega$. In this spirit, the state variables are denoted by $X(\omega,t_{k})$, which is $\operatorname{\mathcal{F}}(t_{k})$-measurable.

To unambiguously describe an American-style option, we have to define two important quantities: $Z(\omega,t_{k})\equiv Z(X(\omega,t_{k}))$ as the exercise payoff at $t_{k}$, and $\tilde{V}(\omega,\tilde{\tau}(t))$ as the sum of discounted (back to $t$) future cash flows from (exclusively) $t$ and up to a stopping time $\tilde{\tau}(t)\in\pi(t)$. Different from $X(\omega,t_{k})$ or $Z(\omega,t_{k})$, $\tilde{V}(\omega,\tilde{\tau}(t))$ is measurable up to $\tilde{\tau}(t)$. For single cash flow derivative, these two quantities are related to each other as $\tilde{V}(\omega,\tilde{\tau}(t))=e^{-r(t_{k}-t)}\tilde{V}(\omega,\tilde{\tau}(t_{k}))=e^{-r[\tilde{\tau}(t_{k})-t]}Z(X(\omega,\tilde{\tau}(t_{k})))$ for $t_{k-1}<t\leq t_{k}$.

Now one can define the expected option value, i.e., present value (PV), as

Similar to previous literature, we limit the attention to square integrable payoff functions, i.e., $\operatorname{\mathbb{E}}[\tilde{V}^{2}(\omega,\cdot)]<\infty$. Let $L^{2}(\Omega,\operatorname{\mathcal{F}},(\operatorname{\mathcal{F}}_{t})_{0\leq t\leq T},\operatorname*{\mathbb{Q}})$ denote the space of square integrable functions. In this article, we take Bermudan options on constant-weighted basket $S_{B}(t)=\sum_{i=1}^{d}w_{i}S_{i}(t)$ with $w_{i}=\frac{1}{d}$ (see Appendix A.1), and worst of issuer callable notes (WIC) on worst-of basket $S_{W}(t)=\min_{i\in\{1,\ldots,d\}}S_{i}(t)$ (see Appendix A.2), respectively. Here the underlying basket is $d$ dimensional and $\vec{S}(t)=(S_{1}(t),\ldots,S_{d}(t))^{T}$ are asset spot prices.

Under discrete-time formulation, a preferable way to solve the Eqn.2.1 is to find out the optimal stopping time strategy, denoted by $\tau(t_{k}):=\tau(\omega,t_{k})$ (removing the tilde compared to $\tilde{\tau}(t_{k})$ for optimality and omitting the dependence on $\omega$ for notational convenience). For further notational clarity, we also introduce $V(\omega,\tau(t_{k}))$ to represent the option path value under optimal stopping times $\tau(t_{k})$, differentiated from $\tilde{V}(\omega,\tilde{\tau}(t))$ by removing the tilde similarly to emphasize that it is calculated under optimal stopping times $\{\tau(t_{k})\}$.

With the help of $V(\omega,\tau(t_{k}))$, we further define the so-called continuation value representing the conditional expected option value upon exercising strictly after, and discounted back to, $t_{k}$ as

Thanks to the $\operatorname{\mathcal{F}}_{t_{k}}$-measurability of $F(\omega,t_{k})$ in the same way as $Z(\omega,t_{k})$, the formal solution of the optimal stopping times can be written as, via a dynamic programming fashion (together with terminal condition $\tau(t_{M})=T$),

where the exercise indicator function is defined as $I(\omega,t_{k}):=\Theta\left(\pm[Z(\omega,t_{k})-F(\omega,t_{k})]\right)$, with Heaviside step function denoted by $\Theta(x)$ and ‘$+$’(‘$-$’) for holder (issuer) exercise. The formal dynamic programming solution of $V(\omega,\tau(t_{k}))$ is also given by

For the ease of illustration without loss generality, we have assumed that the product doesn’t have early termination features, eg. European Option on $S_{B}(t)$. We define the Credit value adjustment (CVA) on a financial transaction as the expected loss resulting from the potential future default of the counterparty, given by:

where $R_{ctpy}$ is the recovery rate and $PD(t)$ denotes the default probability of the counterparty at $t$, which depends on the counterparty hazard rate denoted by $h_{cpty}(t)$. Furthermore, $w(t):=\operatorname{\mathbb{E}}[\tilde{V}(\omega,T)|{\operatorname{\mathcal{F}}}_{t}]$ is the mark-to-market value of the option at $t$.

We define the Expected Positive Exposure (EPE), as well as its discounted version EPE*, as:

Here we re-use the same notation $\tilde{V}(\omega,T)$ defined in previous subsection, to represent the option value before expiry $T$. Efficient computation of $w(t)$ is of prime importance in the calculation of CVA and other XVA’s, where the difficulty in practice arises from the lack of analytical tractability in many cases.

To further take into account the Wrong Way Risk (WWR), we consider a hazard rate with a general functional dependence on $w(t)$(Hull & White, (2012)):

where $a(t)$ and $b$ are WWR model parameters.

As can be seen from above, once the $\operatorname{\mathcal{F}}_{t_{k}}$-measurable conditional expectations in Eqn.2.2 and Eqn.2.6 can be efficiently crystallized into a functional form $F(X(\omega,t_{k}))$, the optimal stopping times in American style option pricing and EPE/CVA calculation, can be determined in a path-by-path manner independently. The core idea of regression based Least-Squares approach (LSM) is to assume this functional form can be decomposed as a linear combination of a set of $R$ basis functions $\vec{\phi}_{R}(x)=\{\phi_{r}(x)|r=1,\ldots,R\}$, where the linear coefficients are obtained via regression using the cross-sectional information computed in the simulation.

For notation convenience, let’s define the least-squares estimator of the continuation value $F_{k}^{P}$ using the $L^{2}$ projection operator $P_{k}$ at $t_{k}$ upon option value $V$ under basis functions $\vec{\phi}_{R}(x)$ in $\operatorname*{\mathbb{R}}^{R}$ space, as

where

It means that $\vec{\beta}_{k}=(\beta_{k,1},\ldots,\beta_{k,R})^{T}$ are to be obtained via regression using the cross-sectional information computed in the simulation as described in Appendix A.3. $\mathcal{L}^{P}(\vec{b};V(\omega,\tau(t_{k+1})))$ in above is the objective function to be minimized, whose lower bound is given by, in the limit of $\vec{b}\cdot\vec{\phi}_{R}(X(\omega,t_{k}))$ being exact:

where $\operatorname{\mathrm{Var}}[\cdot|\operatorname{\mathcal{F}}_{t_{k}}]$ denotes the conditional variance.

In Eqn.2.7, the least-squares estimator $P_{k}[V(\omega,\tau(t_{k+1}))]$ is also called regressed continuation value.

After solving $\vec{\beta}_{k}$ for all $t_{k}$ recursively backward in time from maturity, $F_{k}^{P}(X(\omega,t_{k}))$ are fully determined. Then with a forward pricing Monte Carlo, one can then compute PV as

where the superscript $(i)$ indicates that the variable is under $i^{th}$ of the total $N_{P}$ pricing paths; and the optimal stopping rule is to only exercise the option at

To have a better convergence and calculation efficiency, one usually draws a different set of Monte Carlo paths in this pricing stage more than the previous regression stage to determine the continuation value, i.e., $N_{P}>N_{R}$. This is so-called out-of-sample implementation to eliminate foresight bias without the need of sophisticated regression correctors as mentioned before.

As an analogy of the EPE in the CVA problem is calculated as (by replacing the value function with least-squares estimator)¹¹1Note that we are aware of alternative hybrid approach to calculate EPE, by making use of both the simulated cash flows and regressed function (Joshi & Kwon, (2016)). The drawback of this proposal is that one has to generate all the cash flows at given individual product level, without being decoupled from the path generation at upper portfolio level. It would result in an aggregation of contractual timelines, rather than simple algorithmic operation upon regressed functions at product level when it comes to portfolio aggregation. :

Similarly, the Monte Carlo estimator of CVA can be defined accordingly as:

Alternative approach to solve the American-style option pricing is the finite difference (FD) method. The elegance of FD method is that it marches from maturity to present, which naturally fits into the requirement of backward recursive optimization in stopping time problems. We restrict the discussion within low-dimensional only, i.e., 1D, for maximum efficiency and accuracy. The most generic model reliable in replicating the market characteristics and solvable under 1D FD PDE is the local volatility, under which the option price $V^{FD}(S,t)$ as function of underlying spot $S$ and time $t$ satisfies:

within the time domain $(t_{k},t_{k+1}]$, and subjected to the initial condition $V^{FD}(S,T)=Z(S,T)$ and transition conditions around $t_{k}$ due to obstacle problem as

for holder’s exercise; and

for issuer’s exercise, respectively. As we can see later in Sec.2.6, $q_{FD}(S,t)$ and $\sigma_{FD}(S,t)$ in Eqn.2.13 corresponds to the continuous dividend yield and local volatility for a single asset, respectively.

After solving the 1D backward propagating PDE via Crank-Nicolson method, we will denote by, from now on, $f^{FD}_{k}(S)$ the FD solution of continuation value at $t_{k}$, which is related to backward PDE solution as

The proof of Prop.1 is omitted as it is an elementary consequence of conditional expectation as $L^{2}$ projection that minimizes the objective function of LSM. This conclusion is free from model setting assumption.

For a general problem of interest, one could find an auxiliary variable with existing exact (or quasi-exact) solution as an ansatz to start with as one of the basis functions in the $L^{2}$ projection, with the expectation of dominating regression weight on this ansatz in the spirit of Prop.1. Inspired by the simplicity of FD solution for 1D PDE in Sec.2.4 with efficiency and quasi-exact accuracy, we would intuitively take the option price under 1D process as auxiliary variable, whose FD solution of continuation value $f_{k}^{\text{FD}}(x)$ is given by Eqn.2.14. From now on, we call this new method by incorporating $f_{k}^{\text{FD}}(x)$ into one of the LSM basis functions as FD Solution Ansatz based Least-Squares Monte Carlo (FD-LSM).

Mathematically, FD-LSM is defined by a modified least-squares continuation estimator with a new $L^{2}$ projection operator $P_{k}^{\text{FD-LSM}}$, as opposed to the classic one $P_{k}^{\text{LSM}}$(previously denoted by $P$ defined in Eqn.2.7) operating on $\vec{\phi}_{R}^{\text{LSM}}(x)$(previously denoted by $\vec{\phi}_{R}(x)$). The new basis functions of $P_{k}^{\text{FD-LSM}}$ is given by $\vec{\phi}_{R}^{\text{FD-LSM}}(x):=\vec{\phi}^{\text{FD}}(x)\cup\vec{\phi}_{R-1}^{\text{LSM}}(x)$ as a union of $\vec{\phi}^{\text{FD}}(x)$ and $\vec{\phi}_{R-1}^{\text{LSM}}(x)$ ($R$ is decremented by one to cater for the added new FD ansatz), where $\vec{\phi}^{\text{FD}(x)}=\{1,f_{k}^{\text{FD}}(x)\}$³³3We allow a constant term in the regressors to cater for mismatch of center of mass between independent and dependent variables. under which the associated projection operator $P_{k}^{\text{FD}}$ is also introduced for illustrative purpose later.

Below proposition is established for the illustration of fundamental improvement of $P_{k}^{\text{FD-LSM}}$ over $P_{k}^{\text{LSM}}$ from a theoretical point of view.

Based on Prop.2, proof of convergence on LSM from previous studies (Longstaff & Schwartz, (2001); Stentoft, (2001)) is equally applicable to FD-LSM. As $R\rightarrow\inf$, both $P_{k}^{\text{FD-LSM}}$ and $P_{k}^{\text{LSM}}$ converge to the exact continuation function and are thus unbiased estimators. However, the variance of projection kernel in the former is reduced tremendously in a control variated manner by observing

with

defining the correlation between $V_{k+1}$ and $\hat{V}_{k+1}^{1D}$.

Considering Eqn.2.9 together with above, one can see that least-squares error (lower bound) in FD-LSM is also reduced by the same scaling factor as the variance:

There are several further remarks on this result:

• •

  It is well known that control variates are very effective Monte Carlo variance reduction techniques, see e.g. Chapter 4 in Glasserman, (2004). In the context of American exercise option pricing, up til now, such toolkit has only been applied in the Monte Carlo pricing stage to reduce the final PV variance from literature (see Søndergaard Rasmussen, (2002); Tian & Burrage, (2003)). On the contrary, FD-LSM is designed as equivalently applying the variance reduction technique during the regression stage, which is crucial for optimal stopping strategy determination subsequently in pricing stage. Such improvement cannot be achieved under the usual payoff control variates method.

• •

  As mentioned before, $\hat{V}_{k+1}^{1D}$ can be seen as option value under a 1D simulation process. However, instead of employing the right hand side in Eqn.2.17 directly for least-squares error reduction, Prop.2 allows one to achieve exactly the same effect by just adding the FD ansatz into the projection basis functions. This helps mitigate computational cost as it saves additional tedious Monte Carlo simulation on $\hat{V}^{1D}$, as well as the variance and covariance calculation in Eqn.2.15.

• •

  From the proof of Prop.2, up to first order least-squares errors, we don’t see benefit of introducing cross term between $\vec{\phi}^{\text{FD}}(x)$ and $\vec{\phi}_{R-1}^{\text{LSM}}(x)$. This simplifies the construction of new basis and avoid overfitting due to further introduction of higher order terms.

• •

  The most effective FD auxiliary variable $\hat{V}^{1D}$ is obtained by achieving the largest possible correlation with $V$ based on Eqn.2.18. This provides a important theoretical foundation on the choice of FD ansatz in FD-LSM. For instance, under a multi-asset single-factor setting (e.g. Black Scholes or Local Volatility), we can choose the ansatz from single-factor 1D PDE after inter-asset dimension reduction; for single-asset multi-factor settings (e.g. Heston stochastic volatility), one can construct the equivalent 1D process (Local volatility in general) by matching same marginal distributions, see Gyöngy, (1986); finally, for the most general case of multi-asset multi-factor setting, one could apply the dimension reduction on intra-asset level (reducing the factors other than spot, e.g. stochastic volatility on the same asset) to start with, and then on inter-asset level, eventually arriving at 1D problem which is solvable under Eqn.2.13.

Error accumulation in a multi-period setting

Here we provide quantitative analysis on the error accumulation of under dynamic programming on LSM and FD-LSM. This piece of analysis is largely an extension of previous studies, with a focus on relating the accumulated option pricing errors to the approximation introduced in the regressed continuation value during the backward evolution.

From this point on, for notational simplicity, we drop the spatial dependence of random variables and condense notation by replacing arguments depending on time by subscripts and further replace $t_{k}$ with simply $k$. For example, $V_{k+1}=V(\omega,\tau(t_{k+1}))$, $Z_{k}=Z(\omega,t_{k})$ and $F_{k}=F(\omega,t_{k})$. We restrict the discussion in the context of holder exercisability, but the final conclusion can be applied to issuer exercisability too.

With the discounting factor from period to period $e^{-r\Delta T}$ dropped out for ease of illustration from now on(see Clément et al., (2002)), the dynamic programming solution of option value is re-written here from Eqn.2.4 as

In a similar way, the approximated solution under least-squares Monte Carlo defined by a given projection method $P$ (which can be either within LSM or FD-LSM) follows:

where $P_{k}$ can be either $P_{k}^{\text{LSM}}$ or $P_{k}^{\text{FD-LSM}}$, and $P_{k}[V^{P}_{k+1}]$ is the projection on the option value with error accumulated up to the next step due to approximation from backward.

We are interested in the total expectation value of the difference between $V_{k}^{P}$ and $V_{k}$, and it can be calculated as:

where $\operatorname{\mathbb{E}}(\operatorname{\mathbb{E}}(V_{k+1}^{P}|\operatorname{\mathcal{F}}_{t_{k}}))=\operatorname{\mathbb{E}}(\operatorname{\mathbb{E}}(V_{k+1}^{P}|\operatorname{\mathcal{F}}_{t_{k+1}}))$ and $\operatorname{\mathbb{E}}(\operatorname{\mathbb{E}}(V_{k+1}|\operatorname{\mathcal{F}}_{t_{k}}))=\operatorname{\mathbb{E}}(\operatorname{\mathbb{E}}(V_{k+1}|\operatorname{\mathcal{F}}_{t_{k+1}}))$ have been applied starting with law of total expectation.

We can see easily that in above both terms are non-positive by induction starting from $\operatorname{\mathbb{E}}(V_{M}^{P})=\operatorname{\mathbb{E}}(V_{M})$. Intuitively, the reason is that the least-squares approximation introduces a sub-optimal exercise strategy giving rise to biased low error for holder exercisable option price.

Denoting $e^{P}_{k}=|\operatorname{\mathbb{E}}(V_{k}^{P}-V_{k})|$, we have

where $\mathcal{Y}^{P}_{k}=\operatorname{\mathbb{E}}(\Theta(P_{k}[V^{P}_{k+1}]-Z_{k}))$. For $|\mathcal{B}^{P}_{k}|:=\operatorname{\mathbb{E}}((Z_{k}-F_{k})[\Theta(Z_{k}-P_{k}[V^{P}_{k+1}])-\Theta(Z_{k}-F_{k})])$, we have, by applying Jensen’s inequality,

where we have introduced $\mathcal{X}^{P}_{k}:=\operatorname{\mathbb{E}}(\Theta(|P_{k}[V^{P}_{k+1}]-F_{k}|-|F_{k}-Z_{k}|))$ in the penultimate line, which in practice should be small as $P_{k}[V^{P}_{k+1}]$ is close to $F_{k}$.

Applying basic property of projection to the first term in above (see similar derivation in the proof of Theorem 3.1 in Clément et al., (2002)), we have

where in the last line, the last term is the local regression least-squares error on the projection $P_{k}[V_{k+1}]$, which doesn’t consider accumulated error yet (as opposed to $P_{k}[V^{P}_{k+1}]$). Denoting this term by $\Xi_{k}^{P}[V_{k+1}]$, we arrive at

where we have assumed during error accumulation, $\operatorname{\mathrm{Var}}[V^{P}_{k+1}-V_{k+1}]\ll(e^{P}_{k+1})^{2}$ and thus the former variance term is neglected.

Combining with Eqn.2.22, we have the dynamical programming solution on the upper bound of $e^{P}_{k}$, denoted by $\epsilon^{P}_{k}$, for a given projection $P$, to be

with initial condition $\epsilon^{P}_{M}=0$. This defines the relation between accumulated errors and local regression errors on option prices under least-squares Monte Carlo at each step, although looking simple enough but no closed form solution for further simplification.

To proceed further, we have to make some assumption on

under fixed finite number of cutoff and regression paths within LSM⁴⁴4The asymptotic convergence rate of this term going to zero under big O notation has been analyzed by Stentoft, (2001) in Sec. 3.2.1.. Let’s assume that, $\Xi_{k}^{P^{\text{LSM}}}[V_{k+1}]\sim\alpha\operatorname{\mathbb{E}}[V_{k+1}]$ with a constant $\alpha$, which reflects the consumption of problem-agnostic basis functions. Without loss of generality, let’s further assume $|\operatorname{\mathbb{E}}[F_{k}]|\ll 1$; in practice, this corresponds to a swap form product where the continuation value is closed to zero. Then we have $\Xi_{k}^{P^{\text{LSM}}}[V_{k+1}]\sim\alpha\operatorname{\mathbb{E}}(\operatorname{\mathrm{Var}}[V_{k+1}|\operatorname{\mathcal{F}}_{t_{k}}])$.

According to Eqn. 2.17 and Eqn.2.18, we have

As we can see here, we expect that $\Xi_{k}^{P^{\text{FD-LSM}}}[V_{k+1}]$ can be much smaller than $\Xi_{k}^{P^{\text{LSM}}}[V_{k+1}]$, depending on the correlation between the original and FD auxiliary problems. Let’s assume a high correlation $\tilde{\rho}(V_{k+1},\hat{V}^{1D}_{k+1})=90\%$; and if $\Xi_{k}^{P^{\text{LSM}}}[V_{k+1}]=0.05^{2}$ for all $k$, then $\Xi_{k}^{P^{\text{FD-LSM}}}[V_{k+1}]=(0.05\times\sqrt{1-90\%^{2}})^{2}=0.0218^{2}$. To have an illustrative example, let’s further put $\mathcal{X}^{P}_{k}=0.01$ and $\mathcal{Y}^{P}_{k}=0.8$ for all $k$ in both LSM and FD-LSM as an example, then we have the plots for (upper bound) accumulated option price errors $\epsilon^{P}_{k}$ going backward at each time step in Fig.2.1.