AAD and least-square Monte Carlo: Fast Bermudan-style options and XVA Greeks

Abstract

We show how Adjoint Algorithmic Differentiation (AAD) can be used to calculate price sensitivities in regression-based Monte Carlo methods reliably and orders of magnitude faster than with standard finite-difference approaches. We present the AAD version of the celebrated least-square algorithms of Tsitsiklis and Van Roy (2001) and Longstaff and Schwartz (2001). By discussing in detail examples of practical relevance, we demonstrate how accounting for the contributions associated with the regression functions is crucial to obtain accurate estimates of the Greeks, especially in XVA applications.

Keywords

Adjoint algorithmic differentiation (AAD)Monte Carlo methods Bermudan-style options valuation adjustments (XVA)

1 Introduction

The efficient calculation of the risk factor sensitivities of financial derivatives, also known as the “Greeks”, is an essential component of modern risk management practices. Indeed, the aftermath of the recent financial crisis has seen remarkable changes in the market practice whereby financial institutions need to quantify (and risk-manage) counterparty, funding and capital exposures, collectively known as XVA, in large portfolios, see e.g., Crépey et al. (2014).

The traditional approach for the calculation of the Greeks is the so-called bump and reval or bumping technique. This comes with a significant computational cost as it generally requires repeating the calculation of the P&L of a portfolio under hundreds of market scenarios in order to form finite-difference estimators. As a result, in many cases, even after deploying vast amounts of computer power, these calculations cannot be completed in a practical amount of time.

Conversely, Adjoint Algorithmic Differentiation (AAD), a numerical technique recently introduced to financial engineering (see e.g., Capriotti (2011); Capriotti et al. (2011); Capriotti and Giles (2012, 2010); Henrard (2014)), has proven to be effective for speeding up the calculation of risk factor sensitivities, both for Monte Carlo (MC) and deterministic numerical methods, see Capriotti and Lee (2014); Savickas et al. (2014); Capriotti et al. (2015); Xu et al. (2016); Geeraert et al. (2017).

The main ideas underlying AAD can be illustrated by considering a computer implemented function of the form

$Y = FUNCTION (X)$ (1) mapping a vector $X \in ℝ^{n}$ to a vector $Y \in ℝ^{m}$ through a sequence of intermediate steps $X \to \dots \to U \to V \to \dots \to Y .$ (2)

Here, the real-valued vectors U and V represent intermediate variables utilised in the calculation. Each step may be a distinct high-level function or even a specific instruction.

AAD, sometimes also simply known as adjoint mode of Algorithmic Differentiation (AD), results from “propagating” the derivatives of the final output with respect to all the intermediate variables – the so called adjoints – until the derivatives with respect to the independent variables are formed. Using the standard AD notation, the adjoint of any intermediate variable V_k is defined by ${\bar{V}}_{k} = \sum_{j = 1}^{m} {\bar{Y}}_{j} \frac{\partial Y_{j}}{\partial V_{k}},$ (3) where $\bar{Y}$ is a vector in $ℝ^{m}$ . For each of the intermediate variables U_i, by applying the chain rule, we get ${\bar{U}}_{i} = \sum_{j = 1}^{m} {\bar{Y}}_{j} \frac{\partial Y_{j}}{\partial U_{i}} = \sum_{j = 1}^{m} {\bar{Y}}_{j} \sum_{k} \frac{\partial Y_{j}}{\partial V_{k}} \frac{\partial V_{k}}{\partial U_{i}},$ which corresponds to the adjoint mode equation for the intermediate step represented by the function V = V (U). We thus have a function of the form $\bar{U} = \bar{V} (U, \bar{V})$ where ${\bar{U}}_{i} = \sum_{k} {\bar{V}}_{k} \frac{\partial V_{k}}{\partial U_{i}} .$

Starting from the adjoint of the outputs $\bar{Y}$ , we may apply this rule to each step in the calculation, working from the right to the left, $\bar{X} \leftarrow \dots \leftarrow \bar{U} \leftarrow \bar{V} \leftarrow \dots \leftarrow \bar{Y}$ (4) until we obtain $\bar{X}$ , namely, the linear combination of the rows of the Jacobian of the function X → Y: ${\bar{X}}_{i} = \sum_{j = 1}^{m} {\bar{Y}}_{j} \frac{\partial Y_{j}}{\partial X_{i}}$ (5) for i = 1, …, n.

One particularly important theoretical result is that given a computer program performing some high-level function (1), the execution time of its adjoint counterpart $\bar{X} = FUNCTION_b (X, \bar{Y})$ (6) (with suffix _b for “backward” or “bar”) that computes the linear combination (5), is bounded by three to four times the cost of execution of the original one. That is, $\frac{Cost [FUNCTION_b]}{Cost [FUNCTION]} \leq ω_{A}$ (7) where ω_A ∈ [3, 4], see Griewank and Walther (2008).

In this paper we present the application of AAD to regression-based MC approaches (also known as least-square MC) such as those that are widely utilised for Bermudan-style options, see Carriere (1996); Tsitsiklis and Van Roy (2001); Longstaff and Schwartz (2001), or for XVA applications, see Cesari et al. (2009) and Joshi and Kwon (2016). We develop the AAD implementation of the well-known least-square algorithm for the computation of conditional expectations, and we investigate numerically the impact on the Greeks arising from the sensitivities of the regression functions, a component that is generally ignored for Bermudan-style options by invoking arguments of quasi-optimality of the exercise boundary.

The paper is organised as it follows: In the next section, the regression-based MC algorithm for both Bermudan-style options and XVA is presented. In Section 3 we discusses the AAD algorithms for the regression-based MC method. We give two numerical examples, the best of two stocks Bermudan-style call and its corresponding XVA in Section 4. Here we show how smoothening out discontinuities associated with suboptimal exercise boundaries improves the accuracy of the Greeks of Bermudan-style options, and why the contribution to the sensitivities arising from the regression boundaries is essential for an accurate computation of XVA sensitivities. The efficiency and accuracy of AAD is also compared with bump and reval approaches in the same section.

2 Valuation of Bermudan-style options and XVA by regression-based Monte Carlo

2.1 Bermudan-style options

While European-style options can be exercised only at final maturity, Bermudan-style options can be exercised on multiple dates up to trade expiry.

We denote by T₁, …, T_M the exercise dates of the option and define $D (t) = {T_{m} \geq t}$ . We denote by η (t) the smallest integer such that T_η(t)+1 > t. An exercise policy is represented mathematically by a stopping time taking values in $D (t)$ . We denote by $T (t)$ the set of stopping times taking values in $D (t)$ , see e.g., Andersen and Piterbarg (2010).

A rational investor will exercise the option that he holds in such a way as to maximise its economic value. As a result, the value of a Bermudan-style option is the supremum of the option value over all possible exercise policies. With the notation introduced above, the value of a Bermudan-style option at time t can thus be expressed by $\frac{V (t)}{N (t)} = sup_{τ \in T (t)} 𝔼_{t} [\frac{E (τ)}{N (τ)}],$ (8) where E (t) is the exercise value of the option, and N (t) is the chosen numéraire 1 , Andersen and Piterbarg (2010). In this equation, V (t) is to be interpreted as the value of the option conditional on exercise not having taken place strictly before time t. Here and below we indicate with $𝔼_{t}$ the expectation, under the risk-neutral measure, conditional on the filtration of the random process driving the option price, up to time t. For notations and standard mathematical assumptions of this conventional setting see, e.g., Wilmott et al. (1995).

