Variable annuities in a L\'evy-based hybrid model with surrender risk

TL;DR

This paper introduces a market consistent valuation framework for variable annuities using a hybrid model with Le9vy processes, capturing dependence between financial and surrender risks, and providing explicit formulas and numerical methods.

Contribution

It develops a novel hybrid model incorporating Le9vy processes and surrender risk dependence, with explicit formulas and practical evaluation procedures.

Findings

Surrender behavior significantly impacts annuity pricing.

The model provides explicit analytical formulas for valuation.

Sensitivity analysis highlights the importance of accurate surrender risk modeling.

Abstract

This paper proposes a market consistent valuation framework for variable annuities with guaranteed minimum accumulation benefit, death benefit and surrender benefit features. The setup is based on a hybrid model for the financial market and uses time-inhomogeneous L\'evy processes as risk drivers. Further, we allow for dependence between financial and surrender risks. Our model leads to explicit analytical formulas for the quantities of interest, and practical and efficient numerical procedures for the evaluation of these formulas. We illustrate the tractability of this approach by means of a detailed sensitivity analysis of the price of the variable annuity and its components with respect to the model parameters. The results highlight the role played by the surrender behaviour and the importance of its appropriate modelling.

Tables

Table 1. Table 1. Parameters for the Variable Annuity contract. Reference process: NIG. Financial Market Model parameters - source: Eberlein and Rudmann ( 2018 ) . Mortality Model parameters - source: Escobar et al. ( 2016 ) .

Variable Annuity			Financial Market Model				Surrender Model			Mortality Model
				$L_t^1$	$L_t^2$
$T$	3, 4, 10 years		$\alpha$	4	5.73		$\beta$	0.05		$b$	12.1104
$\delta$	0.01 p.a.		$\beta$	-3.8	-2.13		$C$	0.01		$z$	76.139
$P\left(t_l\right)$	$0.95+0.05t_l/T$		$\delta$	1.34	8.3					$\kappa$	0.4806
$\Delta t_l$	1 year		$a$	0.0020898	-					$\lambda$	0.0195
${\overline{t}}_i-{\overline{t}}_{i-1}$	6 months		${\sigma }_2$	-	0.1818					$\sigma$	0.0254
			$b$	-	0.0065

Table 2. Table 2. Benchmarking Monte Carlo integration with importance sampling - the case of the GMAB. Parameters: Table 1 . ‘Quadrature’: Matlab built-in functions integral, integral2 and integral3. Bias/standard error expressed as percentage of the actual value. Monte Carlo iterations: 100 batches of size 10 6 superscript 10 6 10^{6} . CPU time expressed in seconds and referred to the average time of 1 batch of 10 6 superscript 10 6 10^{6} iterations.

			Quadrature	Monte Carlo integration (Imp. Sampling)
$T$	$K$		Value	Value	Bias (%)	Std. Error (%)
3 years	2	$A_1$	0.9867	0.9867	0.0035	0.0050
		$A_2$	0.1487	0.1482	0.3584	0.2683
		CPU	6.6323	31.2944
			($A_1$: 0.1280)
			($A_2$: 6.5043)
4 years	3	$A_1$	0.9703	0.9702	0.0139	0.0076
		$A_2$	0.1669	0.1669	0.0155	0.0647
		CPU	589.0926	51.6269
			($A_1$: 1.3042)
			($A_2$: 587.7884)

Table 3. Table 3. Variable Annuity by Monte Carlo integration with importance sampling. Parameters: Table 1 . Standard error expressed as percentage of the actual value. Monte Carlo iterations: 100 batches of size 10 6 superscript 10 6 10^{6} . CPU time expressed in seconds and referred to the average time of 1 batch of 10 6 superscript 10 6 10^{6} iterations.

$T$	GMAB	DB	SB	VA
4 years	112.5121	4.9280	2.6997	120.1399
(Std. Error %)	0.0001	0.0256	0.1142
CPU	51.6269	264.3202	31.7341
10 years	93.0783	14.2344	15.4533	122.7661
(Std. Error %)	0.0003	0.0098	0.0534
CPU	179.0842	1971.3826	717.3945

Table 4. Table 4. Benchmarking Monte Carlo integration with importance sampling. Parameters: Table 1 . ‘Quadrature’: Matlab built-in functions integral, integral2 and integral3. Bias/standard error expressed as percentage of the actual value. Monte Carlo iterations: 100 batches of size 10 6 superscript 10 6 10^{6} . CPU time expressed in seconds and referred to the average time of 1 batch of 10 6 superscript 10 6 10^{6} iterations.

GMAB			Quadrature	Monte Carlo integration (Imp. Sampling)
$T$	$K$		Value	Value	Bias (%)	Std. Error (%)
3 years	2	$A_1$	0.9867	0.9867	0.0035	0.0050
		$A_2$	0.1487	0.1482	0.3584	0.2683
		CPU	6.6323	31.2944
			($A_1$: 0.1280)
			($A_2$: 6.5043)
4 years	3	$A_1$	0.9703	0.9702	0.0139	0.0076
		$A_2$	0.1669	0.1669	0.0155	0.0647
		CPU	589.0926	51.6269
			($A_1$: 1.3042)
			($A_2$: 587.7884)
DB			Quadrature	Monte Carlo integration (Imp. Sampling)
$T$	$K$		Value	Value	Bias (%)	Std. Error (%)
4 years	3	$A_{1,1}^1$	0.9866	0.9865	0.0047	0.0051
	$(j,i)=(1,1)$	$A_{1,1}^2$	0.1122	0.1114	0.6452	2.0849
		CPU	9.0793	30.3230
			($A_{1,1}^1$: 0.4369)
			($A_{1,1}^2$: 8.6424)
		$A_{1,2}^1$	0.9866	0.9866	0.0040	0.0047
	$(j,i)=(1,2)$	$A_{1,2}^2$	0.1239	0.1247	0.6851	0.3349
		CPU	9.0793	30.2637
			($A_{1,2}^1$: 0.4369)
			($A_{1,2}^2$: 8.6424)
		$A_{2,1}^1$	0.9703	0.9704	0.0100	0.0074
	$(j,i)=(2,1)$	$A_{2,1}^2$	0.1349	0.1353	0.2873	0.2220
		CPU	834.0646	52.1560
			($A_{2,1}^1$: 1.1743)
			($A_{2,1}^2$: 832.8903)
		$A_{2,2}^1$	0.9703	0.9705	0.0159	0.0079
	$(j,i)=(2,2)$	$A_{2,2}^2$	0.1464	0.1461	0.2016	0.1338
		CPU	834.0646	52.0597
			($A_{2,2}^1$: 1.1743)
			($A_{2,2}^2$: 832.8903)
		$A_{2,3}^1$	0.9703	0.9703	0.0080	0.0069
	$(j,i)=(2,3)$	$A_{2,3}^2$	0.1570	0.1569	0.0378	0.0685
		CPU	834.0646	52.9188
			($A_{2,3}^1$: 1.1743)
			($A_{2,3}^2$: 832.8903)
		$A_{2,4}^1$	0.9703	0.9703	0.0012	0.0075
	$(j,i)=(2,4)$	$A_{2,4}^2$	0.1670	0.1669	0.0275	0.0588
		CPU	834.0646	52.9617
			($A_{2,4}^1$: 1.1743)
			($A_{2,4}^2$: 832.8903)
SB			Quadrature	Monte Carlo integration (Imp. Sampling)
$T$	$K$		Value	Value	Bias (%)	Std. Error (%)
4 years	3	$B_2^1$	0.9871	0.9871	0.0030	0.0029
	$\left(i=2\right)$	$B_2^2$	0.9717	0.9717	0.0066	0.0041
		CPU	1.1314	31.7341
			($B_2^1$: 0.1280)
			($B_2^2$: 1.0034)

Equations

$$\displaystyle E_{Q}\Big{[}e^{iuL^{j}_{t}}\Big{]}=\exp\bigg{(}\int_{0}^{t}\Big{(}iub_{s}^{j}-\frac{1}{2}c_{s}^{j}u^{2}+\int_{\mathbb{R}}\big{(}e^{iux}-1-iux{\mathds{1}}_{\{|x|\leq 1\}}\big{)}F^{j}_{s}(dx)\Big{)}ds\bigg{)}$$

$$\int_{0}^{T^{*}}\bigg{(}|b_{s}^{j}|+c_{s}^{j}+\int_{\mathbb{R}}(\min{\{|x|^{2},1\}})F_{s}^{j}(dx)\bigg{)}ds<\infty.$$

$${\int }_0^{T^{\ast }}{\int }_{\{|x|>1\}}{e}^{ux}{F}_s^j(dx)ds<\infty .$$

$$\displaystyle E_{Q}\Big{[}e^{iuL^{j}_{t}}\Big{]}=\exp\bigg{(}\int_{0}^{t}\Big{(}iub_{s}^{j}-\frac{1}{2}c_{s}^{j}u^{2}+\int_{\mathbb{R}}\big{(}e^{iux}-1-iux\big{)}F^{j}_{s}(dx)\Big{)}ds\bigg{)}.$$

$${\theta }_s^j(z)=b_s^{j}z+\frac{1}{2}{c}_s^j{z}^2+{\int }_{\mathbb{R}}(e^{zx}-1-zx)F_s^j(dx),$$

