A flexible lattice framework for valuing options on assets paying discrete dividends and variable annuities embedding GMWB riders

In a market where a stochastic interest rate component characterizes asset dynamics, we propose a flexible lattice framework to evaluate and manage options on equities paying discrete dividends and variable annuities presenting some provisions, like a guaranteed minimum withdrawal benefit. The framework is flexible in that it allows to combine financial and demographic risk, to embed in the contract early exercise features, and to choose the dynamics for interest rates and traded assets. A computational problem arises when each dividend (when valuing an option) or withdrawal (when valuing a variable annuity) is paid, because the lattice lacks its recombining structure. The proposed model overcomes this problem associating with each node of the lattice a set of representative values of the underlying asset (when valuing an option) or of the personal subaccount (when valuing a variable annuity) chosen among all the possible ones realized at that node. Extensive numerical experiments confirm the model accuracy and efficiency.


Introduction
In this paper, we propose a flexible lattice framework that allows us to obtain accurate evaluation and management of the risk affecting both long-term options written on assets paying discrete dividends, and variable annuities (VAs hereafter) embedding some provisions, like a guaranteed minimum withdrawal benefit (GMWB hereafter).In order to consider a realistic context that permits to achieve consistent risk estimates, the model has been developed in a market where the traded asset dynamics is characterized by stochastic interest rates since, as evidenced by Aas et al. (2018), the choice of an adequate interest rate model represents an additional key component in life insurance and financial evaluations.Under this perspective, the model allows to choose the most appropriate dynamics, among the ones widely used in finance and insurance, for both the interest rate and the underlying asset, by combining financial and demographic risk and embedding in the contract early exercise features.
Pricing options when the underlying asset pays discrete dividends is a nontrivial problem that, in the financial literature, has been faced applying different techniques.Attempts working on European options have been provided by Merton (1973) and Beneder and Vorst (2001).An approximation for American options with a single dividend has been proposed, among others, by Black (1975) and, separately, by Roll (1977), Geske (1979) and Whaley (1981) (the well known Roll-Geske-Whaley method).Unfortunately, in some cases, both the Black and the Roll-Geske-Whaley methods underestimate the call option price.An alternative approach widely used to price both European and American-type options written on assets paying discrete dividends is represented by lattice-based methods.They can be easily applied to take into account single or multiple discrete dividends, even if computational problems arise due to the jumps at the dividend payment dates that cause the lack of the lattice recombining structure.To overcome this problem, Vellekoop and Nieuwenhuis (2006) present a modified Cox et al. (1979) tree (CRR hereafter) based on suitable integration techniques that are introduced at the dividend payment dates.Later, in Vellekoop and Nieuwenhuis (2011), the same authors derived an integral formula to deal, among others, with discrete dividends in a free-boundary problem framework.Dai (2009) introduces a new recombining tree structure, the stair tree, which uses extra nodes only when it needs to simulate the price jumps due to the dividend payments.Indeed, the stair tree follows the CRR tree structure between dividend dates while, for nodes at the dividend dates, trinomial branching schemes are devised to connect the two adjacent CRR tree structures.Finally, it is worth citing the recent work of Costabile et al. (2018), who propose a suitable construction of a tree with a reconnecting property.
A computational problem, similar to the one encountered when evaluating financial options written on assets paying discrete dividends, may be faced in the procedures that use lattice-based methodologies to estimate the premiums of VA embedding a GMWB rider.In its classical form, such a rider returns to the policyholder the entire initial invested premium through periodical withdrawals during the policy lifetime regardless of the subaccount performance, plus the remaining subaccount value balance, if any, at maturity.As a consequence, the effect of the periodical withdrawals on the lattice structure used to discretize the subaccount value is similar to the effect of the paid dividends on the lattice discretizing the underlying asset, i.e., the lattice loses its recombining structure when the withdrawals are paid.Evaluation models for VA with GMWB have been proposed by many stems of the actuarial literature.For instance, when the policyholder receives a constant withdrawn amount at each withdrawal epoch, contributions have been developed both under the assumption of a constant interest rate by Milevsky and Posner (2001), and stochastic interest rate by Lin and Tan (2003), and Peng et al. (2012).On the other hand, when considering the more involved and challenging pricing problem that embeds a surrender option in the policy contract or, alternatively, gives the policyholder the possibility to optimally decide the amount to withdraw at each anniversary in order to maximize the current policy value, various contributions have been proposed by Milevsky and Salisbury (2006), Bauer et al. (2008), Dai et al. (2008), Chen and Forsyth (2008), Yang and Dai (2013), Luo and Shevchenko (2014), Hyndman and Wenger (2014), Bacinello et al. (2016), and Costabile et al. (2020).
The main contribution of this paper is to provide a lattice-based evaluation method suitable to evaluate both European or American options written on assets paying discrete dividends during the option lifetime, and VAs characterized by GMWB riders even when the policy embeds a surrender option.The choice of a lattice methodology has been driven by its flexibility and ease of implementation that makes it very useful for practitioners.The proposed framework is quite general because it may accommodate a wide range of diffusions used in finance and insurance, both for the stochastic interest rate and for the underlying asset, by simply modifying the process parameters.The considered continuous time processes are approximated by means of recombining lattices, which work directly on the original processes.Indeed, the generated lattice values discretize only the diffusion part of the processes, while branching probabilities are computed to assure that the first and the second order local moments of the discrete distribution match the corresponding continuous-time ones, at least within the limit.The main feature of the proposed approach is relative to the fact that it does not make use of any preliminary transformation of the original processes to obtain the approximating recombining lattices.For instance, to support this assertion, we evidence that the straightforward application of the Cox and Rubinstein (1985) model to discretize the spot rate process, which is generally heteroskedastic, leads to a non-recombining bushy tree that makes the evaluation problem computationally unmanageable.As a consequence, to construct approximating recombining lattices, preliminary transformations of the original heteroskedastic processes have been proposed to obtain homoskedastic processes, as it happens in Nelson and Ramaswamy (1990).To sum up, our model is based on the construction of reconnecting lattices for the original processes bypassing such transformations.After discretizing the spot rate process, the algorithm works to discretize the dividend-paying underlying asset process in the case we want to evaluate an option, or the subaccount value in the case we want to evaluate a VA with GMWB.The discretization is still based on a recombining binomial lattice that is similar to the one used to approximate the interest rate dynamics but, in addition, presents a vector of representative values in correspondence with each node, generated by taking into account the payments of the dividends or withdrawals up to that node.Indeed, the main problem to face is relative to the fact that the jumps due to the payments of the dividends or withdrawals cause the loss of the recombining shape of the tree.Consequently, the number of trajectories connecting the root of the tree to a generic node grows up quickly when increasing the number of time steps.Since considering all the paths connecting the root to a generic node is computationally infeasible, we propose to select a pre-specified number of representative trajectories along which we compute the underlying asset or the subaccount values.The resulting two lattices are combined in order to establish a bivariate tree presenting four branches for each node.Joint probabilities are computed for the possible jumps in order to take into account the correlation between the two processes.As the model works on representative values, a linear interpolation technique is used when solving backward through the lattice to compute the option price or the VA present value in terms of the discounted payoffs over the bivariate lattice branches.
The remaining of the paper is organized as follows.In Sect.2, after presenting the continuous-time model, we establish the discretizing lattice for the interest rate and the underlying asset (when valuing an option) or the personal subaccount (when valuing a variable annuity), then introducing the resulting bivariate lattice.Section 3 concerns the application of the algorithm presented in Sect. 2 for the pricing of options written on assets paying discrete dividends and the valuation of VAs with GMWB riders.Section 4 presents some numerical results confirming the accuracy of the proposed model and provides some comparisons with the existing methods.Finally, Sect. 5 concludes the paper.

