1 Introduction↩︎

The world population is ageing fast with life expectancy raising near 90 years in some countries. Longevity risk, i.e., the risk of outliving one’s savings, is becoming a critical concern for retirees. Variable annuities (VA) with living and death benefit guarantees are financial products that can help to manage this risk. These products provide savings protection while taking advantage of market growth at the same time. VAs are typically classified based on the type of guarantees they provide. These include guaranteed minimum death benefit (GMDB), guaranteed minimum withdrawal benefit (GMWB), guaranteed minimum maturity benefit (GMMB), guaranteed minimum accumulation benefit (GMAB), and guaranteed minimum income benefit (GMIB) [1]. A good overview of VA products and the development of their market can be found in [2], [3] and [4]. Insurers started to sell these types of products in the 1990s in United States. Later, these products became popular in Europe, UK and Japan as well. The market for VAs has become very large. In the United States, according to the LIMRA U.S. Individual Annuities Survey, annual sales of VAs between 2003 and 2014 ranged from $129 to $184 billion (peaking in 2007). Sales then declined to around $100 billion per year by 2016 and remained near this level until 2023, before increasing to $126.5 in 2024 and $142.6 billion in 2025. It is also interesting to observe that sales of fixed annuities were below VA sales during 2003-2015, reaching as little as about half of VA sales in some years. However, this trend later reversed, with fixed annuity sales reaching $307.6 in 2024 and $321.5 billions in 2025, more than double the corresponding VA sales).

VAs with guarantees provide flexibility to the policyholder to withdraw funds from the contract, though some penalty may apply for deviating from the contractual rate. If the policyholder behaves passively and makes withdrawals at the contractual rate defined at the beginning of the contract, then the behaviour of the policyholder is called static. In this case the paths of the wealth account can be simulated forward, and a standard Monte Carlo (MC) simulation method can be used for GMxB pricing. On the other hand, if the policyholder optimally decides the amount to withdraw at each withdrawal date, then the behaviour of the policyholder is called dynamic. Under the optimal withdrawal strategy, the fair pricing of variable annuities with GMWB becomes a stochastic optimal control problem.

In the literature, popular numerical methods to price variable annuities include

• Tree/lattice algorithms: GLWB in [5], GMWB in [6]–[8], GMWB and GMDB in [9], GMMB and GMAB in [10],

• Numerical PDE methods: GMMB in [11], GLWB in [5], GMWB in [7], [12], GMWB and GMDB in [9], and

• Numerical integration: GMWB in [13], GMWB with surrender in [14], GLWB in [15], GMAB in [16], [17], GMWDB in [18], GMMB and GMAB in [10], including the Fourier cosine expansion method (GMWB in [19], GMDB in [20], GMMB in [21]–[24]) and the frame project method (GMDB in [25], GMAB and GMDB in [26], GMMB and GMDB in [27]).

All these numerical methods are expected to perform well when the number of stochastic state variables of the problem is low (one or two), and to become quickly infeasible when this number grows. When the problem involves multiple state variables, an attractive numerical method to solve discrete-time stochastic control problems with finite horizon is the Least-Squares Monte Carlo (LSMC, [28], [29]), due to the better convergence guarantees of Monte Carlo methods in high dimension. The LSMC algorithm was introduced to variable annuity pricing in [30] in the case of a semi-static strategy with only one surrender guarantee. This was later extended in [31] in the specific case of a dynamic strategy with bang-bang control (the holder of a GLWB can maximize the contract writer’s losses by only performing non-withdrawal, withdrawal at exactly the contract rate, or full surrender). In such situation, the pricing problem can be formulated as an optimal switching problem, for which the LSMC algorithm is known to be easily applicable [32]–[34]. Similarly, [7] considered dynamic withdrawals under the simplifying assumption that the optimal withdraws should be multiples of the guaranteed rate.

In this paper, we revisit the problem of pricing variable annuities with embedded guarantees by LSMC, and generalize it to the case of fully dynamic optimal strategies with fully controlled underlying stochastic state variables, and no ex ante simplifications of the optimal dynamic strategy.

Using the GMWB example, the standard LSMC algorithm cannot be applied directly, as the paths of the underlying wealth process are altered by the optimal cash withdrawals that should be found from the backward-in-time solution, meaning that the underlying wealth process cannot be initially simulated forward in time as expected from the LSMC procedure. To overcome this problem, we introduce the use of the control randomization technique [35], [36] to VA pricing, to extend LSMC to handle stochastic optimal control problems with controlled Markov processes.

As emphasized in [37] or [34] for example, the LSMC algorithm can perfectly work with any other type of regression than ordinary least squares regression (in which case “Regression Monte Carlo” is probably a more appropriate name for the resulting algorithm). In fact, while [29] used OLS regression, the seminal article [28] used nonparametric regression instead (spline regression and local kernel regression). Since then, many other types of regression have been investigated for LSMC applied to optimal stopping problems (see [38] and references therein).

In the context of discrete-time, finite horizon stochastic control problems, estimating continuation values by neural network regression was investigated in [39]. This approach has been applied to option hedging [40], [41], storage optimization [41], [42], portfolio optimization [43], [44], utility maximization [45], [46], and, recently, to the pricing of convertible bonds with path-dependent provisions [47]. In this article, we introduce this methodology to the pricing of variable annuities with embedded guarantees.

Finally, we also establish how variations of the methods can be implemented to easily obtain confidence intervals for the optimal contract value.

The paper is organized as follows. Section 2 details how to model the pricing of a variable annuity with Guaranteed Minimum Withdrawal Benefit (GMWB) under optimal dynamic withdrawals as a stochastic optimal control problem with controlled underlying stochastic factors. Section 3 describes how to extend the Least Squares Monte Carlo (LSMC) algorithm to solve such problem, and obtain confidence intervals for the optimal contract value. Section 4 provides numerical experiments illustrating the accuracy of the method, and Section 5 concludes the paper.

2 Pricing VA riders as a stochastic control problem↩︎

The specification details of VA riders vary across different companies offering these products and the results for specific GMxB riders presented in the academic literature often refer to different specifications. However, in general, the present value of the overall payoff of the VA contract with a guarantee can be written as

\[\label{total95payoff95eq} H_0({X},{\pi})= \beta_{0,N} R_N({X}_N)+\sum_{n=1}^{N-1} \beta_{0,n} R_n({X}_n,\pi_n),\tag{1}\] where \(R_N(X_N)\) is the cashflow amount received by the policy holder at the contract maturity, \(R_n({X}_n,\pi_n)\) is the cashflow amount received by the policyholder at time \(t_n\), \(\pi_n\) is the withdrawal amount from the contract account, and \(\beta_{i,j}\) is the discounting factor from \(t_j\) to \(t_i\) \[\beta_{i,j}=\exp\left(-\int_{t_i}^{t_j}r(t)dt\right),\;t_j>t_i,\] where \(r(t)\) is the risk-free interest rate (possibly stochastic). That is, we consider the time discretization \(0=t_0<t_1<\cdots<t_N=T\) corresponding to the contract withdrawal dates, where \(t_0=0\) is today and \(T\) is the contract maturity. Here, \({X}=\{X_n\}_{n=0}^N\), and \({X}_n\) is the state variable vector (containing variables such as wealth account, guarantee account, etc.) before the withdrawal \(\pi_n\) at time \(t_n\) evolving as \[\label{eq:transitionfunction} {X}_{n+1} = \mathcal{T}_n \left ({X}_n, \pi_n, {Z}_{n+1} \right ),\tag{2}\] where \({Z}_1,...,{Z}_N\) are independent random disturbances (i.e., the state variable \({X}=\{X_n\}_{n=0}^N\) is a Markov process). Under the assumption of no-arbitrage market with respect to the financial risk, the price of the VA with GMxB can be expressed as an expectation with respect to the risk-neutral probability measure for the underlying risky asset. Let \((\Omega,\mathcal{F},\mathbb{Q})\) be a probability space with sample space \(\Omega\), filtration \(\mathcal{F}=\{\mathcal{F}_t : t\ge 0\}\) and risk-neutral probability measure \(\mathbb{Q}\). Then, the contract fair price under the given withdrawal strategy \({\pi}=(\pi_1,\ldots,\pi_{N-1})\), can be calculated as \[\label{GMWDB95general95eq} V_0\left({X}_0\right)=\mathbb{E}\left[H_0({X},{\pi})|\mathcal{F}_0\right].\tag{3}\] Here, \(\mathbb{E}[\cdot|\mathcal{F}_t]\) denotes an expectation under the risk-neutral probability measure \({\mathbb{Q}}\) with respect to the state vector \({X}\), conditional on information available at time \(t\).