$$E_{Q}\big{[}\exp(zL^{j}_{t})\big{]}=\exp\bigg{(}\int_{0}^{t}\theta^{j}_{s}(z)ds\bigg{)}.$$

$$E_{Q}\bigg{[}\exp\bigg{(}\int_{0}^{t}f(s)dL^{j}_{s}\bigg{)}\bigg{]}=\exp\bigg{(}\int_{0}^{t}\theta^{j}_{s}(f(s))ds\bigg{)},$$

$$f(t,T)=f(0,T)+{\int }_0^t\alpha (s,T)ds-{\int }_0^t{\sigma }_1(s,T)d{L}_s^1,0\le t\le T\le T^{\ast },$$

$$B(t,T)=\exp\bigg{(}-\int_{t}^{T}f(t,s)ds\bigg{)}.$$

$$A(s,T)$$

$$B(t,T)=B(0,T)\exp\bigg{(}\int_{0}^{t}(f(s,s)-A(s,T))ds+\int_{0}^{t}\Sigma(s,T)dL^{1}_{s}\bigg{)},$$

$$\Sigma (s,T)\le \frac{M_1}{3},$$

$$A(s,T)={\theta }_s^1(\Sigma (s,T)),0\le s\le T.$$

$$S_{t}=S_{0}\exp\Big{(}\int_{0}^{t}r(s)ds+\int_{0}^{t}\sigma_{2}(s)dL^{2}_{s}+\int_{0}^{t}\beta(s)dL^{1}_{s}-\omega(t)\Big{)}.$$

$${\sigma }_2(s)\le \frac{M_2}{2},$$

$$\omega (t)={\int }_0^t[{\theta }_s^2({\sigma }_2(s))+{\theta }_s^1(\beta (s))]ds.$$

$$\max (I{S}_T,G(T)),$$

$$I{S}_{t_i}P(t_i),$$

$$\max (I{S}_{{\bar{t}}_i},G({\bar{t}}_i)),$$

$$Q\left({\tau }^m(x)>t\right)={\mathbb{E}}_Q\left(e^{-{\int }_0^t{\lambda }_u^m(x+u)du}\right),$$

$${\lambda }_t^m(x+t)={\lambda }^{m,0}(x+t){\xi }_t,$$

$${\lambda }^{m,0}(x+t)=\frac{1}{b}{e}^{\left(\frac{x+t-z}{b}\right)},$$

$$d\xi (t)=\kappa (\exp (-\lambda t)-\xi (t))dt+\sigma d{W}_t.$$

$$\displaystyle E_{Q}\Big{[}{\mathds{1}}_{\{\tau^{m}(x)>t\}}{\mathds{1}}_{A}\Big{]}=Q(A)E_{Q}\Big{[}e^{-\int_{0}^{t}\lambda^{m}_{u}(x+u)du}\Big{]}.$$

$$D(t)=Y_t-p(t)+{\int }_t^{T}f(t,s)ds-\delta T,0\le t\le T.$$

$${\lambda }^s(t)=\beta D^2(t_i)+C,t_i\le t<t_{i+1},$$

$$\displaystyle Q(\tau^{s}\geq t_{i}|{\mathscr{F}}_{t_{i}}^{L^{1},L^{2}})=\exp\Big{(}-\int_{0}^{t_{i}}\lambda^{s}_{u}du\Big{)},$$

$$Q({\tau }^s\ge t|{\mathcal{F}}_t^{L^1,L^2})$$

$$\mathrm{GMAB}(T)={\text{mathds}1}_{\{{\tau }^m(x)>T\}}{\text{mathds}1}_{\{{\tau }^s>T\}}\max (I{S}_T,G(T)).$$

$$\mathrm{SB}(t_i)={\text{mathds}1}_{\{{\tau }^s=t_i\}}{\text{mathds}1}_{\{{\tau }^s<{\tau }^m(x)\}}I{S}_{t_i}P(t_i),$$

Full text

Keywords: Finance; Variable Annuities; Hybrid models; Lévy processes; Surrender Risk

JEL Classification: G13, G12, G22, C63

(a) Faculty of Finance, Cass Business School, City, University of London, UK

(b) Freiburg Institute for Advanced Studies (FRIAS), Germany

(c) Department of Mathematical Stochastics, University of Freiburg, Germany

(d) University of Strasbourg Institute for Advanced Study (USIAS), France

*∗*Corresponding Author: Thorsten Schmidt; email: thorsten.schmidt@stochastik.uni-freiburg.de

1. Introduction

Variables Annuities (VAs) are unit-linked investment policies providing a post-retirement income, which is generated by the returns on a suitably managed financial portfolio. Various guarantees are applied with the aim of providing protection of the policyholders’ saving accounts. VAs are popular insurance products in the US, Japan, the UK, and are increasingly present in the other European markets as well. According to the Life Insurance and Market Research Association (LIMRA) Secure Retirement Institute and the Insured Retirement Institute (IRI), VAs sales for 2018 in the US were more than $100 billion - a 2% increase compared to 2017.

Common types of guarantees offered by variable annuity contracts are the so-called Guaranteed Minimum Accumulation Benefit (GMAB), the Guaranteed Minimum Death Benefit (henceforth DB) - which applies in case of early death - the Guaranteed Minimum Income Benefit and the Guaranteed Minimum Withdrawal Benefit. The former two offer protection during the accumulation period, i.e. up to the expiration of the contract, whilst the latter two provide payouts after expiration, during the so-called ‘annuitization’ period. For an extensive overview and classification of these products we refer to Bacinello et al. (2011) and references therein.

Due to the construction of these contracts, the underwriting insurance companies are exposed to financial and mortality risk, as well as surrender risk originated by the policyholder behaviour. Indeed, the option to leave the contract prior to maturity is a common additional feature of insurance contracts, which might cause significant cash outflows for the insurer, and negatively impact on the market growth of VAs (LIMRA Secure Retirement Institute).

The study of the pricing of life insurance contracts in presence of financial risk has been pioneered by Brennan and Schwartz (1976) and Boyle and Schwartz (1977); an extensive literature has developed since, from the seminal contributions of Albizzati and Geman (1994), Bacinello and Ortu (1996) and Grosen and Jørgensen (2002), to the more recent ones of Bacinello et al. (2011), Deelstra and Rayée (2013), Giacinto et al. (2014) and Gudkov et al. (2018), to mention a few. These contributions distinguish themselves in terms of the specific product under consideration, although they are all based on diffusion-driven market models. Extensions to financial dynamics driven by Lévy processes are considered in Ballotta (2005, 2009), although in these papers the focus is primarily on equity markets as the rate of interest is assumed constant and the possibility of surrender is not considered.

In light of the above, the aim of this paper is to provide a realistic framework for the modelling of these risks and the market consistent pricing of variable annuity contracts. We focus specifically on the pricing of the GMAB, the DB and the quantification of the surrender risk by means of the so-called Surrender Benefit (SB), contributing to the current state of the literature in a number of ways. Firstly, contrary to the literature mentioned above, we propose a general joint model for financial and insurance risks, which is driven by time-inhomogeneous Lévy processes. Our choice is motivated by the increased distributional flexibility offered by these stochastic processes in portraying observed market trends (see, for example, Eberlein and Keller, 1995, for an extensive empirical analysis of equity markets).

In more details, for the financial market we adopt the hybrid construction of Eberlein and Rudmann (2018) in which the stochastic dependence between interest rate markets and stock markets is taken into account explicitly. This feature is an important aspect in this context given the typically long maturity of insurance contracts.

Further, we subdivide insurance risk into surrender risk and mortality risk. For the surrender risk, we follow market practice by considering two components: one capturing the baseline surrender behaviour due to non-economic factors and personal contingencies, and the second one which instead is responsive to changes in the conditions of the financial markets (see Kolkiewicz and Tan, 2006, Le Courtois and Nakagawa, 2011, Ducuroir et al., 2016, for example). This component in particular includes a function of the spread between the rate of the contract and the rate offered on the market for equivalent products, and therefore incorporates the stochastic inputs from the financial market model. The function of choice is designed as to suitably accommodate for surrender triggering arguments based on the moneyness hypothesis (see Knoller et al., 2016, for example), the interest rate hypothesis and the emergency fund hypothesis (see, for example, Nolte and Schneider, 2017, and references therein).

For mortality risk, we adopt an extended Gompertz-Makeham model with stochastic mortality improvement ratio as in Krayzler et al. (2016) and Escobar et al. (2016); however we make the natural assumption of stochastic independence between the demographic and the financial risks.

Secondly, through the proposed general framework we obtain closed analytical formulas (up to a multidimensional integral) for the price of the guarantees and the SB. The dimensionality of these integrals is dictated by the frequency with which the policyholder is allowed to terminate the contract.

Finally, we also develop a practical and efficient numerical scheme for the evaluation of these formulas by means of Monte Carlo integration with importance sampling, and we illustrate its applicability by performing a sensitivity analysis of the contract’s value with respect to the model parameters. The results from this analysis underline the importance of a correct calibration of the model to the observed surrender rates.

The paper is organized as follows. In section 2 we introduce the model for the financial market; the contract features of the VAs are introduced in section 3 together with the additional modelling assumptions regarding mortality and surrender risk. The closed analytical pricing formulas are derived in section 4; the corresponding numerical scheme is offered in section 5. We present the results of the sensitivity analysis in section 6, whilst section 7 concludes.