The continuous-time framework
To present contemporaneously the framework used to evaluate options written on dividend paying stocks and VAs embedding GMWB riders, we consider a market where, under the risk-neutral probability measure, interest rates fluctuate stochastically according to the following differential equation and affect the values assumed by an underlying asset that pays at times t h , h = 1, . . ., m fixed amounts, D t h , being t 0 = 0 the contract inception and T ≥ t m the contract maturity measured on annual basis.Hence, the asset value jumps down by an amount equal to the paid amount and its dynamics may be described by the following stochastic differential equation where t h−1 < t ≤ t h and S(t h ) = S(t − h ) − D t h , with S(t − 0 ) = S(t 0 ); furthermore, σ S is a positive constant and α is about to be specified below, while W S (t) and W r (t) are standard Brownian motions with constant correlation ρ.Depending on the definition of the drift term μ(r (t)) and the diffusion term σ r (r (t)), the stochastic differential equation (1) may describe several processes like: • a Cox et al. (1985) process when μ(r (t)) = γ (θ − r (t)) and σ r (r (t)) = σ r √ r (t), with the reversion speed γ , the long-term revision target θ , and the risk-free rate volatility σ r being positive constants.The condition 2γ θ > σ2 r ensures that r (t) is never negative; • a Vasicek (1977) process when μ(r (t)) = γ (θ − r (t)) and σ r (r (t)) = σ r ; • a Hull and White (1994) process when μ(r (t)) = γ (β(t) − r (t)), with β(t) being a deterministic function of time, and σ r (r (t)) = σ r .
If we consider the evaluation problem of an option, equation (2) will represent the dynamics of the underlying asset paying a discrete dividend D t h at time t h and the quantity α = 0.1 Indeed, the latter identifies the insurance fee for the GMWB whenever we want to evaluate a VA (see, among others, Milevsky and Salisbury 2006;Dai et al. 2015;Costabile et al. 2020 for further details).In this case, equation (2) will describe the dynamics of the personal subaccount of a single premium VA embedding a GMWB with withdrawals D t h allowed, as usual in actuarial practice, at each contract anniversary t h = h, h = 1, . . ., T .It is worth evidencing that a VA contract is characterized by the presence of two accounts for the policyholder: a personal subaccount invested in a well diversified reference fund, for which we have supposed the dynamics showed in equation ( 2); a guarantee account having value A t at time t, with A t − h denoting its value at time t h before the payment of the withdrawal and A t h being the value at the same time after the withdrawal has been paid.The initial investment S(0) forms the initial balance of both the accounts, while S(t) remains trapped in zero once it has reached zero.To simplify matters, we detail how both the accounts work by considering a static approach according to which the policyholder withdraws a fixed amount D = S(0) T , i.e., D t h = D with t h = h and h = 1, . . ., T , at each policy anniversary.At the first epoch t 1 = 1, an instant before the withdrawal payment, the subaccount has value S(t − 1 ) generated by (2) starting from S(0), while the guarantee account has value A t − 1 = S(0).After the withdrawal, both the accounts decrease by D assuming value S(t 1 ) = max(S(t − 1 ) − D, 0), 2 and A t 1 = A t − 1 − D, and so on in the next payment anniversaries.In general, the guarantee account in such an approach decreases by D at each policy anniversary, that is A t h = S(0) − h D, h = 1, . . ., T , and when the VA contract expires it assumes value A T = A T − − D = 0. Finally, it worth mentioning that the mixed approach is similar to the static one but it adds the possibility for the policyholder to surrender the contract early. 3Fig. 1 The binomial lattice discretizing the r -process

