Impact Mortality Heterogeneity in Pricing Life Annuities in the Dutch population

(1)

Impact Mortality Heterogeneity in Pricing Life

Annuities in the Dutch population

J. Achterhof (s1763490)

May 29, 2015

Abstract

(2)

1 Introduction

The Dutch population is aging very rapidly, see for e.g. Carrey (2002). Car-rey (2002) states that the number of people of age 65 and older will double between 2010 and 2030. Due to the changing age structure in the Nether-lands the demand for annuity products, pensions, will increase. In order to price these annuity products we need to model and forecast the mortality rates of the Netherlands. The aim of this paper is to investigate how to price the annuity products and how to model the mortality rates in the Nether-lands.

From studies, e.g. Keyfitz (1978), it is well-known that individuals differ in their mortality. The difference in mortality, mortality heterogeneity, was studied by Vaupel and Stallard (1979). In the study by Vaupel and Stallard (1979) the heterogeneity of the mortality rates of the population was mod-eled by using a frailty factor. The proposed frailty factor is used to represent the heterogeneity in a population. Vaupel and Stallard (1979) chose for the assumption that the frailty factor is gamma distributed. Vaupel and Stallard (1979) chose the gamma distribution because it is analytically tractable and readily computable. Another property of the gamma distribution is that the variables cannot be negative which is a property that is needed for modeling mortality rates. Manton and Vaupel (1986) and Hougaard (1984) used as an alternative distribution for the distribution of the frailty factor, the inverse Gaussian distribution. The inverse Gaussian distribution has the desired property that the variables cannot be negative.

The heterogeneity of mortality rates has an impact on the pricing of life in-surance products and on the pricing of annuity products, pensions. In the life insurance market some of the heterogeneity of the mortality rates is reduced by underwriting (e.g. BMI, health history, smoker etc.). This phenomenon of underwriting is however not common in the life annuity market.

(3)

(4)

2 Formulation of the Mortality Models

This section reviews the three types of models we will apply to Dutch mor-tality data and compare afterwards. We will apply two frailty models and one phase-type aging model. In the next subsection the two frailty models will be discussed.

2.1 Formulation of Frailty Models

By applying frailty models to mortality data we can model mortality het-erogeneity. In a frailty model we assume that every individual has his own frailty factor which is fixed at birth. The frailty factor is seen as the vulner-ability of a person to death. We consider the following multiplicative model for the force of mortality,

µ(x, z) = µ(x, 1)z (1)

(5)

(6)

2.1.1 Gamma Frailty

In this section we will make the assumption that the frailty factor at birth, z is Gamma distributed with shape parameter δ and scale parameter θ, z ∼

Gamma(δ, θ). Due to this assumption the marginal density, fZ(z), of frailty

at the age of zero is given by,

fZ(z) = fZ|X(z|X = 0) =

θδ

Γ(δ)tz

δ−1

exp(−θz). (3)

The mean of a gamma distribution is given by δ_θ and the variance of a gamma

distribution is given by _θδ2. Thus, the mean and the variance of the gamma

distributed frailty at birth are given by,

E[z] = δ

θ

Var[z] = δ

θ2

We want the mean of the frailty at birth to be equal to one, E[z] = 1. Thus, we let δ = θ.

In order to derive the joint probability density function of frailty z and age x we will, first derive the probability density function of the random lifetime

for a newborn, T0, given frailty z by using the following definition of the force

of mortality, µ(x, z) = F 0 X|Z(x|Z = z) 1 − FX|Z(x|Z = z) = fX|Z(x|Z = z) sX|Z(x|Z = z) . (4)

By definition the conditional survival function given frailty z is given by,

(7)

where H(x) is the cumulative hazard of the standard force of mortality at age x, H(x) =R₀xµ(t, 1)dt . Hence, µ(x, z) = fX|Z(x|Z = z) sX|Z(x|Z = z) z µ(x, 1) = fX|Z(x|Z = z) exp(−z H(x)) ⇒ fX|Z(x|Z = z) = z µ(x, 1) exp(−z H(x)). (6)

fX|Z(x|Z = z) is the probability density function of the random lifetime for

a newborn, T0, given frailty, z. The joint probability density function of age

x and frailty z is the product of the conditional probability density functions

fX|Z(x|Z = z) (6) and fZ(z) (3). Thus, it is given by,

fX,Z(x, z) = fX|Z(x|Z = z) fZ(z) = zµ(x, 1) exp(−z H(x)) θ δ Γ(δ) z δ−1 _{exp(−θ z)} = zδ µ(x, 1) exp(−z (H(x) + θ)) θ δ Γ(δ). (7)

In order to obtain the marginal probability density function of frailty in the

surviving population with survival up to age x, integration of fX,Z(x, z) (7)

with respect to age from x to ∞ is executed. The integration is given by,

fX,Z(x, z|X > x) = Z ∞ x fT ,Z(t, z)dt = Z ∞ x fT |Z(t|Z = z) fZ(z)dt = fZ(z) Z ∞ x fT |Z(t|Z = z)dt = θ δ Γ(δ) z δ−1 _{exp(−θ z) exp(−z H(x))} = θ δ Γ(δ) z δ−1 _{exp(−z(H(x) + θ))} ₍₈₎