Additional material, including some of the proofs, is provided in the appendix.

2. The interest rate and equity market

The aim of this section is to introduce a joint model for the interest rate and the equity markets used to develop the pricing framework for variable annuities. Specifically, we follow Eberlein and Rudmann (2018) and build this joint model on time-inhomogenous Lévy processes. Thus we first offer some necessary preliminary results with the pricing application in view.

2.1. Time-inhomogenous Lévy process and pricing.

Given a finite time horizon $T^{*}>0$, consider a stochastic basis $(\Omega,{\mathscr{F}},\mathbb{F},Q)$ with a filtration $\mathbb{F}=({\mathscr{F}}_{t})_{t\in[0,T^{*}]}$ satisfying the usual conditions. Due to our focus on pricing, the probability measure $Q$ represents a risk-neutral martingale measure. Let $L^{1}$ and $L^{2}$ be two independent time-inhomogeneous Lévy processes, i.e. continuous-time processes with independent increments, and with characteristic function

[TABLE]

for $j=1,2$; the local characteristic $(b_{s}^{j},c_{s}^{j},F_{s}^{j})_{s\in[0,T^{*}]}$ satisfies the integrability condition

[TABLE]

We note that a financial market built on (exponential) Lévy processes is in general incomplete, and consequently the risk neutral martingale measure is not unique. Standard practice in this case is to fix the pricing measure via calibration, using market quotes for traded derivatives contracts written on the quantities of interest, i.e. bonds and stocks in the specific case of our application. As model calibration is not in the scope of this paper, we refer to Eberlein and Rudmann (2018) for a detailed illustration of the calibration procedure for the adopted joint model. Nevertheless, for any risk neutral martingale measure to be well defined, we require that the exponential moments of a certain order exist. To this purpose, the following assumption holds throughout the rest of this paper.

Assumption 2.1** (Exponential moments).**

For $j=1,2$, there exist positive constants $M_{j}$ and $\epsilon_{j}$ such that for each $u\in[-(1+\epsilon_{j})M_{j},(1+\epsilon_{j})M_{j}]$

[TABLE]

All standard processes used in mathematical finance such as hyperbolic, normal Inverse Gaussian, Variance Gamma and CGMY processes satisfy the above condition. Assumption 2.1 implies the existence of the first moment of the process; this allows us to rewrite (1) as

[TABLE]

In this setting, we define the cumulant function of $L^{j}$

[TABLE]

for any $z\in\mathbb{C}$ such that $Re(z)\in[-(1+\epsilon_{j})M_{j},(1+\epsilon_{j})M_{j}]$. Then, $E_{Q}[\exp(zL^{j}_{t})]<\infty$ and

[TABLE]

Further, let $f:\mathbb{R}^{+}\to\mathbb{C}$ be a continuous function with $|Re(f)|\leq M_{j}$, then

[TABLE]

where the integrals are defined component-wise for the real and imaginary part (see Eberlein and Raible, 1999, for full details).

2.2. The fixed income market.

For the modelling of the fixed income market we follow the approach in Eberlein et al. (2005) (see also Eberlein and Raible, 1999), so that the starting point is the definition of the dynamics of the instantaneous forward rates $(f(t,T)_{0\leq t\leq T\leq T^{*}})$. Let us assume that

[TABLE]

for a deterministic and bounded function $f(0,T)$. The drift function $\alpha(\cdot)$ and the volatility function $\sigma_{1}(\cdot)$ are assumed to satisfy the usual conditions of measurability and boundedness (see Eberlein et al., 2005, (2.5)). The price of a zero coupon bond at time $t$ with maturity $T\geq t$ is

[TABLE]

Let us denote

[TABLE]

From Fubini’s theorem and equation (3) it follows that the dynamics of the bond price is

[TABLE]

We remind that the short rate $r_{t}$ is implicitly given by the forward rate dynamics in equation (3) by setting $T=t$. Finally, we assume that

[TABLE]

where $M_{1}$ is the constant from Assumption 2.1; this guarantees that the exponential of the stochastic integral has finite expectation.

For the market to be arbitrage free, we require that $(B_{t}^{-1}B(t,T))_{0\leq t\leq T}$, with B_{t}=\exp\big{(}\int_{0}^{t}r(s)ds\big{)}, is a martingale; it follows from (2) that $Q$ is a martingale measure if

[TABLE]

In the following we will always assume that the drift condition (5) holds.

2.3. The stock market

For the modelling of the equity market, we consider the case of a single asset, be it a single stock or a stock index. More general settings can be obtained in a straightforward manner.

It is well known from a number of empirical studies that equity and fixed income markets influence each other; this interaction is of particular importance in the context of long-dated insurance contracts. Thus, following Eberlein and Rudmann (2018) we choose an approach which allows for stochastic dependence between the two markets. Consequently, we model the price process of the asset as

[TABLE]

In this hybrid approach, the driver of the interest rate dynamics affects explicitly the stock price in reason of $\beta(\cdot)$ through the term $\int_{0}^{t}\beta(s)dL^{1}_{s}$. Further dependence is originated endogenously by the (integrated) short rate, i.e. the classical risk neutral term $\int_{0}^{t}r(s)ds$. $\sigma_{2}(\cdot)$ is a positive function and denotes the volatility of the stock price. Both, $\sigma_{2}(\cdot)$ and $\beta(\cdot)$, could be chosen as random processes, but having numerical aspects in mind, in the following we consider deterministic functions $\sigma_{2}(\cdot)$ and $\beta(\cdot)$. To ensure the existence of exponential moments we require

[TABLE]

with $M_{1},M_{2}$ the constants from Assumption 2.1. The drift term $\omega(t)$ in (6) is chosen such that the discounted stock price $(B_{t}^{-1}S_{t})_{t\in[0,T^{*}]}$ is a $Q$-martingale. Equation (2) and the independence of the two driving processes $L^{1}$ and $L^{2}$ imply that

[TABLE]

Hence, under (5) and (8) the joint market model for bonds and security $S$ is free of arbitrage.

3. The variable annuity contract

A variable annuity (VA) is an insurance contract which gives the holder a variety of benefits depending on the notional $I$ in exchange for an initial premium which we will determine in the following. This premium is paid once at inception of the VA contract. The maturity $T$ of the contract is assumed to satisfy $0<T\leq T^{*}$. In the specification considered here, the VA includes three features: a Guaranteed Minimum Accumulation Benefit (GMAB), a Surrender Benefit (SB), and a Death Benefit (DB).

In details, at maturity $T$ the GMAB pays to the policyholder

[TABLE]

with $G(T)=I\exp(\delta T)$, $\delta>0$. In other words, the GMAB offers the best of the investment of an amount $I$ in either the asset $S$ or in a risk-free account with guaranteed rate $\delta$. Note that, in order to simplify the notation, in the following we assume that the price process of the asset is normalized so that $S_{0}=1.$ This payoff, however, can only be claimed if the policyholder is still alive at time $T$ and did not exercise the surrender option before.

In case of early surrender, the right of refund is restricted to the current value of the fund reduced by a compulsory surrender penalty. Let the penalty $P:[0,T]\to(0,1]$ be an increasing function with $P(T)=1$, and define ${\bf t}:=(t_{0},t_{1},\ldots,t_{K})^{\top}$ with $0=t_{0}<t_{1}<\ldots<t_{K}<T$. We assume that premature surrender is possible at any time point $t_{i}\in{\bf t}$ with $i=1,\ldots,K-1$, in which case the policyholder would receive the amount

[TABLE]

Small values of the penalty $P$ at early dates would allow the insurer to recover any expenses related to the writing of the contract even in case of early surrender. The penalty could also serve to the insurer to hedge against a significant rise in the fund value.

Finally, in case of death before maturity $T$, the death benefit pays (to the beneficiaries)

[TABLE]

for $\bar{t}_{i}\in\bar{\bf t}$, $i=1,\ldots,N$, and $\bar{\bf t}:=(\bar{t}_{0},\bar{t}_{1},\ldots,\bar{t}_{N})^{\top}$ with $0=\bar{t}_{0}<\bar{t}_{1}<\ldots<\bar{t}_{N}=T$. The time points $\bar{t_{i}}$ denote the frequency with which mortality is monitored by the insurer over the lifetime of the contract.

In principle, the time scales ${\bf t}$ and $\bar{\bf t}$ could be arbitrary, but in typical cases they are not: we assume that ${\bf t}\subset\bar{\bf t}$ in the sense that any surrender time $t_{i}$ is contained in $\{\bar{t}_{1},\dots,\bar{t}_{N}\}$. This is a very natural assumption since death of the policyholder might be monitored by the insurer at the end of every month or every quarter, whereas surrender of the contract might be allowed only at the end of each year during the life of the contract, or at the policy anniversary.

The description of the contract benefits given above highlights the need for an accurate financial model, but also an appropriate modelling for mortality risk and surrender risk. This is offered in the rest of this section, which we conclude with the discussion of the market consistent valuation of the VA.

3.1. Mortality model

For the modelling of mortality risk, we follow standard literature in the field and adopt a stochastic intensity-based approach. Pioneered by Milevsky and Promislow (2001), further developed by Dahl (2004), Dahl and Møller (2006), and more recently generalized by Li and Szimayer (2011), this framework builds on a given initial curve for mortality rates by superimposing a stochastic process capturing random improvements and fluctuations.