A useful concept is the hold value of the Bermudan-style option: We denote by H (t) the value of the Bermudan-style option when the exercise dates are restricted to $D (T_{η (t) + 1})$ , that is $\frac{H (t)}{N (t)} = 𝔼_{t} [\frac{V (T_{η (t) + 1})}{N (T_{η (t) + 1})}],$ (9) where we have assumed, for simplicity of exposition, no cashflow between t and T_η(t)+1.

The option holder, following an optimal exercise policy, will exercise his option if the exercise value is larger than the hold value, i.e., $V (T_{η (t)}) = max (E (T_{η (t)}), H (T_{η (t)})) .$ (10)

This, when combined with Equation (9), leads to the so-called dynamic programming formulation: $\frac{H (t)}{N (t)} = 𝔼_{t} [max (\frac{E (T_{η (t) + 1})}{N (T_{η (t) + 1})}, \frac{H (T_{η (t) + 1})}{N (T_{η (t) + 1})})],$ (11) for T_η ≤ t < T_η+1, and η = 1, …, M - 1. Starting from the terminal condition H (T_M) ≡0, Equation (11) defines a backward iteration in time for H (t). By definition, this is also equal to V (t) if t is not an exercise date, i.e., if T_η(t) < t < T_η(t)+1. Conversely, if t is an exercise date, t = T_η(t), then V (T_η(t)) = max(E (T_η(t)) , H (T_η(t))).

The dynamic programming formulation above implies that the stopping time, defining the optional exercise date as seen at time t, is given by $τ^{★} = inf [T_{m} \geq t : E (T_{m}) \geq H (T_{m})] .$ (12)

The optimal exercise strategy defined by Equation (12) requires the computation of the hold value H (t), m = η (t) + 1, …, M - 1. In a setting in which the underlying risk factor process {X (t) } _0≤t≤T is a generic k-dimensional Markov process, the hold value H (t) is a function of the state vector at time t. That is, $H_{t} (x) : = 𝔼 [\frac{N (X (t))}{N (X (T_{m + 1}))} V (X (T_{m + 1})) | X (t) = x] .$ (13)

When the dimension of the Markov process k is small enough, the conditional expectation value in Equation (13) can be computed in a straightforward way by discretising the risk-factor process and performing standard backward induction on a tree or a grid, or by discretising an associated Partial Differential Equation (PDE). Here we refer to, e.g., Wilmott et al. (1995). However, the complexity of grid-based calculations is exponential in the dimension of the Markov process and numerical implementations become infeasible when k ≥ 4.

As we will review in Section 2.3, regression-based MC techniques provide an effective way of computing conditional expectation values of the form (13).

2.2 XVA

We next consider the computation of the Credit Valuation Adjustment (CVA) and the Debt Valuation Adjustment (DVA) as the main measures of a dealer’s counterparty credit risk, see e.g., Cr'epey et al. (2014). For a given portfolio of trades with the same investor or institution, the CVA (resp. DVA) aims to capture the expected loss (resp. gain) associated with the counterparty (resp. dealer) defaulting in a situation in which the position, netted for any collateral posted, has a positive mark-to-market for the dealer (resp. counterparty).

This can be evaluated at time T₀ = 0 by

$\begin{matrix} XVA & = - & 𝔼 [𝕀 (τ_{c} \leq T) \frac{L_{c}}{N (τ_{c})} {(V (τ_{c}))}^{+} \\ + 𝕀 (τ_{d} \leq T) \frac{L_{d}}{N (τ_{d})} {(V (τ_{d}))}^{-}], \end{matrix}$ (14) where τ_c (resp. τ_d) is the default time of the counterparty (resp. the dealer), V (t) is the net present value of the portfolio or netting set at time t from the dealer’s point of view (the so-called conditional future exposure), L_c (resp. L_d) is the loss given default of the counterparty (resp. the dealer), and $𝕀 (τ_{c} \leq T)$ (resp. $𝕀 (τ_{d} \leq T)$ ) is the indicator that the counterparty’s (resp. dealer’s) default happens before the longest deal maturity T in the portfolio. Here, for simplicity of notation, we consider the unilateral CVA and DVA, the generalization to the bilateral formulation, see e.g., Crépey et al. (2014), being straightforward.

Equation (14) is typically computed on a discrete time grid of “horizon dates” 0 = T₀ < T₁ < … < T_M = T. For instance, we may have

$\begin{matrix} XVA \\ ≃ - \sum_{m = 1}^{M} 𝔼 [L_{c} ({SP}_{c} (T_{m - 1}) - {SP}_{c} (T_{m})) \frac{{(V (T_{m}))}^{+}}{N (T_{m})} \\ + L_{d} ({SP}_{d} (T_{m - 1}) - {SP}_{d} (T_{m})) \frac{{(V (T_{m}))}^{-}}{N (T_{m})}], \end{matrix}$ (15) where SP_c (t) (resp. SP_d (t)) is the survival probability of the counterparty (resp. the dealer) up to time t, e.g., conditional on a realization of the default intensity (or hazard rate) process in a Cox framework, see Lando (1998). Here we assume that the default times τ_c and τ_d are independent of the portfolio values V_{τ
_c} and V_{τ
_d}, respectively. In general, the right hand side of Equation (15) depends on several correlated random market factors, including interest rate, recovery amounts, and all the market factors the net conditional future exposure of the portfolio, V (t), depends on. As such, its calculation typically requires a MC simulation.

In the k-dimensional Markov setting introduced above the conditional future exposure V (T_m) is a function V (X (T_m)) of the state vector at time T_m. However, only for vanilla securities and simple models for the evolution of the risk factors, such conditional future exposures can be expressed in closed form, and regression based Monte Carlo is commonly employed to produce approximate estimators see, e.g., Cesari et al. (2009).

2.3 Conditional expectation values and Bermudan-style options by regression

Regression methods are based on the observation, see e.g., Friedman et al. (2001), that given a real-valued random input vector $X \in ℝ^{d}$ and Y a real valued random output, the conditional expectation $𝔼 [Y | X]$ is the function of X that best approximates in the least-square sense the output Y. That is, $𝔼 [Y | X] = \underset{c}{argmin} 𝔼 [(Y - c)^{2}] .$ (16)

In particular, assuming that the conditional expectation is a linear function of some unknown vector of parameters β, $𝔼 [Y | X] = β^{T} X,$ (17)

Equation (16) reduces to the well-known linear regression conditions, giving for the optimal vector of parameters $β = 𝔼 {[X X^{T}]}^{- 1} 𝔼 [XY] .$ (18)

In the context of the valuation of Bermudan-style options, the hold value (13) on an exercise date T_m is assumed to be of the form ${\hat{H}}_{m} (x) : = {\hat{H}}_{T_{m}} (x) = β_{m}^{T} ψ (x)$ (19) where ψ (x) = (ψ₁ (x) , …, ψ_d (x)) ^T is a vector of d basis functions and β_m = (β_1m, …, β_dm) ^T is the vector of coefficients to be determined by regressing N (X (T_m))/N (X (T_m+1)) V (X (T_m+1)) versus ψ (X (T_m+1)). This gives $β_{m} = Ψ_{m}^{- 1} Ω_{m},$ (20) where we define the d × d matrix $Ψ_{m} = 𝔼 [ψ (X (T_{m})) ψ^{T} (X (T_{m}))]$ (21) and the d × 1 vector $Ω_{m} = 𝔼 [\frac{N (X (T_{m})) V (X (T_{m + 1}))}{N (X (T_{m + 1}))} ψ (X (T_{m}))] .$ (22)