By normalizing result (8) we obtain the probability density function of the frailty in the surviving population up to age x. Thus,

fX,Z(x, z|X > x) =

(θ + H(x))δ

Γ(δ) z

(8)

The probability density function of z conditional on X > x, the frailty in the surviving population up to age x, is a Gamma(δ, θ + H(x)). The mean and the variance of the gamma distributed frailty z conditional on X > x, the frailty in the surviving population with survival up to age x are given by,

E[z|X > x] = δ

(θ + H(x))

Var[z|X > x] = δ

(θ + H(x))2

and with θ = δ we have,

E[z|X > x] = θ (θ + H(x)) (10) Var[z|X > x] = θ (θ + H(x))2 (11) The mean frailty (10) with survival up to age x, frailty z conditional on X > x, is decreasing as age increases. This is due to that the members of the cohort with the higher values of frailty have a higher probability of dying and therefore the members of the cohort with a higher frailty level die first. This selection process lowers the mean frailty with survival up to age x, as age increases. Due to the same reason the variance (11) is decreasing with age.

(9)

equal to all the achievable values of z, fX(x) = Z ∞ 0 fX,Z(x, z)dz fX(x) = Z ∞ 0 zδ µ(x, 1) exp(−z (H(x) + θ)) θ δ Γ(δ)dz = µ(x, 1) k θ δ (θ + H(x))δ+1 Z ∞ 0 zδ exp(−z (H(x) + θ)) (θ + H(x))δ+1 Γ(δ + 1) dz (13) = µ(x, 1) δ θ δ (θ + H(x))δ+1 (14)

note that (13)R₀∞zk exp(−z (H(x)+θ)) (θ+H(x))_Γ(δ+1) δ+1dz is the integrand of a Gamma

density. The obtained distribution fX(x) (14) is the unconditional

distribu-tion where age equals the age x. From the uncondidistribu-tional distribudistribu-tion of x, the age, and the joint distribution for the frailty z and the age x, we can derive the conditional distribution of frailty. The distribution of frailty z conditional on age x is the distribution of frailty of those who die at age x. The derivation is given by,

fZ|X(z|X = x) = fX,Z(x, z) fX(x) = z δ _{µ(x, 1) exp(−z (H(x) + θ))} θδ Γ(δ) µ(x,1) δ θδ (θ+H(x))δ+1 = (H(x) + θ) δ+1 Γ(δ + 1) z δ _{exp(−z(H(x) + θ))} ₍₁₅₎

where the last equation (15) is equal to the density function of that of a Gamma density of the form Gamma(δ + 1, θ + H(x)). The mean and the variance of the frailty for the deaths of people at the age of x are given by,

E[z|X = x] = δ + 1

θ + H(x)

Var[z|X = x] = δ + 1

(10)

and with θ = δ we have,

E[z|X = x] = θ + 1

(θ + H(x)) (16)

Var[z|X = x] = θ + 1

(θ + H(x))2 (17)

The mean (16) and the variance (17) of the gamma distributed frailty con-ditional on X = x are decreasing as age increases. Also note that the mean (10) of the gamma distributed frailty E[z|X > x] in the surviving population at age x is smaller than the mean (16) of the frailty,E[z|X = x] for the deaths at the age x.

The gamma distribution has the scaling property. That is, if X ∼ Gamma(δ, θ)

then for any c > 0 cX ∼ Gamma(δ, θ_c). From this scaling property of the

gamma distribution we can easily derive the distribution for the force of mor-tality at age x. The distribution of the gamma distributed frailty z in the sur-viving population at age x is Gamma(δ, H(x)+θ) and since µ(x, z) = µ(x, 1) z

(1) the distribution for the force of mortality at age x is Gamma(δ,H(x)+θ_µ(x,1) ).

And with δ = θ we have that the distribution for the force of mortality at

(11)

2.1.2 Inverse Gaussian Frailty

In this section we will make the assumption that the frailty factor at birth, z, is Inverse Gaussian distributed. We will use the probability density function of the Inverse Gaussian frailty z at birth proposed by Hougaard (1984). The probability density function is given by,

fZ(z) = fZ|X(z|X = 0) = λ π 1/2 exp(4λδ)1/2 z−3/2 exp −δ z −λ z (19) with shape parameter λ > 0 and parameter δ > 0. The mean and the variance of a Inverse Gaussian distributed frailty at birth are given by,

E[z] = λ δ 1/2 Var[z] = 1 2 λ δ3 1/2

Again we want the mean of the frailty at birth to be equal to one, E[z] = 1. Thus, we let λ = δ.

In order to derive the joint probability density function of frailty z and age x we will, first derive the probability density function of the random lifetime

for a newborn, T0, given frailty z by using the following definition of the force

of mortality, µ(x, z) = F 0 X|Z(x|Z = z) 1 − FX|Z(x|Z = z) = fX|Z(x|Z = z) sX|Z(x|Z = z) . (20)

By definition the conditional survival function with the given frailty z is given by,

(12)

where H(x) is the cumulative hazard of the standard force of mortality at

age x, H(x) =R₀xµudu. Thus,

µ(x, z) = fX|Z(x|Z = z)

sX|Z(x|Z = z)

z µ(x, 1) = fX|Z(x|Z = z)