In details. Let $\tau^{m}(x)$ be a random time capturing the remaining lifetime of a $x$ years old individual. The corresponding survival probability with respect to the given probability measure is

[TABLE]

where $\lambda_{t}^{m}(x+t)$, $t>0$, is the stochastic mortality intensity for an individual aged $x+t$ at time $t$. This intensity process is modelled as

[TABLE]

for an initial curve of the mortality intensity, $\lambda^{m,0}$, and a strictly positive process $\xi_{t}$ such that $\xi_{0}=1$, capturing the mortality improvements from time 0 to time $t$ for a person aged $x+t$. Finally, the mortality intensity satisfies the property that $\lambda^{m}_{0}(x)=\lambda^{m,0}(x)$ (see Dahl and Møller, 2006 and references therein for full details).

A popular choice for the initial mortality curve is represented by the Gompertz-Makeham model (Dahl, 2004, Dahl and Møller, 2006, for example). An alternative specification corresponding to the structure for the UK standard tables for annuitant and pension population is adopted in Ballotta and Haberman (2006). Numerous choices for the process of the mortality improvement ratio, $\xi_{t}$, have been put forward as well: Dahl (2004), Biffis (2005), Dahl and Møller (2006) for example focus on affine diffusion specifications with particular emphasis on time-inhomogeneous CIR processes; Ballotta and Haberman (2006) instead extend the generalized linear model by superimposing a standard mean reverting Ornstein-Uhlenbeck process and suitably accommodating longevity effects.

In light of the above, in the following we choose the standard Gompertz-Makeham model for the initial curve, so that

[TABLE]

and an extended Ornstein-Uhlenbeck with mean reversion level $e^{-\lambda t}$ for the mortality improvement ratio as in Krayzler et al. (2016), i.e.

[TABLE]

We assume that $W$ is a Brownian motion independent of $(L^{1},L^{2})$, $z$, $\kappa$ and $\sigma$ are non-negative, $b$ is positive and $\lambda\in\mathbb{R}$. This implies that the mortality intensity is independent of the financial market.

Let us denote by $\mathbf{F}^{L^{1},L^{2}}=({\mathscr{F}}^{L^{1},L^{2}}_{t})_{0\leq t\leq T^{*}}$ the filtration generated by $L^{1}$ and $L^{2}$. Then, for any set $A\in{\mathscr{F}}^{L^{1},L^{2}}_{t}$, it follows that

[TABLE]

A few considerations are in order. Firstly, as $\xi$ is modelled as a Gaussian Ornstein-Uhlenbeck process, the mortality intensity can become negative with positive probability. However, this probability is negligible (less than $10^{-5}$) as shown in Appendix A.2. in Escobar et al. (2016). For the occurring difficulties with negative values in the intensity, see for example Bielecki and Rutkowski (2002), Lemma 9.1.4 and the related remarks. Secondly, the above assumption of independence between the mortality risk driver, $W$, and the financial risk drivers, $L^{1}$ and $L^{2}$, is plausible from the perspective of the insurer: precise modelling would be highly client-specific, and would require information rarely accessible by insurance companies. Moreover, the independence assumption implies a high degree of tractability which is important for complex products such as VAs considered here.

Finally, as a consequence of the independence between demographic and financial risks, the computation of the survival probability is carried out under the risk neutral measure $Q$. For specific considerations about the interaction between mortality risk and market risk, we refer the interested reader to Dahl and Møller (2006) and references therein.

3.2. Surrender model

Surrender risk is notoriously difficult to assess and model due to the nature of the decisions leading to it. Policyholders can indeed surrender both because of rising alternative financial opportunities, and apparently non rational (in the financial sense) behaviours due to personal considerations and contingencies. Nevertheless, surrender represents one of the main risks faced by life insurance companies due to the liquidity issues it can generate, with potential loss of market share (see, for example, Loisel and Milhaud, 2011, and references therein).

A common market and academic practice to surrender modelling (see Kolkiewicz and Tan, 2006, Le Courtois and Nakagawa, 2011, Ducuroir et al., 2016, for example) is to consider two components: a deterministic baseline hazard rate function capturing lapses111‘Lapse’ was originally used to denote termination of an insurance policy and loss of coverage because the policyholder had failed to pay premia, whilst ‘surrender’ denotes termination accompanied by the payout of a surrender benefit. Nowadays ‘lapse’ often denotes both situations. due to non-economic factors, and a stochastic process representing additional shocks to the baseline due to changes in the market.

This random component is usually linked to the spread between the return offered by the contract and the one offered on the market for equivalent products. Indeed, dependence of the surrender on interest rates is relatively intuitive: higher interest rates are strong incentives for the policyholder to switch to higher yield investments, whilst very low interest rates - such as the ones currently observed in all major economies - could represent advantageous opportunities for refinancing. The random component should also be linked to the change in value of the underlying asset, as it directly impacts on the amount received in case of surrender.

Thus, following this line of reasoning, let $\tau^{s}$ denote the random time of the policyholder decision to surrender. As specified above, surrender is allowed at time points $t_{i}$, $i=1,\dots,K-1$; by convention $\{\tau^{s}=\infty\}$ corresponds to no surrender. Let $\lambda^{s}$ denote the corresponding intensity of surrender; then $\lambda^{s}(t)=0$ for $t\in[0,t_{1})\cup[t_{K},T]$. Further, we model the baseline surrender behaviour due to non-economic factors and personal contingencies by a non-negative constant $C$.

Let $D(t)$ be the process driving the dynamic lapse component; consistently with the considerations offered above, we build this process on the spread between the return offered by the surrender benefit (net of any penalty charge) plus the market rates at which this amount can be invested, and the total yield of the policy represented by the value at maturity of the guaranteed amount. Thus, let $Y_{t}=\log S_{t}$, and $p(t)=-\log P(t)$. Then

[TABLE]

The overall surrender intensity consequently is defined as

[TABLE]

so that it is piecewise constant on the interval $[t_{i},t_{i+1})$ for $i=1,\dots,K-1$.

The non-negative constant $\beta$ captures the dependence between the surrender intensity and the market and is a measure of the investors’ rationality (in the pure economic sense), and their response to personal financial motivations. Equation (12) uses the square of the spread function $D(t)$ in order to capture both situations of favourable market conditions offering more remunerative investment opportunities, and financial market turmoils in which policyholders might lack sufficient resources to finance their expenses (emergency fund hypothesis).

Although in spirit similar to others in the literature (see for example Le Courtois and Nakagawa, 2011 and Escobar et al., 2016), our construction distinguishes itself also for the non-Gaussian dynamics of the underpinning risk drivers $L^{1}$ and $L^{2}$.

The resulting probability of no surrender is given by

[TABLE]

for all $1\leq i\leq K-1$, and

[TABLE]

for $t_{K-1}<t$. We note that the last integral equals $\exp(-\int_{0}^{T}\lambda^{s}_{u}du)$. The set-up chosen here can be obtained in a doubly-stochastic model or in a model where immersion holds, see Aksamit and Jeanblanc (2017) for a comprehensive treatment in this regard. An alternative form of the intensity is investigated in the Appendix A.

Finally, we observe that the intensity functions $\lambda^{s}$ and $\lambda^{m}$ are independent due to the assumed independence between demographic and financial risks.

3.3. The price of the variable annuity

Using the notation introduced above, we can now formulate the actual cash-flows associated with the considered variable annuity. Firstly, recall that the GMAB provides a payoff only at the maturity time $T$ if the policyholder is alive (i.e., $\{\tau^{m}(x)>T\}$) and if there was no surrender until this time (i.e., $\{\tau^{s}>T\}$). Consequently, the associated cash-flow at maturity $T$ is

[TABLE]

Secondly, the surrender option can be exercised only once if the policyholder is still alive (i.e., $\{\tau^{s}<\tau^{m}(x)\}$). Therefore, should surrender occur, the surrender benefit at time $t_{i}$ pays

[TABLE]

where $t_{i}$ is one of the possible premature surrender dates. Finally, the death benefit provides a payoff only in case of no early surrender, and is quantified as

[TABLE]

where $\bar{t}_{i}$ is one of the possible payoff dates of the death benefit.

The price ${\rm P}^{\rm VA}$ of the variable annuity at time $t=0$ is equal to the sum of the prices of its constituents, i.e.

[TABLE]

with, by standard risk neutral valuation argument

[TABLE]

[TABLE]

[TABLE]

Tractable pricing formulae for these expressions are provided in the following section.

4. Market-consistent valuation

In this section we derive analytical expressions for the components of the variable annuity contract discussed above in the model setup provided in section 2 in presence of mortality and surrender risk. Useful results and representations needed in the following are provided in the Appendix B.

4.1. Guaranteed minimum accumulation benefit

The independence between $\tau^{m}(x)$ and the financial market and result (14) imply that

[TABLE]

and therefore

[TABLE]

We introduce the $T$-forward measure $Q^{T}$ defined as

[TABLE]

Denoting the expectation with respect to $Q^{T}$ by $E^{T}$, we obtain

[TABLE]

Observe that

[TABLE]

consequently

[TABLE]

with obvious definitions in the last line. Thus, the term $A_{1}$ captures the cost faced by the insurance company of no early surrender; the term $A_{2}$, instead, is the cost of the option embedded in the GMAB, conditional on no early surrender.

Let

[TABLE]

