A unified pricing of variable annuity guarantees under the optimal stochastic control framework

In this paper, we review pricing of variable annuity living and death guarantees offered to retail investors in many countries. Investors purchase these products to take advantage of market growth and protect savings. We present pricing of these products via an optimal stochastic control framework, and review the existing numerical methods. For numerical valuation of these contracts, we develop a direct integration method based on Gauss-Hermite quadrature with a one-dimensional cubic spline for calculation of the expected contract value, and a bi-cubic spline interpolation for applying the jump conditions across the contract cashflow event times. This method is very efficient when compared to the partial differential equation methods if the transition density (or its moments) of the risky asset underlying the contract is known in closed form between the event times. We also present accurate numerical results for pricing of a Guaranteed Minimum Accumulation Benefit (GMAB) guarantee available on the market that can serve as a benchmark for practitioners and researchers developing pricing of variable annuity guarantees.


Introduction
Many wealth management and insurance companies worldwide are offering investment products known as variable annuities (VA) with some guarantees of living and death benefits to assist investors with managing pre-retirement and post-retirement plans. These products take advantage of market growth while provide protection of the savings against market downturns. Insurers started to offer these products from the 1990s in United States. Later, these products became popular in Europe, UK and Japan and more recently in Australia. The VA contract cashflows received by the policyholder are linked to the investment portfolio choice and performance (e.g. the choice of mutual fund and its strategy) while traditional annuities provide a pre-defined income stream in exchange for the lump sum payment. According to LIMRA (Life Insurance and Market Research Association) reports, the VA market is huge: VA sales in United States were $158 billion in 2011, $147 billion in 2012 and $145 billion in 2013.
The types of VA guarantees (referred in the literature as VA riders) offered for investment portfolios are classified as guaranteed minimum withdrawal benefit (GMWB), guaranteed minimum accumulation benefit (GMAB), guaranteed minimum income benefit (GMIB) and guaranteed minimum death benefit (GMDB). These guarantees, generically denoted as GMxB, provide different types of protection against market downturns and policyholder death. GMWB allows withdrawing funds from the VA account up to some pre-defined limit regardless of investment performance during the contract; GMAB and GMIB both provide a guaranteed investment account balance at the contract maturity that can be taken as a lump sum or standard annuity respectively. Guaranteed lifelong withdrawal benefit (GLWB), a specific type of GMWB, allows withdrawing funds at the contractual rate as long as the policyholder is alive. GMDB provides a specified payment if the policyholder dies. Precise specifications of the products within each type can vary across companies and some products may include combinations of these guarantees.
A good overview of VA products and the development of their market can be found in Bauer et al. (2008), Ledlie et al. (2008) and Kalberer and Ravindran (2009). There have been a number of papers in academic literature considering pricing of these products. Most of these are focused on pricing VA riders under the pre-determined (static) policyholder behaviour in withdrawal and surrender. Some studies include pricing under the active (dynamic) strategy when the policyholder 'optimally' decides the amount of withdrawal at each withdrawal date depending on the information available at that date. Standard Monte Carlo (MC) method can easily be used to estimate price in the case of pre-defined withdrawal strategy but handling the dynamic strategy requires backward in time solution that can be done only via the partial differential equation (PDE), direct integration or regression type MC methods.
In brief, pricing under the static and dynamic withdrawal strategies via PDE based methods has been developed in Milevsky and Salisbury (2006), Dai et al. (2008) and . Bauer et al. (2008) develops a unified approach with numerical estimation via MC and direct integration methods. The direct integration method was developed further in Luo and Shevchenko (2015a,b) using Gauss-Hermite quadrature and cubic interpolations. Bacinello et al. (2011) consider many VA riders under stochastic interest rate and stochastic volatility if the policyholder withdraws at the pre-defined contractual rate or completely surrenders the contract. Their pricing is accomplished either by the ordinary MC or Least-Squares MC to account for the optimal surrender. Typically, pricing of VA riders is considered under the assumption of geometric Brownian motion for the risky asset underlying the contract, though a few papers looked at extensions such as stochastic interest rate and/or stochastic volatility, see e.g. Forsyth and Vetzal (2014), Luo and Shevchenko (2016), Bacinello et al. (2011), Huang and Kwok (2015). Azimzadeh and Forsyth (2014) prove the existence of an optimal bang-bang control for GLWB contract when the contract holder can maximize contract writer's losses by only ever performing non-withdrawal, withdrawal at the contract rate or full surrender. However, they also demonstrate that the related GMWB contract does not satisfy the bang-bang principle other than in certain degenerate cases. Huang and Kwok (2015) developed a regression-based MC method for pricing GLWB under the bang-bang strategy in the case of stochastic volatility. GMWB pricing under the bang-bang strategy was studied in Luo and Shevchenko (2015c). The difficulty with applying the well known Least-Squares MC introduced in Longstaff and Schwartz (2001) for pricing VA riders under the optimal strategy is due to the fact that the paths of the underlying VA wealth account are affected by the withdrawals. In principle, one can apply control randomization methods extending Least-Squares MC to handle optimal stochastic control problems with controlled Markov processes recently developed in Kharroubi et al. (2014), but the accuracy and robustness of this method for pricing VA riders have not been studied yet.
One common observation in the above mentioned literature is that pricing under the optimal strategy often leads to prices significantly higher than observed on the market. These studies rely on the option pricing risk-neutral methodology in quantitative finance to find a fair fee. Here, the fundamental idea is to find the cost of a dynamic self-financing replicating portfolio which is designed to provide an amount at least equal to the payoff of the contract. The cost of establishing this hedging strategy is the no-arbitrage price of the contract. This is under the assumption that the contract holder adopts an optimal strategy (exercise strategy maximising the monetary value of the contract). If the purchaser follows any other exercise strategy, the contract writer will generate a guaranteed profit if continuous hedging is performed. Of course the strategy optimal in this sense is not related to the policyholder circumstances. In pricing VA with guarantees, it is reasonable to consider alternative assumptions regarding the investor's withdrawal strategy. This is because an investor may follow what appears to be a sub-optimal strategy that does not maximise the monetary value of the option. This could be due to reasons such as liquidity needs, tax and other personal circumstances. Moreover, mortality risk is diversified by the contract issuer through selling many contracts to many people while the policyholder cannot do it. Also, there might be no liquid secondary market for VAs on which the policy could be sold (or repurchased) at its fair value. The policyholder may act optimally with respect to his preferences and circumstances but it may be different from the optimal strategy that maximises the monetary value of the contract. In this case we calculate a fair fee to be deducted in order to finance a dynamic replicating portfolio for the guarantees (options) embedded in the contract under the assumption of a particular exercise strategy. The replicating portfolio will provide sufficient funds to meet any future payouts that arise from writing the contract.
However, the fair fee obtained under the assumption that investors behave optimally to max-imise the value of the guarantee does offer an important benchmark because it is a worst case scenario for the contract writer. Also, as noted in Hilpert et al. (2014), secondary markets for equity linked insurance products (where the policyholder can sell their contracts) are growing. Thus, third parties can potentially generate guaranteed profit through hedging strategies from financial products such as VA riders which are not priced under the assumption of the optimal withdrawal strategy. Knoller et al. (2015) mentions several companies recently suffering large losses related to increased surrender rates, indicating that either charged fees were not sufficiently large or that hedging program did not perform as expected. One way to analyze the withdrawal behavior of VA holder and evaluate the need of these products is to solve the life-cycle utility model accounting for consumption, housing, bequest and other real life circumstances. Developing a full life-cycle model with all preferences and required parameters is challenging but there are already several contributions reporting some interesting findings in this direction: Moenig (2012); Horneff et al. (2015); Gao and Ulm (2012); Steinorth and Mitchell (2015). This topic will not be considered in this paper. It is also important to note a recent paper by Moenig and Bauer (2015) considering the pricing under the optimal strategy in the presence of taxes via subjective risk-neutral valuation methodology. They demonstrated that including taxes significantly affects the value of the VA withdrawal guarantees producing results in line with empirical market prices.
In this paper we review pricing of living and death benefit guarantees offered with VAs, and present a unified optimal stochastic control framework for pricing these contracts. The main ideas have been developed and appeared in some forms in a number of other papers. However, we believe that our presentation is easier to understand and implement. We also present direct integration method based on computing the expected contract values in a backward timestepping through a high order Gauss-Hermite integration quadrature applied on a cubic spline interpolation. This method can be applied when transition density of the underlying asset between the contract cashflow event dates or its moments are known in closed form. We have used this for pricing specific financial derivatives and some simple versions of VA guarantees in Shevchenko (2014, 2015a). Here, we adapt and extend the method to handle pricing VA riders in general. As a numerical example, we calculate accurate prices of GMAB with possible annual ratchets (reset of the guaranteed capital to the investment portfolio value if the latter is larger on anniversary dates) and allowing optimal withdrawals. The contract that we consider is very similar in specifications to the real product marketed in Australia, see for example MLC (2014) and AMP (2014). Numerical difficulties encountered in pricing this VA rider are common across other VA guarantees and at the same time comprehensive numerical pricing results for this product are not available in the literature. These results (reported for a range of parameters) can serve as a benchmark for practitioners and researchers developing numerical pricing of VA riders.
In the next section, a general specification of VA riders is given. In Section 3 we discuss stochastic models used for pricing these products. Section 4 provides precise specification for some popular VA riders. In Section 5 we outline the calculation of the fair price and fair fee as a solution of an optimal stochastic control problem. Section 6 presents the numerical methods and algorithms for pricing VA riders. In Section 7 we present numerical results for the fair fees of GMAB rider. Concluding remarks are given in Section 8.