These equations provide a straightforward recipe to compute the regression coefficients β_m by substituting Ψ_m and Ω_m with their sample average over N_MC MC replications, namely:

Simulate N_MC independent MC paths $X_{m}^{(n)}$ of X (T_m) by the recursion $X_{m + 1}^{(n)} = F (T_{m}, X_{m}^{(n)}, θ)$ (23) for m = 0, …, M - 1, and n = 1, …, N_MC. Here F is a function based on the chosen models for the risk factors, and $θ \in ℝ^{N_{θ}}$ is a vector of N_θ model parameters.

For n = 1, …, N_MC, compute the terminal payoff of the contract by setting $V_{M}^{(n)} = E_{M}^{(n)},$ (24) where $E_{M} : = E (X_{M}^{(n)})$ is the final exercise value of the option.

Apply the following backward induction steps for m = M - 1, …, 1:

Compute the MC sample average of Ψ_m and Ω_m 2 by $Ψ_{m} = \frac{1}{N_{MC}} \sum_{n = 1}^{N_{MC}} ψ_{m}^{(n)} (ψ_{m}^{(n)})^{T},$ (25) $Ω_{m} = \frac{1}{N_{MC}} \sum_{n = 1}^{N_{MC}} ψ_{m}^{(n)} \frac{N_{m}^{(n)} V_{m + 1}^{(n)}}{N_{m + 1}^{(n)}},$ (26) where $ψ_{m}^{(n)} : = ψ (X_{m}^{(n)})$ and $N_{m}^{(n)} : = N (X_{m}^{(n)})$ .

Compute the regression coefficients β_m by matrix inversion and multiplication: $β_{m} = Ψ_{m}^{- 1} Ω_{m} .$ (27)

For the estimate of the hold value $H_{m}^{(n)} : = H_{m} (X_{m}^{(n)})$ , set $H_{m}^{(n)} = β_{m}^{T} ψ_{m}^{(n)},$ (28) for n = 1, …, N_MC.

For the estimate of the Bermudan-style option value at time T_m, set $V_{m}^{(n)} = max (E_{m}^{(n)}, H_{m}^{(n)}),$ (29) where $E_{m}^{(n)} : = E (X_{m}^{(n)})$ is the exercise value at time T_m, for n = 1, …, N_MC.

Compute the MC estimate of the Bermudan-style option at time T₀ by $V_{0} = \frac{1}{N_{MC}} \sum_{n = 1}^{N_{MC}} \frac{V_{1}^{(n)}}{N_{1}^{(n)}} .$ (30)

This approach was introduced by Tsitsiklis and Van Roy (2001) and they showed that the estimator V₀ converges for n→ ∞ to the true value V (0) provided that the representation (19) holds exactly.