Further, define $\Delta t_{l}:=t_{l}-t_{l-1}$, as well as $R:=(0,\ldots,0,r)\in\mathbb{R}^{K}$, with $1<r<2$, and let for all $0\leq s\leq T$, and $u\in\mathbb{R}^{K-1}$, $v\in\mathbb{C}^{K}$

[TABLE]

[TABLE]

[TABLE]

[TABLE]

The value of the GMAB is then given in the following.

Theorem 4.1**.**

The price ${\rm P}^{\rm GMAB}$ is given by

[TABLE]

with

[TABLE]

with $A(s,T)$ as in (5).

Proof.

We start by calculating $A_{1}$ and $A_{2}$, then at the end of the proof we present an explicit formula for $Q(\tau^{m}(x)>T)$. By the definition of $\lambda^{s}$ in (12), we have

[TABLE]

where f(x_{1},\ldots,x_{K-1}):=\prod_{l=2}^{K}\big{(}e^{-\beta\Delta t_{l}x_{l-1}^{2}}\big{)}=\prod_{l=2}^{K}f_{l}(x_{l-1}). For a generic function $f$ we denote by $\hat{f}$ its Fourier transform. Then, for any $y\in\mathbb{C}$,

[TABLE]

This implies that

[TABLE]

and we observe that $\hat{f}\in L^{1}(\mathbb{R}^{K-1})$. By Theorem 3.2 in Eberlein et al. (2010)

[TABLE]

where for any $u=(u_{1},\ldots,u_{K-1})$, $\tilde{M}({\rm i}u)$ is defined as follows

[TABLE]

Using representation (47) of $D$ given in the appendix together with equations (3) and (11), we obtain

[TABLE]

Moreover, in virtue of representation (49), the density of $Q^{T}$ given in (19) can be written as

[TABLE]

Consequently

[TABLE]

with $E$ and $F$ as in (4.1). The last equality follows from the independence of $L^{1}$ and $L^{2}$ and equation (2). Therefore, with $D(u,T)$ defined in (4.1) we have

[TABLE]

and the representation for $A_{1}$ follows from (27).

We continue by calculating

[TABLE]

With the definitions $S_{T}=\exp(Y_{T})$, $G(T)=I\exp(\delta T)$, and (47), we obtain that

[TABLE]

Therefore, by the same arguments that have been used in the calculation of $A_{1}$, we can prove that

[TABLE]

for

[TABLE]

To ensure integrability, we define $g(x_{1},\ldots,x_{K}):=F(x_{1},\ldots,x_{K})e^{-rx_{K}}$, with $1<r<2$. Further, let

[TABLE]

Then, $g_{K}\in L^{1}(\mathbb{R})$, such that $g\in L^{1}(\mathbb{R}^{K})$. Moreover, elementary integration shows that for all $y\in\mathbb{R}$

[TABLE]

Observe that $|\hat{g}_{K}(y)|_{\mathbb{C}}=(((1-r)^{2}+y^{2})(r^{2}+y^{2}))^{-1/2}$, thus, $\hat{g}_{K}\in L^{1}(\mathbb{R})$. Therefore, combining the last result with (26), we deduce that $\hat{g}\in L^{1}(\mathbb{R}^{K})$, and

[TABLE]

Since $g,\hat{g}\in L^{1}(\mathbb{R}^{K})$, we can apply Theorem 3.2 in Eberlein et al. (2010) and obtain

[TABLE]

with $R:=(0,\ldots,0,r)\in\mathbb{R}^{K}$, $1<r<2$, and $\tilde{N}(R+{\rm i}u)$ defined as

[TABLE]

As above, using the notation from (4.1), we derive

[TABLE]

Observe that as $1<r<2$, and in virtue of (4) and (7), $r\sigma_{2}(s)\leq M_{2}$ and $|r\beta(s)+(1-r)\Sigma(s,T)|\leq(2r-1)\frac{M_{1}}{3}\leq M_{1}$. Thus, the above expectation exists. Using independence of $L^{1}$ and $L^{2}$ and (2), we obtain that

[TABLE]

On the other hand, observe that for any $u\in\mathbb{R}^{K}$,

[TABLE]

Thus, we deduce that

[TABLE]

Plugging (30) and (31) into (29), the claim follows.

It remains to compute $Q(\tau^{m}(x)>T)$. Given the setup in section 3.1, Itô’s lemma and Fubini’s theorem (see Escobar et al., 2016, as well for an alternative argument) imply that for any $0\leq t\leq T$

[TABLE]

with

[TABLE]

and

[TABLE]

with $\kappa,b,z,\lambda$ and $\sigma$ as in section 3.1. ∎

4.2. Death benefit

In this section we compute the value of the death benefit previously defined as

[TABLE]

Recall the definition of $w_{l}$ from Equations (22) and (23). Further, for $i\in\{1,\ldots,N\}$, we define

[TABLE]

We recall that $\Delta t_{l}:=t_{l}-t_{l-1}$. For all $j\in\{1,\ldots,K-1\}$, let $R:=(0,\ldots,0,r)\in\mathbb{R}^{j+1}$ with $1<r<2$, $u\in\mathbb{R}^{j}$, and $v\in\mathbb{C}^{j+1}$. For all $0\leq s\leq T$, for all $j\in\{1,\ldots,K-2\}$ and all $i\in\{1,\ldots,N\}$ such that $t_{j}<\bar{t}_{i}\leq t_{j+1}$, and for $j=K-1$ and all $i$ such that $t_{K-1}<\bar{t}_{i}\leq T$, we define

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Finally, we define for $0\leq s\leq T$, $u\in\mathbb{R}$ and $i\in\{1,\ldots,N\}$

[TABLE]

The following result for the DB value holds.

Theorem 4.2**.**

The price ${\rm P}^{\rm DB}$ is given by

[TABLE]

where, using the notation from (34), for $j\in\{1,\ldots,K-1\}$

[TABLE]

The proof of Theorem 4.2 is similar to the proof of Theorem 4.1. For this reasons we defer it to the Appendix C.

Remark 4.1**.**

We have

[TABLE]

As shown in the previous theorem this last expectation can be computed directly. However, this expression can also be traced back to ${\rm P}^{\rm GMAB}$. Note, that by definition,

[TABLE]

To facilitate the computation of ${\rm P}^{\rm GMAB}(\bar{t}_{i})$ we suggest the following approximation: first, for $j$ such that $t_{j}<\bar{t}_{i}\leq t_{j+1}$, we obtain that

[TABLE]

Now, we define $\hat{D}(\bar{t}_{i})=Y(\bar{t}_{i})-\delta\bar{t}_{i}$, then

[TABLE]

The suggested approximation is

[TABLE]

which can be evaluated using the results from Theorem 4.1, replacing $T$ by $\bar{t}_{i}$ and $K$ by the $j+1$ for which $t_{j}<\bar{t}_{i}\leq t_{j+1}$.

4.3. Surrender benefit

Finally, we compute the value of the surrender benefit

[TABLE]

We use the spot measure, i.e. the measure with the stock price chosen as numeraire. This allows to exploit the dependence between the surrender intensity and the stock price.

Recall the definition of $w_{l}$ from Equations (22) and (23). Further, for $0\leq s\leq T$, $i\in\{2,\ldots,K-1\}$, $u\in\mathbb{R}^{i-1}$, and $v\in\mathbb{R}^{i}$, we define

[TABLE]

[TABLE]

[TABLE]

[TABLE]

Theorem 4.3**.**

The price ${\rm P}^{\rm SB}$ is given by

[TABLE]

where $B_{1}^{1}=1$ and for $i\in\{2,\ldots,K-1\}$,

[TABLE]

and, for $i\in\{1,\ldots,K-1\},$

[TABLE]

Proof.

By construction {\rm P}^{\rm SB}=\sum_{i=1}^{K-1}E_{Q}\Big{[}e^{-\int_{0}^{t_{i}}r(u)du}{\rm SB}(t_{i})\Big{]}. As

[TABLE]

we are interested in computing the expression

[TABLE]

It follows from (13) that, for $A\in\mathscr{F}_{t_{i}}^{L^{1},L^{2}}\subset\mathscr{F}_{t_{i+1}}^{L^{1},L^{2}}$,

[TABLE]

Consequently, equation (36) can be written as

[TABLE]

$Q(\tau^{m}(x)>t_{i})$ is as in equation (32). Let us write $B^{1}_{i}$ for the first expectation and $B^{2}_{i}$ for the second one. In order to compute these two expectations we introduce the spot probability measure $Q^{S,i},\ i=1,\dots,K-1$ defined by its Radon-Nikodym derivative

[TABLE]

The above defines indeed a density process as the discounted stock price \big{(}e^{-\int_{0}^{t}r_{u}du}S(t)\big{)}_{t} is a $Q$-martingale. We can use the new measure to simplify $B^{1}_{i}$ and $B^{2}_{i}$ as follows

[TABLE]

Note that $\lambda^{s}(u)$ for $t_{i}\leq u<t_{i+1}$ is defined by $D(t_{i})$ and consequently it is $\mathscr{F}_{t_{i}}^{L^{1},L^{2}}$-measurable. Focusing first on $B^{1}_{i}$, by construction $\lambda^{s}(u)=0$ for $u\in[0,t_{1})$, and therefore $B_{1}^{1}=1$. For $i\in\{2,\ldots,K-1\}$, we follow the same strategy as in the proof of Theorem 4.1 for the computation of $A_{1}$ and consider