The risk-free rate process discretization
The first step of the proposed model is to discretize the stochastic interest rate process (1) on the entire interval [0, T ] by a recombining binomial lattice.Following the Costabile and Massabo' (2010) algorithm that allows to force a heteroskedastic process to be homoskedastic, we obtain a recombinant and computationally simple lattice in which the number of nodes grows up linearly with the number of time steps.
As usual in lattice methods, we start by dividing the interval [0, T ] into n subintervals of equal length t = T /n, and choose n as a multiple of T in order to assure that each contract anniversary coincides with a layer of lattice nodes.In this way, we avoid the biases that may occur whenever, for instance, the VA withdrawal epochs do not coincide with discrete observation times in the lattice approximation.For i = 0, . . ., n, we denote the lowest node at time i t by (i, 0), the second lowest one by (i, 1), and so on.Hence, we denote by r (i, k), k = 0, . . ., i, the value of the interest rate discrete approximating process at node (i, k).Rooting the tree at r (0, 0) = r (0), at the generic i-th time interval, with i = 1, . . ., n, we have the values: , the value at node (i, k) is set equal to the value associated with node (i − 2, k − 1).
In other words, as depicted in Fig. 1, once the interest rate values on the lattice edges have been computed, the values for the inner nodes are simply determined by generating horizontal layers of nodes starting from the ones located on the two edges.At this point, we are left to compute transition probabilities.Starting from a generic value r (i, k), upward, p r (i, k), and downward, q r (i, k) = 1 − p r (i, k), transition probabilities are computed by imposing that the two successor points for the state variable, r (i + 1, k + 1) and r (i + 1, k), bracket the expected value of the interest rate r(0, 0) 2, 1)   r(3, 2)   r(3, 1)   r(4, 2) r(4, 3) r(4, 1) •r(3, 0) + µ(r(3, 0))Δt Fig. 2 An example of multiple jumps in the r -process in the next time interval, that is Whenever r (i, k) + μ(r (i, k)) t does not lie between r (i + 1, k + 1) and r (i + 1, k), p r (i, k) is not a legitimate probability, i.e., a number between 0 and 1, so that multiple upward or downward jumps are required.An example, starting from node (3, 0), is given in Fig. 2.
In general, a multiple jump for the process located at node (i, k) is achieved by defining k m as otherwise, and q r (i, k) = 1 − p r (i, k), respectively.