The insurer is collecting the fee \(\alpha\) for the guarantee affecting the evolution of \(X_t\) and the contract value \(V_0(X_0;\alpha)\) depends on this fee rate. Then, the fair value of the fee \(\alpha\) corresponds to the solution of the equation \(V_0(X_0;\alpha)=W_0\), where \(W_0\) is the full amount invested into the contract at \(t_0\).

The withdrawal strategy \(\pi\) can depend on the time and state variables and is assumed to be given when the price of the contract is calculated in 3 . Different types of withdrawal strategies exist. They may be classified as either static, optimal, or sub-optimal.

• Static strategy. Under this strategy, the policyholder decisions are deterministically determined at the beginning of the contract and do not depend on the evolution of the wealth and benefit base accounts. For example, the policyholder may withdraw at the contractual rate only.

• Optimal strategy. Under the optimal withdrawal strategy, the decision on the withdrawal amount \(\pi_n\) depends on the information available at time \(t_n\), i.e. depends on the state variable \(X_n\). The optimal strategy is calculated as \[\label{optimalstrategy95eq} {\pi}^\ast({X})=\underset{{\pi\in\mathcal{A}}}{\mathrm{\arg\sup}}\;\mathbb{E}\!\left[H_0({X},{\pi})|\mathcal{F}_0\right],\tag{4}\] where the supremum is taken over all admissible strategies \(\pi\) in space \(\mathcal{A}\). Any other strategy \({\pi}({X})\) different from \({\pi}^\ast({X})\) is called sub-optimal and leads to a smaller price.

Evaluating the contract 3 under the optimal withdrawal strategy 4 is a standard optimal stochastic control problem for a controlled Markov process. Note that, the control variable \(\pi_n\) affects the transition law of the underlying state variable \({X}_t\) from \(t_n\) to \(t_{n+1}\) and thus the process is controlled. For a good textbook treatment of stochastic control problems in finance, see [48]. This type of problems can be solved recursively to find the contract value \(V_{n}(x)\) at \(t_n\) when \(X_n=x\) for \(n=N-1,\ldots,0\) via the backward induction Bellman equation \[\label{Bellman95eq} V_{n}(x)=\sup_{\pi_n\in \mathcal{A}_n}\left(R_n(x,\pi_n)+ \mathbb{E}\left[\beta_{n,n+1} V_{{n+1}}(X_{n+1})\bigg| X_{n}=x;\pi_n \right]\right),\tag{5}\] starting from the final condition \(V_N(x)=R_N(x)\). Note that the transition probability to reach state \(X_{n+1}=x^\prime\) at time \(t_{n+1}\) if the withdrawal (action) \(\pi_n\) is applied in the state \(x\) at time \(t_n\) depends on \(\pi_n\). Obviously, the above backward induction can also be used to calculate the fair contract price in the case of a static strategy \(\pi\); in this case the space of admissible strategies \(\mathcal{A}_n\) contains only one pre-defined value and \(\sup(\cdot)\) becomes redundant.

3 Numerical valuation of GMWB via LSMC↩︎

In general, solving the dynamic programming problem 5 needs to be done numerically, and can be very computationally intensive if a quadrature-based method is used for evaluating the expectation in 5 as the number of states, stochastic and control variables increases. The idea behind utilitising the LSMC method is to approximate the conditional expectation in 5 , \[\label{eq:conditionalexp} \Phi_n(X_{n},\pi_n) = \mathbb{E}\left[ \beta_{n,n+1} V_{n+1}(X_{n+1}) | X_{n} ; \pi_n \right],\tag{6}\] by a regression scheme with independent variables \(X_n\) and randomised \(\pi_n\), and response variable \(\beta V_{n+1}(X_{n+1})\). The approximation of the function is denoted as \(\widehat{\Phi}_n\). Hereafter, for notational convenience we drop indexes for discounting factors.

3.1 LSMC algorithms↩︎

The LSMC proceeds first by simulation of the random state, control and disturbance variables for \(X_t^m\), \(\pi_t^m\), \(m=1,...,M\), \(t=0,...,T\) (forward simulation) as in Algorithm 1, where \(Rand\) corresponds to random sampling from some distribution that could be designed for the specific problem. Then, the problem is solved with backward in time induction that can be accomplished using various algorithms. The most popular algorithms are: Realised Value Algorithm 2 and Regression Surface Algorithm 3 implemented in Regress Now fashion. Both can be implemented in Regress Later fashion; see Algorithms 4 and 5. Typical implementations of LSMC are based on the ordinary least squares (OLS) regression, which can be replaced with nonlinear regression such as, e.g. neural network regression used in our study.

For Regress Now algorithms, at each point in time \(t<T\), the value function is approximated using the least-squares regression \[\begin{align} &\beta V_{n+1}(X^{m}_{n+1})=f_{\theta_n}(X^{m}_n,\pi^m_n)+\epsilon^{m}_n, \\ &\epsilon^{m}_n\overset{iid}{\sim} F_n(\cdot),\;\;\mathbb{E}[\epsilon^{m}_n]=0,\;\;\mathbb{V}\mathrm{ar}[\epsilon^{m}_n]=\sigma_n^2,\;\;m = 1,...,M, \end{align}\] where \(f_{\theta_n}(X^{m}_n,\pi^m_n)\) is a parameterised regression function (parameterised by a vector \(\theta_n\)). In the case of the OLS regression, then \(f_{\theta_n}(X^{m}_n,\pi^m_n)=\boldsymbol{\Lambda}_n' \mathbf{L}(X^{m}_n,\pi^m_n)\), where \(\mathbf{L}(X^m_n,\pi^m_n)\) is a vector of basis functions and \(\mathbf{\Lambda}_n\) the regression coefficient vector. In our study we consider polynomial basis functions and a deep neural network parameterisation for \(f_\theta(X^{m}_n,\pi^m_n)\). Then, after the parameters \(\hat{\theta}_n\) have been estimated, \[\label{unbiasedestimate95eq} \widehat{{\Phi}}_n(X_n,\pi_n)={f}_{\hat{\theta}_n}(X_n,\pi_n).\tag{7}\] For Regress Later, the value function is approximated via the least-squares regression \[\beta V_{n+1}(X^{m}_{n+1})=f_{\theta_{n+1}}(X^{m}_{n+1})+\epsilon^{m}_{n+1},\] where \(f_{\theta_{n+1}}(X^{m}_{n+1})=\boldsymbol{\Lambda}_{n+1}' \mathbf{L}(X^{m}_{n+1})\) in the case of OLS, and then the expectation \[\Phi_n(X_{n},\pi_n) = \mathbb{E}\left[ \beta_{n,n+1} V_{n+1}(X_{n+1}) | X_{n} ; \pi_n \right]\] is calculated. It can be calculated numerically via quadrature methods if the number of stochastic disturbances is one or two. If a polynomial basis is used and the moments of \(X_{n+1}|X_n,\pi_n\) are known in closed form, then this expectation can be computed in closed form (see 6.2 for an example).

Sometimes, to avoid difficulties in the approximation of the value function, a transformation that has a similar shape as the value function can be applied on the value function in the above regression. Also, heteroscedasticity can be present in residuals, and then conditional variance can be modelled using another regression. These improvements of the standard LSMC algorithm are described in [50].

Two types of LSMC implementations are Regression Surface and Realised Value. In the case of Regression Surface, the value function is calculated in the backward recursion as \[\hat{V}_n(X_n)=R_n(X_n,\pi^\ast(X_n))+\Phi_n(X_n,\pi^\ast_n(X_n)).\] The Realised Value method calculates the value function as \[\hat{V}_n(X_n)=R_n(X_n,\pi^\ast(X_n))+\beta\hat{V}_{n+1}(X_{n+1}),\] and requires the recomputation of the sample paths from \(t + 1\) to \(T\) after each backwards-in-time iteration, as the optimal control affects the future state variables, hence changes the simulated paths [51]. Due to recalculation along the simulated path for each iteration, the computational complexity of the realised value method grows quadratically in the number of time steps compared with the linear growth for the regression surface method. However, the realised value method tends to accumulate much less regression errors over time and, from experience, this method is much more suitable for problems prone to regression errors when the number of time periods increases.