[TABLE]

where $f(x_{1},\ldots,x_{i-1}):=\prod_{l=2}^{i}e^{-\beta\Delta t_{l}x_{l-1}^{2}}$ with $\Delta t_{l}=t_{l}-t_{l-1}$. As in (27) we obtain that the last expectation is

[TABLE]

with

[TABLE]

and $\tilde{M}^{i-1}({\rm i}u)$ defined as

[TABLE]

From (47) it follows that

[TABLE]

By equation (2) it follows that

[TABLE]

Therefore, given the definition of $D^{i}(u,T)$ in (35) we have

[TABLE]

and the representation of $B^{1}_{i}$ follows from (41).

Finally, we note that

[TABLE]

Repeating mutatis mutandis the above arguments, the expression for $B^{2}_{i}$ follows. ∎

5. Numerical implementation and testing

This section discusses the numerical computation of the pricing equations in Theorems 4.1, 4.2 and 4.3, based on the model features introduced in section 2. For the purpose of the numerical analysis, we choose as a relevant Lévy process the Normal Inverse Gaussian (NIG) process introduced by Barndorff-Nielsen (1997) with cumulant function

[TABLE]

for $\mu\in\mathbb{R}$, $\delta>0$, $0\leq|\beta|<\alpha$. The parameter $\alpha$ controls the steepness of the density (and therefore its tail behaviour), $\beta$ primarily controls the (sign of the) skewness of the distribution, whilst $\delta$ is the scale parameter; the location parameter $\mu$ is instead set to zero, without loss of generality.

Further, we assume a simplified Vasiček structure for the function $\sigma_{1}(s,T)$ so that for $a>0$

[TABLE]

and

[TABLE]

For the equity part, we assume $\sigma_{2}(s)=\sigma_{2}>0$ and $\beta(s)=b\in\mathbb{R}$.

All numerical experiments refer to contracts with a parameter setting as in Table 1; the parameters of the financial model are taken from the calibration exercise of Eberlein and Rudmann (2018). The short maturity contracts are used for benchmarking purposes, whilst the 10 year maturity contract represents a realistic specification for practical purposes. Concerning the possible termination dates, for illustration purposes we use a 1 year frequency, i.e. $\Delta t_{l}=1$, whilst the time grid for mortality, which is assumed to have finer granularity, is built with a frequency of 6 months. We note though that the results obtained in Theorems 4.1, 4.2 and 4.3 hold for any choice of the relevant time steps. The numerical schemes are implemented in Matlab R2018a and run on a computer with an Intel i5-6500, 3.20 gigahertz CPU, and 8 gigabytes of RAM.

5.1. Implementation

The multidimensional integrals in Theorems 4.1, 4.2 and 4.3 are computed by means of Monte Carlo integration (see for example Pharr and Humphreys, 2010, and references therein). Thus, the Monte Carlo estimate of a high-dimensional integral over the domain $\Omega$

[TABLE]

is obtained by evaluating the function $f(x)$ at $M$ points $x$, drawn randomly in $\Omega$ with a given probability density $p(x)$, so that

[TABLE]

The error is measured by means of the (unbiased) sample variance

[TABLE]

as in standard Monte Carlo simulation; this point shows that the rate of convergence is independent of the dimensionality of the integrand, which is the advantage of Monte Carlo integration compared to standard numerical quadrature techniques in computing high-dimensional integrals.

Classical Monte Carlo integration uses the uniform probability density. However, visual inspection of the integrand functions $M(u,T)$ and $N(u,T)$ of Theorem 4.1 reveals that $M(u,T)$ is strongly peaked around the origin, whilst $N(u,T)$ is characterised by alternations of peaks and troughs of significant magnitude - as shown in Figure 1 for $T=3$ ($K=2$). Similar considerations hold for the integrand functions in Theorems 4.2 and 4.3, due to the similarities between the corresponding payoff functions. These features, in general, can induce large variances in the Monte Carlo estimate of the corresponding integral.

Consequently, for variance reduction purposes and in order to speed up convergence, we implement importance sampling; in details, given that the integrand functions are strongly peaked around the origin, the importance distribution of choice is the multivariate Gaussian distribution with zero mean (i.e. centered around the peak), independent components, and a given variance matrix, which is treated as a parameter. Full deployment of the pricing algorithm in the setting of VEGAS and MISER Monte Carlo integration is left to future research.

5.2. Benchmarking and Testing

Numerical results for the pricing functions considered in this paper are reported in Table 2 and 3. In order to provide a reliable benchmark for the Monte Carlo integration procedure introduced above, we consider first some simple examples for which the relevant multidimensional integrals can be tackled with standard quadrature packages. We illustrate the results obtained for the case of the GMAB. Similar performances are obtained for the DB and the SB as well, and we refer the interested reader to the Appendix D for more details.

In details. Table 2 reports values obtained with deterministic quadrature methods, the corresponding estimate from Monte Carlo integration, and the CPU time. Together with the Monte Carlo estimate, $I_{MC}$, we also report measures of its accuracy in terms of the absolute value of the bias of the estimator expressed as percentage of the value obtained by quadrature $I_{Q}$, i.e.

[TABLE]

In addition, we also report the percentage standard error of the Monte Carlo estimate

[TABLE]

For increased reliability of the estimate of the CPU time, we base the numerical experiment on 100 repetitions of the Monte Carlo algorithm with sample size $10^{6}$.

We consider two examples for a contract with maturity $T=3$ years (i.e. $K=2$) and $T=4$ years (i.e. $K=3$) respectively. With annually spaced termination dates $t_{i}$, these choices imply that $A_{1}$ and $A_{2}$ are respectively one- and two-dimensional integrals in the first case, and two- and three-dimensional integrals in the second case, which can also be computed using packages for deterministic quadrature methods.

For importance sampling, numerical experiments show that relatively small biases and standard errors can be obtained by using the same variance fixed at 0.25 for the first $K-1$ dimensions, and increasing this value to 1 for the last $K^{th}$ dimension, as to cater for the higher variability of the integrand function $N(u,T)$. The results in Table 2 confirm the quality of the estimates as all biases and standard errors are below 0.5%.

The value of the VA together with its components are reported in Table 3. We observe that the GMAB value decreases with increasing maturities, due to the larger number of dates at which the contract can be terminated by either death or surrender of the policyholder. Consistently with these findings, we also observe the increase in the value of both the DB and the SB for contracts with longer maturities. The overall value of the VA increases as well with longer maturities, denoting the dominant impact of the DB and the SB components.

Finally, we observe the higher computational cost of the DB and SB parts of the full variable annuity; this is due to the large number of integrals of the cumulant functions of interest (see definitions (34) and (35)) be computed. This number depends on the granularity of both the $\mathbf{t}$ and $\bar{\mathbf{t}}$ grids: the finer these grids, the more computational demanding these components of the contract become.

6. Results: Sensitivity Analysis

The sensitivity analysis is carried out by perturbing the parameters of interest one at a time, ceteris paribus. The benchmark case is given by a 10 years contract with the parameters set as in Table 1.

In Figure 2 we show the sensitivity of the VA and its components to the parameters $\beta$ and $C$, controlling the surrender intensity. Consistently with intuition, the higher the value of $\beta$, the stronger the influence of the financial market on the policyholder decision to surrender. This is captured in particular by the value of the surrender benefit, SB, which increases with $\beta$. Due to the higher surrender probability, the values of the GMAB and the DB reduce, and so does the resulting value of the VA.

The impact of the parameter $\beta$ though becomes less relevant in correspondence of higher values of the constant $C$ which represents the baseline surrender behaviour. In this case, the non-economic factors leading to the decision of surrender dominate the influence of the financial market to the point that the VA as well as its components are almost insensitive to $\beta$.

These results highlight that a correct quantification of the baseline surrender parameter $C$ is of paramount importance given the significant impact on the value of the VA and its components. In this respect, it is crucial for insurance companies to correctly calibrate the surrender model to the information that they could collect regarding the policyholder behaviour towards lapses.

Figure 3 shows that the impact of the parameter $\beta$ is also affected by the dependence between the equity and the fixed income markets, here captured by the parameter $b$. Indeed, in presence of explicit dependence between the financial markets (i.e. $b\neq 0$), the behaviour of the policyholders becomes more responsive to changes in the market conditions. This is reflected in the rate of change of the contract’s value.

Finally, in Figure 4 we illustrate the behaviour of the contract with respect to the guaranteed rate $\delta$. A higher value of this rate corresponds to a more valuable VA, as one would expect; the SB though is relatively insensitive to $\delta$, also due to the low values of the latter used in this analysis.

7. Conclusion

We proposed a general framework for the valuation of a number of guarantees commonly included in variable annuity products, such as guaranteed minimum accumulation benefit, death benefit and surrender benefit, in the setting of a hybrid model based on multivariate Lévy processes, with surrender risk captured endogenously.

This framework proves to be tractable, and allows for the deployment of efficient numerical schemes based on Monte Carlo integration. We used this setting to gain insights into the contract sensitivity with respect to the model parameters, with special focus on surrender risk. The emphasis on surrender in particular finds its motivation in the large percentage of early terminations experienced by the insurance companies issuing variable annuities, and the resulting net cash outflow they sustain as a consequence (as also pointed out by LIMRA Secure Retirement Institute and IRI). The results obtained in our analysis confirm the importance of the surrender behaviour and its appropriate modelling.