A modification of this algorithm was proposed by Longstaff and Schwartz (2001) and it entails replacing Equation (29) in Step R3 (d) with $V_{m}^{(n)} = {\begin{matrix} E_{m}^{(n)} & if E_{m}^{(n)} > H_{m}^{(n)}, \\ N_{m}^{(n)} V_{m + 1}^{(n)} / N_{m + 1}^{(n)} & otherwise, \end{matrix}$ (31) which, in the examples considered, was shown to lead to more accurate results. In the following, however, for simplicity of exposition, we will consider the estimator in Equation (29).

2.4 Lower bound algorithm for Bermudan-style options

The hold value obtained by regression as described in the previous section defines an exercise policy whereby on each exercise date T_m the option is exercised if $E (X (T_{m})) > β_{m}^{T} ψ (X (T_{m})) .$ (32)

Such policy, being an approximation of the solution of the dynamic programming Equation (11), will in general correspond to a suboptimal stopping time. As a result, when utilised in a second, independent, MC simulation, the exercise policy obtained by regression, will result in a lower-bound estimator for the Bermudan-style option value. The corresponding algorithm can be schematically described as it follows.

For each MC replication indexed by n = 1, …, N_MC perform steps (L1) to (L4) below:

Simulate the path $X_{m}^{(n)}$ of the risk factor vector X (T_m) as in (R1).

For m = 1, 2, . . . , M - 1, compute the approximate hold value of the option at time T_m using the associated regression vector β_m, and regression functions ψ, by $H_{m}^{(n)} = β_{m}^{T} ψ_{m}^{(n)}$ (33) with the hold value at expiry T_M set to zero.

Compute the path-wise estimator for the discounted cash-flows of the option $P^{(n)} = \sum_{m = 1}^{M} [1^{(n)} (t_{1}, t_{m}) 1 (E_{m}^{(n)} > H_{m}^{(n)}) \frac{E_{m}^{(n)}}{N_{m}^{(n)}}],$ (34) where $1^{(n)} (t_{1}, t_{m}) = (\prod_{i = 1}^{m - 1} 1 (H_{i}^{(n)} > E_{i}^{(n)})),$ (35) for m > 1, and the convention 1⁽ⁿ⁾ (t₁, t₁) =1.

Compute the MC estimate of the Bermudan-style option at time T₀ = 0 by $V_{0} = \frac{1}{N_{MC}} \sum_{n = 1}^{N_{MC}} P^{(n)} .$ (36)

2.5 XVA by regression

As described in Section 2.2, the calculation of the XVA in Equation (15) requires the conditional future exposure V (t) on a set of dates determined by a discretisation time grid T₁ … T_M. The regression algorithm described in the previous section can be easily adapted to compute such quantity. Indeed, the conditional value of each of the options contained in the netting set can be obtained using the same least-square procedure.

Once the regression algorithm is completed, we can use the regression functions to compute the hold value of each option in the portfolio on the discretisation time grid by $H_{m} = β_{m}^{T} ψ_{m}$ . If the discretisation time T_m is not an exercise opportunity for the option under consideration, then this is also its conditional future exposure. Conversely, the conditional future exposure is obtained by comparing the hold value to the exercise value as in Equations (29) and (31). These observations translate in the following algorithm.

For each MC replication indexed by n = 1, …, N_MC perform steps (X1) to (X3) below:

Simulate the path $X_{m}^{(n)}$ of the risk factor vector by the recursion: $X_{m + 1}^{(n)} = F (T_{m}, X_{m}^{(n)}, θ),$ (37) for m = 0, …, M - 1. Simulate the path of the counterparty’s and the dealer’s default intensity, $λ_{m}^{d, c} = λ^{d, c} (T_{m})$ by the recursions $λ_{m + 1}^{c, (n)} = G_{c} (T_{m}, λ_{m}^{c, (n)}, θ),$ (38) $λ_{m + 1}^{d, (n)} = G_{d} (T_{m}, λ_{m}^{d, (n)}, θ),$ (39) for m = 0, …, M - 1, where G_c (resp. G_d) is the function describing the dynamics of the counterparty’s (resp. the dealer’s) hazard rate.

Compute the (discretised) path-wise survival probabilities for the counterparty and the dealer by ${SP}_{m}^{c, (n)} = exp [- \sum_{j = 0}^{m - 1} λ_{j}^{c, (n)} (T_{j + 1} - T_{j})],$ (40) ${SP}_{m}^{d, (n)} = exp [- \sum_{j = 0}^{m - 1} λ_{j}^{d, (n)} (T_{j + 1} - T_{j})],$ (41) for m = 1, 2, . . . , M.

For m = 1, 2, . . . , M - 1, approximate the hold value of the p-th option in the portfolio at time T_m using the associated regression vector β_p,m and regression functions ψ_p by $H_{p, m}^{(n)} = β_{p, m}^{T} ψ_{p, m},$ (42) with p = 1, …, P. The hold value at the expiry date T_M is set to zero. The conditional expectation value of the portfolio is given by $V_{m}^{(n)}$ = $\sum_{p = 1}^{P} V_{p, m}^{(n)}$ , where $V_{p, m}^{(n)} = {\begin{matrix} max {H_{p, m}^{(n)}, E_{p, m}^{(n)}}, & if T_{m} is an exercise date for the p - th option \\ H_{p, m}^{(n)}, & otherwise, \end{matrix}$ (43) for m = 1, 2, . . . , M, where $E_{p, m}^{(n)}$ is the exercise value of the p-th option at time T_m on the n-th path.

Compute the path-wise XVA by ${XVA}^{(n)} = - \sum_{m = 1}^{M} [L_{c} ({SP}_{m - 1}^{c, (n)} - {SP}_{m}^{c, (n)}) \frac{{(V_{m}^{(n)})}^{+}}{N_{m}^{(n)}} + L_{d} ({SP}_{m - 1}^{d, (n)} - {SP}_{m}^{d, (n)}) \frac{{(V_{m}^{(n)})}^{-}}{N_{m}^{(n)}}] .$ (44)

Form the MC estimator $XVA = \frac{1}{N_{MC}} \sum_{n = 1}^{N_{MC}} {XVA}^{(n)} .$ (45)

3 The AAD algorithm for regression-based Monte Carlo

3.1 Function regularizations

In the context of MC methods, Capriotti and Giles (2012) show that AAD allows to calculate the sensitivities by differentiating the relevant estimator on a path by path basis. As a path-wise method, the MC estimators must satisfy specific regularity conditions, see Glasserman (2003). For instance, all the functions appearing in each step leading to the computation of the payout estimator must be Lipschitz continuous. A practical way of addressing non-Lipschitz estimators is to smoothen out the singularities they contain. This can be achieved by observing that in most cases the singularities in the payout functions, although not necessarily implemented as such, can be expressed in terms of Heaviside functions. For instance, the payoff function of a digital option (for a one-dimensional underlying asset) is, $P (X (T)) = 1 (X (T) > K),$ (46) while the payoff of a knock-out, path-dependent option with barrier monitored at the time T₁, …, T_M is of the form $P (X (T_{1}), \dots, X (T_{M})) = \prod_{m = 1}^{M} 1 (X (T_{m}) > B_{m})$ (47) where B_m is the the barrier level at time T_m.

The singularities in such payoff functions can be regularized by replacing the indicator function with one of its smoothened counterparts. A very common choice, for instance, is to approximate the step function with a “call spread” payoff functions,

$\begin{matrix} 1 (x > K) \approx H_{δ} (x - K) \\ = {(min (\frac{x - (K - δ)}{2 δ}, 1))}^{+}, \end{matrix}$ (48) where 0 < δ ⪡ K. This is a standard choice for digital options, because it has a useful interpretation in terms of the hedging portfolio of a long and a short position in two calls with strike price K - δ and K + δ. This regularization gives rise to functions that are Lipschitz continuous with respect to x and can be differentiated in a straightforward manner. In particular, the adjoint regularized Heaviside function (48) is given by ${\bar{H}}_{δ} (x - K, \bar{x}) = {\begin{matrix} \bar{x} / 2 δ & if K - δ \leq x \leq K + δ \\ 0 & otherwise \end{matrix}$ (49) where $\bar{x}$ is the adjoint of the input variable. When δ → 0, (48) gives the correct derivative of the Heaviside function in the distributional sense, i.e., the Dirac delta function. However, while approaching this limit, its derivative is zero or vanishingly small, apart from a very small portion of the sample space where instead it is very large. This leads to exceedingly large variances in the MC sampling of the estimators expressed in terms of such functions, signalling the breakdown of the Lipschitz continuity condition. Hence, the choice of the smoothening parameter δ is necessarily a tradeoff between the bias, vanishing for δ → 0, and the statistical errors of the MC sampling, diverging in the same limit.

In general, the payoff estimator for Bermudan-style options (34) is not differentiable with respect to the pathwise value of the approximate exercise boundary $H_{m}^{(n)}$ , and it requires the regularization described above. A common approximation among practitioners, see e.g., Leclerc et al. (2009), is to assume that that exercise boundary implied by the rule (32) is close to optimality so that the value of the contract is approximatively continuous across the exercise boundary and no regularization is required. Under this assumption, no contribution to the sensitivities is associated with the perturbations of the exercise boundary, and one can therefore keep the regression coefficients fixed while calculating the sensitivities, see Andersen and Piterbarg (2010). As discussed in Section 4, depending on the accuracy of the basis functions in representing the exercise boundary, this may or may not be an accurate approximation.

3.2 AAD for the lower bound algorithm for Bermudan-style options

The AAD implementation of the lower bound algorithm for Bermudan-style options described in Section 2.4, producing the MC estimators for the sensitivities of the estimator (36) with respect to a set of model parameters θ_k, k = 1, …, N_θ, given by ${\bar{θ}}_{k} = \frac{\partial V (θ)}{\partial θ_{k}},$ (50) consists of the adjoint of steps (L1) to (L4) executed backwards for each MC replication n = 1, …, N_MC. After setting the adjoint of the option value ${\bar{V}}_{0} = 1$ , the adjoint of the model parameters $\bar{θ} = 0$ , and the adjoint of the regression coefficients ${\bar{β}}_{m} = 0$ , the steps to perform are:

Set ${\bar{P}}^{(n)} = {\bar{V}}_{0} \frac{1}{N_{MC}} .$ (51)

Assuming the indicator functions in the estimator in (34) have been regularized as discussed in Section 3.1, compute $\begin{matrix} {\bar{H}}_{m}^{(n)} & = & {\bar{H}}_{δ} (H_{m}^{(n)} - E_{m}^{(n)}, {\bar{P}}^{(n)}) (J_{m}^{(n)} - Q_{m}^{(n)}), \\ {\bar{E}}_{m}^{(n)} & = & {\bar{P}}^{(n)} R_{m}^{(n)} + {\bar{H}}_{δ} (E_{m}^{(n)} - H_{m}^{(n)}, {\bar{P}}^{(n)}) \\ (J_{m}^{(n)} - Q_{m}^{(n)}), \end{matrix}$ (52) ${\bar{N}}_{m}^{(n)} = - {\bar{P}}^{(n)} \frac{Q_{m}^{(n)} H_{δ} (E_{m}^{(n)} - H_{m}^{(n)})}{N_{m}^{(n)}}$ for m = M, …, 1 and where $\begin{matrix} J_{m}^{(n)} = \sum_{k = m + 1}^{M} \frac{E_{k}^{(n)}}{N_{k}^{(n)}} 1 (E_{k}^{(n)} > H_{k}^{(n)}) \\ \prod_{j = 2, j \neq m}^{k - 1} 1 (H_{j}^{(n)} > E_{j}^{(n)}), \\ Q_{m}^{(n)} = \frac{E_{m}^{(n)}}{N_{m}^{(n)}} 1^{(n)} (t_{1}, t_{m}), \\ R_{m}^{(n)} = \frac{1}{N_{m}^{(n)}} 1 (E_{m}^{(n)} > H_{m}^{(n)}) 1^{(n)} (t_{1}, t_{m}) . \end{matrix}$

Here we also adopt the conventions $\prod_{j = 2, j \neq m}^{1} \dots = 1$ , $\sum_{M + 1}^{M} \dots = 0$ . Initialise the adjoints of the risk factors so that ${\bar{X}}_{m}^{(n)} = {\bar{E}}_{m}^{(n)} \frac{\partial E_{m}^{(n)}}{\partial X_{m}^{(n)}} + {\bar{N}}_{m}^{(n)} \frac{\partial N_{m}^{(n)}}{\partial X_{m}^{(n)}} .$ (53)

For m = M, …, 1, update the adjoint of the regression coefficients β_m to give ${\bar{β}}_{m} + = ψ_{m}^{(n)} {\bar{H}}_{m}^{(n)},$ (54) as well as the adjoints of the basis functions $ψ_{m}^{(n)}$ , ${\bar{ψ}}_{m}^{(n)} = β_{m} {\bar{H}}_{m}^{(n)},$ (55) and update the adjoints of the state vector ${\bar{X}}_{m}^{(n)} + = ({\bar{ψ}}_{m}^{(n)})^{T} \frac{\partial ψ_{m}^{(n)}}{\partial X_{m}^{(n)}},$ (56) where we use the standard notation “+ =” for the addition assignment operator.

For m = M, …, 0 compute the adjoint of the risk factor evolution such that $\begin{matrix} {\bar{X}}_{m}^{(n)} + = {\bar{X}}_{m + 1}^{(n)} \frac{\partial F}{\partial X_{m}^{(n)}} (T_{m}, X_{m}^{(n)}; θ), \\ \bar{θ} + = {\bar{X}}_{m + 1}^{(n)} \frac{\partial F}{\partial θ} (T_{m}, X_{m}^{(n)}, θ), \end{matrix}$ (57) where the gradients are computed by applying the rules of adjoint differentiation following the instructions that implement the function F. Finally, the adjoint ${\bar{X}}_{0}$ is utilised to populate the component of $\bar{θ}$ corresponding to the adjoint of the model parameter X₀.

3.3 AAD for XVA by regression

The AAD implementation of the algorithm for the calculat ion of XVA described in Section 2.5, producing the MC estimators for the sensitivities of the estimator (45) with respect to a set of model parameters θ_k, k = 1, …, N_θ, given by ${\bar{θ}}_{k} = \frac{\partial XVA (θ)}{\partial θ_{k}},$ (58) consists of the adjoint of steps (X1)-(X5) executed backwards for each MC replication n = 1, …, N_MC:

Set the adjoint of the XVA value, $\bar{XVA} = 1$ , the adjoint of the model parameters $\bar{θ} = 0$ and ${\bar{XVA}}^{(n)} = \bar{XVA} \frac{1}{N_{MC}} .$ (59)

For m = M, …, 1 compute: $\begin{matrix} {\bar{V}}_{m}^{(n)} & = & - \frac{{\bar{XVA}}^{(n)}}{N_{m}^{(n)}} [L_{c} ({SP}_{m - 1}^{c, (n)} - {SP}_{m}^{c, (n)}) 1 (V_{m}^{(n)} > 0) \\ + L_{d} ({SP}_{m}^{d, (n)} - {SP}_{m}^{d, (n)}) 1 (V_{m}^{(n)} < 0)], \end{matrix}$ (60) $\begin{matrix} {\bar{N}}_{m}^{(n)} & = & \frac{{\bar{XVA}}^{(n)}}{{(N_{m}^{(n)})}^{2}} [L_{c} ({SP}_{m - 1}^{c, (n)} - {SP}_{m}^{c, (i)}) {(V_{m}^{(n)})}^{+} \\ + L_{d} ({SP}_{m - 1}^{d, (n)} - {SP}_{m}^{d, (n)}) {(V_{m}^{(n)})}^{-}], \end{matrix}$ (61) ${\bar{SP}}_{m}^{c, (n)} = {\bar{XVA}}^{(n)} [\frac{V_{m}^{(n)}}{N_{m}^{(n)}} (1 - δ_{m, 0}) - \frac{V_{m + 1}^{(n)}}{N_{m + 1}^{(n)}}] 1 (V_{m}^{(n)} > 0),$ (62) ${\bar{SP}}_{m}^{d, (n)} = {\bar{XVA}}^{(n)} [\frac{V_{m}^{(n)}}{N_{m}^{(n)}} (1 - δ_{m, 0}) - \frac{V_{m + 1}^{(n)}}{N_{m + 1}^{(n)}}] 1 (V_{m}^{(n)} < 0),$ (63) where δ_m,n is the Kronecker symbol, and $V_{M + 1}^{(n)} / N_{M + 1}^{(n)} = 0$ .

For m = M, …, 1, set ${\bar{V}}_{p, m}^{(n)} = {\bar{V}}_{m}^{(n)}$ and compute the adjoint of Equation (43) by $\begin{matrix} {\bar{H}}_{p, m}^{(n)} & = & {\bar{V}}_{p, m}^{(n)} 1 (H_{p, m}^{(n)} > E_{p, m}^{(n)}), \\ {\bar{E}}_{p, m}^{(n)} & = & {\bar{V}}_{p, m}^{(n)} 1 (H_{p, m}^{(n)} < E_{p, m}^{(n)}), \end{matrix}$ (64) if T_m is an exercise date for the p-th option, and ${\bar{H}}_{p, m}^{(n)} = {\bar{V}}_{p, m}^{(n)}, {\bar{E}}_{p, m}^{(n)} = 0,$ (65) otherwise. Initialise the adjoints of the risk factors, ${\bar{X}}_{m}^{(n)} = {\bar{E}}_{p, m}^{(n)} \frac{\partial E_{p, m}^{(n)}}{\partial X_{m}^{(n)}} + {\bar{N}}_{m}^{(n)} \frac{\partial N_{m}^{(n)}}{\partial X_{m}^{(n)}} .$ (66) Initialise the adjoint of the regression coefficients β_p,m, and update them together with the adjoints of the basis functions ψ_p,m, and the adjoints of the state vector as in Equations (54) to (56).

Update the adjoint of the simulated hazard rates according to ${\bar{λ}}_{m}^{c, (n)} + = - (T_{m + 1} - T_{m}) \sum_{j = m + 1}^{M} {\bar{SP}}_{m}^{c, (n)} {SP}_{m}^{c, (n)},$ (67) ${\bar{λ}}_{m}^{d, (n)} + = - (T_{m + 1} - T_{m}) \sum_{j = m + 1}^{M} {\bar{SP}}_{m}^{d, (n)} {SP}_{m}^{d, (n)} .$ (68)

For m = M, …, 0 compute the adjoint of the hazard rate evolution $\begin{matrix} {\bar{λ}}_{m}^{c, (n)} + = {\bar{λ}}_{m + 1}^{c, (n)} \frac{\partial G_{c}}{\partial λ_{m}^{c, (n)}} (T_{m}, λ_{m}^{c, (n)}; θ), \\ {\bar{λ}}_{m}^{d, (n)} + = {\bar{λ}}_{m + 1}^{d, (n)} \frac{\partial G_{d}}{\partial λ_{m}^{d, (n)}} (T_{m}, λ_{m}^{d, (n)}; θ), \end{matrix}$ (69) $\begin{matrix} \bar{θ} + = {\bar{λ}}_{m + 1}^{c, (n)} \frac{\partial G_{c}}{\partial θ} (T_{m}, λ_{m}^{c, (n)}; θ), \\ \bar{θ} + = {\bar{λ}}_{m + 1}^{d, (n)} \frac{\partial G_{d}}{\partial θ} (T_{m}, λ_{m}^{d, (n)}; θ), \end{matrix}$ (70) where the gradients can be computed through the AAD implementation of the function G and of the risk factor evolution as in step ( $\bar{L 1}$ ).

3.4 AAD for the regression algorithm

In the AAD implementations presented in Sections 3.2 and 3.3 the adjoint of the regression coefficients in Equation (54) do not contribute to the calculation of sensitivities. As mentioned in Section 3.1, and illustrated with numerical examples in Section 4, this can be justified as a reasonable approximation in the case of Bermudan-style options when the exercise boundary approximated by regression is close to the optimal one. However, such arguments are generally approximations for Bermudan-style options, and cannot be invoked at all when regression is utilised for XVA. In this case, the contributions associated with the sensitivities of the regression coefficients must be taken into account in order to obtain accurate estimates of the model parameters sensitivities. In this section, we discuss how these contributions can be computed by the AAD implementation of the least-square MC algorithm described in Section 2.3.

Skip this step if the regression algorithm is utilised in conjunction with a second independent simulation for Bermudan-style options or in the context of XVA. Initialise the adjoint of the option value ${\bar{V}}_{0} = 1$ , the adjoint of the model parameters $\bar{θ} = 0$ , the adjoint of the regression coefficients ${\bar{β}}_{m} = 0$ , m = 1, …, M - 1 and set $\begin{matrix} {\bar{V}}_{1}^{(n)} = \frac{{\bar{V}}_{0}}{N_{MC}} \frac{1}{N_{1}^{(n)}}, \\ {\bar{N}}_{1}^{(n)} = - \frac{{\bar{V}}_{0}}{N_{MC}} \frac{V_{1}^{(n)}}{{(N_{1}^{(n)})}^{2}}, \end{matrix}$ (71) for n = 1, …, N_MC.

For m = 1 to M - 1, we have:

For n = 1, …, N_MC, compute $\begin{matrix} {\bar{E}}_{m}^{(n)} = {\bar{V}}_{m}^{(n)} 1 (E_{m}^{(n)} > H_{m}^{(n)}), \\ {\bar{H}}_{m}^{(n)} = {\bar{V}}_{m}^{(n)} 1 (E_{m}^{(n)} < H_{m}^{(n)}), \end{matrix}$ (72) and initialise the adjoint of the risk factor path value $X_{m}^{(n)}$ by ${\bar{X}}_{m}^{(n)} = {\bar{E}}_{m}^{(n)} \frac{\partial E_{m}^{(n)}}{\partial X_{m}^{(n)}} .$ (73)

Set: ${\bar{β}}_{m} + = \sum_{n = 1}^{N_{MC}} ψ_{m}^{(n)} {\bar{H}}_{m}^{(n)},$ (74) and initialise the adjoint of the basis functions ψ as in Equation (55)

Compute the adjoint of the two intermediate variables Ω_m and Ψ_m in Equation (27) using the results in Giles (2008) by ${\bar{Ω}}_{m} = Ψ_{m}^{- T} {\bar{β}}_{m}, {\bar{Ψ}}_{m} = - {\bar{Ω}}_{m} β_{m}^{T} .$ (75)

For n = 1, …, N_MC compute the adjoint of Equation (26) by $\begin{matrix} {\bar{ψ}}_{m}^{(n)} + = \frac{1}{N_{MC}} N_{m}^{(n)} V_{m + 1}^{(n)} {\bar{Ω}}_{m}, \\ {\bar{V}}_{m + 1}^{(n)} = \frac{1}{N_{MC}} N_{m}^{(n)} (ψ_{m}^{(n)})^{T} {\bar{Ω}}_{m}, \end{matrix}$ (76) ${\bar{N}}_{m}^{(n)} = \frac{1}{N_{MC}} V_{m + 1}^{(n)} (ψ_{m}^{(n)})^{T} {\bar{Ω}}_{m},$ where $N_{m}^{(n)} = N_{m}^{(n)} / N_{m + 1}^{(n)}$ , and compute the adjoint of Equation (25) by ${\bar{ψ}}_{m}^{(n)} + = \frac{1}{2 N_{MC}} ({\bar{Ψ}}_{m}^{(n)} + {({\bar{Ψ}}_{m}^{(n)})}^{T}) ψ_{m}^{(n)} .$ (77)

Then update the adjoint of the risk factor vectors by $\begin{matrix} {\bar{X}}_{m}^{(n)} + = \frac{{\bar{N}}_{m}^{(n)}}{N_{m + 1}^{(n)}} \frac{\partial N_{m}^{(n)}}{\partial X_{m}^{(n)}} + {\bar{ψ}}^{T} \frac{\partial \bar{ψ}}{\partial X_{m}^{(n)}}, \\ {\bar{X}}_{m + 1}^{(n)} + = - \frac{{\bar{N}}_{m}^{(n)} N_{m}^{(n)}}{N_{m + 1}^{(n)}} \frac{\partial N_{m + 1}^{(n)}}{\partial X_{T_{m + 1}}^{(n)}} . \end{matrix}$ (78)

Compute the adjoint of the risk factor vector at expiry by ${\bar{X}}_{M}^{(n)} = {\bar{V}}_{M}^{(n)} \frac{\partial E_{M}^{(n)}}{\partial X_{M}^{(n)}},$ (79) for n = 1, …, N_MC.

For m = M, …, 0, compute the adjoint of the risk factor evolution as in step ( $\bar{L 1}$ ).

4 Numerical results

4.1 Best of two assets Bermudan-style option

As a first example, we consider the classical case of a Bermudan-style option on the maximum of two assets under a standard lognormal model for the asset price processes {X₁ (t)} and {X₂ (t)}. The payoff function at exercise time t is given by ${(max {X_{1} (t), X_{2} (t)} - K)}^{+},$ (80) where K is the strike price. A similar example is also studied in Glasserman (2003); Broadie and Glasserman (1997) and Broadie and Glasserman (2004). Andersen and Broadie (2004) obtain nearly exact estimates for this example using twelve basis functions, including the value of a European-style option. We assume the option can be exercised every three months up to 3 years. We assume the underlying assets are independent geometric Brownian motion processes with the same initial values X (0), volatilities σ, and dividend rates d. In particular, we choose X₁ (0) = X₂ (0) =1, σ₁ = σ₂ = 0.2 and d₁ = d₂ = 0.1 and the risk-free rate r = 0.05. We compute the Bermudan-style option prices, the Deltas (i.e., the sensitivity with respect to the spot values, X_i (0)) and the Vegas (i.e., the sensitivity with respect to the volatilities σ_i) for K = 0.9, 1.0 and 1.1 with AAD and bumping. We also compare these with their “exact” value, which is obtained by a PDE approach (using finite differences to estimate sensitivities). The numerical results are found in Table 1 and in Fig. 1 demonstrating an excellent agreement between the AAD and PDE results. The regression basis functions we adopt are the terms of the third-order polynomial ${1, x_{1}, x_{2}, x_{1}^{2}, x_{2}^{2}, x_{1} x_{2}, x_{1}^{3}, x_{1}^{2} x_{2}, x_{1} x_{2}^{2}, x_{2}^{3}}$ (81) and the payoff functions up to power three, i.e.,

$\begin{matrix} {{(max (x_{1}, x_{2}) - K)}^{+}, {({(max (x_{1}, x_{2}) - K)}^{+})}^{2}, \\ {({(max (x_{1}, x_{2}) - K)}^{+})}^{3}} . \end{matrix}$ (82)

Fig.1

Deltas (left panel) and Vegas (right panel) for the Bermudan-style option in (80) as a function of strike. The smoothening parameter in the call spread regularization (48) is δ = 0.005. The number of MC paths is 400,000. The results are obtained with the setting (81) and (82). The MC uncertainty (in parenthesis) is computed using the binning technique with 20 bins for each set of simulations.

Table 1

Prices, Deltas and Vegas for the Bermudan-style option in (80) with three different strike values. The smoothening parameter in the call spread regularization (48) is δ = 0.005. The number of MC paths employed in the AAD approach is 400,000. The results are obtained with the setting (81) and (82). The MC uncertainty (in parenthesis) is computed using the binning technique of Capriotti and Giles (2010) with 20 bins for each set of simulations

K = 0.9	Price	Delta	Vega
AAD	0.2012(2)	0.415(3)	0.456(2)
PDE	0.20107	0.41423	0.45740
K = 1.0
AAD	0.1396(1)	0.337(2)	0.486(2)
PDE	0.13959	0.33588	0.48440
K = 1.1
AAD	0.0941(2)	0.258(1)	0.463(2)
PDE	0.09431	0.25635	0.46253

Here the smoothening parameter for the calculation of the Greeks, discussed in Section 3.1, δ = 0.005, was chosen as a reasonable compromise between variance and bias of the estimator. This is illustrated in Fig. 2, showing how for δ = 0.005 the bias introduced by the finite δ is of the same order of magnitude of the statistical uncertainty for the chosen computational budget, and is negligible for any practical application.

Fig.2

AAD Delta (left panel) and Vega (right panel) of the Bermudan-style option in (80) for K = 1 vs the smoothening parameter δ for the call spread regularization (48). The number of simulated paths is 3,000,000 for δ = 0.01 and is increased as δ is decreased in order to keep statistical uncertainties roughly constant. The results are obtained with the setting (81) and (82). The values in the graphs are fitted based on a quadratic polynomial function (green lines).

As discussed in Section 3.1, neglecting to smoothen out the exercise boundary, although common in the financial practice, introduces a bias in the computation of sensitivities because the exercise boundary is in general not optimal. This is illustrated in Fig. 3, in which we compare the Delta, with and without smoothening, for different choices of the basis functions. These results are obtained with the inputs (81) and (82). Here for Delta, the smoothened estimator turns out to be more accurate, especially when the exercise boundary obtained by regression is a less accurate approximation of the real one. However, it is difficult to establish a priori whether the unsmoothened estimator provides a smaller or a larger bias than the smoothened one. This is because the bias introduced by the lack of smoothening may be offset by the bias arising from the sub-optimality of the exercise boundary. This is illustrated in the right panel of Fig. 3 showing that for Vega the smoothened and unsmoothened estimators have a comparable accuracy. In any case, a consideration to keep in mind is that smoothening the exercise boundary is generally required to obtain stable second-order risk values.

Fig.3

AAD Delta (left panel) and Vega (right panel) of the Bermudan-style option in (80) for K = 1 as obtained with the unsmoothened and the smoothened estimators with the call-spread regularization (48) for five choices of the regression basis functions. The number of simulated paths is 1,000,000 with 20 bins and the smoothening parameter is δ = 0.005.

The quasi-optimality of the exercise boundary is also generally invoked among industry practitioners as a justification for neglecting the contributions to the sensitivities arising from the exercise boundary. Clearly, the quality of this approximation is dependent on the accuracy of the regression functions in reproducing the actual exercise boundary. This is illustrated in Fig. 4 where we plot Delta and Vega of the Bermudan-style option (80) for different choices of the regression basis functions. Here we compare the results obtained by the AAD algorithm, as described in Section 3.4, in the case where a) the exercise boundary is kept fixed, and b) when accounting instead for its contributions to the sensitivities. As expected, the difference between the two approaches vanishes as the regression functions become more accurate. However, as it is shown for the Vega, it can lead to a significant bias if a simple (e.g., linear or quadratic) representation of the exercise boundary is adopted.

Fig.4

AAD Delta (left panel) and Vega (right panel) of the Bermudan-style option in (80) for K = 1 as obtained by keeping the exercise boundary fixed (Fixed) and accounting instead for its contributions to the sensitivities (Flexible). The number of simulated paths is 1,000,000 with 20 bins and the call-spread smoothening parameter is δ = 0.005.

4.2 XVA sensitivities

As another example, we compute the sensitivities of XVA (14) for the same option defined in Equation (80). Here, for simplicity, we assume that the hazard rate and the volatility are piecewise constants. The following results are obtained with the specifications (81) and (82). The XVA sensitivities with respect to some of the model parameters, namely the term structure of hazard rates of the counterparty and of volatilities of the underlying, obtained with the AAD algorithm described in Section 2.5 are compared with the ones obtained by standard bumping in Tables 2 and 3. As expected, the AAD sensitivities are in excellent agreement with those obtained by bumping with any discrepancies attributable to the bias of the finite-difference approach completely masked in this example by the MC uncertainties.

Table 2
The XVA sensitivities with respect to the piecewise volatility. The results are obtained with the setting (81) and (82). Both sets of the results are computed with 1,000,000 MC paths and 25 bins

σ ₀ σ ₁ σ ₂ σ ₃ σ ₄ σ ₅

AAD 0.00298(1) 0.00297(2) 0.00294(1) 0.00283(1) 0.00271(1) 0.00255(1)

Bumping 0.00298(1) 0.00297(2) 0.00294(1) 0.00283(1) 0.00271(1) 0.00255(1)

σ ₆ σ ₇ σ ₈ σ ₉ σ ₁₀ σ ₁₁

AAD 0.00236(1) 0.00218(1) 0.00201(1) 0.00183(1) 0.00165(1) 0.001450(8)

Bumping 0.00236(1) 0.00218(1) 0.00201(1) 0.00183(1) 0.00165(1) 0.001450(8)

	σ ₀	σ ₁	σ ₂	σ ₃	σ ₄	σ ₅
AAD	0.00298(1)	0.00297(2)	0.00294(1)	0.00283(1)	0.00271(1)	0.00255(1)
Bumping	0.00298(1)	0.00297(2)	0.00294(1)	0.00283(1)	0.00271(1)	0.00255(1)
	σ ₆	σ ₇	σ ₈	σ ₉	σ ₁₀	σ ₁₁
AAD	0.00236(1)	0.00218(1)	0.00201(1)	0.00183(1)	0.00165(1)	0.001450(8)
Bumping	0.00236(1)	0.00218(1)	0.00201(1)	0.00183(1)	0.00165(1)	0.001450(8)

Table 3

The XVA sensitivities with respect to the piecewise hazard rates. The results are obtained with the setting (81) and (82). Both sets of the results are computed with 1,000,000 MC paths and 25 bins

	λ ₀	λ ₁	λ ₂	λ ₃	λ ₄	λ ₅
AAD	0.03337(1)	0.03336(2)	0.03329(2)	0.03306(3)	0.03269(3)	0.03218(3)
Bumping	0.03337(1)	0.03336(2)	0.03329(1)	0.03306(3)	0.03269(3)	0.03218(3)
	λ ₆	λ ₇	λ ₈	λ ₉	λ ₁₀	λ ₁₁
AAD	0.03151(3)	0.03074(4)	0.02989(4)	0.02893(4)	0.02784(4)	0.02669(4)
Bumping	0.03151(3)	0.03074(3)	0.02989(4)	0.02893(4)	0.02784(4)	0.02669(4)

Similar to the case of Bermudan-style option Greeks, keeping the regression boundary fixed while computing the sensitivities, introduces a bias. However, for XVA, this issue is much more severe than in the case of the Bermudan-style option Greeks because no quasi-optimality argument can be invoked. As shown in Fig. 5, the volatilities sensitivities obtained by bumping and AAD keeping into account of the regression contributions are almost identical. Instead, the results of AAD keeping the regression coefficients fixed are remarkably different and, if utilised for risk management, would lead to significant mis-hedging. On the other hand, as expected, the XVA sensitivities with respect to the hazard rates are not affected by the contribution arising from the regression functions. This is because the hazard rates do not enter the computation of the portfolio value conditional on default, and hence do not appear in the regression parameters.

Fig.5

XVA sensitivities with respect to the piecewise volatility (left panel) and hazard rate (right panel), computed by AAD keeping into account of the contributions of the regression coefficients (Flexible), by AAD keeping the regression coefficients fixed (Fixed), and by finite-differences (Bumping). The number of MC paths is 1,000,000 and the number of bins is 25. The results are obtained with the setting (81) and (82).

Finally, the remarkable computational efficiency of the AAD approach is illustrated in Fig. 6. Here we plot the cost of computing the XVA sensitivities with respect to the term structure of the counterparty hazard rate and the underlying volatility, relative to the cost of performing a single valuation. The calculation of the sensitivities by means of AAD can be performed in about three times the cost of computing the XVA value, that is, well within the theoretical bound (7). In contrast, the cost of bumping, for one-sided finite-difference estimators, is in general (1 + N) times the cost of as single valuation, where N is the number of model parameters, i.e., in this case over 20 times the cost of computing the value of the XVA.

Fig.6

Ratio of the CPU time required for the calculation of the XVA, and its sensitivities, and the CPU time spent for the computation of the XVA alone.

5 Conclusions

We have shown how Adjoint Algorithmic Differentiation (AAD) can be utilised to implement efficiently the computation of sensitivities in regression-based MC methods. We have derived the AAD version of the celebrated least-square algorithms of Tsitsiklis and Van Roy (2001) and Longstaff and Schwartz (2001) and, by discussing in detail examples of practical relevance, we have demonstrated how accounting for the sensitivities contributions associated with the regression functions is crucial to obtain accurate estimates of the Greeks in XVA applications and for Bermudan-style options, especially when the exercise boundary is not particularly accurate. We have also shown how smoothening out the discontinuities associated with suboptimal exercise boundaries can lead in some situations to more accurate estimates of the sensitivities.

From a computational stand point, similarly to what was previously found in other MC and PDE settings, see e.g., Capriotti and Giles (2012); Capriotti et al. (2015), the proposed method allows the computation of the complete first-order risk at a cost which is at most four times the cost of calculating the value of the portfolio itself. This typically results in orders of magnitude of savings in computational time with respect to standard finite-difference approaches, thus making AAD an ever so indispensable technique in the toolbox of modern financial and insurance engineering.

Footnotes

In the following, for simplicity of notation, we will set N(0) =1.

Here and in the following, to keep the notation simple, we do not introduce different symbols for expectations and the respective sample averages.

Acknowledgments

It is a pleasure to thank Jacky Lee, for numerous useful discussions, and for a careful reading of the manuscript. We are also grateful to the anonymous referee for the feedback and comments. After completion of this manuscript we became aware of the work performed simultaneously by Antonov on a similar application of AAD to Callable Exotics, see . The opinions and views expressed in this paper are uniquely those of the authors, and do not necessarily represent those of Credit Suisse Group.

References

Andersen

, Broadie

2004. Primal-dual simulation algorithm for pricing multidimen-sional American options, Management Science 50(9), 1222–1234.

Andersen

, Piterbarg

2010. Interest Rate Modeling, Volume I: Foundations and Vanilla Models, Atlantic Financial Press.

Antonov

2016. Algorithmic Differentiation for Callable Exotics. SSRN: https://ssrn.com/abstract=2881992

Broadie

, Glasserman

1997. Pricing American-style securities using simulation, Journal of Economic Dynamics and Control 21(8), 1323–1352.

Broadie

, Glasserman

2004. A stochastic mesh method for pricing high-dimensional American options, Journal of Computational Finance 7(4), 35–72.

Capriotti

2011. Fast greeks by algorithmic differentiation, Journal of Computational Finance 14(3), 3–35.

Capriotti

, Giles

2010. Fast correlation Greeks by adjoint algorithmic differentiation, Risk 23, 79–85.

Capriotti

, Giles

2012. Algorithmic differentiation: Adjoint greeks made easy, Risk 25, 92–98.

Capriotti

, Jiang

, Macrina

2015. Real-time risk management: An AAD-PDE approach, International Journal of Financial Engineering 2(04), 1550039–1550070.

10.

Capriotti

, Lee

S. J.

2014. Adjoint credit risk management, Risk 27, 92–98.

11.

Capriotti

, Peacock

, Lee

2011. Real time counterparty credit risk management in Monte Carlo, Risk 24, 86–90.

12.

Carriere

J.F.

1996. Valuation of the early-exercise price for options using simulations and non-parametric regression, Insurance: Mathematics and Economics 19(1), 19–30.

13.

Cesari

, Aquilina

, Charpillon

, Filipovic

, Lee

, Manda

2009. Modelling, pricing, and hedging counterparty credit exposure: A technical guide. Springer Science & Business Media.

14.

Crépey

, Bielecki

T.R.

, Brigo

2014. Counterparty risk and funding: A tale of two puzzles. CRC Press.

15.

Friedman

, Hastie

, Tibshirani

2001. The elements of statistical learning (Vol. 1). Springer, Berlin: Springer series in statistics.

16.

Geeraert

, Lehalle

C.A.

, Pearlmutter

, Pironneau

, Reghai

2017. Mini-symposium on automatic differentiation and its applications in the financial industry. https://arxiv.org/abs/1703.02311

17.

Giles

M.B.

2008. Collected matrix derivative results for forward and reverse mode algorithmic differentiation. In Advances in Automatic Differentiation, 35–44. Springer Berlin Heidelberg.

18.

Glasserman

2003. Monte Carlo methods in financial engineering (Vol. 53). Springer Science & Business Media.

19.

Griewank

and Walther

2008. Evaluating derivatives: Principles and techniques of algorithmic differentiation, Society for Industrial and Applied Mathematics.

20.

Henrard

2014. Adjoint algorithmic differentiation: Calibration and implicit function theorem, Journal of Computational Finance 17(4), 37–47.

21.

Joshi

, Kwon

O.K.

2016. Least squares Monte Carlo credit value adjustment with small and unidirectional bias, International Journal of Theoretical and Applied Finance 19(08), 1650048–1650064.

22.

Lando

1998. On Cox processes and credit risky securities, Review of Derivatives Research 2(2-3), 99–120.

23.

Leclerc

, Liang

, Schneider

2009. Fast Monte Carlo Bermudan greeks, Risk 22(7), 84–88.

24.

Longstaff

F.A.

, Schwartz

E.S.

2001. Valuing American options by simulation: A simple least-squares approach, Review of Financial Studies 14(1), 113–147.

25.

Savickas

, Hari

, Wood

, Kandhai

2014. Super fast greeks: An application to counterparty valuation adjustments, Wilmott 69, 76–81.

26.

Tsitsiklis

J.N.

, Van Roy

2001. Regression methods for pricing complex American-style options, IEEE Transactions on Neural Networks 12(4), 694–703.

27.

Wilmott

, Howison

, Dewynne

1995. The mathematics of financial derivatives: A student introduction. Cambridge University Press.

28.

, Chen

, Coleman

T.F.

2016. The effcient application of automatic differentiation for computing gradients in financial applications, Journal of Computational Finance 19(3), 71–96.