exp(−z H(x))

⇒ fX|Z(x|Z = z) = z µ(x, 1) exp(−z H(x)). (21)

The probability density function (21) is the probability density function of the the random lifetime for a newborn given frailty z. The joint

probabil-ity densprobabil-ity function fX,Z(x, z) of age x and frailty z is the product of the

conditional probability density functions fX|Z(x|Z = z) (21) and fZ(z) (19).

Thus, it is given by,

fX,Z(x, z) =fX|Z(x|Z = z) fZ(z) =z µ(x, 1) exp(−z H(x)) λ π 1/2 exp (4λδ)1/2 z−3/2 exp −δ z −λ z . (22)

By integrating fX,Z(x, z) (22) with range from x to ∞ for the age we obtain

distribution of frailty z conditional on X > x or stated different the frailty given survival up to age x. The integration is given by,

(13)

By normalizing result (23) we obtain the probability density function of z in the surviving population up to age x,

fX,Z(x, z|X > x) = λ π 1/2 exp(4λ(δ + H(x)))1/2 z−3/2 exp −(δ + H(x))) z −λ z . (24)

The probability density function (24) of z in the surviving population at

age x is a Inverse Gaussian density. The mean and the variance of the

Inverse Gaussian distributed frailty z conditional on X > x, the frailty in the surviving population at age x are given by,