The underlying asset process discretization
The second step of the algorithm is to establish a recombining lattice discretizing the underlying asset dynamics (2).We start from its initial value S(0, 0) = S( 0) and compute at each node (i, j), i = 0, . . ., n, j = 0, . . ., i, the quantity S(i, j) = S(0, 0)u j d i− j , with u = e σ S √ t and d = 1 u .To define the transition probability assigned to each lattice node (i, j), again we impose that the successor points for S(i, j), S(i + 1, j + 1) and S(i + 1, j), bracket the expected value of the process in the next time interval.Observing the drift term in (2), we have to take into account all the possible values assumed by the risk-free rate at nodes (i, k) lain at each i-th time step.This is done by defining a triplet (i, j, k) in correspondence of which the risk-free rate assumes value at node (i, k) and the asset assumes value at node (i, j).In this way, upward transition probability p S (i, j, k) associated with node (i, j, k) is computed by The corresponding probability of a downward jump is q S (i, so that the successors of S(i, j) are S(i + 1, j m + 1) and S(i + 1, j m ) with transition probabilities and q S (i, j, k) = 1 − p S (i, j, k).

The algorithm to manage dividends or withdrawals
At this point we have to take into account the effect on the lattice of the periodical dividends if we want to evaluate an option, or withdrawals if we want to evaluate a VA with a GMWB.Indeed, the lattice loses its recombining structure as shown in Fig. 3 where a six time-step asset evolution is depicted and it is supposed, for simplicity, that the time horizon between two consecutive payments is constant.
The tree originates at time t = 0 and ends at time t = 6 t, while at time t 1 = 2 t the asset pays the first dividend/withdrawal D t 1 and at time t 2 = 4 t the asset pays the second dividend/withdrawal D t 2 .Figure 3 shows clearly how the presence of the Fig. 3 The effect on the lattice structure when two discrete dividends/withdrawals equal to D t 1 and D t 2 are paid at time t 1 = 2 t and t 2 = 4 t, respectively dividend/withdrawal payments causes a huge increment in the number of possible asset values since the number of nodes grows deeply when the number of time steps increases.It is easy to understand that, proceeding forward along the tree in this way, the problem becomes computationally difficult to manage and this is a key point to look at when developing a reliable and efficient lattice-based evaluation method.
In order to overcome this obstacle, we propose an algorithm that saves the lattice recombinant shape even in the case of dividend or withdrawal payments.The idea is to build a tree for which, at each node (i, j), we do not consider a single asset value but a set of possible values capturing the asset evolution along each trajectory connecting the root of the tree with the considered node, taking properly into account all the dividends/withdrawals paid up to time i t.
To present the procedure, as usual in financial practice when managing options written on dividend paying stocks, we suppose that the first dividend may be paid at a generic date t 1 after the contract inception, while the subsequent dividends at time t h , h = 2, . . ., m, are yearly paid starting from t 1 .To capture this aspect, we denote by = n T the number of tree steps falling in an entire contract year, and with λ 1 = t 1 the number of time steps falling between the contact inception and the first payment date t 1 .The original choice of n as a multiple of T assures that is an integer but, to avoid biases, we have to choose n in a away that also λ 1 is an integer.Satisfying both such condition, we assure that all the dividend paying dates coincide with one of the lattice layers of nodes.Indeed, the first dividend payment date coincides with the λ 1 -th step, the generic h-th dividend payment date t h , h = 2, . . ., m, coincides with step λ 1 +(h −1) (because the subsequent dividends at time t h , h = 2, . . ., m, are yearly paid starting from t 1 ) and, finally, the contract maturity T coincides with the n-th step of the lattice.When evaluating a VA with a GMWB, the notation simplifies a bit because, generally, withdrawals are paid in correspondence with the policy anniversaries so that t h = h, h = 1, . . ., T , and obviously also λ 1 = .To capture the effects of the dividends/withdrawals on the underlying asset value, their amounts are represented in terms of the number of asset shares transferred to the beneficiary multiplied by the asset value at the payment date.To detail how the dividends/withdrawals are managed, let us consider a generic path starting from inception and reaching node (i, j) characterized by the asset values S(l, j l ) and denoted by τ (i, j) = {(l, j l ), l = 0, . . ., i; j i = j}.Following the considered trajectory, at time t 1 , i.e., after λ 1 time steps, with the first dividend/withdrawal payment the asset value decreases by D t 1 = n S (λ 1 , j λ 1 )S(λ 1 , j λ 1 ), where n S (λ 1 , j λ 1 ) represents the number of shares deducted from the underlying asset with the aim of paying the dividend/withdrawal D t 1 .Similarly, at time t 2 , i.e., after λ 1 + time steps, we may define the second dividend/withdrawal payment In general, at time t h , h = 1, . . ., m, with t h ≤ i t, i.e., after λ 1 + (h − 1) time steps, with the h-th dividend/withdrawal payment the asset value decrements by ).Hence, if the asset value follows the trajectory τ (i, j) to reach node (i, j), the total number of shares deducted to pay the dividends/withdrawals up to the i-th step is given by where H τ (i, j) = h : (λ 1 + (h − 1) , j λ 1 +(h−1) ) ∈ τ (i, j), h = 1, . . ., m , so that the underlying asset value arising following trajectory τ (i, j) is computed by It is worth noting that the computed number of shares in equation ( 5) represents a map that does not change its value between two consecutive payment dates, but its value is strictly dependent upon τ (i, j).It means that, generally, different paths reaching node (i, j) lead to different values for the total number of shares deducted to pay the dividends/withdrawals but, clearly, considering all the i j possible paths reaching node (i, j) is computational infeasible.To reduce the computational complexity of the evaluation problem, we propose to select among all the possible paths reaching node (i, j) a set made up of η(i, j) representative ones used to compute the resulting asset values, where η(i, j) is given by the explicit formula reported in the following proposition.
Proposition 1 In a binomial lattice characterized by n time steps, the number of representative underlying asset values, η(i, j), associated with a generic node (i, j), i = 0, . . ., n, j = 0, . . ., i, is given by: Proof See the Appendix.
To give an idea of the computational advantages that arise when using the proposition, consider for instance node (25, 7) in a lattice used to discretize an asset paying the first dividend/withdrawal at time t 1 = 2 t so that λ 1 = 2, the second dividend/withdrawal at time t 2 = 6 t so that = 4, and the general h-th dividend/withdrawals at time t h = (λ 1 + (h − 1) ) t.We should detect 25 7 = 480700 different trajectories, but we compute only η(25, 7) = 33 representative values for the underlying asset.To detail how we compute the different η(i, j) representative values for node (i, j), we present the following iterative procedure: Step a) start from the lowest trajectory τ min (i, j) reaching node (i, j) made up by i − j down steps followed by j up steps.An example of such trajectory is depicted in Fig. 4 for (i, j) = (4, 1).The first representative asset value, RS(i, j; 1), is obtained by applying formula (6) upon trajectory τ min (i, j), that is S(i, j); (7) Step b) the remaining subset elements RS(i, j; l), l = 2, . . ., η(i, j), are generated iteratively starting from τ min (i, j).Consider all the quantities ), and choose the one that corresponds to the maximum value assumed by h, namely h max , at node If such a h max does not exist, τ min (i, j) is the only representative trajectory and there is only one representative asset value (computed by ( 7)) associated with node (i, j); otherwise, the second representative asset value RS(i, j; 2) is computed as in (7) by substituting τ min (i, j) with τ new (i, j).Repeating this procedure, starting from τ new (i, j) until such a h max does not exist, allows us to obtain exactly η(i, j) representative asset values associated with node (i, j).

The resulting bivariate lattice
Now, we are in the position to establish a bivariate approximating lattice that takes into account the joint evolution of the underlying asset and the stochastic interest rate.At each time step, we consider all the values RS(i, j; l) and r (i, k), where i = 0, . . ., n, j = 0, . . ., i, k = 0, . . ., i, and l = 1, . . ., η(i, j).The index i represents the time step, j and k the positions at the i-th time step of the asset value and the risk-free rate, respectively, and l is the position index in the subset of the representative values associated with node (i, j) of the asset lattice.Consequently, denoting by (i, j, k; l) each possible state of nature, the joint process presents four possible scenarios because each discrete process can show an upward movement, denoted by "u," or a downward movement, denoted by "d."Each scenario is labeled with an ordered pair where the first element indicates the S-movement and the second element refers to the r -movement: the first scenario is associated to the case of an upward movement both in S and r discrete processes with probability p uu ; the second one is associated with the case of an upward movement in S and a downward movement in r with probability p ud ; the third one is associated with a downward movement in S and an upward movement in r with probability p du ; the latter one is associated with a downward movement both in S and r with probability p dd .The probability of reaching each scenario is computed by adjusting the product of the marginal probabilities associated with the corresponding movements in the lattice approximating r (t) and S(t) to take into account their correlation, and are computed by solving the linear system induced by the following equations: 1. the probabilities sum to 1, 2. the constraint upon the marginal probability of the S-process is set up as 3. the constraint upon the marginal probability of the r -process is set up as 4. the covariance between the discretized r -process and S-process must equal the covariance between the continuous time ones, 4 that is Solving simultaneously equations ( 8)-( 11) leads to the following solutions: 4 Equation ( 11) can be obtained as follows.Starting form a generic state (i, j, k), the discrete S-process deviates from its mean by ±σ S S(i, j) √ t ("+" in case of an upward movement and "-" in case of a downward movement), while the r -process deviates from its mean by ±σ r (r (i, k)) √ t.The covariance of the joint discrete bivariate distribution is computed by the difference between the expected value of the product of the two variables, that is and the product between the expected value of each variable, that is Hence, the obtained quantity must replicate the covariance of the corresponding continuous time distributions, that is ρσ S S(i, j)