We envisage the applicability of the model setup offered in this paper, and the proposed numerical scheme, also for other benefits offered by VAs, such as the post-retirement guaranteed minimum income benefit. This would represent a promising route for further research.

Acknowledgements

The authors gratefully acknowledge the financial support from the Freiburg Institute for Advanced Studies (FRIAS) within the FRIAS Research Project ‘Linking Finance and Insurance: Theory and Applications (2017–2019)’.

Appendix A An alternative specification for the surrender intensity

The surrender intensity $\lambda^{s}$ is defined in (12) by

[TABLE]

where $W_{1}(x):=\beta x^{2}+C$. Consequently $W_{1}(x)\to\infty$ as $|x|\to\infty$. One might be interested in the behaviour of the contract price for a surrender intensity whose values do not become arbitrarily large. Therefore, in the following we study the valuation of the variable annuity’s components under an alternative specification which keeps the surrender intensity bounded.

We start with an observation concerning the behaviour of $D(t_{i})$, $i=1,\dots,K-1$. The following lemma shows that these quantities stay in a bounded interval with high probability. Based on this fact, Proposition A.2 and its corollary show that the contract values obtained with the alternative surrender intensity approximate the values under the original assumption.

Lemma A.1**.**

For any $\epsilon>0$, there exists $L>0$ such that

[TABLE]

where the last result holds for $j\in\{1,\ldots,K-2\}$ and any $i$ such that $t_{j}<\bar{t}_{i}\leq t_{j+1}$, and for $j=K-1$ and any $i$ such that $t_{K-1}<\bar{t}_{i}\leq T$. $Q^{\bar{t}_{i}}$ is defined in (51), $Q^{S,j}$ is defined in (37), and $Q^{T}$ is defined in (19).

Proof.

Let $L_{1}$ denote a positive constant. Applying the Markov inequality and the inequality $e^{|x|}\leq e^{x}+e^{-x}$, we get

[TABLE]

Now we use the representation given in (47), the definition of $Q^{T}$ in (19) and (49), to derive

[TABLE]

Observe that, by independence of $L^{1}$ and $L^{2}$, we have

[TABLE]

where the last quantities are finite due to (4) and (7). Therefore, $C_{i}:=E^{T}(e^{D(t_{i})})$ is finite, and by similar arguments $C_{i}^{\prime}:=E^{T}(e^{-D(t_{i})})$ is finite as well. Observe that $e^{-L_{1}}\sum_{i=1}^{K-1}(C_{i}+C_{i}^{\prime})\leq\epsilon$ is equivalent to $L_{1}\geq-\log(\epsilon)+\log(\sum_{i=1}^{K-1}(C_{i}+C_{i}^{\prime}))$. As a consequence of (45) and the last argument, we deduce that for $L_{1}\geq-\log(\epsilon)+\log(\sum_{i=1}^{K-1}(C_{i}+C_{i}^{\prime}))$, $Q^{T}(\bigcup_{i=1}^{K-1}\{|D(t_{i})|>L_{1}\})\leq\epsilon$.

By similar arguments, we can prove that for any $j\in\{1,\ldots,K-1\}$, there exists $\hat{L}_{j}>0$, such that $Q^{S,j}(\bigcup_{i=1}^{j}\{|D(t_{i})|>\hat{L}_{j}\})\leq\epsilon$. We can prove also that for $j\in\{1,\ldots,K-2\}$ and any $i$ such that $t_{j}<\bar{t}_{i}\leq t_{j+1}$, and for $j=K-1$ and any $i$ such that $t_{K-1}<\bar{t}_{i}\leq T$, there exists $L_{i,j}$ such that $Q^{\bar{t}_{i}}(\bigcup_{l=1}^{j}\{|D(t_{l})|>L_{i,j}\})\leq\epsilon$. Therefore, by defining $\bar{L}:=\max_{(i,j)}L_{i,j}$ and $L:=\max(L_{1},\hat{L}_{1},\ldots,\hat{L}_{K-1},\bar{L})$, we deduce that (44) holds true. ∎

Remark A.1**.**

The constant $L$ in Lemma A.1 depends on $\epsilon$ and should be denoted $L_{\epsilon}$ but we omit this notation for simplicity.

Based on the above, let us replace $W_{1}$ by another positive and continuous function $W_{2}$ which coincides with $W_{1}$ on $[-L,L]$, $L$ as in Lemma A.1, is constant outside this compact set on the positive half line and converges to 0 as $x\to-\infty$. Specifically, we define $W_{2}$ as

[TABLE]

with the constants $\beta$ and $C$ as in (12). The new surrender intensity $\bar{\lambda}^{s}$ is defined as

[TABLE]

for $i\in\{1,\ldots K-1\}$, and $\bar{\lambda}^{s}(t)=0$ for $t\in[0,t_{1})\cup[t_{K},T]$. Then, as we will see in Corollary A.3, the corresponding prices, denoted by ${\rm\bar{P}}^{\rm GMAB},\bar{{\rm P}}^{\rm DB}$ and $\bar{{\rm P}}^{\rm SB}$, approximate ${\rm P}^{\rm GMAB},{\rm P}^{\rm DB}$ and ${\rm P}^{\rm SB}$.

Starting with the value of the GMAB, we rewrite equation (21) as

[TABLE]

with

[TABLE]

Following the same calculations as in Theorem 4.2, we deduce that

[TABLE]

where, $A_{0,i}$ is the same as in Theorem 4.2, and for $j\in\{1,\ldots,K-2\}$ and any $i$ such that $t_{j}<\bar{t}_{i}\leq t_{j+1}$, and for $j=K-1$ and any $i$ such that $t_{K-1}<\bar{t}_{i}\leq T$

[TABLE]

Finally, by the same argument as in Theorem 4.3, we deduce that

[TABLE]

with $\bar{B}_{1}^{1}=1$, due to $\bar{\lambda}^{s}(u)=0$, for $u\in[0,t_{1})$ by construction, and

[TABLE]

Given the above, the following result holds.

Proposition A.2**.**

We have

• (1)

$|\bar{A}_{1}-A_{1}|\leq 2\epsilon$,

• (2)

$|\bar{A}_{2}-A_{2}|\leq 2C_{2}\epsilon^{1/2}$, with C_{2}:=E^{T}\big{[}(e^{D(T)}-1)^{2}\big{]}^{1/2},

• (3)

$|\bar{A}^{1}_{j,i}-A^{1}_{j,i}|\leq 2\epsilon$,

• (4)

$|\bar{A}^{2}_{j,i}-A^{2}_{j,i}|\leq 2C_{2,i}\epsilon^{1/2}$, with C_{2,i}:=E^{\bar{t}_{i}}\big{[}(e^{D_{\bar{t}_{i},\bar{t}_{i}}+p(\bar{t}_{i})}-1)^{2}\big{]}^{1/2},

• where in (3) and (4) we have either $j\in\{1,\ldots,K-2\}$ and $i$ such that $t_{j}<\bar{t}_{i}\leq t_{j+1}$, or $j=K-1$ and $i$ such that $t_{K-1}<\bar{t}_{i}\leq T$,
• (5)

$|\bar{B}_{j}^{1}-B_{j}^{1}|\leq 2\epsilon$, for $j\in\{2,\ldots,K-1\}$,

• (6)

$|\bar{B}_{j}^{2}-B_{j}^{2}|\leq 2\epsilon$, for $j\in\{1,\ldots,K-1\}$.

$A_{1}$, $A_{2}$ are given in (21), $A^{1}_{j,i}$, $A^{2}_{j,i}$ in (52), and $B_{j}^{1}$, $B_{j}^{2}$ are given in (38) and (39).

Proof.

It suffices to prove (2) which is representative for the degree of sophistication. The proofs of (1), (3), (4), (5) and (6) follow mutatis mutandis. From the definition of $W_{1}$ given in (43), we derive

[TABLE]

with the obvious definition of $a_{1}$ in the last line. The third equality holds as $W_{1}$ and $W_{2}$ coincide on $[-L,L]$. Observe that

[TABLE]

with the obvious notation in the last line. We deduce that

[TABLE]

As $W_{1}$ is a positive function, the following holds

[TABLE]

where the third line is a consequence of the Cauchy-Schwarz inequality. For the last inequality we use (44) together with the fact that C_{2}:=E^{T}\bigg{[}(e^{D(T)}-1)^{2}\bigg{]}^{1/2}<\infty. This last fact is easy to verify from the definition of $D(T)$ in (47). A similar bound can be achieved for $a_{1}$, therefore

[TABLE]

which proves (2). ∎

Based on the last proposition, the interested reader can prove the following corollary.

Corollary A.3**.**

We have

• (1)

$|{\rm\bar{P}}^{\rm GMAB}-{\rm P}^{\rm GMAB}|\leq 2Q(\tau^{m}(x)>T)B(0,T)G(T)[\epsilon+C_{2}\epsilon^{1/2}],$* *

*with *C_{2}:=E^{T}\big{[}(e^{D(T)}-1)^{2}\big{]}^{1/2}.

• (2)