VA rider contract specification
Consider a VA contract with some guarantees for living and death benefits purchased by an x-year old individual at time t 0 = 0 with the up-front premium invested in a risky asset (e.g. a mutual fund), denoted as S(t) at time t ≥ 0. The VA rider specification includes dates when events such as withdrawal, ratchet (step-up), bonus (roll-up), death benefit payment, etc. may occur. Precise definitions of these events depend on the contract and corresponding examples will be provided in Section 4. We assume that the withdrawal can only take place on the set of the ordered event times T = {t 1 , . . . , t N }, where T = t N is the contract maturity. Also, the set of policy anniversaries when the ratchet may occur is denoted as T r and is assumed to be a subset of T . For simplicity on notation we assume that all other events may only occur on the withdrawal dates. The value of VA contract with guarantees at time t is determined by the three main state variables.
• Wealth account W (t), value of the investment account which is linked to the risky asset S(t) and modelled as stochastic process.
• Guarantee account A(t), also referred in the literature as benefit base. It is not changing between event times but can be stochastic via stochasticity in W (t) at the event times depending on the contract features.
• Discrete state variable I n ∈ {1, 0, −1} corresponding to the states of policyholder is being alive at t n , died during (t n−1 , t n ], or died before or at t n−1 correspondingly. Denote the death probability during (t n−1 , t n ] as q n = Pr[I n = 0|I n−1 = 1], i.e. Pr[I n = 1|I n−1 = 1] = 1 − q n . Note that q n depends on the age of the contract holder at t n and thus depends on the age x at t 0 = 0.
Other state variables are needed if the interest rate and/or volatility are stochastic but these are not affected at the contract event times and typically do not enter formulas for the contract cashflows; these will not be considered explicitly. Extra state variable is also required to track a tax free base to account for taxes; this will be considered in Section 5.4. In principle, different guarantees included in VA may have different benefit base state variables. For notational simplicity and also from practical perspective, we assume that all guarantees in VA are linked to the same benefit base account.
Initially, W (0) and A(0) are set equal to the upfront premium. The contract holder is allowed to take withdrawal γ n at time t n , n = 1, . . . , N − 1. Denote the values of the benefit base just before and just after t n as A(t − n ) and A(t + n ) respectively, and similarly for the wealth account W (t − n ) and W (t + n ). The contract product specification determines: • The contractual (guaranteed) withdrawal amount G n for the period (t n−1 , t n ] that may depend on the benefit base A(t − n ) and/or W (t − n ).
• Jump conditions at the event times relating state variables before and after the event, subject to withdrawals γ n belonging to an admissible space A n : where h W n (·) and h A n (·) are some functions that may also depend on the fee, penalty and annual step-up parameters. For example, if only a ratchet is possible at t n ∈ T r and no other contract events, then In practice, several events such as withdrawal, ratchet, bonus, etc. may occur at the same time t n , and the contract specification determines the order of this events.
• The payout P T (W, A) at the contract maturity if policyholder is alive at t = T .
• The payout D n (W, A) to the beneficiary at t n in the case of the policyholder death during (t n−1 , t n ], n = 1, . . . , N . • The cashflow received by the policyholder f n (W (t − n ), A(t − n ), γ n ) at the event times t n , n = 1, . . . , N − 1, that might be different from γ n due to penalties.
The specification details typically vary across different companies and are difficult to extract from the very long product specification documents. Moreover, results for specific GMxB riders presented in academic literature often refer to different specifications.
Once the above conditions, i.e. functions h W n (·), h A n (·), P T (·), P D (·), f n (γ n ) and admissible range for withdrawal A n are specified by the contract design, and a specific stochastic evolution of the state variables is assumed within (t n−1 , t n ), n = 1, . . . , N , then pricing of the contract can be accomplished by numerical methods. In particular, if withdrawals are optimal then pricing can be accomplished by PDE, direct integration or regression based MC methods. If withdrawals are deterministic, then standard MC along with PDE and direct integration methods can be used. The use of a particular numerical technique is determined by the complexity of the underlying stochastic model.