Method applications
In this section, we detail how the proposed method may be easily applied to evaluate both European and American-type options written on dividend paying stocks and VAs embedding GMWB riders and surrender options.Starting from the constructed lattice structure, we report the procedures needed to proceed backward in order to obtain at inception the option prices, the single premiums of VA policies or, alternatively, the value for the provision α continuously paid for the GMWB riders.

Options on dividend paying stocks
The first application of the presented discretization regards to European call option contracts written on underlying assets paying discrete dividends.We denote by V (i, j, k; l) the option value in state of nature (i, j, k; l), with i = 0, . . ., n, j = 0, . . ., i, k = 0, . . ., i, and l = 1, . . ., η(i, j).To price such a derivative at the contract inception, we proceed backward starting from the n-th step of the bivariate tree coinciding with the option maturity T .Here, the option payoff is computed by V (n, j, k; l) = max (S(n, j; l) − K , 0), with K being the option strike price.
Proceeding backward, the generic continuation value at node (i, j, k; l) is given by where arising when the r -process shows an upward movement from r (i, k) to r (i + 1, k m + 1) and the asset value moves upward from S(i, j) to S(i + 1, j m + 1); arising when the r -process shows a downward movement from r (i, k) to r (i + 1, k m ) and the asset value moves upward from S(i, j) to S(i + 1, j m + 1); • V (i + 1, j m , k m + 1; l du ) is the option value in correspondence of the asset value [RS(i, j; l) − D t h I {t h ,h=1,...,m} ] S(i+1, j m )

S(i, j)
arising when the r -process shows an upward movement from r (i, k) to r (i + 1, k m + 1) and the asset value moves downward from S(i, j) to S(i + 1, j m ); • V (i + 1, j m , k m ; l dd ) is the option value in correspondence of the account value [RS(i, j; l) − D t h I {i t=t h ,h=1,...,m} ] S(i+1, j m ) S(i, j) arising when the r -process shows a downward movement from r (i, k) to r (i + 1, k m ) and the account value moves downward from S(i, j) to S(i + 1, j m ).

S(i, j)
belongs to the set of the representative asset values associated with node (i + 1, j m + 1), V (i + 1, j m + 1, k m +1; l uu ) is already available; otherwise, V (i +1, j m +1, k m +1; l uu ) is computed by linear interpolation as in De Angelis et al. (2016).The same happens for all the other quantities involved in equation ( 12).Whenever we want to evaluate an American option, as usual in financial practice, we have to consider in the backward formula (12) the early exercise value of the option contract, that is Clearly, the pricing procedure may be easily adapted to price European or American put options written on assets paying discrete dividends.