We note that in the forward simulation Algorithm 1, the state variables and controls at time \(t\) are sampled from problem specific distribution and then the state variables at \(t+1\) are simulated according to the transition function (2 ). That is, one does not need to simulate full trajectories of state variables according to the transition function (2 ) for \(n=0,...,N\). Instead, only samples from \(t\) to \(t+1\) are needed for regression in the corresponding backward solution. In some applications such as [45], [46], it appears critical; although in our illustrative GMWB examples, we did simulate full trajectories for \(n=0,...,N\).

Finally, we want to mention that the use of post-decision variable may improve efficiency of the method if the transition equation 2 can be separated into two transitions: \[\label{equation-post-decision} X_{t+1} = \widetilde{\mathcal{T}}_t(F(X_t, \pi_t), Z_{t+1}),\tag{8}\] where the deterministic transition to the post-decision variable \(\widehat{X}_{t} = F(X_{t}, \pi_{t})\) precedes the transition \(X_{t+1} = \widetilde{\mathcal{T}}_t(\hat{X}_{t}, Z_{t+1})\). This allows the conditional expectation 6 to be simplified as: \[\mathbb{E}\big[\beta_{n,n+1} V_{t+1}(X_{t+1})\, \big|\, X_{t}; \pi_{t} \big] = \mathbb{E}\big[\beta_{n,n+1} V_{t+1}(X_{t+1})\, \big|\, \widehat{X}_{t} \big].\] This method offers two main advantages: (1) dimension reduction in the covariates needed for the least-squares approximation of the conditional expectation, and (2) an increase in sampling efficiency by sampling only the post-decision states \(\hat{X}_{t}\) rather than both \(X_{t}\) and \(\pi_{t}\).

3.2 Upper and lower estimators↩︎

All LSMC methods presented above, Algorithms 2-5, estimate the optimal policy \(\widehat\pi_n^\ast(x)\), \(n=1,...,N-1\). The contract value for a given policy is given by 3 . If we replace the optimal policy by its estimate and calculate the expectation by averaging \(M\) contract payoffs over independent Monte Carlo trajectories then, by definition of supremum, the obtained contract value estimator will be a lower biased estimator \[\label{MC95lower95estimator95eq} \widehat{V}_0^L\left({X}_0\right)=\frac{1}{M}\sum_{m=1}^M H_0({X^m},\hat{{\pi}}^\ast).\tag{9}\] Realised Value Algorithms 2 and 4 already calculate this lower estimator if we set \(X_0^m=X_0, m=1,...,M\) and find \(\widehat{V}^L_0=\frac{1}{M}\sum_m \hat{V}_0(\hat{X}_0^m)\). In the case of Regression Surface Algorithms 3 and 5, formula 9 should be implemented once the optimal policy has been estimated. This regression surface estimator \(\widehat{V}_0(X_0)\) calculated by Algorithms 3 and 5 can be shown to be an approximate upper estimator \(\hat{V}_0^U(X_0)\), by using Jensen’s inequality to swap the order of the expectation and supremum operators iteratively, see [52], [33].

4 Numerical study↩︎

In this section we present LSMC results for GMWB prices under the optimal policyholder withdrawal strategies, compared with some established accurate finite-difference results obtained with a fine mesh and small time steps. We show results for Regress Now LSMC methods only. It is possible to obtain results for Regress Later LSMC methods for models considered in our study, but these methods cannot be generalised to more complicated models where numerical integration is required for conditional expectations in which more than two stochastic state variables are involved. Also, in the case of neural network regression, conditional expectations would have to be calculated numerically regardless of the model complexity. For the benefit of the readers, 6.2 presents closed-form formulas for the moments of the state variables up to order three in the case of the Vasicek stochastic interest rate model that can be used for Regress Later LSMC, but we do not pursue numerical results here. We mention that we tried Regress Later LSMC with polynomial regression for a simple model with constant interest rate, but we found that it did not perform as well as Regress Now methods.

4.1 LSMC settings↩︎

For LSMC calculations we use the following settings for polynomial and neural network regressions.

Polynomial regression

1. For forward simulation we just use a standard forward simulation in Algorithm 1, i.e. line 3 is replaced with \(X_t^m=\widetilde{X}_t^m\) and \(X_0=\widetilde{X}_0=x_0\).

2. We do not perform any transformation of value functions before regression.

3. Basis functions \(\mathbf{L}(.)\) for Algorithms 2 and 3:
   In cases of constant interest rate, polynomials up to the 3rd order of \(W, A, \pi\) are used, i.e. every term in the polynomial function takes the form \(W^{j_w}A^{j_a}\pi^{j_\pi}\), where \(j_w,j_a,j_\pi\) are all non-negative integers and \(j_w+j_a+j_\pi \leq 3\). Initially all possible combinations of \(j_w,\;j_a,\;j_\pi\) were used, i.e. the complete 3rd order three-dimensional polynomial function was used. Later, after some numerical experiments the terms \(W^2\pi\), \(A^2\pi\) and \(\pi^3\) were removed from the set, because they did not improve the accuracy of results. For cases with stochastic interest rate, the basis functions are described in Section 4.4.

Neural network regression settings
The neural network uses a fully connected feed forward architecture with input dimension \(D_{\text{in}} = 2\) (post-decision state variables) for constant interest rates and \(D_{\text{in}} = 3\) for stochastic interest rates, three hidden layers each of width \(H = 128\) with SiLU activation (Sigmoid Linear Unit, \(SiLU(x)=x/(1+e^{-x})\)), and output dimension \(D_{\text{out}} = 1\). Training is performed for \(2{,}000\) epochs using the Adam optimizer with learning rate \(0.001\), weight decay \(10^{-5}\), and cosine annealing schedule. Optimal withdrawals are computed by grid search with \(1{,}000\) discretization points.

4.2 Model and VA rider specification↩︎

For numerical illustration we consider the joint dynamics for the reference portfolio of assets \(S(t)\), e.g. a mutual fund, underlying the contract and the stochastic interest rate \(r(t)\), under the risk-neutral probability measure \(\mathbb{Q}\), governed by \[\begin{align} \label{referenceportfolio95eq} \begin{aligned} \frac{dS(t)}{S(t)}&=r(t) dt +\sigma_S \left(\rho d\mathcal{B}_1(t) +\sqrt{1-\rho^2} d\mathcal{B}_2(t)\right),\\ dr(t)&=\kappa(\theta-r(t))dt+\sigma_r d\mathcal{B}_1(t). \end{aligned} \end{align}\tag{10}\] Here, \(\mathcal{B}_1(t)\) and \(\mathcal{B}_2(t)\) are independent standard Wiener processes, \(\rho\) is the correlation coefficient between \(S(t)\) and \(r(t)\) processes, and \(\sigma_S\) is the asset volatility parameter. The process for the interest rate \(r(t)\) is the well-known Vasicek model with constant parameters \(\kappa\), \(\theta\) and \(\sigma_r\). For simplicity of notation, we assume that model parameters are constant in time though the results can be generalized to the case of time dependent parameters.

For this stochastic interest rate model, the price of a zero coupon bond \(P(t,T)\) at time \(t\) with maturity \(T\), can be found in closed-form and then a change of numéraire can simplify the calculation of expectations; see [49]. We will not pursue this in this paper but rather consider a brute force LSMC because this specific change of numeraire trick may not work for other interest rate models.

Consider the following VA contract with a basic GMWB often used in research studies, which is convenient for benchmarking. The actual products may have extra features, but these can be easily incorporated in the model and numerical algorithms developed in this paper.

• The premium paid by the policyholder upfront at \(t_0\) is invested into the reference portfolio/risky asset \(S(t)\). The value of this portfolio (hereafter referred to as wealth account) at time \(t\) is denoted as \(W(t)\), so that the upfront premium paid by the policyholder is \(W(0)\). GMWB guarantees the return of the premium via the withdrawals \(\pi_n\ge 0\) allowed at times \(t_n\), \(n=1,2,\ldots,N\). The total amount of withdrawals cannot exceed the guarantee \(W(0)\), and withdrawals can be different from the contractual (guaranteed) withdrawal \(G_n=W(0)(t_n-t_{n-1})/T\), with penalties imposed if \(\pi_n>G_n\). Denote the annual contractual rate as \(g:=1/T\). Then, the wealth account \(W(t)\) evolves as \[\label{eq95Wt} W(t_{n+1})=\max\left(W(t_n)-\pi_n,0\right)\frac{S(t_{n+1})}{S(t_{n})} e^{-\alpha \Delta_{n+1}},\\\tag{11}\]

  where \(\Delta_n=t_n-t_{n-1}\) and \(\alpha\) is the annual fee continuously charged by the contract issuer. If the account balance becomes zero or negative, then it will stay zero till maturity.