Stochastic Model
Commonly in the literature, stochastic models for the financial risky asset S(t) underlying the VA rider assume that there is no arbitrage in the financial market which means that there is a risk-neutral measure Q under which payment streams can be valued as expected discounted values. Moreover, this means that the cost of portfolio replicating the contract is given by its expected discounted value under Q. Hence, the fair price of the contract can be expressed as an expectation of the contract discounted cashflows with respect to Q. Some models considered in the literature assume that the financial market is complete which means that the risk-neutral measure Q is unique. It is also assumed that market has a risk-free asset that accumulates continuously at risk free interest rate. These are typical assumptions in the academic research literature on pricing financial derivatives, for a good textbook in this area we refer the reader to e.g. Björk (2004).
Regarding the mortality risk, it is assumed that it is fully diversified via selling the contract to many policyholders. In the case of systemic (undiversified) mortality risk, the risk-neutral fair value can be adjusted using an actuarial premium principle, see e.g. Gaillardetz and Lakhmiri (2011). Another common assumption is that mortality and financial risks are independent.
A benchmark model commonly considered in the literature on pricing VA riders is the well-known Black-Scholes dynamics for the reference portfolio of assets S(t) that under the risk-neutral measure Q is known to be Here, B(t) is the standard Wiener process, r(t) is the risk free interest rate and σ(t) is the volatility. Under this model the financial market is complete. Without loss of generality, the model parameters can be assumed to be piecewise constant functions of time for time discretization 0 = t 0 < t 1 < · · · < t N = T . Denote corresponding asset values as S(t 0 ), . . . , S(t N ) and risk free interest rate and volatility as r 1 , . . . , r N and σ 1 , . . . , σ N respectively. That is, σ 1 is the volatility for (t 0 , t 1 ]; σ 2 is the volatility for (t 1 , t 2 ], etc. and similarly for the interest rate. Pricing VA riders in the case of extensions of the above model to the stochastic interest rate and/or stochastic volatility have been developed in e.g. Forsyth and Vetzal (2014), Luo and Shevchenko (2016), Bacinello et al. (2011), Huang and Kwok (2015).
Regarding mortality modelling, the standard way is to use official Life Tables to estimate the death probability q n = Pr[I n = 0|I n−1 = 1] during (t n−1 , t n ]. Life Tables provide annual death probabilities for each age and gender in a given country; probabilities for time periods within a year can be found by e.g. linear interpolation, see Luo and Shevchenko (2015b). Instead of a Life Table, stochastic mortality models such as the benchmark Lee-Carter model introduced in Lee and Carter (1992) can also be used to forecasts the required death probabilities (accounting for systematic mortality risk).
For a given process of risky asset S(t), t ≥ 0, the value of the wealth account W (t) evolves as where dt n = t n − t n−1 and α is the annual fee continuously charged by contract issuer for the provided guarantee. In the case of S(t) following the geometric Brownian motion process (4), we have S(t n ) = S(t n−1 )e (rn− 1 where z 1 , . . . , z N are independent and identically distributed standard Normal random variables. In practice, the guarantee fee is charged discretely and proportional to the wealth account that can easily be incorporated into the wealth process (5). Denoting the discretely charged fee with the annual basis as α, the wealth process becomes Typically, the difference between continuously and discretely charged fees is not material as observed in our numerical results given in Section 7. Another popular fee structure corresponds to fees charged as a proportion of the benefit base, so that Here, it is assumed that discrete fees are deducted before withdrawal but it can be vice versa depending on the contract specifications. For simplicity, we do not consider management fees α m charged by a mutual fund for managing the investment portfolio. If management fees α m is given exogenously, then it will have an impact on the fair fee α that should by charged by the VA guarantee issuer. This can be accomplished as described in e.g. Forsyth and Vetzal (2014) and can be easily incorporated in the framework outlined in our paper. Obviously, α will be larger for given α m > 0 comparing to the case α m = 0. The management fees reduce the performance of the investment account thus increasing the value of the guarantee as reported in e.g.  for GMWB or Forsyth and Vetzal (2014) for GLWB. They commented that insurers wishing to provide the cheapest guarantee could provide the guarantee on the corresponding inexpensive exchange traded index fund rather than on a managed mutual fund account with extra fees.

VA riders
There are many different specifications for GMWB, GLWB, GMAB, GMIB and GMDB in the industry and academic literature. In this section we provide a mathematical formulation for some standard VA rider setups. We assume that the guarantee fee α is charged continuously. If the fee is charged discretely (and before withdrawal and other contract events), then one should make the following adjustment to the formulas in this section: if the fee is proportional to the wealth account and if the fee is proportional to the benefit base.

GMWB
A VA contract with GMWB promises to return at least the entire initial investment through cash withdrawals during the policy life plus the remaining account balance at maturity, regardless of the portfolio performance. Often in academic literature, the studied GMWB type has a very simple structure, where the penalty is applied to the cashflow paid to the contract holder, while the benefit base is reduced by the full withdrawal amount. Specifically, with γ n ∈ A n , A n = [0, A(t − n )]; and cashflow paid to the contract holder is where β ∈ [0, 1] is the penalty parameter for excess withdrawal. The contractual amount is defined as G n = W (0)(t n − t n−1 )/T and the maturity condition is Note that the above specification does not allow early surrender which can be included via extending the withdrawal space A n . Also, there is no death benefit; it is assumed that beneficiary will maintain the contract if the case of policyholder death. This contract has only basic features facilitating comparison of results from different academic studies, such as , Dai et al. (2008), Luo and Shevchenko (2015c), Luo and Shevchenko (2015a). Specifications common in the industry include cases where the contractual amount G n is specified to be different from G n = W (0)(t n − t n−1 )/T and a penalty is applied to both the withdrawn amount and the benefit base. For example, specifications used in Moenig and Bauer (2015) to compare with the industry products include: where x is the age of the policyholder in years at t 0 = 0, β e n and β g n are excess withdrawal and early withdrawal penalty parameters that can change with time, and Moenig and Bauer (2015) also considered several specifications for the benefit base jump conditions.
• Specification 1: • Specification 2: • Specification 3: In addition, a ratchet (reset of the benefit base to the wealth account if the latter is higher) can apply at anniversary dates. If it occurs before the withdrawal, then in the above formulas one should make the following adjustment If the reset is taking place after the withdrawal, then one should have

GLWB
GLWB is similar to GMWB but provides guaranteed withdrawal for life; upon death the remaining wealth account value is paid to the beneficiary. The contractual withdrawal amount G n is typically based on a fixed proportion g of the benefit base A(t), i.e. G n = g ×A(t − n )(t n −t n−1 ). The benefit base can increase via ratchet (step-up) or bonus (roll-up) features. Bonus feature provides an increase of the benefit base if no withdrawal is made on a withdrawal date. Complete surrender refers to the withdrawal of the whole policy account. The withdrawal can exceed the contractual amount and in this case the net amount received by the policyholder is subject to a penalty. Under the typical specification considered e.g. in Huang and Kwok (2015), the cashflow received by the policyholder is , where β is the penalty parameter for excess withdrawal. The benefit base jump condition, including ratchets and bonus features, is given by where b n is the bonus rate parameter that may change in time. Finally, if the policyholder dies during (t n−1 , t n ], the beneficiary receives a death benefit payment D n (W (t − n ), A(t − n )) = W (t − n ) and t N corresponds to the maximum age beyond which survival is deemed impossible.

GMAB
GMAB rider provides a certainty of capital till some maturity (e.g. 10 or 20 years) and the potential for a capital growth. Typical GMAB products sold on the market do not impose penalty on the policyholder withdrawal amount but can penalise the benefit base (protected capital balance) under some conditions. It is also common to have a ratchet feature, where the protected capital balance increases to the wealth account if the latter is higher on an anniversary date. Withdrawals from the account are allowed subject to a penalty. For example, specifications of the product marketed by MLC (2014) and AMP (2014) in Australia are very close to the following formulation: where C n (γ n ) is a penalty function that can be larger than γ n as defined below, and γ n ∈ A n = [0, W (t − n )]. The product is offered for the super and pension account types. The super account is designed for an investor being in an accumulation phase, while the pension account is for a retired investor in an annuitization phase. The difference between the accounts in terms of technical details is only in the penalty applied to the protected capital after withdrawals; the super account discourages withdrawals more than the pension account. In both cases the penalty is in the form of a reduction of the protected capital (benefit base) larger than the withdrawn amount. The penalty only applies if the wealth account balance is below the protected capital amount. A super account penalizes any amount of withdrawals, while the pension account only penalize excessive withdrawals.
Specifically, for a super account, the function C n (γ n ) is given by and for a pension account, the penalty is That is, the penalty for the pension account applies only if the wealth account balance is below the protected capital amount and the withdrawal is above a pre-determined amount G n .
Finally, the terminal condition is given by A total withdrawal of the wealth account balance effectively terminates the contract, as the penalty mechanism ensures the protected capital is always exhausted to zero by a complete withdrawal.

GMIB
At maturity, the holder of GMIB rider can select to take a lump sum of the wealth account W (T ) or annuitise this amount at annuitization rateä T current at maturity or annuitize the benefit base A(T − ) at pre-specified annuitization rateä g . Annuitization rate is defined as the price of an annuity paying one dollar each year. If the account value is below the benefit base, then the customer cannot take A(T − ) as a lump sum but only as an annuity at pre-specified rate. Thus, the payoff of VA with GMIB at time T is The benefit base may include roll-ups and ratchets. Again, this rider can be offered jointly with other riders. For example, it can be part of GMWB or GMAB contract maturity conditions. For discussion and pricing of GMIB in academic literature, see Marshall et al. (2010) and Bauer et al. (2008).

GMDB
GMDB rider provides a death benefit if the policy holder death occurs before or at the contract maturity. Assuming that if the policyholder dies during (t n−1 , t n ], then the beneficiary will be paid an amount D n (·) at t n , where some of the common death benefit types are: Some providers adjust the initial premium W (0) for inflation in the death benefit. For some policies, the death benefit type may change at some age, e.g. death benefit type 1 or type 2 may change to type 0, effectively making the death benefit expiring at some age (e.g. at the age of 75 years). The death benefit can be provided on top of some other guarantees and the contract may provide a spousal continuation option that allows a surviving spouse to continue the contract. The contract may have accumulation phase where the death benefit may increase, and continuation phase where the death benefit remains constant. Pricing GMDB has been considered in e.g. Milevsky and Posner (2001)

Fair Pricing
Denote the state vector at time t n before the withdrawal as X n = (W (t − n ), A(t − n ), I n ) and X = (X 1 , . . . , X N ). Given the withdrawal strategy γ = (γ 1 , . . . , γ N −1 ), the present value of the overall payoff of the VA contract with a guarantee is a function of the state vector Here, is the cashflow at the contract maturity and is the cashflow at time t n . Also, B i,j is the discounting factor from t j to t i

Pricing as Stochastic Control Problem
Let Q t (W, A) be the price of the VA contract with a guarantee at time t, when W (t) = W , A(t) = A and policyholder is alive. For simplicity of notation, if the policyholder is alive, we drop mortality state variable I n = 1 in the function arguments. Assume that financial risk can be eliminated via continuous hedging. Also assume that mortality risk is fully diversified via selling the contract to many people of the same age, i.e. the average of the contract payoffs where I is the real probability measure corresponding to the mortality process I 1 , I 2 , . . .. Then the contract price under the given withdrawal strategy γ can be calculated as Here, E Q,I t [·] denotes an expectation with respect to the state vector X, conditional on information available at time t, i.e. with respect to the financial risky asset process under the risk-neutral probability measure Q and with respect to the mortality process under the real probability measure I. Then the fair fee value of α to be charged for VA guarantee corresponds to Q 0 (W (0), A(0)) = W (0). That is, once a pricing of Q 0 (W (0), A(0)) for a given α is developed, then a numerical root search algorithm is required to find the fair fee.
The withdrawal strategy γ can depend on time and state variables and is assumed to be given when price of the contract is calculated in (25). The withdrawal strategies are classified as static, optimal, and suboptimal.
• 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, policyholder withdraws at the contractual rate only.
• Optimal strategy. Under the optimal withdrawal strategy, the decision on the withdrawal amount γ 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 where the supremum is taken over all admissible strategies γ. Any other strategy γ(X) different from γ * (X) is called suboptimal.
Given that the state variable X = (X 1 , . . . , X N ) is a Markov process and the contract payoff is represented by the general formula (21), calculation of the contract value (25) under the optimal withdrawal strategy (26) is a standard optimal stochastic control problem for a controlled Markov process. Note that, the control variable γ n affects the transition law of the underlying wealth W (t) process 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 Bäuerle and Rieder (2011). This type of problems can be solved recursively to find the contract value Q tn (x) at t n when X n = x for n = N − 1, . . . , 0 via the backward induction Bellman equation starting from the final condition Q T (x) = H N (x). Here, K tn (dx |x, γ n ) is the stochastic kernel representing probability to reach state in dx at time t n+1 if the withdrawal (action) γ n is applied in the state x at time t n . Obviously, the above backward induction can also be used to calculate the fair contract price in the case of static strategy γ; in this case the space of admissible strategies A n consists only one pre-defined value and sup(·) becomes redundant. For clarity, denote Q t − n (·) and Q t + n (·) the contract values just before and just after the event time t n respectively. Then, after calculating expectation with respect to the mortality state variable I n+1 in (27), the required backward recursion can be rewritten explicitly as with the jump condition This recursion is solved for n = N − 1, N − 2, . . . , 0, starting from the maturity condition

Alternative Solution
Given that the mortality and financial asset processes are assumed independent, and the withdrawal decision does not affect mortality process, one can calculate the expected value of the payoff (21) with respect to the mortality process, , and then calculate the price under the optimal strategy sup γ E Q t 0 [ H 0 (W , A)] or under the given strategy It is easy to find that where p n = Pr[τ > t n |τ > t 0 ] and q n p n−1 = Pr[t n−1 < τ ≤ t n |τ > t 0 ] for random death time τ , i.e. p n = p n−1 (1 − q n ). Note that, previously we defined q n = Pr[t n−1 < τ ≤ t n |τ > t n−1 ]. The payoff (30) has the same general form as the payoff (21). Thus, the optimal stochastic control problem Ψ t 0 (W (0), A(0)) = sup γ E Q t 0 [ H 0 (W , A)] can be solved using Bellman equation (27) leading to the following explicit recursion . It is easy to verify that this recursion leads to the same solution Ψ t 0 (W, A) = Q t 0 (W, A) and the same optimal strategy for γ as obtained from the recursion (28-29), noting that Note that, . That is, one cannot find the price under the optimal strategy conditional on the death time and then average over random death times, that would lead to the result larger than Q t 0 (W, A), see Luo and Shevchenko (2015b).

Remarks on Withdrawal Strategy
The guarantee fare fee based on the optimal policyholder withdrawal is the worst case scenario for the issuer, i.e. if the guarantee is hedged then this fee will ensure no losses for the issuer (in other words full protection against policyholder strategy and market uncertainty). Of course this is under the given assumptions about stochastic model for the underlying risky asset. If the issuer hedges continuously but investors deviate from the optimal strategy, then the issuer will receive a guaranteed profit.
Any strategy different from the optimal is sup-optimal and will lead to smaller fair fees. Of course the strategy optimal in this sense is not related to the policyholder circumstances. The policyholder may act optimally with respect to his preferences and circumstances but it may be different from the optimal strategy calculated in (29). On the other hand, as noted in Hilpert et al. (2014), secondary markets for equity linked insurance products (where the policyholder can sell their contracts) are growing. Thus, financial third parties can potentially generate guaranteed profit through hedging strategies from financial products such as VA riders which are not priced according to the worst case assumption of the optimal withdrawal strategy. Thus the development of secondary markets for VA riders would lead to an increase in the fees charged by the issuing companies. Knoller et al. (2015) undertakes an empirical study of policyholders behavior in Japanese VA market and they show that the moneyness of the guarantee has the largest explanatory power for the surrender rates.
One way to introduce a reasonable suboptimal withdrawal model is to assume that the policyholder follows a default strategy withdrawing a contractual amount G n at each event time t n unless the extra value from optimal withdrawal is greater than θ × G n , θ ≥ 0. Setting θ = 0 corresponds to the optimal strategy, while θ 1 leads to the strategy of withdrawals at the contract rate. This is the approach considered e.g. in Forsyth and Vetzal (2014) and . More complicated approach would specify a life-cycle utility model to determine the strategy optimal for the policyholder with respect to his circumstances and preferences, this is the approach studied in Moenig (2012); Horneff et al. (2015); Gao and Ulm (2012); Steinorth and Mitchell (2015). In any case, once the strategy is specified (estimated empirically or by another model), one can use (29) to calculate the fair price and fair fee with the admissible strategy space A n restricted to the specified strategy.

Tax consideration
Withdrawals from the VA type contracts may attract country and individual specific government taxes. Moenig and Bauer (2015) demonstrated that including taxes significantly affects the value of VA withdrawal guarantees. They developed a subjective risk-neutral valuation methodology and produced results in line with empirical market prices. Following closely to Moenig and Bauer (2015), we introduce an extra state variable R(t) to present the tax base which is the amount that may still be drawn tax-free, and assume that all event times t n ∈ T are the policy anniversary dates. The initial premium is assumed to be post-tax and taxes are applied to future investment gains (not the initial investment).
Denote a marginal income tax rate as κ and marginal capital gain tax from investment outside of VA contract as κ. It is assumed that earnings from VA are treated as ordinary income and withdrawals are taxed on a last-in first-out basis. Thus if the wealth account W (t − n ) exceeds the tax base R(t − n ), any withdrawal up to W (t − n ) − R(t − n ) will be taxed at the rate κ and will not affect the tax base; larger withdraws will not be subject to tax but will reduce the tax base. Specifically, the tax base will be changed at withdrawal time t n as The cashflow received by the policyholder will be reduced by taxes i.e. one has to make the following change in the contract specifications listed in Section 4 Using arguments for replicating pre-tax cashlows at t n with post-tax cashflows at t n+1 , it was shown in Moenig and Bauer (2015) that Q t + n (W, A, R) should be found not as the direct expectation (28) but should be found as the solution of the following nonlinear equation where This is referred to as subjective valuation from the policyholder perspective and depends on the investor current position (including possible offset tax responsibilities) and tax rates. Numerical examples in Moenig and Bauer (2015) show that the VA guarantee prices accounting for taxes in the above way are lower than ignoring the taxes (not surprisingly, because it lead to the suboptimal strategy), making the prices overall more aligned with those observed in the market.

Numerical Valuation of VA riders
In the case of realistic VA riders with discrete events such as ratchets and optimal withdrawals, there are no closed form solutions and fair price has to calculated numerically, even in the case of simple geometric Brownian motion process for the risky asset. In general, one can use PDE, direct integration or regression type MC methods, where the backward recursion (28-29) is solved numerically. Of course, if the withdrawal strategy is known, then one can always use standard MC to simulate state variables forward in time till the contract maturity or policyholder death and average the payoff cashflows over many independent realizations. This standard procedure is well known and no further discussion is needed.
In this section, we give a brief review of different numerical methods that can be used for valuation of VA riders. Then, we provide detailed description of the direct integration method that can be very efficient and simple to implement, when the transition density of the underlying asset or it's moments between the event times are known in closed form. Finally, in Section 6.5 we present calculation of hedging parameters (referred in the literature as Greeks). (2001) is designed for uncontrolled Markov process problems and can be used to account for the contract early surrender, as e.g. in Bacinello et al. (2011). However, it cannot be used if an optimal withdrawal strategy is involved. This is because dynamic withdrawals affect the paths of the underlying wealth account and one cannot carry out a forward simulation step required for the subsequent regression in the backward induction. However, it should be possible to apply control randomization methods extending Least-Squares MC to handle the optimal stochastic control problems with controlled Markov processes, as was recently developed in Kharroubi et al. (2014). The idea is to first simulate the control (withdrawals) and the state variables forward in time, where the control is simulated independently from other variables. Then, use regression on the simulated state variables and control to estimate expected value (28) and find the optimal withdrawal using (29). However, the accuracy and robustness of this method for pricing withdrawal benefit type products have not been studied yet. As usual, it is expected that the choice of the basis functions for the required regression step will have significant impact on the performance. We also note that in some simple cases of the withdrawal strategy admissible space such as bang-bang (no withdrawal, withdrawal at the contractual rate, or full surrender), it is possible to develop other modifications of Least-Squares MC such as in Huang and Kwok (2015) for pricing of the GLWB rider.

Simulation based Least-Squares MC method introduced in Longstaff and Schwartz
The expected value (28) can also be calculated using PDE or direct integration methods. In both cases, the modeller discretizes the space of the state variables and then calculates the contract value for each grid point. The PDE for calculation of expected value (28) under the assumed risk-neutral process for the risky asset S(t) is easily derived using Feynman-Kac theorem; for a good textbook treatment of this topic, see e.g. Björk (2004). However, the obtained PDE can be difficult or even not practical to solve in the high-dimensional case. In particular, in the case of geometric Brownian motion process for the risky asset (4), the governing PDE in the period between the event times is the one-dimensional Black-Scholes PDE, with jump conditions at each event time to link the prices at the adjacent periods. Since the benefit base state variable A(t) remains unchanged within the interval (t i−1 , t i ), i = 1, 2, . . . , N , the contract value Q t (W, A) satisfies the following PDE with no explicit dependence on A, This PDE can be solved numerically using e.g. Crank-Nicholson finite difference scheme for each A backward in time with the jump condition (29) applied at the contract event times. This has been done e.g. in Dai et al. (2008) and  for pricing GMWB with discrete optimal withdrawals. Of course, if the volatility or/and interest rate are stochastic, then extra dimensions will be added to the PDE making it more difficult to solve. Forsyth and Vetzal (2014) used PDE approach to calculate VA rider prices in the case of stochastic regime-switching volatility and interest rate.
Under the direct integration approach, the expected value (28) is calculated as an integral approximated by summation over the space grid points, see e.g. Bauer et al. (2008). More efficient quadrature methods (requiring less points to approximate the integral) exist. In particular, in the case of a geometric Brownian motion process for the risky asset, it is very efficient to use the Gauss-Hermite quadrature as developed in Luo and Shevchenko (2014) and applied for GMWB pricing in Luo and Shevchenko (2015a). Section 6.3 provides detailed description of the method for pricing VA riders in general. This method can be applied when the transition density of the underlying asset between the event times or it's moments are known in closed form. It is relatively easy to implement and computationally faster than PDE method because the latter requires many time steps between the event times. In Luo and Shevchenko (2016), this method was also used to calculate GMWB in the case of stochastic interest rate under the Vasicek model.
In both PDE and direct integration approaches, one needs some interpolation scheme to implement the jump condition (29), because state variables located at the grid points of discretized space do not appear on the grid points after the jump event. This will be discussed in detail in Section 6.4. Of course, if the underlying stochastic process is more complicated than geometric Brownian motion (4) and does not allow efficient calculation of the transition density or its moments, one can always resort to PDE method.
In our numerical examples of GMAB pricing in Section 7, we adapt a direct intergation method based on the Gauss-Hermite integration quadrature applied on a cubic spline interpolation, hereafter referred to as GHQC. For testing purposes, we also implemented Crank-Nicholson finite difference (FD) scheme solving PDE (36) with the jump condition (29).

Overall algorithm description
Both PDE and direct integration numerical schemes start from a final condition for the contract value at t = T − . Then, a backward time stepping using (28) or solving corresponding PDE gives solution for the contract value at t = t + N −1 . Application of the jump condition (29) to the solution at t = t + N −1 gives the solution at t = t − N −1 from which further backward in time recursion gives solution at t 0 . For simplicity assume that there are only W (t) and A(t) state variables. The numerical algorithm then takes the following key steps.

Algorithm 6.1 (Direct Integration or PDE method)
• Step 1. Generate an auxiliary finite grid 0 = A 1 < A 2 < · · · < A J to track the benefit base balance A. • Step 2. Discretize wealth account balance W space as W 0 < W 1 < · · · < W M to generate the grid for computing the expectation (28). • Step 4. Evaluate integration (28) for each A j , j = 0, . . . , J, to obtain Q t + N −1 (W, A) either using direct integration or solving PDE. In the case of direct integration method, this involves one-dimensional interpolation in W space to find values of Q t − N (W, A) at the guadrature points different from the grid points.
• Step 5. Apply the jump condition (29) to obtain Q t − N −1 (W, A) for all possible values of γ N −1 and find γ N −1 that maximizes Q t − N −1 (W, A). In general, this involves a two-dimensional interpolation in (W, A) space.

•
Step 7. Evaluate integration (28) for the backward time step from t 1 to t 0 to obtain solution Q 0 (W, A) at W = W (0) and A(0), or may be at several points if these are needed for calculation of some hedging sensitivities such as Delta and Gamma discussed in Section 6.5.
In our implementation of the direct integration method based on the Gauss-Hermite quadrature for numerical examples in Section 7, we use a one-dimensional cubic spline interpolation required to handle integration in Step 4 and bi-cubic spline interpolation to handle jump condition Step 5.
If the model has other stochastic state variables (similar to W ) changing stochastically between the contract event times, such as stochastic volatility and/or stochastic interest rate, then grids for these extra dimensions should be generated and the required integration or PDE to evaluate (28) will have extra dimensions. Also, extra auxiliary state variables (similar to A) unchanged between the contract event times, such as tax base and/or extra benefit base, will require extra dimensions in the grid and interpolation for the jump condition at the event times.
We have to consider the possibility of W (t) goes to zero due to withdrawal and market movement, thus one has to use the lower bound W 0 = 0. The upper bound W M should be set sufficiently far from the initial wealth at time zero W (0). A good choice of such a boundary could be based on the high quantiles of distribution of S(T ). For example, in the case of geometric Brownian motion process (5), one can set conservatively W M = W (0)e |mean(ln(S(T )/S(0)))|+5×stdev(ln(S(T )/S(0))) .
Often, it is more efficient to use equally spaced grid in ln W space. In this case, W 0 cannot be set to zero and instead should be set to a very small value (e.g. W 0 = 10 −10 ). Also, for some VA riders, using equally spaced grid in ln A space is also more efficient.

Direct integration method
To compute Q 0 (W (0), A(0)), we have to evaluate the expectations in the recursion (28). Assuming the conditional probability density of W (t − n ) given W (t + n−1 ) is known in closed form p n (w|W (t + n−1 )), the required expectation (28) can be calculated as where Q t − n (w, A) = B n−1,n (1 − q n )Q t − n (w, A) + q n D n (w, A) . The above integral can be estimated using various numerical integration (quadrature) methods. Note that, one can always find W (t − n ) as a transformation of the standard normal random variable Z as is the standard normal distribution, and F n (·) and F −1 n (·) are the distribution and its inverse of W (t − n ). Then, the integral (37) can be rewritten as This type of integrand is very well suited for the Gauss-Hermite quadrature that for an arbitrary function f (x) gives the following approximation Here, q is the order of the Hermite polynomial, ξ (q) i , i = 1, 2, . . . , q are the roots of the Hermite polynomial H q (x), and the associated weights λ (q) i are given by This approximate integration works very well if function f (x) is without singularities and it calculates the integral exactly if f (x) is represented by a polynomial of degree 2q − 1 or less. Note that Q t (w, ·) is known only at the grid points W m , m = 0, 1, . . . , M and interpolation is required to estimate Q t (w, ·) at the quadrature points. From our experience with pricing different VA guarantees, we recommend the use of the natural cubic spline interpolation which is smooth in the first derivative and continuous in the second derivative; and the second derivative is assumed zero for the extrapolation region above the upper bound.
Of course it can be difficult to find the distribution F n (·) and its inverse F −1 n (·) in general. In the case of geometric Brownian motion process (5), the transition density p n (·|·) is just a lognormal density and Then, a straightforward application of the Gauss-Hermite quadrature for the evaluation of integral (37) gives that should be calculated for each grid point W (t + n−1 ) = W m , m = 0, 1, . . . , M . Often, a small number of quadrature points is required to achieve very good accuracy; in our numerical examples in the next section we use q = 9 but very good results are also obtained with q = 5.
If the transition density function from W (t + n−1 ) to W (t − n ) is not known in closed form but one can find its moments, then the integration can also be done with similar efficiency and accuracy by method of matching moments as described in Shevchenko (2014, 2015a). The method also works very well in the two-dimensional case, see e.g. Luo and Shevchenko (2016) where it was applied for GMWB pricing in the case of stochastic interest rate.

Jump condition application
Either in PDE or direct integration method, one has to apply the jump condition (29) at the event times to obtain Q t − n (W, A). For the optimal strategy, we chose a value of withdrawal γ n ∈ A n maximizing the value Q t − n (W, A). To apply the jump conditions, an auxiliary finite grid 0 = A 1 < A 2 < · · · < A J = W M is used to track the remaining benefit base state variable A. For each A j , we associate a continuous solution using (40) and interpolation. In general, as can be seen from (29), the jump condition makes it impractical, if not impossible, to ensure the values of W and A after the jump to always fall on a grid point. Thus a two-dimensional interpolation is required. In this work we Figure 1: Illustration of the application of jump condition. The value Q t (W i , A j ) at t = t − n and at node point (W i , A j ) equals to Q t (W, A) at t = t + n with W = W i − γ n and A = A j − C(γ n ). The point (W, A) is located inside the grid bounded by (W m , W m+1 ) and (A k , A k+1 ).
adopted the bi-cubic spline interpolation for accuracy and efficiency. Figure 1 illustrates the application of jump conditions. It is natural to form a uniform grid in A so that optimal withdrawal strategies can be tested on a constant increment δA = A j+1 − A j , as has been done successfully in Luo and Shevchenko (2015a) for pricing of a basic GMWB specified by (8-9). However, extensive numerical tests show that if a uniform grid in A is used for pricing GMAB with ratchets and optimal withdrawals (our numerical example in Section 7), then neither linear interpolation nor cubic interpolation in A can achieve an efficient convergence in pricing results. A very fine mesh has to be used before we see a stable solution, which can take up to several hours to obtain a fair fee, in sharp contrast to basic GMWB where less than one minute computer time is required. On the other hand, if we make the grid in A uniform in Y = ln A and use linear or cubic interpolation based on variable Y , then we obtain a very good convergence on a moderately fine grid and the CPU time for a fair fee is about 30 minutes (a few minutes for a fair price). The CPU used for all the calculations in this study is Intel(R) Core(TM) i5-2400 @3.1GHz.
As we have already mentioned, a two-dimensional interpolation has to be used for applying the jump condition. We suggest to use either a bi-linear interpolation or a bi-cubic spline interpolation, e.g. see (Press et al., 1992, section 3.6), in both cases applied on the log-transformed state variables X = ln W and Y = ln A. For numerical examples in this paper, we have adapted the more accurate bi-cubic spline interpolation for all the numerical results.
For uniform grids, the bi-cubic spline is about five times as expensive in terms of computing time as the one-dimensional cubic spline. Suppose the jump condition requires the value Q(W, A) at the point (W, A) located inside a grid: W i ≤ W ≤ W i+1 and A j ≤ A ≤ A j+1 . Equivalently, the point X = ln W and Y = ln A is inside the grid: ln W i = x i ≤ X ≤ x i+1 = ln W i+1 and ln A j = y j ≤ Y ≤ y j+1 = ln A j+1 . Because the grid is uniform in both X and Y variables, the second derivatives ∂ 2 Q/∂X 2 and ∂ 2 Q/∂Y 2 can be accurately approximated by the three-point central difference, and consequently the one-dimensional cubic spline on a uniform grid involves only four neighboring grid points for any single interpolation. For the bi-cubic spline, we can first obtain Q(·, ·) at four points Q(W, A j−1 ), Q(W, A j ), Q(W, A j+1 ), Q(W, A j+2 ) by applying the one-dimensional cubic spline on the dimension X = log W for each point and then we can use these four values to obtain Q(W, A) through a one-dimensional cubic spline in Y = log A. Thus five one-dimensional cubic spline interpolations are required for a single bi-cubic spline interpolation, which involves sixteen grid points neighboring (W, A) point.

Calculating Greeks for hedging
Calculation of the contract price in (28) under the risk neutral probability measure Q means that one can find a portfolio replicating the VA guarantee, i.e. perform hedging eliminating the financial risk. Finding correct hedging depends on the underlying stochastic model for the risky asset. The basic hedging is the so-called delta hedging eliminating randomness due to stochasticity in the underlying risky asset S(t). Here, we use S(t) as a tradable asset to hedge the exposure of the guarantee to the wealth account W (t). One can construct a portfolio consisting of the money market account and ∆ S (t) units of S(t), so that ∆ S (t)S(t) = ∆ W (t)W (t), where ∆ W (t) is the number of units of the wealth account referred as Delta. Denote the value of the VA guarantee as U t (W, A) which is just a difference between the contract value with the guarantee Q t (W, A) and the value of the wealth account W , i.e.
Then, under the delta hedging strategy, one has to select for time t between the contract event times. Of course if there are extra stochasticities in the model such as stochastic interest rate and/or stochastic volatility, delta hedging will not eliminate risk completely and hedging with extra assets will be required which is model specific. See e.g. Forsyth and Vetzal (2014), for constructing hedging in the case of regime switching stochastic volatility and interest rate. A popular active hedging strategy in the case of extra stochastic factors is the minimum variance hedging strategy, where ∆ W (t) is selected to minimize the variance of portfolio's instantaneous changes, e.g. applied in Huang and Kwok (2015) for hedging GLWM in the case of stochastic volatility model. Practitioners also calculate other sensitivities (partial derivatives) of the contract with respect to the interest rate and volatility (referred to as Rho and V ega) and even second partial derivatives such as Gamma = ∂ 2 U t (W, A)/∂W 2 to improve hedging strategies. Here, we refer the reader to the standard textbooks in the area of pricing financial derivatives such as Wilmott (2013) or Hull (2006). Numerical estimation of the contract sensitivities (referred to as Greeks) is more difficult than estimation of the contract price. A general standard approach to calculate Greeks is to perturb the relevant parameter and re-calculate the price. Then one can use a two-point central difference to estimate the first derivatives and a three-point central difference for the second derivatives. In general, the finite difference PDE (or direct integration) methods generally produce superior accuracy in calculating Greeks when compared to Monte Carlo method (at least for low dimensions when finite difference method is practical or direct integration is possible). For Delta and Gamma, the finite difference method (or direct integration method) yields second order accurate values without re-calculating price using prices already calculated at the uniform grid points.
More accurate calculation of the main Greeks, Delta and Gamma, can be achieved using the so-called likelihood method as follows. The contract price at t 0 is calculated in the last time step (t 0 , t 1 ) in backward induction as an integral (37). Differentiating (37) with respect to W (0) = w 0 , Delta can be found as Thus it can be calculated using the same direct integration method as used for Q + t 0 (s 0 , A) with the factor ∂ ln p 1 (W |w 0 )/∂w 0 added to the integrand. Similarly, the required derivative to calculate Gamma can be found as Note, the above integrations for Delta and Gamma are only required for the (t 0 , t 1 ) time step and for a single grid point W (0) = w 0 . Here, t 0 should be understood as the current contract valuation time rather than time when the contract was initiated. Equivalently, for PDE approach using finite difference method, one can sometimes derive the corresponding PDEs for the Greeks and solve these PDEs for the last time step, see e.g. Tavella and Randall (2000). Similarly for Monte Carlo method, simulations used to calculate the price can be used to calculate Delta and Gamma weighted with extra factors under the expectations in (42) and (43).
It is illustrative to show how to derive hedging portfolio under the basic Black-Scholes model. Here, we assume that the underlying risky asset S(t) follows where B P (t) is the standard Brownian motion under the physical (real ) probability measure, and µ(t) is the real drift. Then the wealth account evolution is Here, we assume a continuously charged fee proportional to the wealth account but it is not difficult to deal with the case of discretely charged fees.
To hedge, the guarantee writer takes a long position in ∆ S units of S(t), i.e. forms a portfolio By Ito's lemma, the changes of portfolio within (t n−1 , t n ), n = 1, . . . , N are where the last term αW dt is the fee amount collected over dt. Setting eliminates all the random terms in (45), making the portfolio locally riskless. This means that the portfolio earns at risk free interest rate r(t), i.e. dΠ t = rΠ t dt, leading to the PDE Substituting U t (W, A) = Q t (W, A)−W in the above gives the PDE for Q t (W, A), the total value for the contract with the guarantee, i.e. the same as (36). Recalling Feynman-Kac theorem, it is easy to see that the stochastic process for W corresponding to this PDE is the risk-neutral process (4).

Numerical Example: GMAB pricing
Numerical solutions for pricing VA riders involve many complicated numerical procedures and features. These are more involved when compared to pricing of most exotic derivatives in financial markets. It is important that these solutions are properly tested and validated. As a numerical example for illustration, using direct integration method (GHQC), we calculate accurate prices of GMAB with possible annual ratchets and allowing optimal withdrawals as specified in Section 4.3. With these features, the GMAB is very close to the real product marketed in Australia by e.g. MLC (2014) and AMP (2014). We assume geometric Brownian motion model for the risky asset (4). When applicable, we compare results with MC and finite difference PDE methods. Numerical difficulties encountered in pricing this GMAB rider are common across other VA guarantees. Also, comprehensive numerical pricing results for this particular product are not available in the literature. These validated results (reported for a range of parameters) can serve as a benchmark for practitioners and researchers developing numerical pricing of VAs with guarantees. We consider four GMAB types: 1. GMAB with the annual ratchets but no withdrawals. In this case, a standard MC method can be used to compare with GHQC results -this is a good validation of the implemented numerical functions related to the ratchet feature, in addition to validating the numerical integration by Gauss-Hermite quadrature.
2. GMAB with the annual ratchets and a regular withdrawals of a fixed percentage of the wealth account. In this case a standard MC method can also be used to compare with GHQC results -this is a good validation of implemented numerical functions related to jump conditions due to both ratchet and withdrawal features. In addition, in order to test the numerical functions related to the application of penalties, we assume a pension account where the static withdrawal rate is set above the penalty threshold.
3. GMAB with the optimal quarterly withdrawals and the annual ratchets for a super account, where the penalty (18) is applied on any withdrawal γ n when W (t − n ) < A(t − n ).
4. GMAB with the optimal quarterly withdrawals and the annual ratchets for a pension account, where the penalty (19) is applied if the withdrawal γ n is above the penalty threshold G n and if W (t − n ) < A(t − n ).
As a comparison, results from our PDE finite difference implementation will also be shown for Case 4, the most complicated example among the four listed above. In addition, we will also calculate results for Case 4 in the case of the guarantee fee charged discretely (quarterly). All reported GHQC results are based on q = 9 quadrature points. We did not observe any material difference in results if q is increased further. Results based on q = 5 are also very accurate.

GMAB with ratchet and no withdrawal
Consider a GMAB rider with the annual ratchet and no withdrawals. In this case a standard MC method can be used to compare with GHQC results which is a good validation of implemented numerical functions related to the ratchet feature. Table 1 compares GHQC and MC results for the fare fee α of GMAB with the annual ratchet for the interest rate r ranging from 1% to 7% and the volatility σ = 10% and 20%. The maximum relative difference between the two methods is 0.76% at interest rate r = 5% and σ = 10%. The maximum absolute difference between the two methods is one basis point at the lowest interest rate r = 1% where the fee is the highest. On average, the relative difference is 0.52% and the absolute difference is 0.5 basis point, which is 5 cents per year on a one thousand dollar account. The GHQC results are obtained with the mesh size M = 400 and J = 200, and on the average it takes 22 seconds per price (calculation of the fare fee requires iterations over several prices). The MC results are obtained using 20 million simulations and it takes about 62 seconds per price. The agreement between the two methods at σ = 20% is also very good. In absolute terms, the maximum difference between the two methods is 1.1 basis point at the lowest interest rate r = 1%. In relative terms, the maximum difference between the two methods is 0.34% at the highest interest rate r = 7%. On average, the relative difference is 0.17% which is significantly smaller than the corresponding case at σ = 10%.

GMAB with ratchet and static withdrawal
Consider a GMAB rider with the annual ratchet and a regular quarterly withdrawals of a fixed percentage of the wealth account. In this case a standard MC method can also be used to compare with the GHQC results which is a good validation of implemented numerical functions related to jump conditions due to both ratchet and withdrawal features. Here, we consider a pension type account with the penalty given by (19). In this case, regular withdrawals at a pre-determined percentage level are allowed. In order to test the numerical functions related to the application of penalty, we also consider the static withdrawal above a pre-determined threshold level that will attract a penalty. We set the withdrawal threshold at 15% of the wealth account per annum, and the withdrawal frequency is quarterly, i.e. the quarterly withdrawal threshold is G n = 3.75% of the wealth account.
In the first test we allow a regular quarterly withdrawal of 3.75% of the wealth account balance, i.e. γ n = G n and there is no penalty on the withdrawals. Table 2 compares GHQC and MC results for the fare fee α. In relative terms, the maximum difference between the two methods is 0.08% at interest rate r = 5%. On average, the relative difference is 0.06% and the absolute difference is 0.3 basis point. The GHQC results are obtained with the mesh size M = 400 and J = 400, and the MC results are obtained with 20 million simulations per price.
Comparing with Table 1, the fair fee for the static withdrawal is about 8% higher than the corresponding no-withdrawal case at the lowest interest rate r = 1%, but it is about 13% lower than the corresponding no-withdrawal case at r = 7%. We have also tested static 2.5%  Table 2: Fair fee α in bp as a function of interest rate r for the GMAB with the annual ratchet and static quarterly withdrawal of 3.75% and 4% of the wealth account (15% and 16% annually respectively). The penalty threshold (pension type account) is set at 15% annually. The contract maturity is T = 10 years and volatility is σ = 20%. δ is the relative difference between Monte Carlo (MC) and GHQC results.
quarterly withdrawals and obtained the same pattern: at lower interest rate the fair fee of static withdrawal (which is also below the penalty threshold) is higher than the corresponding no-withdrawal case, and at higher interest it is the opposite. These differences in the fair fees at relatively low and high interest rates can be broadly interpreted as follows. At lower interest rates, where the expected capital growth is relatively slow, it is better to perform a regular withdrawal at or below the penalty threshold and take the protected capital at the maturity. However, at higher interest rates, where the expected capital growth rate is also high, it is beneficial not to carry out a regular withdrawal and keep the capital to grow. The above test also demonstrates that the MC and GHQC methods agree very well for pricing GMAB with a static withdrawal not exceeding the penalty threshold. This confirms the accuracy and efficiency of our numerical implementation of the jump condition using a bi-cubic interpolation in GHQC method.
In the second test of static withdrawal, we allow a regular quarterly withdrawal of 4% (16% per annum), i.e. the annual withdrawal rate is slightly higher than the penalty threshold of 15% per annum and there is a penalty applied for each withdrawal. The GHQC and MC results for this test are also presented in Table 2. In absolute terms, the maximum difference between the two methods is only 0.07 basis point at interest rate r = 6%, which is less than 1 cent per year for a one thousand dollar account. In relative terms, the maximum difference between the two methods is 0.1% at interest rate r = 6%. On average the relative difference is only 0.03%, which shows the two methods also agree very well in the case of excessive static withdrawals, where the penalty is applied for each withdrawal. This is a very convincing validation that our GHQC implementation of all numerical functions associated with the jump conditions, including the bi-cubic spline interpolation, is correct and accurate. The above tests are very close to validation of the entire algorithm in the case of optimal withdrawals. This is because in pricing the optimal withdrawal case, exactly the same integration and interpolation functions are used, and the only extra step required is to find the withdrawal rate maximizing the price. Nevertheless, for optimal strategy cases we will carry out some further validations by comparing results between GHQC and finite difference PDE methods.
Comparing Table 2 with Table 1 for the no-withdrawal case, the fair fee for the static withdrawal in excess of the penalty threshold is dramatically reduced: it is reduced by about 80% from the corresponding no-withdrawal case at r = 1%, and it is reduced by about 65% at r = 7%. Thus, a regular withdrawal above the penalty threshold is a very bad strategy regardless the interest rate level. In this instance, the penalty takes away the protected capital on a regular basis. Note that, the penalty is applied in terms of the whole withdrawal amount, not just on the exceeded part, see penalty function (19). Thus a slight excess over the penalty threshold can cause a large change in the price or in the fair fee, as observed in the second test.
The above two tests show that a regular static withdrawal is only slightly beneficial at very low interest rate and only when the withdrawal rate does not exceed the penalty threshold. In the next section, it will be demonstrated by numerical results that an optimal withdrawal is always beneficial regardless the interest rate level and penalties.

Optimal withdrawal -super account
Consider a GMAB for a super account with the annual ratchet and assume that the policyholder can exercise an optimal withdrawal strategy quarterly. For super account, any withdrawal will penalize the protected capital amount (benefit base) if the wealth account is below the benefit base according to (17) and (18).  Table 3: Fair fee α in bp (1 bp=0.01%) as a function of interest rate r for the GMAB on a super account when withdrawals are optimal. The contract maturity is T = 10 years. ε is the percentage difference between the fair fee in the case optimal withdrawal and the fair fee in the static case (no withdrawal) from Table 1. Table 3 shows the fair fee for a super account as a function of the interest rate at σ = 10% and σ = 20%. The columns under ε show an extra percentage value in the fee due to optimal withdrawal when compared to the static case of no withdrawal in Table 1. The results show, the extra fee is only about one or two percent for most cases, except at the low interest rate and high volatility. This extra fee due to optimal withdrawal is insignificant for the super account in most cases, mainly due to the heavy penalties applied. As will be shown in the next section, if the penalty is less severe as in the case of pension account, the extra fee becomes much more significant. If the penalty is completely removed, then numerical experiments show that the extra fee will be several times, e.g. 300%, larger than the static case, demonstrating the full value of the optimal strategy.

Optimal withdrawal -pension account
Consider a GMAB for a pension account with the annual ratchet and assume that the policyholder can exercise an optimal withdrawal strategy quarterly. For a pension account, any withdrawal above a pre-defined withdrawal level G n will penalize the protected capital amount A(t) if the wealth account W (t) is below A(t) according to (17) and (19). Here, we set the annual withdrawal limit at 15% of the wealth account, i.e. G n = 0.25 × 15% = 3.75% is quarterly withdrawal threshold.  Table 4: Fair fee α in bp (1 bp=0.01%) as a function of interest rate r for the GMAB for a pension account when withdrawals are optimal. The contract maturity is T = 10 years and quarterly withdrawal limit is G n = 0.25 × 15%. ε is the percentage difference between the fair fee in the case optimal withdrawal and the fair fee in the static case (no withdrawal) from Table 1. Table 4 shows the fair fee of GMAB for a pension account as a function of the interest rate r when σ = 10% and 20%.
The columns under ε show an extra percentage value in the fee due to optimal withdrawal when compared to the static case of no withdrawal in Table 1. The results show, the extra fee ranges from about 9% at the highest interest rate r = 7% to about 39% at the lowest interest rate r = 1%. This extra fee is much more significant than in the case of super account, see Table 3, apparently due to reduced penalties. At σ = 20%, the extra fee ranges from about 9% at the highest interest rate r = 7% to about 46% at the lowest interest rate r = 1%. This extra fee is higher than in the case of lower volatility σ = 10%, in both percentage and absolute terms.
Also, in Table 4, the numbers in the parentheses next to the continuous fair fee values α are the GHQC results for the discretely charged fair fee α d = − log(1 − αdt n )/dt n , where at the end of each quarter t n , the policyholder wealth account is charged a fee proportional to the account value αdt n W (t − n ), see the wealth process (6). Results show only little difference between the continuous fee α and the discrete fee α d . On average, the relative difference is 0.15%. Of course, at a higher frequency of charging fee (e.g. a monthly fee), the difference will be even less.
The last column in Table 4 shows the continuous fee calculated using our finite difference PDE method. The agreement between the GHQC (quadrature method) and finite difference PDE method is very good. The average relative difference in the fair fees between the two methods is 0.20%. Figure 2 show the curves of the fair fee for a pension account as a function of interest rate r using results from Table 4, in comparison with the static case (no withdrawal) from Table 1.
For some comparison, the market fees offered by MLC (2014) are 1.75% for a 10 year capital protection of a "balanced portfolio" and 0.95% for a "conservative growth portfolio"; and fees offered by AMP (2014) are 1.3% for a 10 year capital protection of a "balanced strategy" portfolio and 0.95% for a "moderately defensive strategy" portfolio. Though the values of volatility are not known for these market portfolios, it seems that market prices are significantly lower than the fair fee, which is also observed in the literature before; e.g. see Milevsky and Salisbury (2006), Bauer et al. (2008) and .  Table 4 and Table 1, respectively.

Conclusion
In this paper we have reviewed pricing VA riders and presented a unified pricing approach via an optimal stochastic control framework. We discussed different models and numerical procedures applicable in general to most of the VA riders with various contractual specifications. To price these VA riders under the geometric Brownian motion model for the risky asset, often assumed in practice, we have extended and generalized the direct integration method based on the Gauss-Hermite quadrature, introduced earlier in Luo and Shevchenko (2015a) for some specific and simpler product specifications.
As an example, we presented a numerical valuation of capital protection guarantees (GMAB riders), with specifications matching closely the real market products offered in Australia by e.g. MLC (2014) and AMP (2014). Numerical valuation of this guarantee involves all the main numerical difficulties encountered in pricing other VA riders, such as ratchets and optimal withdrawals. Numerical results have been validated by MC and finite difference PDE methods and can serve as a benchmark for practitioners and researchers developing numerical pricing of VA riders. As expected, we observed that the extra fee that has to be charged to counter the optimal policyholder behavior is most significant at lower interest rate and higher volatility levels, and it is very sensitive to the penalty withdrawal threshold.
As we have already discussed in Section 5.3, the fee based on the optimal policyholder withdrawal is the worst case scenario for the issuer, i.e. if the guarantee is hedged then this fee will ensure no losses for the issuer (in other words full protection against policyholder strategy and market uncertainty). If the issuer hedges continuously but investors deviate from the optimal strategy, then the issuer will receive a guaranteed profit. Any strategy different from the optimal is sup-optimal and will lead to smaller fair fees. Of course the strategy optimal in this sense is not related to the policyholder circumstances. The policyholder may act optimally with respect to his preferences and circumstances which may be different from the optimal strategy maximizing losses for the policy issuer. Life-cycle modelling can be undertaken to analyze and estimate sub-optimality of policyholder behavior. However, development of secondary markets for insurance products may expose the policy issuers to some significant risk if a fee for the guarantees is not charged to cover the worst case scenario.
It is important to note that the guarantee could be written on more that one asset (several mutual funds). In this case it is still common for practitioners to use a single-asset proxi model to calculate the price and hedging parameters. Obviously such approach has significant drawbacks (e.g. the sum of geometric Brownian motions is not a geometric Brownian motion). PDE and direct integration methods are not practical in high-dimensions and thus one has to rely on the MC methods to treat multi-asset case accurately. In the case of static withdrawal, it is not difficult to consider full multi-asset model and calculate the price using standard MC as in Ng and Li (2013). However, in the case of optimal withdrawal strategies, numerical valuation in the multi-asset case will require development of regression type MC solving backward recursion for processes affected by the withdrawals. One could apply control randomization methods extending Least-Squares MC developed in Kharroubi et al. (2014), but the accuracy robustness of this method for pricing VA riders have not been studied yet.
The specification details of VA riders typically vary across different companies and are difficult to extract and compare from the very long product specification documents. Moreover, results for specific GMxB riders presented in academic literature often refer to different specifications. As a result, cross-validation and benchmark research studies are rare. Given that numerical solutions used for pricing of VA riders are complex, it is important that these solutions are properly tested and validated. Moreover, new products are appearing in VA market regularly with increasing complexity that raises an important question, as discussed in Carlin and Manso (2011) and mentioned in Moenig and Bauer (2015), whether new complex products are designed to suite the policyholder needs better or introduced for the purpose of obfuscation.