VAs embedding GMWB riders
The second application of the proposed methodology regards to the evaluation of VA policies embedding GMWB riders, which allow withdrawals D t h at each contract anniversary t h = h, h = 1, . . ., T , as usual in the actuarial practice.As a consequence, we have λ 1 = in the binomial discretization.The personal subaccount value dynamics follows equation (2) and has initial value S(0), exactly as the initial guarantee account, i.e., A 0 = S(0).In our framework, we can easily manage both the static and the mixed approach in which the policyholder withdraws a fixed amount D t h = D = S(0) T , at each policy anniversary.The policy evaluation procedure starts again from the policy maturity T where, in correspondence of each state of nature (n, j, k; l), the policy value is computed as where we recall that A T − = D. Proceeding backward, the generic contract value at node (i, j, k; l) if early redemption is not allowed, i.e., in the static approach, is given by where the quantities V (i + 1, j m + 1, k m + 1; l uu ), V (i + 1, j m + 1, k m ; l ud ), V (i + 1, j m , k m + 1; l du ), and V (i + 1, j m , k m ; l dd ) are computed following the procedure detailed in Sect.3.1.Whenever the VA contract includes a surrender option allowing the policyholder to redeem the contract early an instant before each withdrawal date, we fall in the mixed approach where we define the surrender value coherently with the actuarial literature,5 as where φ is the contractually fixed penalty paid by the policyholder for the policy early redemption.Greater is φ, more discouraged is the policyholder to surrender the contract.As a consequence, when considering the mixed approach, equation (13) modifies as At this point, we are left to detail the procedure for computing the insurance fee α paid for the GMWB rider that appears in the subaccount stochastic dynamics (2).Such a fee is generally assumed to be continuously payable as a constant fraction of the subaccount value.The backward procedure described above represents the tool used to evaluate α because it is included in all the computations reported in equations ( 13) and (15).Indeed, the quantity α appears in the formula for computing the probabilities of the S-movements, as reported in equations ( 3) and ( 4), so that it is embedded in equations ( 13) and ( 15) through the probabilities p uu , p ud , p du , and p dd .As a consequence, the policy value at the contract inception results as a function of the provision α, i.e., V (0, 0, 0; 1) = V (α), so that to compute the fair value α paid at the contract inception for the GMWB rider we have to solve numerically the nonlinear equation V (α) − S(0) = 0.The evaluation procedure detailed above does not take into account the effect of mortality risk that, on the contrary, represents a crucial aspect to look at when valuing insurance policies.However, the method may be easily extended to include the mortality effect as detailed hereafter.Suppose to consider an insured aged z in calendar time t, i.e., at the contract inception, and denote by p t z = e −μ z (t) the one-year survival probability, where μ z (t) is the real world mortality rate assumed to be constant in year t, so that the w-year survival probability is given by w p t z = w−1 l=0 e −μ z+l (t+l) .The following Wang (2000) method is used to convert the real world mortality rates into risk-neutral ones: • the real world probability for an individual of age z who dies before age z + w is given by w q t z = 1 − w−1 l=0 e −μ z+l (t+l) ; • the corresponding risk-neutral world distribution function is w q t z = ( −1 ( w q t z )− π), where is the distribution function of the standard normal distribution, and π is the market price of risk; • the risk-neutral mortality rate for an individual of age z + w in calendar year t + w is computed recursively by μ z+w (t + w) = − ln(1 − w+1 q t z ) − w−1 l=0 μ z+l (t +l).
The obtained risk-neutral mortality rates are included in the backward induction scheme as follows.The policy value in equation ( 13) represents the amount received by the policyholder if he remains alive during the time interval Since the mortality rate is constant in each calendar year, we have to consider that at time (i − 1) t the age of the individual is z + (i − 1) t , and the calendar year is t + (i − 1) t .As a consequence, the survival probability is given by e − tμ z+ (i−1) t (t+ (i−1) t ) .On the contrary, if the policyholder dies during the same interval, with probability 1 − e − tμ z+ (i−1) t (t+ (i−1) t ) , the subaccount value is immediately returned to the beneficiary and the GMWB is cancelled.For illustrative purposes, we assume that the real world mortality rate follows the GAR-94 table where the effect of mortality improvement is captured as follows: the mortality rate for an individual aged z in year 1994 + w can be estimated as q 1994+w z = q 1994 z (1 − A A z ) w , where A A z is the annual improvement factor in the mortality rate for age z.Both q 1994 z and A A z are provided by Society of Actuaries Group (1995).To sum up, the contract value is computed backward by summing up two components: • the quantity paid in case of death in state of nature (i, j, k; l) multiplied by the death probability in the interval [(i − 1) t, i t), defined as where the subaccount values RS(i + 1, j m + 1; l uu ), RS(i + 1, j m + 1; l ud ), RS(i + 1, j m ; l du ), and RS(i + 1, j m ; l dd ) are already available in the proposed discretization; • the quantity characterizing the policy value in the same state of nature whenever the insured remains alive during the considered period multiplied by the corresponding survival probability, defined as where the quantities V (i +1, j m , k m +1; l du ), and V (i +1, j m , k m ; l dd ) are again computed following the procedure detailed in Sect.3.1.
The entire contract value in each state of nature (i, j, k; l) is finally given by in the static approach, or by in the mixed approach.Finally, we remark that the procedure to compute the provision α when mortality risk is included in the VA policy contract is the same as illustrated above.

Numerical results
In this section, we present some numerical investigations to assess the accuracy of the proposed algorithm.We analyze both the cases of European and American-style options written on assets paying discrete dividends, and VA embedding GMWB riders.
In addition, we provide comparison with the existing evaluation models in financial and actuarial literature or with benchmark values created "ad hoc" whenever the existing literature does not provide useful results.

Options on dividend paying stocks
The first numerical experiments are relative to options written on assets paying discrete dividends.We start by assessing the model performance providing a comparison with the existing pricing models.Unfortunately, such models do not allow to embed a stochastic dynamics for the interest rate so that we simplify matters by considering a constant interest rate in place of the described stochastic dynamics.Then, we introduce a stochastic dynamics for the interest rate and provide a comparison between the proposed algorithm and a Monte Carlo method chosen as the benchmark.
In detail, in Table 1, we start by considering the simple case of an American call option with time to maturity T = 1 year written on an underlying asset paying a single dividend of amount D t 1 = 7.0 at time t 1 , with t 1 being equal to 0.1, 0.5, and 0.9, respectively.The other parameter values are: the initial asset value S(0) = 100; the asset volatility σ S = 0.30; the strike price K may assume value 70, 100, and 130, respectively.In addition, in order to provide comparisons of our results (denoted by the acronym DDMR followed by the number of the considered time steps, n = 500 and n = 1000) with the ones reported in Vellekoop and Nieuwenhuis (2006) (VN) and in Dai (2009) (Dai), we degenerate the interest rate process in (1) in a constant value r = 0.05.It is worth observing that, in all the analyzed cases the considered models provide very close results that allow us to assess the accuracy of the proposed model.
To further show the model accuracy, in Table 2, we consider the case of a European call option written on a multiple dividend paying stock.The option has maturity T = 7 years, while the underlying asset pays seven different discrete dividends equal to D t 1 = 6.0,D t 2 = 6.5, D t 3 = 7.0, D t 4 = 7.5, D t 5 = 8.0, D t 6 = 8.0 and D t 7 = 8.0.The epoch of the first dividend payment t 1 is reported in the first column of Table 2, while the remaining dividend epochs are given by t h+1 = t 1 + h, h = 1, . . ., 6.
The other parameters are fixed to S(0) = 100 and σ S = 0.25, while we still consider a constant value for the interest rate r = 0.06 in order to provide comparison with the existing models.In detail, the table exhibits the comparison of the option price obtained by means of the proposed algorithm (DDMR) when n = 210, with the values reported in Vellekoop and Nieuwenhuis (2006) (VN) where n assumes value 250, 500, and 1000, in Haug et al. (2003) (HHL), in Bos et al. (2003) (BGS), in Bos and Vandermark (2002) (BvdM), and in Beneder and Vorst (2001) (BV).The benchmark value is assumed to be the price computed by Vellekoop and Nieuwenhuis (2006) (VN) when applying a Richard extrapolation technique (VNRE).As it is shown in Table 2, our results are the closest to the benchmark with respect to the other considered models, despite we compute the option prices with only n = 210 time steps.For instance, the VN model requires a greater and greater number of time steps to be consistent with the benchmark.In order to apply the evaluation model when the interest rate is stochastic, in Table 3, we consider the case of an American call option when the underlying asset pays a single or multiple discrete dividends and the interest rate process is described by a Cox et al. (1985) dynamics.We impose S(0) = K = 100, γ = 0.5, θ = r 0 , while the other parameter values are reported in the table, being m the number of dividend payments.We compute option prices by the proposed methodology (DDMR) with n = 500 and compute also least-squares Monte Carlo (MC) estimates using 10 million trials with the corresponding 95% confidence interval (CI).As it may be observed, our results lay in the corresponding Monte Carlo confidence interval thus evidencing the model accuracy also when a stochastic dynamics for the interest rate is considered.