• Denote the value of the contract guarantee at time \(t\) as \(A(t)\), hereafter referred to as guarantee account, with \(A(0)=W(0)\). The guarantee balance evolves as \[\label{accountbalance95eq} A(t_{n+1})=A(t_n)-\pi_n\;\; n=1,2,\ldots,N-1\tag{12}\] with \(\pi_n\le A(t_{n})\)

• The cashflow received by the policyholder at the withdrawal time \(t_n\) is given by \[C_n(\pi_n)=\left\{\begin{array}{ll} \pi_n, & \mathrm{if}\; 0\le \pi_n\le G_n, \\ G_n+(1-\beta)(\pi_n-G_n), & \mathrm{if}\; \pi_n>G_n, \end{array} \right.\] where \(G_n\) is the contractual withdrawal and \(\beta\in [0,1]\) is the penalty coefficient applied to the portion of withdrawal above \(G_n\).

• Let \(V_n(W,r,A)\) be a price of the VA contract with GMWB at time \(t_n\), when \(W(t)=W\), \(r(t)=r\), \(A(t)=A\). At maturity, the policyholder takes the maximum between the remaining guarantee account net of penalty charge and the remaining balance of the wealth account, i.e. the final payoff is \[\label{finalcond95eq} V_{N}(W,r, A)=\max\left(W,C_N(A)\right).\tag{13}\]

During the contract, the policyholder receives the cashflows \(C_n(\pi_n)\), \(n=1,2,\ldots,N-1\) and the final payoff at maturity. Given the withdrawal strategy \(\boldsymbol{\pi}=(\pi_1,\ldots,\pi_{N-1})\), the present value of the total contract payoff is given by 1 with settings \(X_t=(W(t),r(t),A(t))\), \(R_N(X_n):=\max\left(W,C_N(A)\right)\), and \(R_n(X_n,\pi_n):=C_n(\pi_n)\).

4.3 GMWB pricing results for constant interest rate↩︎

In this section, we consider the case of constant interest rate and GMWB with parameter settings: \(S(0)=1.0\), \(r = 5\%\), \(g=10\%\) (\(T:=1/g=10\) years), \(\beta=10\%\), \(N_w=1\), \(\alpha = 0.0135\). All LSMC price results are averages over 20 independent LSMC runs, and the corresponding standard errors are given in brackets next to the price estimates. Table 1 shows price estimates for a 10-year GMWB contract with annual withdraw frequency at different volatility levels in the case of optimal withdrawal strategy. The reported results are from LSMC Regress Now/Realised Value Algorithm 2 and Regress Now/Regression Surface Algorithm 3, both in the case of polynomial OLS regression and neural network regression.

For the benchmark results in Table 1, we use results from the finite difference method (column “FD”), where we have used 800 node points for \(W\), 200 node points for \(A\), and the time step size \(dt=0.0125\), which is the typical setting to get very accurate results, as shown, e.g., in [13]. Because we do not know the exact solution, FD results will be used as the “benchmark” solution with which all the other numerical results will be compared to estimate their accuracy. To estimate the accuracy of FD results, we calculated FD results when the number of nodes for both \(W\) and \(A\) is half of the fine mesh, i.e. 400 node points for \(W\), 100 node points for \(A\), and the time step is halved to \(dt=0.025\). It can be quantitatively shown (as described in, e.g., [53]) that in general the difference between FD results based on the coarse mesh and the results based on the finer mesh with doubled number of nodes, can be used as an accuracy measure of the FD results of the coarser mesh. Thus, the reported FD error is conservative as it is reported for the fine mesh in our table.

In all LSMC methods, the randomized distribution of withdraw amount in the forward simulation is a combination of discrete uniform and continuous uniform. Specifically, we give \(25\%\) probability to each of no withdraw (\(dW=0\)) and withdraw at the contractual rate (\(dW=G\)), and \(50\%\) probability to the withdraw amount \(0 < dW < A(t)\) (uniform and continuous). This randomization gave the best results in our experiments.

As can be seen, comparing LSMC results with the FD results, the regression surface Algorithm 3 upper estimator \(V^U\) obtained using a 3rd order polynomial regression performs a bit better than in the case of NN regression. Lower estimators from the NN and OLS regressions are very close.

Table 2 shows LSMC results as the number of trajectories \(M\) increases from \(10^4\) to \(4\times 10^6\) indicating that one needs \(M>10^5\) trajectories to get at least 1% accuracy. It also shows that realised value LSMC works better when the number of trajectories is small, and is outperformed by regression surface LSMC when \(M\ge 10^6\). We also note that in this table we show results both for \(V^L\) and \(V^U\) estimators from Algorithm 3. Comparing to regression surface, \(V^L\) for realised value performs better for small \(M\) and the same for large \(M\) and thus in all other tables we show \(V^L\) from realised value LSMC Algorithm 2 only.

Overall, in the reported and other results, we observed that \(V^L\) is lower and \(V^U\) is larger than the “exact” result, and these bounds converge towards the “exact” value as \(M\) increases. As expected, the standard errors of the estimators decrease as \(M\) increases. Also, as expected, the value of the contract and its standard error increase as volatility \(\sigma\) increases.

We also tried post decision variables 8 OLS regressions with polynomials up to order 3, but the obtained accuracy was not as good as for OLS regressions reported in Table 1. This is another highlight of the NN advantage that it does not require such feature engineering.

Finally we note that the computing time in the case of the polynomial regression results reported in Table 1 is approximately 40 minutes over 20 runs (i.e. about 2 minutes per price calculation) on a modern laptop (Intel(R) Core(TM) Ultra 7 255U, 2.00 GHz) without parallel threading using Intel Fortran. Of course, the actual total computational effort is substantially larger, as it also involves experimentation with different basis functions. The computing time in the case of the neural network regression results reported in Table 1 is approximately 3 hours over 20 runs (i.e. about 9 minutes per price calculation) using PyTorch with GPU acceleration running on an NVIDIA H20 GPU.

4.4 GMWB pricing results for stochastic interest rate↩︎

In this section we show GMWB pricing results in the case of the Vasicek stochastic interest rate model. The simulation of the underlying model 10 is done using the exact simulation scheme described in 6.1. Results obtained by LSMC are compared with those obtained by GHQC (Gauss-Hermite quadrature on cubic spline) developed in [54] and used in [49] for pricing GMWB in the case of stochastic interest rate². Unlike the previous section where some very fine finite-difference results were used as a benchmark, here we do not have high accuracy finite-difference results with very fine meshes due to the higher dimension of the problem.

In Table [tab95vasicek1], LSMC pricing results (column 3 to column 6) for a 10-year GMWB contract with annual withdraw frequency for different values of volatility \(\sigma\) under optimal dynamic withdrawals are shown for the OLS and NN regressions, in comparison with GHQC results (column 2). For LSMC, we show numerical results for price estimates (obtained as an average over 100 independent runs of the LSMC calculations) and their corresponding standard errors. All the input parameters are given in the table caption. The numerical setup (e.g. number of quadrature points) for GHQC is exactly the same as in [49]. For OLS LSMC calculations we used the complete set of 2nd order polynomial basis functions, that is, all the terms in the form of \(W^{j_w}A^{j_a}\pi^{j_\pi} r^{j_r}\) are included, where \(j_w,j_a,j_\pi,j_r\) are all non-negative integers and \(j_w+j_a+j_\pi +j_r \leq 2\), totalling 15 terms altogether including the constant term.

We also tried other sets of polynomial basis functions. For example, we considered the set with terms in the form \(W^{j_w}A^{j_a}\pi^{j_\pi}\) (without containing interest rate \(r\)), the same as in the case of constant interest rate used in the previous section (e.g. Table 1), and all the extra terms involving \(r\) have the form \(W^{j_w}A^{j_a}\pi^{j_\pi} r^{j_r}\) with \(1\leq j_w+j_a+j_\pi +j_r \leq 2\) and \(1\leq j_r \leq 2\). In other words, the extra terms involving \(r\) are \(rW\), \(rA\),\(r\pi\), \(r\) and \(r^2\). Altogether in this set there were 23 terms. Results obtained with this set are not materially different from the reported ones (i.e. the difference is well within standard errors). We also tried to add all the cubic terms involving \(r\) to set B, i.e. adding all possible terms of the form \(W^{j_w}A^{j_a}\pi^{j_\pi} r^{j_r}\), with the constraints \(1\leq j_r \leq 3\) and \(j_w+j_a+j_\pi +j_r =3\). This set has altogether 30 polynomial terms. Again, we did not observe material difference in results.

Taking the GHQC results as fairly accurate with relative error of the order of \(10^{-4}\) (see [18] and [49]), we can measure the performance of LSMC. From Table [tab95vasicek1], we can see that the OLS Algorithm 3 upper estimator \(V^U\) appears above the ‘true’ value and the OLS Algorithm 2 lower estimator \(V^L\) appears below the ‘true’ value. The upper and lower NN LSMC estimates also appear to be consistent with ‘true’ values, though OLS estimators appear slightly more accurate. Also, we can see that the standard errors of OLS LSMC estimators are now one order of magnitude larger than in the cases with constant interest rate shown in the previous section. By contrast, the standard errors of NN estimators are similar to the constant interest rate case. The main reason for this phenomenon is that the neural network architecture scales gracefully with the dimensionality of the state space: extending the model from the constant interest rate case (\(D_{\text{in}} = 2\)) to the stochastic interest rate case (\(D_{\text{in}} = 3\)) increases the total number of parameters only marginally, from \(33{,}537\) to \(33{,}665\) (an increase of \(128\) parameters, or about \(0.4\%\)), since the additional input dimension only affects the weight matrix of the first hidden layer (\(D_{\text{in}} \times H\)). This is in sharp contrast to polynomial regression, where the number of basis functions grows combinatorially with the number of state variables, as illustrated above by the expansion from the original basis set to the 15-term set when including interest rates. The NN approach therefore generalizes to higher-dimensional state spaces with essentially no additional model complexity, leading to similar standard errors, and also without the need for manual feature engineering.

5 Conclusion↩︎

In this paper, we demonstrated that the numerical evaluation of VAs under optimal withdrawal strategies can be successfully accomplished using the LSMC method. The main advantages of the method are its flexibility with respect to the underlying stochastic processes, its straightforward extension to additional state and control variables, and its ability to handle complex contract payoff structures. In our experiments, polynomial regression performed slightly better than NN regression. However, polynomial regression required substantially more preliminary analysis to identify suitable basis functions. Overall, we recommend the use of NN regression because, due to its “non-parametric” nature, it can adapt to the shape of the target function without requiring the functional form to be specified in advance. A similar observation was made in macroeconomic models studied by [55] where NN regression was advocated over polynomial regression for approximating conditional expectations in macro-finance models. Moreover, making interest rates stochastic did not cause any trouble to the deep LSMC method, while this caused the standard errors of classical polynomial LSMC to increase by one order of magnitude.

We presented and discussed different variants of LSMC: regress now/regress later, regression surface/realised value. It is important to note that LSMC estimators are downward or upward biased depending on the method, and it is prudent to report both the upper and lower estimators. The difference between the upper and lower estimators is significantly larger than standard errors of the estimators obtained over independent runs.

Regression-based approximations of conditional expectations have long been used in computational economics through forward-in-time simulation combined with regression to evaluate expectations arising in equilibrium conditions [56]. This approach was later extended through the use of neural networks by [57], while [55] revisited related ideas from a modern machine-learning perspective in macro-finance models. In the case of LSMC, we solve a recursive stochastic control problem backward in time by approximating the conditional expectation in the Bellman equation via regression, and then recover the optimal control though a separate numerical maximization. This separation is particularly convenient for VAs because the controls are constrained and the state dynamics is endogenous.

Finally, we mention that we attempted to price GMWB using a global neural network parameterisation of controls with time as an additional input variable and then estimating the NN parameters by maximizing the empirical objective function directly by gradient ascent. While easier to implement, this approach did not reach the same accuracy as LSMC. We therefore leave the exploration of this approach for further research.

Acknowledgement↩︎

Nicolas Langrené acknowledges the support of the Guangdong Provincial/Zhuhai Key Laboratory of IRADS (2022B1212010006).

6 Closed-form formulas↩︎

We consider the model \[\begin{align} dS(t) & =r(t)S(t)dt+\sigma_{S}S(t)(\rho d\mathcal{B}_{1}(t)+\sqrt{1-\rho^{2}}d\mathcal{B}_{2}(t))\tag{14}\\ dr(t) & =\kappa(\theta-r(t))dt+\sigma_{r}d\mathcal{B}_{1}(t)\tag{15} \end{align}\]

6.1 Exact simulation scheme↩︎

The strong solutions of \(S(t)\) and \(r(t)\) are given by: \[\begin{align} S(t) & =S(0)\exp\left(\int^{t}_{0}\left(r(\tau)-\frac{1}{2}\sigma^{2}_{S}\right)d\tau+\int^{t}_{0}\sigma_{S}\left(\rho d\mathcal{B}_{1}(\tau)+\sqrt{1-\rho^{2}}d\mathcal{B}_{2}(\tau)\right)\right)\\ r(t) & =e^{-\kappa t}\left(r(0)+\int^{t}_{0}e^{\kappa\tau}\kappa\theta d\tau+\int^{t}_{0}e^{\kappa\tau}\sigma_{r}d\mathcal{B}_{1}(\tau)\right) \end{align}\] Similarly, between times \(t_{i}\) and \(t_{i+1}\) with \(0\leq t_{i}<t_{i+1}\): \[\begin{align} S(t_{i+1}) & =S(t_{i})\exp\left(\int^{t_{i+1}}_{t_{i}}r(\tau)d\tau-\frac{1}{2}\sigma^{2}_{S}\Delta_{t_{i}}+\int^{t_{i+1}}_{t_{i}}\sigma_{S}\left(\rho d\mathcal{B}_{1}(\tau)+\sqrt{1-\rho^{2}}d\mathcal{B}_{2}(\tau)\right)\right)\\ r(t_{i+1}) & =e^{-\kappa\Delta_{t_{i}}}\left(r(t_{i})+\int^{t_{i+1}}_{t_{i}}e^{\kappa(\tau-t_{i})}\kappa\theta d\tau+\int^{t_{i+1}}_{t_{i}}e^{\kappa(\tau-t_{i})}\sigma_{r}d\mathcal{B}_{1}(\tau)\right) \end{align}\] where \(\Delta_{t_{i}}:=t_{i+1}-t_{i}\). Introduce the following three random variables: \[\begin{align} R:= & r(t_{i+1})-e^{-\kappa\Delta_{t_{i}}}r(t_{i})=e^{-\kappa\Delta_{t_{i}}}\left(\int^{t_{i+1}}_{t_{i}}e^{\kappa(\tau-t_{i})}\kappa\theta d\tau+\int^{t_{i+1}}_{t_{i}}e^{\kappa(\tau-t_{i})}\sigma_{r}d\mathcal{B}_{1}(\tau)\right)\tag{16}\\ Y:= & \int^{t_{i+1}}_{t_{i}}\!(r(\tau)-e^{-\kappa\Delta_{t_{i}}}r(t_{i}))d\tau=\int^{t_{i+1}}_{t_{i}}\!\!e^{-\kappa(u-t_{i})}\left(\int^{u}_{t_{i}}\!e^{\kappa(\tau-t_{i})}\kappa\theta d\tau+\!\int^{u}_{t_{i}}\!e^{\kappa(\tau-t_{i})}\sigma_{r}d\mathcal{B}_{1}(\tau)\right)du\tag{17}\\ E:= & -\frac{1}{2}\sigma^{2}_{S}\Delta_{t_{i}}+\int^{t_{i+1}}_{t_{i}}\sigma_{S}\left(\rho d\mathcal{B}_{1}(\tau)+\sqrt{1-\rho^{2}}d\mathcal{B}_{2}(\tau)\right)\tag{18} \end{align}\] These three random variables are normally distributed and are independent of \(\mathcal{F}_{t_{i}}\). Using these notations, we can write \(S(t_{i+1})\) and \(r(t_{i+1})\) as \[\begin{align} S(t_{i+1}) & =S(t_{i})\exp\left(e^{-\kappa\Delta_{t_{i}}}r(t_{i})\Delta_{t_{i}}+Y+E\right)\tag{19}\\ r(t_{i+1}) & =e^{-\kappa\Delta_{t_{i}}}r(t_{i})+R\tag{20} \end{align}\] where \(\Delta_{t_{i}}:=t_{i+1}-t_{i}\). To complete the definition of the exact simulation scheme 19 20 , we need to compute the mean, variance and correlations of the three Gaussian random variables \(R\), \(Y\) and \(E\). We will make use of the following classical lemma.

We can now compute the mean and variance of \(R\), \(Y\) and \(E\). The following formulas hold: \[\begin{align} R & =e^{-\kappa\Delta_{t_{i}}}\left(\int^{t_{i+1}}_{t_{i}}e^{\kappa(\tau-t_{i})}\kappa\theta d\tau+\int^{t_{i+1}}_{t_{i}}e^{\kappa(\tau-t_{i})}\sigma_{r}d\mathcal{B}_{1}(\tau)\right)=\mu_{R}+\varsigma_{R}N_{R}\label{eq:R95decomposition}\\ \mu_{R} & :=\mathbb{E}[R]=e^{-\kappa\Delta_{t_{i}}}\int^{t_{i+1}}_{t_{i}}e^{\kappa(\tau-t_{i})}\kappa\theta d\tau=\theta\left(1-e^{-\kappa\Delta_{t_{i}}}\right)\nonumber \\ \varsigma^{2}_{R} & :=\mathbb{V}\mathrm{ar}[R]=\frac{\sigma^{2}_{r}}{2\kappa}\left[1-e^{-2\kappa\Delta_{t_{i}}}\right]\nonumber \\ N_{R} & :=\varsigma^{-1}_{R}e^{-\kappa\Delta_{t_{i}}}\int^{t_{i+1}}_{t_{i}}e^{\kappa(\tau-t_{i})}\sigma_{r}d\mathcal{B}_{1}(\tau)\;\sim\;\mathcal{N}(0,1)\nonumber \end{align}\tag{21}\]

\[\begin{align} Y & =\int^{t_{i+1}}_{t_{i}}e^{-\kappa(u-t_{i})}\left(\int^{u}_{t_{i}}e^{\kappa(\tau-t_{i})}\kappa\theta d\tau+\int^{u}_{t_{i}}e^{\kappa(\tau-t_{i})}\sigma_{r}d\mathcal{B}_{1}(\tau)\right)du=\mu_{Y}+\varsigma_{Y}N_{Y}\label{eq:Y95decomposition}\\ \mu_{Y} & :=\mathbb{E}[Y]=\int^{t_{i+1}}_{t_{i}}e^{-\kappa(u-t_{i})}\left(\int^{u}_{t_{i}}e^{\kappa(\tau-t_{i})}\kappa\theta d\tau\right)du\nonumber \\ & =\theta\left(\Delta_{t_{i}}-\frac{1-e^{-\kappa\Delta_{t_{i}}}}{\kappa}\right)\nonumber \\ \varsigma^{2}_{Y} & :=\mathbb{V}\mathrm{ar}[Y]=\frac{\sigma^{2}_{r}}{\kappa^{2}}\int^{t_{i+1}}_{t_{i}}\left(1-e^{-\kappa(\tau-t_{i})}\right)^{2}d\tau\nonumber \\ & =\frac{\sigma^{2}_{r}}{\kappa^{2}}\left[\Delta_{t_{i}}-\frac{2}{\kappa}\left(1-e^{-\kappa\Delta_{t_{i}}}\right)+\frac{1}{2\kappa}\left(1-e^{-2\kappa\Delta_{t_{i}}}\right)\right]\nonumber \\ N_{Y} & :=\varsigma^{-1}_{Y}\int^{t_{i+1}}_{t_{i}}e^{-\kappa(u-t_{i})}\int^{u}_{t_{i}}e^{\kappa(\tau-t_{i})}\sigma_{r}d\mathcal{B}_{1}(\tau)du\;\sim\;\mathcal{N}(0,1)\nonumber \end{align}\tag{22}\]

\[\begin{align} E & =-\frac{1}{2}\sigma^{2}_{S}\Delta_{t_{i}}+\int^{t_{i+1}}_{t_{i}}\sigma_{S}\left(\rho d\mathcal{B}_{1}(\tau)+\sqrt{1-\rho^{2}}d\mathcal{B}_{2}(\tau)\right)=\mu_{E}+\varsigma_{E}N_{E}\label{eq:E95decomposition}\\ \mu_{E} & :=\mathbb{E}[E]=-\frac{1}{2}\sigma^{2}_{S}\Delta_{t_{i}}\nonumber \\ \varsigma^{2}_{E} & :=\mathbb{V}\mathrm{ar}[E]=\sigma^{2}_{S}\Delta_{t_{i}}\nonumber \\ N_{E} & =\varsigma^{-1}_{E}\int^{t_{i+1}}_{t_{i}}\sigma_{S}\left(\rho d\mathcal{B}_{1}(\tau)+\sqrt{1-\rho^{2}}d\mathcal{B}_{2}(\tau)\right)\;\sim\;\mathcal{N}(0,1)\nonumber \end{align}\tag{23}\] So far, we established these three Gaussian decompositions 21 22 23 : \[\begin{align} R & =\mu_{R}+\varsigma_{R}N_{R}\tag{24}\\ Y & =\mu_{Y}+\varsigma_{Y}N_{Y}\tag{25}\\ E & =\mu_{E}+\varsigma_{E}N_{E}\tag{26} \end{align}\] The correlations between the three standard Gaussian random variables \(N_{R}\), \(N_{Y}\) and \(N_{E}\) are given by: \[\begin{align} \rho_{RY}:=\mathbb{E}[N_{R}N_{Y}] & =\frac{\sigma^{2}_{r}}{\varsigma_{r}\varsigma_{Y}}\frac{\left(1-e^{-\kappa\Delta_{t_{i}}}\right)^{2}}{2\kappa^{2}}\tag{27}\\ \rho_{RE}:=\mathbb{E}[N_{R}N_{E}] & =\frac{\rho\sigma_{r}\sigma_{S}}{\kappa\varsigma_{r}\varsigma_{E}}\left(1-e^{-\kappa\Delta_{t_{i}}}\right)\tag{28}\\ \rho_{YE}:=\mathbb{E}[N_{Y}N_{E}] & =\frac{\rho\sigma_{r}\sigma_{S}}{\kappa\varsigma_{Y}\varsigma_{E}}\left[\Delta_{t_{i}}-\frac{1}{\kappa}\left(1-e^{-\kappa\Delta_{t_{i}}}\right)\right]\tag{29} \end{align}\] To sum up, the model 14 15 can be simulated exactly on a discrete-time grid \(0=t_{0}\leq t_{1}<\ldots<t_{N}=T\) by \[\begin{align} S(t_{i+1}) & =S(t_{i})\exp\left(e^{-\kappa\Delta_{t_{i}}}r(t_{i})\Delta_{t_{i}}+\mu_{Y}+\mu_{E}+\varsigma_{Y}N_{Y}+\varsigma_{E}N_{E}\right)\tag{30}\\ r(t_{i+1}) & =e^{-\kappa\Delta_{t_{i}}}r(t_{i})+\mu_{R}+\varsigma_{R}N_{R}\tag{31} \end{align}\] where the standard Gaussian random variables \(N_{Y}\), \(N_{E}\) and \(N_{R}\) are correlated according to the formulas 27 28 29 . This can be further simplified into \[\begin{align} S(t_{i+1}) & =S(t_{i})\exp\left(e^{-\kappa\Delta_{t_{i}}}r(t_{i})\Delta_{t_{i}}+\mu_{Y}+\mu_{E}+\sqrt{\varsigma^{2}_{Y}+\varsigma^{2}_{E}+2\varsigma_{Y}\varsigma_{E}\rho_{YE}}N_{S}\right)\tag{32}\\ r(t_{i+1}) & =e^{-\kappa\Delta_{t_{i}}}r(t_{i})+\mu_{R}+\varsigma_{R}N_{R}\tag{33}\\ \rho_{RS} & =\mathbb{E}[N_{R}N_{S}]=\frac{\varsigma_{Y}\rho_{RY}+\varsigma_{E}\rho_{RE}}{\sqrt{\varsigma^{2}_{Y}+\varsigma^{2}_{E}+2\varsigma_{Y}\varsigma_{E}\rho_{YE}}}\tag{34} \end{align}\] where \(N_{R}\) and \(N_{S}\) are two standard Gaussian random variables with correlation \(\rho_{RS}\).

6.2 Regress-later formulas↩︎

6.2.1 Formulas for general time step↩︎

For the Regress-Later implementation with third-order polynomial basis, we need to compute the discounted conditional moments \(\mathbb{E}[\exp\left(-\int^{t_{i+1}}_{t_{i}}r(\tau)d\tau\right)S(t_{i+1})^{p}r(t_{i+1})^{q}\left|\mathcal{F}_{t_{i}}\right.]\) with \(0\leq p\leq3\), \(0\leq q\leq3\), \(0<p+q\leq3\). We introduce \(\ell\in\mathbb{R}\) and consider the slightly more general problem of computing \[B_{\ell,p,q}:=\mathbb{E}[\exp\left(-\ell\int^{t_{i+1}}_{t_{i}}r(\tau)d\tau\right)S(t_{i+1})^{p}r(t_{i+1})^{q}\left|\mathcal{F}_{t_{i}}\right.]\label{eq:B}\tag{35}\] where we are interested in the case \(\ell=1\), but the computations easily accommodate any other value of \(\ell\), including \(\ell=0\) (no discounting). Using equations 16 17 18 , \(B_{\ell,p,q}\) can be written as \[B_{\ell,p,q}=S(t_{i})^{p}\exp\left((p-\ell)e^{-\kappa\Delta_{t_{i}}}r(t_{i})\Delta_{t_{i}}\right)\mathbb{E}\left[\exp\left((p-\ell)Y+pE\right)\left(e^{-\kappa\Delta_{t_{i}}}r(t_{i})+R\right)^{q}\left|\mathcal{F}_{t_{i}}\right.\right]\label{eq:B95decomposition}\tag{36}\] To compute the remaining conditional expectation in equation 36 , we need to use the formulas 24 , 25 , 26 , 27 , 28 , and 29 for \(R\), \(Y\), and \(E\), as well as the following lemma:

\[\begin{align} C_{\ell,p,q} & =S(t_{i})^{p}\exp\left((p-\ell)(\mu_{Y}+e^{-\kappa\Delta_{t_{i}}}r(t_{i})\Delta_{t_{i}})+p\mu_{E}\right)\nonumber \\ & \times\exp\left(\frac{1}{2}\left((p-\ell)^{2}\varsigma^{2}_{Y}\left(1-\rho^{2}_{RY}\right)+p^{2}\varsigma^{2}_{E}\left(1-\rho^{2}_{RE}\right)+2p(p-\ell)\varsigma_{Y}\varsigma_{E}\left(\rho_{YE}-\rho_{RY}\rho_{RE}\right)\right)\right)\nonumber \\ & \times\mathcal{E}((p-\ell)\varsigma_{Y}\rho_{RY}+p\varsigma_{E}\rho_{RE}\,,\,q)\label{eq:Clpq95explicit} \end{align}\tag{37}\] where \(\mathcal{E}((p-\ell)\varsigma_{Y}\rho_{RY}+p\varsigma_{E}\rho_{RE}\,,\,q)\) is given explicitly in Lemma 2. This can be further simplified by noticing that \(\mathcal{E}(c,q)\) is always proportional to \(e^{\frac{c^{2}}{2}}\), and that \(e^{-\frac{c^{2}}{2}}\mathcal{E}(c,q)\) is a polynomial in \(c\). Define \[\mathcal{P}(c,q):=e^{-\frac{c^{2}}{2}}\mathcal{E}(c,q)\label{eq:Pcq}\tag{38}\] From Lemma 2, we know that \[\begin{align} \mathcal{P}(c,0) & =1\\ \mathcal{P}(c,1) & =c\\ \mathcal{P}(c,2) & =c^{2}+1\\ \mathcal{P}(c,3) & =c^{3}+3c \end{align}\] Using the notation 38 , \(C_{\ell,p,q}\) is given by \[\begin{align} C_{\ell,p,q} & =S(t_{i})^{p}\exp\left((p-\ell)(\mu_{Y}+e^{-\kappa\Delta_{t_{i}}}r(t_{i})\Delta_{t_{i}})+p\mu_{E}\right)\nonumber \\ & \times\exp\left(\frac{1}{2}\left((p-\ell)^{2}\varsigma^{2}_{Y}+p^{2}\varsigma^{2}_{E}+2p(p-\ell)\varsigma_{Y}\varsigma_{E}\rho_{YE}\right)\right)\nonumber \\ & \times\mathcal{P}\left((p-\ell)\varsigma_{Y}\rho_{RY}+p\varsigma_{E}\rho_{RE}\,,\,q\right)\label{eq:Clpq95explicit952} \end{align}\tag{39}\] By plugging this formula into equations ?? ?? ?? ?? , we obtain that the conditional expectations \(B_{\ell,p,q}\) (equation 35 ) in the model 14 15 are given explicitly by \[\begin{align} B_{\ell,p,q} & =\mathbb{E}[\exp\left(-\ell\int^{t_{i+1}}_{t_{i}}r(\tau)d\tau\right)S(t_{i+1})^{p}r(t_{i+1})^{q}\left|\mathcal{F}_{t_{i}}\right.]\nonumber \\ & =S(t_{i})^{p}\exp\left((p-\ell)(\mu_{Y}+e^{-\kappa\Delta_{t_{i}}}r(t_{i})\Delta_{t_{i}})+p\mu_{E}\right)\nonumber \\ & \times\exp\left(\frac{1}{2}\left((p-\ell)^{2}\varsigma^{2}_{Y}+p^{2}\varsigma^{2}_{E}+2p(p-\ell)\varsigma_{Y}\varsigma_{E}\rho_{YE}\right)\right)\times D_{\ell,p,q}\label{eq:Blpq95final} \end{align}\tag{40}\] where \(D_{\ell,p,q}\) is given explicitly by \[\begin{align} D_{\ell,p,0} & =1\tag{41}\\ D_{\ell,p,1} & =e^{-\kappa\Delta_{t_{i}}}r(t_{i})+\mu_{R}+(p-\ell)\varsigma_{R}\varsigma_{Y}\rho_{RY}+p\varsigma_{R}\varsigma_{E}\rho_{RE}\tag{42}\\ D_{\ell,p,2} & =\left(e^{-\kappa\Delta_{t_{i}}}r(t_{i})+\mu_{R}+(p-\ell)\varsigma_{R}\varsigma_{Y}\rho_{RY}+p\varsigma_{R}\varsigma_{E}\rho_{RE}\right)^{2}+\varsigma^{2}_{R}\tag{43}\\ D_{\ell,p,3} & =\left(e^{-\kappa\Delta_{t_{i}}}r(t_{i})+\mu_{R}+(p-\ell)\varsigma_{R}\varsigma_{Y}\rho_{RY}+p\varsigma_{R}\varsigma_{E}\rho_{RE}\right)^{3}\nonumber \\ & +3\varsigma^{2}_{R}\left(e^{-\kappa\Delta_{t_{i}}}r(t_{i})+\mu_{R}+(p-\ell)\varsigma_{R}\varsigma_{Y}\rho_{RY}+p\varsigma_{R}\varsigma_{E}\rho_{RE}\right)\tag{44} \end{align}\] for the first few values of \(q\).

6.2.2 Closed-form formula for last time step↩︎

Consider the final time step \([t_{N-1},t_{N}]\). Let \(\Delta_{t_{N-1}}=t_{N}-t_{N-1}\). We have \[W(t_{N})=\max\left(W(t_{N-1})-\pi_{N-1},0\right)\frac{S(t_{N})}{S(t_{N-1})}e^{-\alpha\Delta_{t_{N-1}}}\] \[A(t_{N})=A(t_{N-1})-\pi_{N-1}\] \[C_{N}(A(t_{N}))=\begin{cases} A(t_{N}) & \mathrm{if\;}0\leq A(t_{N})\leq G_{N}\\ G_{N}+(1-\beta)(A(t_{N})-G_{N}) & \mathrm{if\;}A(t_{N})>G_{N} \end{cases}\] The final payoff is \[\max(W(t_{N}),C_{N}(A(t_{N})))\] The conditional expectation of the discounted final payoff with respect to \(\mathcal{F}_{t_{N-1}}\) is given by \[\begin{align} & \mathbb{E}\left[\exp\left(-\ell\int^{t_{N}}_{t_{N-1}}r(\tau)d\tau\right)\max(W(t_{N}),C_{N}(A(t_{N})))\left|\mathcal{F}_{t_{i}}\right.\right]\\ = & \mathbb{E}\left[\exp\left(-\ell\int^{t_{N}}_{t_{N-1}}r(\tau)d\tau\right)\max\left(\max\left(W(t_{N-1})-\pi_{N-1},0\right)\frac{S(t_{N})}{S(t_{N-1})}e^{-\alpha\Delta_{t_{N-1}}},C_{N}(A(t_{N}))\right)\left|\mathcal{F}_{t_{i}}\right.\right]\\ = & \mathbb{E}\left[\exp\left(-\ell\left(Y+e^{-\kappa\Delta_{t_{N-1}}}r(t_{N-1})\Delta_{t_{N-1}}\right)\right)\right.\\ & \left.\times\max\left(\max\left(W(t_{N-1})-\pi_{N-1},0\right)e^{-\alpha\Delta_{t_{N-1}}}\exp\left(e^{-\kappa\Delta_{t_{N-1}}}r(t_{N-1})\Delta_{t_{N-1}}+Y+E\right),C_{N}(A(t_{N}))\right)\left|\mathcal{F}_{t_{i}}\right.\right] \end{align}\] In other words, it is of the form \[a\mathbb{E}\left[\exp\left(-\ell Y\right)\max\left(b\exp\left(Y+E\right),c\right)\right]\] with the \(\mathcal{F}_{t_{i}}\)- measurable parameters \[\begin{align} a & =\exp\left(-\ell e^{-\kappa\Delta_{t_{N-1}}}r(t_{N-1})\Delta_{t_{N-1}}\right)>0\\ b & =\max\left(W(t_{N-1})-\pi_{N-1},0\right)e^{-\alpha\Delta_{t_{N-1}}}\exp\left(e^{-\kappa\Delta_{t_{N-1}}}r(t_{N-1})\Delta_{t_{N-1}}\right)\geq0\\ c & =C_{N}(A(t_{N}))\geq0 \end{align}\]

6.2.3 Case \(b=0\)↩︎

If \(b=0\), then the final conditional expectation is given by \[\begin{align} & a\mathbb{E}\left[\max\left(b\exp\left((1-\ell)Y+E\right),c\exp\left(-\ell Y\right)\right)\right]\\ = & ac\mathbb{E}\left[\exp\left(-\ell(\mu_{Y}+\varsigma_{Y}N_{Y})\right)\right]\\ = & ac\exp\left(-\ell\mu_{Y}+\ell^{2}\varsigma^{2}_{Y}/2\right) \end{align}\]

6.2.4 Case \(c=0\)↩︎

If \(c=0\), then the final conditional expectation is given by \[\begin{align} & a\mathbb{E}\left[\max\left(b\exp\left((1-\ell)Y+E\right),c\exp\left(-\ell Y\right)\right)\right]\\ = & ab\mathbb{E}\left[\exp\left((1-\ell)Y+E\right)\right]\\ = & ab\exp\left((1-\ell)\mu_{Y}+\mu_{E}\right)\mathbb{E}\left[\exp\left((1-\ell)\varsigma_{Y}N_{Y}+\varsigma_{E}N_{E}\right)\right]\\ = & ab\exp\left((1-\ell)\mu_{Y}+\mu_{E}\right)\mathbb{E}\left[\exp\left(\sqrt{(1-\ell)^{2}\varsigma^{2}_{Y}+\varsigma^{2}_{E}+2(1-\ell)\rho_{YE}\varsigma_{Y}\varsigma_{E}}N_{YE}\right)\right]\\ = & ab\exp\left((1-\ell)\mu_{Y}+\mu_{E}+\left((1-\ell)^{2}\varsigma^{2}_{Y}+\varsigma^{2}_{E}+2(1-\ell)\rho_{YE}\varsigma_{Y}\varsigma_{E}\right)/2\right) \end{align}\] where we used equation ?? with \(p=1\).

6.2.5 Case \(b>0\) and \(c>0\)↩︎

In this case \[\begin{align} & a\mathbb{E}\left[\max\left(b\exp\left((1-\ell)Y+E\right),c\exp\left(-\ell Y\right)\right)\right]\\ = & a\mathbb{E}\left[\exp(\max\left(\log(b)+(1-\ell)Y+E\;,\;\log(c)-\ell Y\right))\right] \end{align}\] Define \[\begin{align} X_{1} & :=\log(b)+(1-\ell)Y+E\\ & =\log(b)+(1-\ell)\mu_{Y}+\mu_{E}+\sqrt{(1-\ell)^{2}\varsigma^{2}_{Y}+\varsigma^{2}_{E}+2(1-\ell)\rho_{YE}\varsigma_{Y}\varsigma_{E}}N_{YE}\\ & =m_{1}+\sigma_{1}N_{YE}\\ m_{1} & :=\log(b)+(1-\ell)\mu_{Y}+\mu_{E}\\ \sigma_{1} & :=\sqrt{(1-\ell)^{2}\varsigma^{2}_{Y}+\varsigma^{2}_{E}+2(1-\ell)\rho_{YE}\varsigma_{Y}\varsigma_{E}} \end{align}\] and \[\begin{align} X_{2} & :=\log(c)-\ell Y\\ & =\log(c)-\ell\mu_{Y}+\ell\varsigma_{Y}(-N_{Y})\\ & =m_{2}+\sigma_{2}(-N_{Y})\\ m_{2} & :=\log(c)-\ell\mu_{Y}\\ \sigma_{2} & :=\ell\varsigma_{Y} \end{align}\] Using equation ?? with \(p=1\), the correlation between the two standard Gaussian random variables \(N_{YE}\) and \(-N_{Y}\) is given by \[\begin{align} \rho_{12}:= & \mathbb{E}\left[N_{YE}(-N_{Y})\right]\\ & =\mathbb{E}\left[\frac{(1-\ell)\varsigma_{Y}N_{Y}+\varsigma_{E}N_{E}}{\sqrt{(1-\ell)^{2}\varsigma^{2}_{Y}+\varsigma^{2}_{E}+2(1-\ell)\rho_{YE}\varsigma_{Y}\varsigma_{E}}}(-N_{Y})\right]\\ & =\frac{-(1-\ell)\varsigma_{Y}\mathbb{E}\left[N_{Y}N_{Y}\right]-\varsigma_{E}\mathbb{E}\left[N_{E}N_{Y}\right]}{\sqrt{(1-\ell)^{2}\varsigma^{2}_{Y}+\varsigma^{2}_{E}+2(1-\ell)\rho_{YE}\varsigma_{Y}\varsigma_{E}}}\\ & =\frac{-(1-\ell)\varsigma_{Y}-\varsigma_{E}\rho_{YE}}{\sqrt{(1-\ell)^{2}\varsigma^{2}_{Y}+\varsigma^{2}_{E}+2(1-\ell)\rho_{YE}\varsigma_{Y}\varsigma_{E}}} \end{align}\] To sum up, \[\begin{align} & a\mathbb{E}\left[\max\left(b\exp\left((1-\ell)Y+E\right),c\exp\left(-\ell Y\right)\right)\right]\\ = & a\mathbb{E}\left[\exp\left(\max\left(X_{1},X_{2}\right)\right)\right] \end{align}\] where \(X_{1}\sim\mathcal{N}(m_{1},\sigma_{1})\), \(X_{2}\sim\mathcal{N}(m_{2},\sigma_{2})\) and \(\mathbb{C}\mathrm{orr}(X_{1,}X_{2})=\rho_{12}\). To conclude, we use the closed-form formula for the moment generating function of the maximum between two correlated Gaussian random variables, which can be found for example in [58]: \[\begin{align} & \mathbb{E}\left[\exp\left(-\ell\int^{t_{N}}_{t_{N-1}}r(\tau)d\tau\right)\max(W(t_{N}),C_{N}(A(t_{N})))\left|\mathcal{F}_{t_{i}}\right.\right]\nonumber \\ = & a\mathbb{E}\left[\exp\left(-\ell Y\right)\max\left(b\exp\left(Y+E\right),c\right)\right]\nonumber \\ = & a\mathbb{E}\left[\exp\left(\max\left(X_{1},X_{2}\right)\right)\right]\nonumber \\ = & a\exp\!\left(m_{1}\!+\!\frac{\sigma^{2}_{1}}{2}\right)\!\times\!\Phi\!\left(\frac{m_{1}-m_{2}+\sigma^{2}_{1}-\rho_{12}\sigma_{1}\sigma_{2}}{\sqrt{\sigma^{2}_{1}+\sigma^{2}_{2}-2\rho_{12}\sigma_{1}\sigma_{2}}}\right)\!+a\exp\!\left(m_{2}\!+\!\frac{\sigma^{2}_{2}}{2}\right)\!\times\!\Phi\!\left(\frac{m_{2}-m_{1}+\sigma^{2}_{2}-\rho_{12}\sigma_{1}\sigma_{2}}{\sqrt{\sigma^{2}_{1}+\sigma^{2}_{2}-2\rho_{12}\sigma_{1}\sigma_{2}}}\right)\label{eq:final95payoff95conditional95expectation} \end{align}\tag{45}\] where \(\Phi\) is the cumulative distribution function of the standard Gaussian distribution.

References↩︎