E[z|X > x] = λ (δ + H(x) 1₂ Var[z|X > x] = 1 2 s λ (δ + H(x))3

and with δ = λ the mean, variance and the coefficient of variation are given by, E[z|X > x] = λ λ + H(x) 1₂ (25) Var[z|X > x] = 1 2 s λ (λ + H(x))3 (26)

The mean (25) and the variance (26) are decreasing when the age x in-creases. The more heterogeneous the population is, the faster the decrease of the mean (25) and the variance (26) is. Similar to the gamma distribu-tion the inverse Gaussian distribudistribu-tion has the scaling property. From the scaling property of the inverse Gaussian distribution we can easily derive the distribution for the force of mortality at age x. The distribution of the inverse Gaussian distributed frailty z in the surviving population at age x is IG(λ, H(x) + δ) and since (1) the distribution for the force of mortality at

age x is IG(λ µ(x, 1),H(x)+δ_µ(x,1) ). And with δ = λ we obtain the distribution of

(14)

2.2 Phase-Type Aging Model

In this section we will introduce the phase-type aging model, which is a finite-state continuous-time Markov process, that is going to be used to model the death rate. This phase-type aging model was introduced by Lin and Liu (2007). Lin and Liu (2007) introduced a finite-state continuous Markov process in order to model the aging process, a phase-type aging model. In the phase-type aging model the different states represent a physiological age state or a so called health condition state. The process of aging in the phase-type model is the transition from one state to the next state. Death of an individual is modeled as the transition from any state to the absorbing state, the death state. If we let n be the number of phases of the model,

λi, i = 1, ..., n be the transition rates and qi, i = 1, ..., n are the transition

rates to the death state the model can be displayed by figure 1.

1

2

3 n

λ 1 λ2 λ3 λn-1

...

q₃ q₂ q₁ q_n

Figure 1: Representation phase-type aging process

In the phase-type aging model people enter the process in state 1, thus the

initial state vector is given by, α = 1 0 ... 0. The generator of this

described phase-type aging model, or Markov process, is,

Q =∆ q

>

0 0

where we have that q = q1 q2 ... qn and the transition rate matrix for

the transient states is,

_−(λ

(15)

For these model assumptions the time until the event of death follows the phase-type distribution (α, ∆). For the phase-type distribution the distri-bution function and the density function for X are given by,

f (x) = α exp (∆ x)∆0

F (x) = 1 − α exp (∆ x)ι

where ι =1 1 ... 1>, ∆0 = −∆ι and exp(X) is the matrix exponential.

The matrix exponential is given by the power series exp(X) =P∞

k 1 k!X

k_{. By}

definition the survival function of time to the event of death X = x is given by, s(x) = P ({X > x}) = Z ∞ x f (u)du = 1 − F (x) = α exp (∆ x)ι. (29)

From the definition of the survival function 29 we can obtain the death probability, which is given by,

ˆ

q(x) = s(x) − s(x + 1)

(16)

3 Model Estimation

In this section we will describe the estimation procedures applied to the different types of models.

3.1 Frailty Model Estimation

For the estimation of the parameters of the gamma- and Inverse Gaussian distributed force of mortality based on the observed force of mortality the proposed method by Su and Sherris (2012) is applied. The proposed method takes the variability of the observed mortality data of the cohorts into account and the method requires assumptions in order to estimate the parameters of the gamma- and Inverse Gaussian distributed force of mortality. The assumption made is that the observed number of deaths in the cohorts is a sample drawn from the underlying population which has the size of the cohort

exposure to risk, Ex. It is assumed that the individual forces of mortality are

distributed randomly with the mean being equal to, E(µ), and the variance being equal to, Var(µ). Due to the standard force of mortality being known and the mortality rates being independent we can apply the classical Central Limit Theorem. Thus, the distribution for the sample mean with the size of

the sample being Ex is close to the normal distribution with mean E(µ) and

variance Var(µ)_E

x .

(17)

Based on the assumptions, the mean and variance of the sample mean, ˆ

µ(x, z), are given by,

E(ˆµ(x, z)) = Ex E(µ(x, z)) Ex = µ(x, 1)δ δ + H(x) Var(ˆµ(x, z)) = Ex Var(µ(x, z)) E2 x = µ(x, 1) 2_δ Ex (δ + H(x))2

Thus, the log likelihood function that needs to be maximized is given by,

L(ˆµ(x, z)|Ex, x, δ, α, β) = X x −0.5 log(2π) + log(σ2_{(x, δ, α, β)} − (ˆµ(x, z) − µ(x, δ, α, β)) 2 2σ2_{(x, δ, α, β)} . (31)

with µ(x, δ, α, β) and σ2_{(x, δ, α, β) being the mean and variance of the gamma}

distributed frailty.

Next we assume that the frailty follows a Inverse Gaussian distribution, the mean and the variance of the individual force of mortality at age x are then given by, E(µ(x, z)) = µ(x, 1) λ λ + H(x) 1₂ Var(µ(x, z)) = µ(x, 1) 2 2 λ (λ + H(x))3 1₂

Based on the assumptions, the mean and the variance for the sample mean, ˆ

µ(x, z), are given by,

E(ˆµ(x, z)) = Ex E(µ(x, z)) Ex = µ(x, 1) λ λ + H(x) 1₂ Var(ˆµ(x, z)) = Ex Var(µ(x, z)) E2 x = µ(x, 1) 2 2 Ex λ (λ + H(x))3 1₂

Thus, the log likelihood function that needs to be maximized is given by,

(18)

with µ(x, λ, α, β) and σ2_{(x, λ, α, β) being the mean and variance of the}

(19)

3.2 Phase Type Model Estimation

Since many physiological functions decrease after a certain age Su and Sher-ris (2012) make the assumption, which we will also adopt, that the transition

rate, λi, is constant after a certain age.

λi =

(

λi if i ≤ k

λ otherwise

where k = 4. For the death rate, qi, the assumption is made that it is an

increasing function with age. From the data it can be observed that the Dutch death rate grows exponentially with age, see figures 2 and 3 . The following assumption is made,

qi = γ + α exp(β i)

The parameters a, b, γ, λi for i ≤ k (k = 4) and λ for are all parameters of

the model. The intention of the formulation of the model is to mark out a de-velopmental period of the phases for very young ages, by the different values

of λi. Different estimation methods can be applied for estimating the

(20)

squared errors for the death probability, SSE = ω X x=0 (q(x) − ˆq(x))2 s(x) (33)

(21)

4 Data

The Dutch mortality data was collected from the Human Mortality Database University of California, Berkeley (USA), and Max Planck Institute for De-mographic Research (Germany) (HMD). The data obtained consists of cen-tral rates of death for males and females from the cohorts 1945 until 1950. The obtained central rates of death are treated as an estimation for the co-hort forces of mortality. The reason we select the coco-horts from 1945 until 1950 is because the people in these cohorts represent the population that are retired recently or will retire soon. The selected population is the population that has the highest demand for pension, annuity contracts. In the following figures we present the log transformed cohort forces of mortality for males and females, where the age range starts at 30.

30 40 50 60 70 2.5 3.0 3.5 4.0 4.5 Male Age log( m *10,000) ₁₉₄₅ 1946 1947 1948 1949 1950

(22)

30 40 50 60 70 2.0 2.5 3.0 3.5 4.0 Female Age log( m *10,000) 1945 1946 1947 1948 1949 1950

Figure 3: Log transformed female forces of mortality

(23)

0 20 40 60 80 100 1 2 3 4 5 6 7 8 Age log( m *10,000) 1880 1900 1920 1940 1960

Figure 4: Log mortality rates, cohorts 1880-1960

1850 1900 1950 2000 30 40 50 60 70 80 Year lif e e xpectancy

Figure 5: Life expectancy at birth for males and females, 1850-2009

(24)

expectancy increases and the mortality rate decreases over time.

For the frailty models the data needed as input, is the central exposure to risk of the cohort and the force of mortality of the cohort. The central rates of death are treated as an estimation for the cohort forces and the cohort central exposure to risk are collected directly from the Human Mortality Database.

The required format of data for the estimation of the phase-type aging model are the death rate and the survival probability. The death rate can be col-lected from the relation shown below and assuming an uniform distribution of deaths,

q(x) = m(x)

1 + 0.5 m(x)

From the obtained death rates the survival probability can be collected.

(25)

0 10 20 30 40 50 60 70 1 2 3 4 5 6 7 8 Females Age log( m *10,000) 1945 1946 1947 1948 1949 1950

Figure 7: Log transformed female death rates

(26)

5 Results

5.1 Results Model Estimation

5.1.1 Gamma Frailty Model

The estimated values of the parameters of the gamma frailty model estimated from the maximization of the log likelihood function, 31, are given in table 2 and 3. Next to the estimated values of the parameters of the gamma frailty model, the maximum likelihood and the coefficient of determination are shown in table 2 and 3.

The high values of the coefficients of determination, R2, indicate a good

(27)

30 40 50 60 70 2.5 3.0 3.5 4.0 4.5 Male 1945 Gamma Age ln(10,000 u) Observed Gamma Fit 30 40 50 60 70 2.0 2.5 3.0 3.5 4.0 4.5 Male 1946 Gamma Age ln(10,000 u) Observed Gamma Fit 30 40 50 60 70 2.0 2.5 3.0 3.5 4.0 4.5 Male 1947 Gamma Age ln(10,000 u) Observed Gamma Fit 30 40 50 60 70 2.0 2.5 3.0 3.5 4.0 4.5 Male 1948 Gamma Age ln(10,000 u) Observed Gamma Fit 30 40 50 60 70 2.5 3.0 3.5 4.0 Male 1949 Gamma Age ln(10,000 u) Observed Gamma Fit 30 40 50 60 70 2.5 3.0 3.5 4.0 Male 1950 Gamma Age ln(10,000 u) Observed Gamma Fit

(28)

30 40 50 60 70 2.0 2.5 3.0 3.5 4.0 Female 1945 Gamma Age ln(10,000 u) Observed Gamma Fit 30 40 50 60 70 1.5 2.0 2.5 3.0 3.5 4.0 Female 1946 Gamma Age ln(10,000 u) Observed Gamma Fit 30 40 50 60 70 1.5 2.0 2.5 3.0 3.5 4.0 Female 1947 Gamma Age ln(10,000 u) Observed Gamma Fit 30 40 50 60 70 1.5 2.0 2.5 3.0 3.5 4.0 Female 1948 Gamma Age ln(10,000 u) Observed Gamma Fit 30 40 50 60 70 1.5 2.0 2.5 3.0 3.5 4.0 Female 1949 Gamma Age ln(10,000 u) Observed Gamma Fit 30 40 50 60 70 2.0 2.5 3.0 3.5 Female 1950 Gamma Age ln(10,000 u) Observed Gamma Fit

(29)

30 40 50 60 70 2.5 3.0 3.5 4.0 4.5 Male 1945 Gamma Age ln(10,000 u) Observed Gamma Fit No Frailty 30 40 50 60 70 2.0 2.5 3.0 3.5 4.0 4.5 Male 1946 Gamma Age ln(10,000 u) Observed Gamma Fit No Frailty 30 40 50 60 70 2.0 2.5 3.0 3.5 4.0 4.5 Male 1947 Gamma Age ln(10,000 u) Observed Gamma Fit No Frailty 30 40 50 60 70 2.0 2.5 3.0 3.5 4.0 4.5 Male 1948 Gamma Age ln(10,000 u) Observed Gamma Fit No Frailty 30 40 50 60 70 2.5 3.0 3.5 4.0 Male 1949 Gamma Age ln(10,000 u) Observed Gamma Fit No Frailty 30 40 50 60 70 2.5 3.0 3.5 4.0 Male 1950 Gamma Age ln(10,000 u) Observed Gamma Fit No Frailty

(30)

30 40 50 60 70 2.0 2.5 3.0 3.5 4.0 Female 1945 Gamma Age ln(10,000 u) Observed Gamma Fit No Frailty 30 40 50 60 70 1.5 2.0 2.5 3.0 3.5 4.0 Female 1946 Gamma Age ln(10,000 u) Observed Gamma Fit No Frailty 30 40 50 60 70 1.5 2.0 2.5 3.0 3.5 4.0 Female 1947 Gamma Age ln(10,000 u) Observed Gamma Fit No Frailty 30 40 50 60 70 1.5 2.0 2.5 3.0 3.5 4.0 Female 1948 Gamma Age ln(10,000 u) Observed Gamma Fit No Frailty 30 40 50 60 70 1.5 2.0 2.5 3.0 3.5 4.0 Female 1949 Gamma Age ln(10,000 u) Observed Gamma Fit No Frailty 30 40 50 60 70 2.0 2.5 3.0 3.5 Female 1950 Gamma Age ln(10,000 u) Observed Gamma Fit No Frailty

(31)

5.1.2 Inverse Gaussian Frailty Model

The estimated values of the parameters of the inverse Gaussian frailty model estimated from the maximization of the log likelihood function 32 are given in table 4 and 5. Next to the estimated values of the parameters of the inverse Gaussian frailty model, the maximum likelihood and the coefficient of determination are given in table 4 and 5.

(32)

30 40 50 60 70 2.5 3.0 3.5 4.0 4.5 Male 1945 IG Age ln(10,000 u) Observed IG Fit 30 40 50 60 70 2.0 2.5 3.0 3.5 4.0 4.5 Male 1946 IG Age ln(10,000 u) Observed IG Fit 30 40 50 60 70 2.0 2.5 3.0 3.5 4.0 4.5 Male 1947 IG Age ln(10,000 u) Observed IG Fit 30 40 50 60 70 2.0 2.5 3.0 3.5 4.0 4.5 Male 1948 IG Age ln(10,000 u) Observed IG Fit 30 40 50 60 70 2.5 3.0 3.5 4.0 Male 1949 IG Age ln(10,000 u) Observed IG Fit 30 40 50 60 70 2.5 3.0 3.5 4.0 Male 1950 IG Age ln(10,000 u) Observed IGFit

(33)

30 40 50 60 70 2.0 2.5 3.0 3.5 4.0 Female 1945 IG Age ln(10,000 u) Observed IG Fit 30 40 50 60 70 1.5 2.0 2.5 3.0 3.5 4.0 Female 1946 IG Age ln(10,000 u) Observed IG Fit 30 40 50 60 70 1.5 2.0 2.5 3.0 3.5 4.0 Female 1947 IG Age ln(10,000 u) Observed IG Fit 30 40 50 60 70 1.5 2.0 2.5 3.0 3.5 4.0 Female 1948 IG Age ln(10,000 u) Observed IG Fit 30 40 50 60 70 1.5 2.0 2.5 3.0 3.5 4.0 Female 1949 IG Age ln(10,000 u) Observed IG Fit 30 40 50 60 70 2.0 2.5 3.0 3.5 Female 1950 IG Age ln(10,000 u) Observed IGFit

(34)

30 40 50 60 70 2.5 3.0 3.5 4.0 4.5 Male 1945 IG Age ln(10,000 u) Observed IG Fit No Frailty 30 40 50 60 70 2.0 2.5 3.0 3.5 4.0 4.5 Male 1946 IG Age ln(10,000 u) Observed IG Fit No Frailty 30 40 50 60 70 2.0 2.5 3.0 3.5 4.0 4.5 Male 1947 IG Age ln(10,000 u) Observed IG Fit No Frailty 30 40 50 60 70 2.0 2.5 3.0 3.5 4.0 4.5 Male 1948 IG Age ln(10,000 u) Observed IG Fit No Frailty 30 40 50 60 70 2.5 3.0 3.5 4.0 Male 1949 IG Age ln(10,000 u) Observed IG Fit No Frailty 30 40 50 60 70 2.5 3.0 3.5 4.0 Male 1950 IG Age ln(10,000 u) Observed IG Fit No Frailty

(35)

1945-30 40 50 60 70 2.0 2.5 3.0 3.5 4.0 Female 1945 IG Age ln(10,000 u) Observed IG Fit No Frailty 30 40 50 60 70 1.5 2.0 2.5 3.0 3.5 4.0 Female 1946 IG Age ln(10,000 u) Observed IG Fit No Frailty 30 40 50 60 70 1.5 2.0 2.5 3.0 3.5 4.0 Female 1947 IG Age ln(10,000 u) Observed IG Fit No Frailty 30 40 50 60 70 1.5 2.0 2.5 3.0 3.5 4.0 Female 1948 IG Age ln(10,000 u) Observed IG Fit No Frailty 30 40 50 60 70 1.5 2.0 2.5 3.0 3.5 4.0 Female 1949 IG Age ln(10,000 u) Observed IG Fit No Frailty 30 40 50 60 70 2.0 2.5 3.0 3.5 Female 1950 IG Age ln(10,000 u) Observed IG Fit No Frailty

(36)

5.1.3 Phase-type aging Model

The estimated parameters for the phase-type aging model and the coefficient of determination are given in table 6 for male cohort and in table 7 for the female cohort.

(37)

(38)

(39)

(40)

5.2 Heterogeneity in the models

5.2.1 Heterogeneity in Frailty models

(41)

0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0

Density gamma frailty Male 1945

Frailty Density Age 0 Age 30 Age 60 Age 90 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0

Frailty Density Age 0 Age 30 Age 60 Age 90

(42)

0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0

Density gamma frailty Female 1945

Frailty Density Age 0 Age 30 Age 60 Age 90

(43)

0 20 40 60 80 100 0.0 0.2 0.4 0.6 0.8 1.0

Mean gamma frailty male 1945

age mean fr ailty 0 20 40 60 80 100 0.0 0.2 0.4 0.6 0.8 1.0

Mean gamma frailty female 1945

age mean fr ailty 0 20 40 60 80 100 0 2 4 6 8 10

Variance gamma frailty male 1945

age v ar fr ailty 0 20 40 60 80 100 0 5 10 15

Variance gamma frailty female 1945

age

v

ar fr

ailty

Figure 20: The mean and variance of the gamma distributed frailty for the male and female cohort 1945

In figure 21 the death probabilities with the different values of frailty are shown. The death probability is obtained via the individual force mortality,

q(x, z) = 1 − exp(−µ(x, z))

(44)

vary significanlty and survival prospects of individuals with different frailties are substantially different.

0 20 40 60 80 100 0.0 0.2 0.4 0.6 0.8 1.0

Death probability IG Male 1945

Age Death Prob z=1 z=0.01 z=0.001 z=0.0005 z=0.00001 cohort 0 20 40 60 80 100 0.0 0.2 0.4 0.6 0.8 1.0

Age Death Prob z=1 z=0.01 z=0.001 z=0.0005 z=0.00001 cohort

(45)

5.2.2 Heterogeneity in the Phase-type model

Under a phase-type distribution we assume that every individual starts at state 1, the initial state. As the individuals age they go from state to state. The probability of an individual to be in state i at age x under the phase-type distribution is given by,

P (i, x) = [α exp(∆x)]_i

The conditional probability of an individual that is in state i, conditional on survival of that individual to age x is given by,

Π(i, x) = P (i, x) s(x) = α exp(∆x) α exp(∆x)ι i .

(46)

0 20 40 60 80 100 0.00 0.05 0.10 0.15 0.20

Distibution physiological Male cohort 1945

Age Density age 20 age 40 age 60 age 80

Figure 22: Distribution of different physiological ages at different ages for the male 1945 cohort

(47)

6 Models Applied to Life Annuities

In this section we will apply both the models of mortality heterogeneity to pricing whole life annuities at age 65. If the annuity rates differ among differ-ent levels of heterogeneity than the heterogeneity is significant. If the annuity rates do not differ that much among the different levels of heterogeneity and are close to the cohort annuity rate than the the heterogeneity is not signif-cant. First we consider the frailty model. For the pricing of the whole life annuities we will use different values for the frailty parameter. This way we can evaluate the impact of the different levels we have of mortality on the whole life annuities. The whole life annuities at age 65 are calculated by,

∞

X

k=0

(1 + i)−k kp65.

where have that

kpx = s(x + t) s(x) = exp(− Rx+t 0 µ(y, z)dy)

exp(−R₀xµ(y, z)dy)

= exp(−

Rx+t

0 α exp(βy)zdy)

exp(−Rx

0 α exp(βy)zdy)

(48)

0.000 0.004 0.008 0 5 10 15 20

Annuity values for different fralties

frailty

Ann

uity v

alue

Figure 23: The annuity values for different frailty levels for the male cohort 1945

From figure 23 we see that the when the health condition of an individual is getting worse, the frailty increases, and the annuity rates drop significantly. In table 8 the annuity rates together with the corresponding value of the cumulative distribution function are shown.

Second we consider the phase type aging model. In figure 24 the distribution of physiological ages at age 65 for the male 1945 cohort is displayed. The whole life annuities at age 65 are again calculated by,

(49)

where have that kpx= s(x + t) s(x) = α exp(∆(x + t)) α exp(∆x)

In table 9 the annuity rates for the male cohort 1945 are shown.

0 20 40 60 80 100 0.00 0.05 0.10 0.15 0.20

Distibution physiological Male cohort 1945

Age

Density

age 65

Figure 24: distribution of physiological ages at age 65 for the male 1945 cohort

(50)

the annuity rates of the phase type aging model are higher compared to the annuity rates of the frailty model. The difference in the annuity rates of the two models can be explained by the different assumptions the two models are based on. Under the frailty model the assumption of fixed frailty at birth is made. This is in contrast with the phase type aging model where the dis-tribution by physiological age is not being fixed at birth. The disdis-tribution by physiological age changes over time.

Both the models show that heterogeneity is significant by the different values of the annuity rates and therefore there is a substantial level of heterogeneity in the population.

In table 1 an overview of the annuity rates for individuals with different phys-iological age and different frailty is shown. We can again observe that the range of the annuity rates of the frailty model is much wider compared to the range of the annuity rates of the phase type aging model and the annuity rates of the phase type aging model are higher compared to the annuity rates of the frailty model.

Phys. Age Annuity value % Popu.

80 11.95 99.7307 %

75 13.47 98.82551%

65 16.78 80.799388%

60 18.45 57.75189%

50 21.59 12.40602%

Frailty Annuity value % Popu.

0.01 0.43 99%

0.001 3.91 97 %

0.0001 10.76 72%

0.00005 12.94 56%

0.00001 17.53 13 %

(51)

7 Conclusion

(52)

8 Appendix

Cohort α β k, frailty parameter R2 _{Max LL}

1945 male 2.467e-05 1.132e-01 1.001e-01 0.998 -2010.644

1946 male 3.085e-05 1.047e-01 1.256e-01 0.999 -2150.144

1947 male 2.436e-05 1.102e-01 1.114e-01 0.999 -2341.506

1948 male 2.455e-05 1.107e-01 1.031e-01 0.999 -1878.179

1949 male 3.030e-05 1.078e-01 8.548e-02 0.998 -2308.683

1950 male 2.971e-05 1.074e-01 8.461e-02 0.998 -1558.217

(53)

Cohort α β k, frailty parameter R2 _{Max LL}

1945 female 1.293e-05 1.190e-01 5.858e-02 0.996 -1932.956

1946 female 1.148e-05 1.223e-01 5.063e-02 0.997 -2178.405

1947 female 1.396e-05 0.117 0.059 0.998 -1758.175

1948 female 1.152e-05 1.200e-01 5.911e-02 0.998 -1398.599

1949 female 1.167e-05 1.196e-01 5.931e-02 0.999 -630.6041

1950 female 1.258e-05 1.208e-01 4.636e-02 0.998 -1042.629

(54)

Cohort α β λ, frailty parameter R2 _{Max LL}

1945 male 0.00128 0.1735 0.0000127 0.9992 63.39181

1946 male 1.149e-03 1.644e-01 2.061e-05 0.9995 66.59183

1947 male 8.824e-04 1.699e-01 1.967e-05 0.9994 -57.56566

1948 male 1.149e-03 1.645e-01 2.009e-05 0.9993 29.65012

1949 male 1.699-03 1.582e-01 1.854e-05 0.9994 63.76586

1950 male 1.255e-03 1.592e-01 2.307e-05 0.9994 75.01698

(55)

Cohort α β λ, frailty parameter R2 _{Max LL}

1945 female 1.761e-03 1.749e-01 3.685e-06 0.9986 58.21946

1946 female 5.302e-04 1.681e-01 1.685e-05 0.9982 -160.424

1947 female 2.239e-04 1.667e-01 4.407-05 0.9986 -6.791753

1948 female 2.397e-04 1.714e-01 3.077e-05 0.9986 -39.73842

1949 female 5.951e-06 1.638e-01 1.992e-03 0.9991 70.00895

1950 female 0.0000498241 0.165 0.00021 0.9986 34.18398

(56)

par 1945 1946 1947 1948 1949 1950 λ1 4.023455 4.166523 3.760061 3.794597 3.993929 3.972738 λ2 1.674297 3.083609 2.132490 2.764864 2.224616 2.660362 λ3 0.5349795 2.344817 1.009833 0.8246569 0.7050273 0.7358404 λ4 0.4506401 0.102221 1.398826 1.184276 0.8408396 0.817814 q1 0.7799214 0.5449797 0.5898627 0.3408845 0.311359 0.4189027 q2 1 0.8404001 0.4624179 0.7693929 0.372837 0.5753705 q3 0 0.1676188 0 0 0 0.005121341 q4 0.002287283 0.001259815 0.004034511 0.007931748 0.005635723 0.003808611 λ .9138704 0.9607926 0.2833856 0.4250113 0.8356901 0.4019085 γ 0.0006787966 0.0005849612 0.0006498605 0.0005837152 0.0005965649 0.0006075781

a 2.813006e-06 5.487379e-06 1.162409e-06 7.093437e-06 6.591696e-06 5.572684e-06

b 0.1398511 0.1296796 0.3634568 0.2174078 0.1294664 0.2347134

R2 _0.9990464 _0.9979122 _0.9914491 _0.996336 _0.9977468 _0.9967373

(57)

par 1945 1946 1947 1948 1949 1950 λ1 3.574106 3.839413 3.767153 3.545113 4.003097 3.958333 λ2 2.076802 1.739036 2.110291 2.691664 2.2656 1.808674 λ3 0.1680747 0.1 0.2058198 0.3793745 1.199628 0.116744 λ4 0.8494958 0.80174 0.8991174 0.9325953 0.7120547 0.8801047 q1 0.1070345 0.7964649 0.5154226 0.5769437 0.3694142 0.3310571 q2 0.8596015 0.268925 0.315971 0.4750308 0.2685284 0.1754383

q3 0 0.000003018779 3.594927e-04 1.413280e-04 0 8.195661e-04

q4 0.008709431 0.004822894 0.005842015 0.005787393 5.476177e-03 0.002252384

λ 0.8346289 0.2146633 0.2614269 0.9518814 0.9015491 0.8630845

γ 0.00003929205 1e-06 1.839094e-04 3.148544e-04 3.213094e-04 3.196015e-05

a 0.00002610645 2.474928e-05 4.334547e-06 1.156928e-06 1.908727e-06 3.167886e-05

b 0.1045025 0.3068625 0.3333582 0.1467397 0.1406395 0.1008609

R2 _0.9996442 _0.9944895 _0.9983522 _0.9986491 _0.9990683 _0.9982658

(58)

Frailty Annuity value F (z) % Population 0.01 0.43 99% 0.001 3.91 97 % 0.0001 10.76 72% 0.00005 12.94 56% 0.00001 17.53 13 %

(59)

Physiological Age q65 Annuity value % Population 50 0.003822012 21.59 12.40602% 55 0.00666484 20.07 31.64566% 60 0.01166427 18.45 57.75189% 65 0.01977099 16.78 80.799388% 70 0.0316786 15.10 94.01152% 75 0.04745467 13.47 98.82551% 80 0.06646672 11.95 99.7307 %

(60)

References

Human mortality database.

Adele H. Marshall, Mariangela Zenga (2009). The XIIIth International Con-ference Applied Stochastic Models and Data Analysis. Vilnius Gediminas Technical University Publishing House Technika.

Carrey, D. (2002). Coping with population ageing in the netherlands. ECD Publishing (325).

Gompertz, Benjamin (1825). On the nature of the function expressive of the law of human mortality, and on a new mode of determining the value of life contingencies. Philosophical transactions of the Royal Society of London, 513–583.

Hougaard, Philip (1984). Life table methods for heterogeneous populations: Distributions describing the heterogeneity. Biometrika 71 (1), 75–83. Keyfitz, N. (1978). Improving life expectancy: an uphill road ahead.

Amer-ican Journal of Public Health 68 (10), 654–656.

Lang A., Arthur J.L. (1997). Parameter approximation for phase-type distri-butions in chakravarthy, s.r. and alfa a.s. editors. Matrix-analytic methods in stochastic models 189, 151–206.

Lin, X. Sheldon and Xiaoming Liu (2007). Markov aging process and phase-type law of mortality. North American Actuarial Journal 11 (4), 92–109. Manton, Kenneth G., Eric Stallard and James W. Vaupel (1986). Alternative

models for the heterogeneity of mortality risks among the aged. Journal of the American Statistical Association 81 (395), 635–644.

Neuts (1981). Matrix-geometric solutions in stochastic models - an algorith-mic approach. John Hopkins University Press.

(61)

mortality and implications for a viable life annuity market. Insurance: Mathematics and Economics 51 (2), 322–332.