VA with a GMWB rider
The second algorithm application is relative to VAs with GMWB riders.As usual in the actuarial practice, we suppose that the withdrawals are paid at each policy anniversary measured on yearly basis, i.e., denoting by T the policy maturity the withdrawal epochs are t h = h, h = 1, . . ., T .We present numerical experiments both with or without considering mortality risk and with or without embedding a surrender option.In absence of such an option, we fall in the static approach where the policyholder withdraws a fixed amount D t h = D = S(0)/T at each anniversary while, whenever the surrender option is embedded in the contract, we fall in the mixed approach and suppose that the early redemption is allowed only an instant before each withdrawal epoch.Furthermore, when it occurs, a penalty φ = 0.1 is charged to the part of the withdrawn amount exceeding D as specified in equation ( 14).
At first, to assess the goodness of the proposed model, we provide comparisons with the existing evaluation methods so that we degenerate the stochastic interest rate process being constant at fixed levels.We start by presenting Table 4, where we consider VA policies with GMWBs in absence of surrender options and without considering mortality risk, when the constant interest rate assumes value r = 0.0325, the initial investment in the subaccount is S(0) = 100, while the other parameter values are reported in the table.We exhibit the policy single premiums evaluated by means of the proposed algorithm (DDMR) in comparison with the ones computed by Yang and Dai (2013) (YD) and the Monte Carlo estimates (MC) with relative standard errors (SE) reported in Yang and Dai (2013).The DDMR single premiums, coherently with the ones obtained by YD, increase when σ S increases and converge smoothly to the MC benchmark values.
In Table 5, we evaluate the insurance fee α solving, by means of the Brent's method,6 the following nonlinear equation where we recall that V (α) = V (0, 0, 0; 1), i.e., the policy value computed at the contract inception following the proposed method.For different values of σ S and time to maturity T equal to 20 and 25 years, we evaluate the provision α for the GMWB rider by considering different cases characterized by the presence or absence of the surrender option and mortality risk.In this last case, an individual of age z = 65 is considered and the Wang transform is applied with the market price of risk set to π = 0.565, as it has been already estimated in Costabile et al. (2021).The other parameters are set to S(0) = 100 and, in order to provide comparison with the existing models, we again degenerate the interest rate being it constant at level r = 0.0325.The values computed with the proposed algorithm (DDMR) have been obtained by considering n = 300 while the comparison is done with the values reported in Costabile et al. (2020) (CMR) and Yang and Dai (2013) (YD) when their algorithm is based on n = 1000 time steps.All the results are coherent with the ones obtained with CMR and YD algorithms despite they have been computed with the DDMR algorithm based on only n = 300 steps with respect to n = 1000 steps used in CMR and YD.In Table 6, we introduce a stochastic dynamics for the interest rate and report a comparison of the VA single premiums evaluated by the proposed bivariate tree with the ones obtained by means of the three-dimensional algorithm provided by Dai et al. (2015) (DYL).The risk-free rate dynamics is described by the Hull and White (1994) model.To provide benchmark values, we develop a crude Monte Carlo (MC) method with 10 million samples and the table reports both the MC estimates and the corresponding 95% confidence interval (CI).No surrender option and mortality risk is here considered and the parameters are set as: S 0 = 100, r (0) = 0.0325, σ r = 0.01, γ = 0.1, β(t) = 0.0325, and ρ = −0.25.The other parameter values are reported in the table, which evidences how the results generated by the proposed method are very close to the benchmark values and lay in the corresponding confidence interval, while this circumstance is not satisfied by the results provided by Dai et al. (2015).
To introduce the effect of the mortality risk on the VA valuations when the interest rate is stochastic, we conduct some experiments varying the insured age z and report the computed VA single premium in Table 7 and 8.We consider VA policies embedding GMWB riders with or without considering the surrender option, while the stochastic dynamics for the interest rate is described by means of the Cox et al. (1985) (CIR) model.In detail, in Table 7 we considers a policyholder with age z = 65 at the contract inception, while in Table 8 we impose z = 40 and compute the policy single premiums for different values of the number time steps n, maturity T , and personal subaccount volatility σ S .The other parameters are set as follow: r 0 = 0.0325, ρ = 0.25, γ = 0.5, θ = r 0 , σ r = 0.10, and π = 0.565.For each case, we report the policy single premium without considering the surrender option and, in round brackets, the corresponding premium when embedding the surrender option in order to show the increment in the policy premium impressed by the presence of the surrender option.
To validate the reported results, we have also developed a Monte Carlo approach in absence of the surrender option assumed to be the benchmark.The VA premiums are obtained with the Monte Carlo (MC) approach by considering 1 million trials, and we report the corresponding 95% confidence interval (CI).In all the considered cases, the experiments show that, fixing the age, when σ S increases the VA premiums increase, and when the maturity decreases the VA premiums increase, as expected.Furthermore, when the insured age decreases the VA premiums increase.It is also worth noting that the results (without considering the surrender option) lay in the corresponding Monte Carlo confidence interval in all the cases, thus assessing the model accuracy.
Table 9 shows the provision values α for the GMWB riders, reported in basis points (b.p.), for all the test cases reported in Tables 7 and 8.They are computed as specified in Sect.4.1 by fixing the number of time steps at n = 200.Coherently with the results reported in Tables 7 and 8, Table 9 exhibits that the inclusion of a surrender option in the contract causes the increment of α.The same occurs, for instance, when the insured age at the contract inception or the policy maturity decreases.
To conclude the numerical experiments, in Fig. 5 we show the behavior of the provision α for different values of σ r and σ S when it is computed by the proposed algorithm with n = 100 and considering a CIR dynamics for the risk-free interest rate.More in detail, 0 ≤ σ r ≤ 0.5 and 0.2 ≤ σ S ≤ 0.5.The remaining parameters are set to: S(0) = 100, r (0) = 0.0325, γ = 0.5, θ = 0.05, ρ = 0.25, T = 20, D = 5, and φ = 0.1.The figure confirms, as expected, that when σ r or σ S increases, α moves in the same direction.