|{\rm\bar{P}}^{\rm DB}-{\rm P}^{\rm DB}|\leq 2\sum_{i:\>\bar{t}_{i}>t_{1}}Q\big{(}\tau^{m}(x)\in[\bar{t}_{i-1},\bar{t}_{i})\big{)}B(0,\bar{t}_{i})G(\bar{t}_{i})[\epsilon+C_{2,i}\epsilon^{1/2}],

with C_{2,i}:=E^{\bar{t}_{i}}\big{[}(e^{D_{\bar{t}_{i},\bar{t}_{i}}+p(\bar{t}_{i})}-1)^{2}\big{]}^{1/2}.

• (3)

$|{\rm\bar{P}}^{\rm SB}-{\rm P}^{\rm SB}|\leq 2IP(t_{1})Q(\tau^{m}(x)>t_{1})\epsilon+4I\epsilon\sum_{j=2}^{K-1}P(t_{j})Q(\tau^{m}(x)>t_{j}).$**

Appendix B Some useful results and representations

Let us recall that $A(u,T)=\int_{u}^{T}\alpha(u,s)ds$ and $\Sigma(u,T)=\int_{u}^{T}\sigma_{1}(u,s)ds$. We derive here a representation result in the Lévy forward rate framework.

Lemma B.1**.**

For any $0\leq t\leq T$, we have that

[TABLE]

where $f(t,s)$ is defined in (3).

Proof.

Using Fubini’s theorem for stochastic integrals, we deduce that

[TABLE]

and the claim follows. ∎

The above lemma allows the following representation of $D(t)$ defined in (11).

[TABLE]

In addition, setting $t=T$ in the above lemma, we obtain

[TABLE]

and hence, the bank account has the following well-known representation, see for example Eberlein and Raible (1999) or (13) in Eberlein and Rudmann (2018),

[TABLE]

Appendix C Proof of Theorem 4.2

Proof.

By definition,

[TABLE]

From equation (10), we obtain that

[TABLE]

We distinguish two cases: when $\bar{t}_{i}\leq t_{1}$ and when $t_{1}<\bar{t}_{i}$. We start with the detailed description of the second case which is the more involving one. The first case is treated at the end.

For $j\in\{1,\ldots,K-2\}$ and $i$ such that $t_{j}<\bar{t}_{i}\leq t_{j+1}$, as well as for $j=K-1$ and $i$ such that $t_{K-1}<\bar{t}_{i}\leq T$, we work along the same line as in the proof of Theorem 4.1, and get

[TABLE]

We introduce the $\bar{t}_{i}$-forward measure $Q^{\bar{t}_{i}}$ defined by its Radon-Nikodym density

[TABLE]

Denoting the expectation with respect to $Q^{\bar{t}_{i}}$ by $E^{\bar{t}_{i}}$ the quantity in equation (50) is

[TABLE]

with an obvious notation in the last line. We note that

[TABLE]

with $f(x_{1},\ldots,x_{j})=\prod_{l=2}^{j+1}e^{-\beta\Delta t_{l}x_{l-1}^{2}}$, and $\Delta t_{l}=t_{l}-t_{l-1}.$ As in (27) we obtain that the expectation in (53) can be represented as

[TABLE]

with

[TABLE]

and $\tilde{M}^{j}_{i}({\rm i}u)$ defined as

[TABLE]

By the representation of the bank account in (49), the Radon-Nikodym density (51) can be written as

[TABLE]

Consequently, from the representation of $D(t)$ in (47), it follows

[TABLE]

This last expectation is finite in virtue of (4), and returns

[TABLE]

in virtue of the definitions in (34) and equation (2). Therefore, with $D^{j,i}(u,T)$ defined in (34) we have

[TABLE]

Finally, combining (53) with (54) and the definition of $M^{j,i}(u,T)$ in (34), we deduce that

[TABLE]

For the purpose of the computation of $A_{j,i}^{2}$, we replace $D(t)$ defined in (11) with the following quantity

[TABLE]

so that

[TABLE]

We note that

[TABLE]

for $h(x_{1},\ldots,x_{j+1}):=f(x_{1},\dots,x_{j})(e^{x_{j+1}+p(\bar{t}_{i})}-1)^{+},$ with $f$ given in (53). In order to ensure integrability, let us define $H(x_{1},\ldots,x_{j+1}):=h(x_{1},\ldots,x_{j+1})e^{-rx_{j+1}}$, for some $1<r<2$, and

[TABLE]

Then, $H_{j+1}\in L^{1}(\mathbb{R})$, and $H\in L^{1}(\mathbb{R}^{j+1})$. Moreover, elementary integration shows that for all $y\in\mathbb{R}$

[TABLE]

Observe that $|\hat{H}_{j+1}(y)|_{\mathbb{C}}=e^{rp(\bar{t}_{i})}(((1-r)^{2}+y^{2})(r^{2}+y^{2}))^{-1/2}$, thus, $\hat{H}_{j+1}\in L^{1}(\mathbb{R})$. Therefore, combining the last result with (26), we deduce that $\hat{H}\in L^{1}(\mathbb{R}^{j+1})$, and

[TABLE]

As $H,\hat{H}\in L^{1}(\mathbb{R}^{j+1})$, it follows from Theorem 3.2 in Eberlein et al. (2010) that

[TABLE]

for $R=(0,\ldots,0,r)\in\mathbb{R}^{j+1}$, $1<r<2$, and $\tilde{N}_{i}^{j+1}(R+{\rm i}u)$ defined as

[TABLE]

Using (56) and (48), we get

[TABLE]

Plugging the last quantity in (60), in virtue of (47) and (55), we deduce that

[TABLE]

Using the definitions from (34), the above can be rewritten as

[TABLE]

Observe that, due to $1<r<2$, as well as (4) and (7), $r\sigma_{2}(s)\leq M_{2}$ and $|r\beta(s)+(1-r)\Sigma(s,T)|\leq(2r-1)\frac{M_{1}}{3}\leq M_{1}$. Thus, the above expectation exists. Using the independence of $L^{1}$ and $L^{2}$ and (2), we obtain that

[TABLE]

On the other hand, we observe that for any $u\in\mathbb{R}^{j+1}$,

[TABLE]

Consequently, we deduce that

[TABLE]

Plugging (61) and (62) in (59), and using the definition of $N^{j,i}(u,T)$ in (34), it follows that

[TABLE]

We now consider the case $\bar{t}_{i}\leq t_{1}$ for $i\in\{1,\ldots,N\}$. Using the same arguments as above, we derive

[TABLE]

where $h_{1}(x):=(e^{x+p(\bar{t}_{i})}-1)^{+}$. For some $1<r<2$, we define the function $H_{1}$ as $H_{1}(x):=(e^{x+p(\bar{t}_{i})}-1)^{+}e^{-rx}.$ By Theorem 3.2 in Eberlein et al. (2010), we get

[TABLE]

with

[TABLE]

Using the same arguments as above, and the definitions in (34), we deduce that

[TABLE]

On the other hand, we have

[TABLE]

Plugging (64) and (65) in (63), we get

[TABLE]

Therefore,

[TABLE]

Finally, we mention that for any $1\leq i\leq N$, we can compute $Q(\tau^{m}(x)>\bar{t}_{i})$ (and consequently $Q(\tau^{m}(x)\in[\bar{t}_{i-1},\bar{t}_{i}))$) in virtue of (32).

∎

Appendix D Further numerical results: benchmarking

In this Appendix, we provide full results from the benchmarking exercise of the Monte Carlo integration pricing procedure.

From section 5, we recall that we consider contracts with short maturity, and surrender frequency $\Delta t_{l}=1$ year. In addition, we assume half-annually spaced $\bar{t}_{i}$ for the mortality monitoring.

Consequently, a sensible benchmarking exercise for the value of both the DB and the SB with deterministic quadrature methods can be achieved by considering the 4 year maturity contract, i.e. $K=3$.

Starting with the value of the DB, the term $A_{0,i}$ in Theorem 4.2 is a 1-dimensional integral for $i=1,\dots,N$, which can be obtained by direct quadrature, and therefore is not considered by this benchmarking exercise.

Monte Carlo integration is instead used to compute the remaining terms appearing in Theorem 4.2, i.e. $A_{1,i}^{1}$, $A_{1,i}^{2}$, $A_{2,i}^{1}$, which are two-dimensional integrals, and $A_{2,i}^{2}$ which is a three-dimensional integral, for all $i=1,\dots,N$. These are benchmarked against the values obtained by deterministic quadrature procedures in Matlab.

For importance sampling, we choose the same values of the variance of the importance sampling distribution as in the GMAB case, i.e. 0.25 for the first $K-1$ dimensions and 1 for the final $K^{th}$ dimension. The quality of the estimate is confirmed by the negligible bias and standard errors reported in Table 4.

Concerning the value of the SB, we notice that the first term in the sum defining $P^{SB}$ in Theorem 4.3 is composed by a constant ($B_{1}^{1}$) and a one-dimensional integral ($B_{1}^{2}$), which is obtained by deterministic quadrature; similarly to the previous case this term is not considered in this analysis. The second term in this sum is instead formed by a one- and a two-dimensional integral ($B_{2}^{1}$ and $B_{2}^{2}$ respectively) for which we deploy Monte Carlo integration. Benchmarking is performed against the corresponding quadrature routines in Matlab.

For importance sampling, given the relatively simple forms of the integrand functions, we use the same variance fixed at 0.16 across all the $K$ dimensions. The goodness of the estimate is confirmed by the negligible bias and standard errors shown in Table 4.