Conclusions
In a framework characterized by a stochastic interest rate, we have proposed a model useful for computing both the price of European and American-style options written on assets paying discrete dividends, and the fair single premiums of VAs embedding GMWB riders and surrender options.The presence of multiple discrete dividends, in the option case, or withdrawals, in the VA case, makes the lattice not recombining and, In what follows, we give the details of "Step 2," "Step 3," and ""Step 4," as "Step 1" is easy to understand.
Step 2 To count the nodes in AB L M, it is divided into three parts: a triangle AE F, and two polygon of vertices E BG H F and G L M H.In the triangle AE F, in correspondence of each payment epoch, we count the nodes enclosed between AE and F A, including the nodes lain on F A and excluding those ones lain on AE.To this end, we have to know the number of payment dates falling before the min( j, i − j)-th epoch starting Again, the nodes on the left border H G of the polygon G L M H are excluded because already counted in the polygon E BG H F.
As a consequence, the total number of nodes considered in the quadrilateral AB L M is obtained by summing up the quantities appearing in formula (A.1), (A.2), and (A.3).Whenever min( j, i − j) < λ 1 , the number of nodes considered in AB L M is given by min( j, i − j) i − λ 1 − min( j, i − j) − 1 .
Step 3 At this point, we need to compute the number of nodes belonging to the quadrilateral DC L M. We start from node D coinciding with the max( j, i − j)-th time step starting from inception.The first payment epoch immediately after this node is the max( j,i− j)−λ 1 + 2 -th.To count the number of nodes belonging to DC L M at this epoch, we have to consider two cases: • the first epoch after node D falls before the highest vertices of the quadrilateral, M. In this case the number of nodes considered at that date coincides with the number of time steps needed to reach the first payment epoch starting from node D, that is max( j, i − j) − λ 1 + 1 + λ 1 − max( j, i − j); • such an epoch falls after the highest vertex of the quadrilateral, M. In this case the number of nodes is given by i − λ 1 + λ 1 − i, (A.4) which represents also the number of steps to reach the next payment date starting from the i-th step of the lattice (in Fig. 6, it is the length of the sides C L and M D).
To sum up, the number of nodes belonging to the quadrilateral DC L M at the first possible payment date is computed as min max( j, i − j) where the second factor contemplates the possibility that the max( j,i− j)−λ 1 + 2th epoch coincides with the i−λ 1 -th one, implying that the quadrilateral DC L M has 123 no node to be considered.Now, starting from the max( j,i− j)−λ 1 + 2 -th epoch, the number of payment dates remaining up to the i−λ 1 -th epoch is max i − λ 1 − max( j, i − j) − λ 1 − 2, 0 , but, to take into account the possibility that the max( j,i− j)−λ 1 + 2 -th epoch coincides with the i−λ 1 -th epoch and no node belongs to the quadrilateral DC L M, we have to consider Consequently, the total number of nodes belonging to DC L M is obtained by summing up the quantity in (A.5) to the product between (A.4) and (A.6), that is In the case illustrated in Fig. 6, there are two nodes considered in DC L M and they are indicated by white circles.

LFig. 6
Fig.6The number of representative asset values for node(i, j)

Table 1
American call option on a single dividend paying stock

Table 2
European call option on a multiple dividend paying stock

Table 3
American call prices when the interest rate is stochastic

Table 5
Values of the provision α for the GMWB rider

Table 6
VAs embedding GMWBs under the Hull-White model with T = 25 and T = 20

Table 7
VAs embedding GMWB riders when insured age at the inception is z = 65

Table 8
VAs embedding GMWB riders when insured age at the inception is z = 40