On Insuring Mortality Risk in a Variable Pension
Agreement: A Study of a Stochastic Mortality
Model
Daniel Peters
S2377063
Master’s Thesis Actuarial Studies
Supervisor University of Groningen: Prof. Dr. T.K. Dijkstra Supervisors TKP Pensioen: Dirk-Jan Been & Ren´e de Vries AAG Second Assessor University of Groningen: Prof. Dr. P.A. Bekker
On Insuring Mortality Risk in a Variable Pension
Agreement: A Study of a Stochastic Mortality
Model
Daniel Peters
S2377063
Supervisors:
Prof. Dr. T.K. Dijkstra (University of Groningen)
Dirk-Jan Been (TKP Pensioen)
Ren´
e de Vries AAG (TKP Pensioen)
Second Assessor:
Prof. Dr. P.A. Bekker (University of Groningen)
March, 2018
Abstract
Contents
1 Introduction 1
2 Theoretical Framework 4
2.1 Actuarial notation . . . 4
2.1.1 Whole life annuity-due . . . 5
2.1.2 Joint life and reversionary annuity . . . 6
2.1.3 Yield curve . . . 7
2.1.4 Expected Future Lifetime . . . 9
2.1.5 Force of Mortality . . . 11
2.2 The Pension Plan . . . 11
2.2.1 Variable pension agreement . . . 11
2.2.2 Spouse Pension . . . 13
2.2.3 Characteristics of the specific fund . . . 14
2.2.4 Gender neutrality . . . 16
2.3 Reinsurance Contract . . . 18
3 Models and Simulations 20 3.1 Population Data . . . 21
3.2 Mortality Model Specifications . . . 22
3.2.1 Lee-Carter model . . . 22
3.2.2 Li-Lee model . . . 23
3.2.3 Poisson log-bilinear model . . . 25
3.2.4 AG2016 model . . . 26
3.3 The Bootstrap Approach . . . 29
3.3.1 Residual Bootstrap . . . 30
3.3.2 Poisson Bootstrap . . . 31
3.4 The Simulation Model . . . 32
3.4.1 Pension Fund Data . . . 32
3.4.2 Simulation Setup . . . 34
4 Results 39 4.1 Residual Bootstrap . . . 39
4.2 Poisson Bootstrap . . . 43
Bibliography 55
Appendices 57
A Tables 57
B Graphs and Figures 59
1
Introduction
The negative impact of the latest financial crisis on the Defined Benefit (DB) pension schemes has received extensive coverage in the Dutch media. Pension funds are struggling to make good on their promises: most do not compensate for inflation and some even lower nominal benefits. But Defined Contribution (DC) schemes have had a hard time as well. Since the interest rates are at historically low levels, the annuities that participants in a DC scheme are obliged to buy when they retire will generate significantly lower benefits than anticipated and hoped for. The matter was taken up by Dutch parliament and as a result on 1 September 2016 the ‘Premium Schemes Improvements Act’ (in Dutch: Wet Verbeterde Premieregeling) came into effect.
The act aims to improve the performance of DC plans as it permits the retiree to keep investing his retirement account after the retirement date. This should result in a higher expected pension benefit compared to a life annuity bought at retirement. However, it also gives more uncertainty about the height of the pension benefit: investment returns could disappoint, decreasing the amount of retirement capital and subsequently monthly retirement benefits. Apart from the exposure to investment risk, a retiree with a variable retirement benefit runs the risk that his retirement capital will not suffice if he would live longer than expected. It is therefore crucial to include a mechanism that protects the individual against outliving his retirement capital.
life-long retirement benefit. However, in a variable pension agreement there is no such arrangement. As a consequence, the Premium Schemes Improvement Act provides clear rules which state that an individual’s retirement benefit cannot fluctuate due to the individual surviving longer than expected. The individual is not allowed to bear this risk and it should therefore be shared collectively or be reinsured. Since the group of participants that is required to properly share this type of mortality risk should be significantly large and at the start of a variable pension agreement this group will not be, reinsurance is necessary.
in the DC plan insure the mortality risk in the variable retirement benefit at the fund’s DB population.1
We aim to assess the risk that a reinsurer takes when it takes over mortality risk from individual retirees in a variable pension agreement. Before we discuss the reinsurance contract, we will study the stochastic mortality model that is commonly used by pension providers. We will discuss the model’s properties and with the aid of a bootstrap technique, we test its robustness. The mortality model produces forecasts of future mortality rates using standard time series models. The randomness around future mortality rates comes solely from the simulation of the model’s time-varying parameters. We enhance the model by incorporating uncertainty in the model’s other parameters.
The purpose of this thesis is therefore twofold. We will focus on the following research questions: “What is an actuarially fair price for insuring mortality risk in a variable pension agreement?” and “What is the effect of parameter uncertainty on forecasted life expectancies in the AG2016 stochastic mortality model?”. Fur-thermore, the hybrid pension fund in question decided to initially not ask a risk premium. However, financial institutions are required to hold capital reserves and we will therefore discuss what an appropriate risk premium that covers these cap-ital reserves should be. As the majority of pension providers and life insurers in the Netherlands base the prices of their products on the stochastic mortality model of the Royal Actuarial Association (AG), we begin our study there. The mortality
model will be examined and a practical investigation into the pricing of a real-world problem is conducted.
This thesis is structured as follows. Section 2 provides the theoretical framework with the actuarial notation, the characteristics of the variable pension agreement and the corresponding reinsurance contract. In Section 3 we discuss the stochastic mor-tality models that are used, we explain the bootstrap technique that will be applied and we define the Monte Carlo simulation model that we apply to investigate the reinsurance contract. The results are presented in Section 4 and Section 5 discusses, concludes and provides suggestions for further research.
2
Theoretical Framework
2.1
Actuarial notation
Actuarial science has its own notation, International Actuarial Notation, which cap-tures the probabilities and functions that are of interest and are useful to actuaries. The next section will explain the basic concepts of survival and mortality probabili-ties, and several type of annuities as presented in Dickson et al. (2013).
Let (x) denote a life aged x, with x ≥ 0. The death of (x) can occur at any age greater than x, and let Tx be the future lifetime of (x). We define tpx as the
probability that (x) survives to at least age x + t, andtqx as the probability that (x)
dies before age x + t. Subscript t may be dropped if its value is 1, so pxrepresents the
probability that (x) survives to at least age x + 1. Similarly, qx is the probability that
life contingent annuities. An annuity is a series of payments to an individual as long as the individual is alive on the payment date. The series of payments can occur at different intervals, such as annually, monthly or (theoretically) continuously, but we will only consider annual life annuities. Each year 1 unit will be paid out, conditional on the survival of one or more lives at the payment date.
2.1.1 Whole life annuity-due
We will start with the simplest situation, a whole life annuity-due, which is an annuity paid in advance conditional on the survival of one life, that of (x), with the first payment occurring immediately. The second payment will occur one year later, provided that (x) is alive, and so on. For the valuation of the series of payments, we use a fixed discount factor v = (1 + i)−1 where i is the annual interest rate. Figure 1 shows the series of payments and the associated discount factors and probabilities.
0 1 2 3 Time Amount Discount Probability 1 1 1 1 v px 1 v2 2px 1 v3 3px
Figure 1: Time-line diagram for whole life annuity-due
We can derive the present value of the random series of payments with the help of an indicator function
Note that tpx = E [I(Tx > t)], i.e. the probability that the death of (x) happens at
a later time than time t. We can then construct the random variable
Y = I(Tx > 0) + vI(Tx > 1) + v2I(Tx > 2) + v3I(Tx > 3) + . . . (2)
Since the expected value of the sum of random variables is equal to the sum of their expected values we have
E[Y ] = 1 + vpx+ v22px+ v33px+ . . . (3)
and we can write
¨ ax = E[Y ] = ∞ X k=0 vkkpx. (4)
where we denote ¨ax as the expected present value of a whole life annuity-due.
2.1.2 Joint life and reversionary annuity
Two more types of annuities of interest are the joint life annuity and the reversionary annuity. Since we will be dealing with two lives, we will introduce the notation (y) for the second individual aged y. When dealing with pension plans the worker is often called (x) and his or her spouse is called (y). A frequently made assumption is that the male in a relationship is three years older than the female, and throughout our calculations we too will make this assumption. A joint life annuity ¨axy is a series of
payments of 1 per year while both (x) and (y) are still alive. A reversionary annuity ¨
Assuming the probabilities of death of (x) and (y) are independent we can write
tpxy = P [Tx > t and Ty > t] =tpx tpy. (5)
Using a similar derivation as for ¨ax, we find for the joint life annuity
¨ axy = ∞ X k=0 vkkpxy, (6)
and for the reversionary annuity
¨ ax|y = ∞ X k=0 vk(1 −kpx)kpy. (7) 2.1.3 Yield curve
In practice, the value of annuities are not derived using a fixed interest rate as dis-cussed above. As the annual interest earned on a risk-free investment for a period of 1 year is typically different from the annual interest of a 20-year risk-free investment, it is therefore crucial to discount the expected series of the payments of an annuity with different discount factors. Three key quantities of interest are the yield curve, discount curve and forward rate curve. Let yt denote the t-year spot rate of interest:
the yield per year on a t-year zero-coupon bond. If one would plot yt over time, one
would get a graphical expression of the yield curve. We then define v(t) = (1 + yt)−t
To derive the expression for the one-year forward rate curve, consider the following two-year investments. 1) Investing 1 unit at a pre-determined 2-year rate versus 2) investing 1 unit at a pre-determined 1-year rate for one year and then investing the result at the end of the first year for another year at a pre-determined rate. These two investments should yield the same return in two years. For investment 1) we have (1 + y2)2 at the end of the second year. Investment 2) will yield (1 + y1)(1 + x).
x is actually called the one-year forward rate in year 1 (f1,2) and we can construct it
by (1 + y2)2 = (1 + y1)(1 + f1,2), and hence f1,2 = (1 + y2)2 (1 + y1) − 1.
The one-year forward rates for higher years are constructed analogously, and for the first year we have f0,1 = y1. A graphical illustration of the relation between spot and
one-year forward rates is given in Figure 2.
0 (1 + f0,1) 1 (1 + f1,2) 2 (1 + f2,3) 3 (1 + y1) (1 + y2)2 (1 + y3)3
Figure 2: Relation between spot and one-year forward rates
The same holds for the years 30 to 40, 40 to 50 and so on. The Dutch regulator determines the yield curve with which a pension provider has to discount its future cash flows. For our calculations we will use the yield curve provided by the Dutch Central Bank on 31-12-2017. A plot of the yield curve is provided in Appendix B.
Following the guidelines from the Dutch regulator we will extract the future yield curve from the given forward rates. The method is described as follows. At t = 0 we observe the spot yields of interest y1, y2, . . . , y100and corresponding one-year forward
rates f0,1, f1,2, . . . , f99,100. For t = 1 we assume that the one-year forward rates do
not change and therefore the next one-year forward rate we will observe still equals f1,2. We then derive the one-year spot rate at t = 1 by y1 = f1,2; and for the two-year
spot rate we find (1 + y2)2 = (1 + f1,2)(1 + f2,3), et cetera. This means that for the
valuation of an annuity at t = 1 we use a (slightly) different discount curve than at t = 0.
2.1.4 Expected Future Lifetime
We have already defined Tx as the future lifetime of (x). Another quantity of interest
is the expected future lifetime of (x): ex = E[Tx]. If we assume that deaths happen
on average halfway a calender year, we can define the expected future lifetime as
ex = 1 2+ ∞ X k=1 kpx. (9)
life expectancy does not take future improvement (or deterioration) of survival prob-abilities into account. For example, p67 (= 1 − q67) in 2017 is (slightly) different from
p67 in 2018 (e.g. due to improvements in health care). We can construct the period
life expectancy of (67) in 2017 with the probabilities of death q67, . . . , q120 only for the
year 2017. However, (67) will actually face another q68 in 2018, and another q69 in
2019, and so on. The cohort life expectancy does take this into account. An example of an application of the period life expectancy is the determination of the State Pen-sion age2 in the Netherlands. In most applications however, calculations are based
on death probabilities of a projection table where future improvements are taken into account. As an illustration of the vast differences between the two quantities: based on the projection table of the Royal Actuarial Association (AG2016) we find a period life expectancy of a 65 year old in 2016 of 18.4 years and the corresponding cohort life expectancy is 20.0 years.
To be precise, we should use a notation that distinguishes different (one year) death probabilities for different years, such as q67(2016) and q67(2017). We choose not
to (except for Section 3) for the simple reasons that it would make notation somewhat confusing and it should be evident which probabilities are used. We furthermore note that Dutch pension law states that pension providers have to use a projection table that incorporates a time trend.3
2In Dutch: AOW-leeftijd
3
2.1.5 Force of Mortality
Another fundamental concept in modelling future lifetime is the force of mortality. We have defined qx as the probability that (x) dies before age x + 1, that is qx =
P [Tx ≤ 1]. The force of mortality at age x is denoted by µx and is defined as
µx = lim dx→0+
1
dxP [Tx ≤ dx]. (10)
We can interpret µx dx as the probability that (x) dies before age x + dx. dx can
for example be interpreted as a single day. The central rate of mortality mx is then
the average force of mortality, averaged over the year of age. If we assume that the force of mortality remains constant over each year of age x we get the equality (see Cairns et al., 2009)
mx = µx,
and we can therefore write
qx = 1 − exp(−µx).
2.2
The Pension Plan
2.2.1 Variable pension agreement
lower pension benefit for the retiree. To illustrate the effect of interest rate changes on the value of an annuity: using the mortality rates from the AG2016 projection table, with an accrued capital of 1,000,000 Euro and a fixed annual interest rate of 4%, a male aged 67 in 2016 can buy an annual income until death of 84,650 Euro. When valuing annuities using the risk free rate provided by the Dutch Central Bank (on 30 September 2017), the annual pension benefit diminishes by 24% to 64,250 Euro.
As in the variable pension agreement the retiree is allowed to keep investing his retirement account in the capital market, every year the height of the monthly pension benefit will be determined by the pension provider and paid out of the retiree’s account. The retiree is therefore able to bear investment risk with which he could earn a positive return. The probabilities of survival tpx that are used for
valuing annuities are values between 0 and 1, but the survival of an actual person is binary: one survives to age x + t or one dies between ages x and x + t. The group of retirees who choose for a varying pension agree that in the event of death, their retirement account is distributed among the survivors in the group. So from the perspective of the survivor, one receives a positive biometric return: a bonus for staying alive. The biometric return assures that the individual retirement account does not run low and is able to provide a sufficient pension benefit. In the event of death, one’s biometric return is -1.
The height of the variable pension benefits is thus based on projected mortality rates. In a pool of n retirees of age x, n × px are projected to survive to age
and hence the biometric return will not be sufficient to complement annual pension benefits. For a large population, actual deaths will be close to the projected amount of deaths, but for small populations it would be a high risk that subsequent years there will be no (or far less) deaths than expected. As a large population might be able to properly share this mortality risk, a small group of people would risk severe reductions in retirement benefits. Hence, the risk that a small population actually behaves unlike expected should be reinsured in the variable pension agreement.
2.2.2 Spouse Pension
When choosing between a fixed and a variable pension benefit, the retiree is also able to choose whether in the event of his death, part of his pension will still be paid to his living spouse. It is custom to set spouse pension on 70% of one’s old age pension. Clearly, with the same amount of accrued capital, a person choosing for a spouse pension gets a lower annual pension benefit than a person who only receives an old age pension. One has to choose between a whole life annuity-due or a combination of a whole life annuity-due with a reversionary annuity. The latter is expected to pay out for a longer period and therefore its annual payment should be lower. Let X be the percentage of the old age pension of a whole life annuity-due compared with the combination of a whole life annuity-due with a reversionary annuity, then the following equation holds:
¨
Hence,
X = ¨ax
¨
ax+ 0.7¨ax|y
. (12)
In other words, if we assume X = 80%, then if (x) could receive an old age pension (OP) of 10,000 Euro per year conditional on his own life, he also has a choice for an OP of 8,000 Euro conditional on his life, and in the event of death his living spouse will receive 70% of 8,000 Euro (5,600 Euro).
2.2.3 Characteristics of the specific fund
For the specific pension fund, we will focus on the following hybrid pension agree-ment: Annual wages below 45,000 Euro are considered in the DB plan and wages between 45,000 Euro and the fiscal maximum (for the year 2018: 105,075 Euro) are in a DC plan. Dutch pension law states that a wage above the fiscal maximum is not considered in tax-favourable pension plans; this maximum is adjusted annually. Until recently, the DC participant in this particular fund was forced to purchase a life annuity at the DB fund. At the moment of retirement, the DB fund received the pension capital and invested in the capital market to earn a return with which it would finance and provide indexation on the DC participant’s pension. In the variable pension agreement the DB fund will not receive the pension capital and hence will not be able to earn a (possibly) positive return. The DC fund will invest the retirement account in the capital market on behalf of the retiree and therefore investment risk stays with the individual.
interest and mortality rates. The pension capital is invested in the capital market and each month the benefit is paid out from the capital. When the retiree is still alive at the end of the year (t = 1), next year’s benefit will be determined. The difference between the two benefits as a result of investment returns and changed interest rates are borne by the retiree. On top of the return on investments and changes in interest rates, the retiree will receive a biometric return. This ex ante fixed bonus is based on projected mortality rates of the past year and paid by actual deaths of the past year, i.e. the biometric return per person is fixed but the financing is stochastic. Changes in biometric return that could happen due to actual death rates (the retirement population lives longer than expected) will be reinsured in the variable annuity contract. At t = 2, when new mortality forecasts are published4, the new monthly benefits for the coming year are determined. Moreover, the pension fund has the right to adjust projected mortality rates for their population in the event of a significant deviation in actual death rates compared to projected death rates. As projected mortality rates are adjusted, ones next year’s annual pension benefit will also be adjusted. Changes in pension benefits due to the newly forecasted mortality rates are again borne by the retiree (if the whole Dutch population is expected to live longer than previously estimated, this will have an effect on the individual’s retirement benefit). Table 1 summarises how changes in the variable pension payments after the retirement date are borne.5
4New mortality forecasts are published by the Royal Actuarial Association every two years. 5In line with art.
Table 1: Premium Schemes Improvement Act
Type of Risk How borne?
Investment Risk Individually
Actual Mortality Results Collectively Changes in Life Expectancy Individually
2.2.4 Gender neutrality
According to European law, men and women have to be treated equally. The effect on pensions is that even though there is strong evidence that men and women have a different life expectancy, their pension benefit cannot differ due to a difference in gender. In other words, a man and woman of the same age, with the same accrued capital should receive the same pension. The height of the pension will be based on the ratio of males and females in a pension fund. One way to determine a gender neutral annuity factor is as follows:
Assume that at t = 0, 80% of the workers (x) are male and the remaining 20% are female. All participants will retire coming year (at age 67) and therefore the fund decides to set a sex ratio of 80/20. The annuity factor for males is denoted ¨aM
x ,
the corresponding factor for females is ¨aF
x. Then the gender neutral annuity factor
is then given by
¨
ax = 0.8¨aMx + 0.2¨aFx,
i.e. we have a weighted average of the gender-dependent annuity factors ¨aM
x and ¨aFx.
not forget that sex ratios change due to the fact that the probabilities of survival for males are (slightly) lower than for females of the same age.
We started with a ratio of 80/20 and using the projected survival probabilities we can see that in t years, the percentage of men will be
0.8tpMx
0.8tpMx + 0.2tpFx
.
Figure 3 shows the dynamics of the sex-ratio for the years 2018 to 2058 using the mortality rates provided by the AG2016 projection table. We see that after 20 years, approximately 75 percent of our fictional retirees will be men. As there is a rather small difference between the annuity factors for men and women at high ages, the effect of working with the old 80/20 ratio instead of the true ratio is minimal.
It gets more interesting when determining the corresponding annuity factors for the spouse (y) and working with a fixed age difference. Using the convention that the male spouse is 3 years older and a female spouse is 3 years younger we find at time t = 0:
¨
ay = 0.8¨aMy + 0.2¨a F y,
2020 2030 2040 2050
0.65
0.70
0.75
0.80
Sex Ratio using AG2016
Time
P
ercentage of Men
Figure 3: Percentage of Males
2.3
Reinsurance Contract
(with a high accrued capital) staying alive is larger than for non-wealthy.
Even if the reinsurer correctly projects mortality rates, actual deviations from the probabilities of death can still happen. Assume we have a population of n=10,000 with a probability of death qx of 1%. Assuming the total number of deaths are
binomially distributed with probability qx, then the expected number of deaths is
100 and a 95% confidence interval for the total number of deaths is [81, 120]. Hence, additionally to the risk that the reinsurer does not correctly estimates its population’s true mortality rate, there is the risk that (especially for a small population) actual deaths and forecasted deaths do not match. We can therefore identify two types of risk borne by the reinsurer: the risk of incorrectly estimating the true probabilities of deaths and the risk of deviations from the actual death rate.
3
Models and Simulations
We will ultimately investigate the mortality rates and life expectancies that are obtained from the mortality model that was published by the Royal Actuarial Asso-ciation (2016). A vast amount of literature on (stochastic) mortality models exists with the most prominent publications being the Lee-Carter model (Lee and Carter, 1992) and the Cairns-Blake-Dowd model (Cairns et al., 2006). Renshaw and Haber-man (2006) extend and generalise the Lee-Carter model by including age-specific cohort effects and Currie (2006) presents an alternative Age-Period-Cohort model. For the interested reader, a comparison of the above mentioned models is presented in Cairns et al. (2009).
As an extension to the Lee-Carter model, Li and Lee (2005) proposed a model to forecast mortality for a group of populations. They intentionally left the definition of the term ‘group’ vague as they leave it to the forecaster what to define as groups. The two applications that Li and Lee discuss in their paper are: 1) two-sex mortality and 2) mortality for 15 similar (low-mortality) countries. The idea is that in the long-run mortality of the separate populations converge to a common trend.
AG model is used by most pension funds in the Netherlands due to its simplicity and transparency. Moreover, the AG model is one of the few models accepted by the Dutch regulator for determining the provision for pension liabilities. In order to incorporate all sources of variability (not only the uncertainty in the time-varying parameters) in our mortality projections we also resort to the technique of bootstrap (see Efron, 1992).
3.1
Population Data
We will use population data from the Netherlands and similar European countries to estimate age-specific death rates. The AG define similar European countries as countries with a Gross Domestic Product (GDP) per capita above the European average. These countries are: Austria, Belgium, Denmark, Finland, France, Ger-many, Iceland, Ireland, Luxembourg, Norway, Sweden, Switzerland and the United Kingdom. The total number of deaths (Dx,t) and exposures-to-risk (the number of
person years from which Dx,t occurred) (Ex,t) for ages 0 to 90 are obtained for the
3.2
Mortality Model Specifications
Before we will specify the final mortality model that will be investigated, we will briefly discuss the specifications of the Lee-Carter, Li-Lee and Poisson log-bilinear model that form the basis of the AG2016 model.
3.2.1 Lee-Carter model
Lee and Carter can be regarded as the pioneers of stochastic mortality models. In 1992 they proposed a model that is based on standard time series methods to make long-run forecasts of age-specific mortality of the total U.S. population. The Lee-Carter (LC) method is therefore not based on (subjective) expert opinion but instead extrapolates historical trends to forecast age-specific death rates.
Let µx(t) denote the central death rate6 of age x at time t, for x = 0, 1, . . . , w
and t = 0, 1, . . . , T and define αx = T1
PT
t=0log µx(t), i.e. the average of log µx(t)
over time. Lee and Carter aim to find the least-squares solution of
log µx(t) = αx+ βxκt+ εx,t. (13)
Since there are no regressors on the right hand side of the equation, the solution cannot be found by ordinary regression methods. They therefore applied singular-value decomposition on log µx(t) − αx to find estimates for the parameters βx and
κt. ˆβx equals the first left-singular vector and ˆκt the product of the leading value
and the first right vector.7 We can interpret the time-varying index κt as the overall
6The time-varying equivalent of the assumed constant (over the year of age) of µ
xin (10). 7Let Z
x(t) = log µx(t) − αx, then SVD[Zx(t)] =P r
mortality index, and βx as the rate in which mortality rates respond to changes in κt.
αx is then regarded as the general shape of the mortality schedule across age. The
homoskedastic error terms εx,t are normally distributed with mean 0 and variance
σ2ε.
Since the model is underdetermined, as for any constant c 6= 0 the transformation (αx, βx/c, c · κt) or (αx − c · βx, βx, κt+ c) yields the same log µx(t) as (αx, βx, κt),
the parameters are normalised to ensure model identification by P
xβx = 1 and
P
tκt = 0. Before modelling κt as a time series process it should be adjusted such
that it equates the estimated number of deaths to the actual observed deaths. κt
can then be modelled and forecasted as an ARIMA time series model using standard Box-Jenkins methodology (Box and Jenkins, 1976). In the original publication, and in many other applications of the LC model, a random walk with drift,
κt= κt−1+ θ + et, (14)
with (a usually negative) drift term θ and et ∼ N (0, σ2) (an ARIMA(0,1,0) model)
is found to describe κt well.
3.2.2 Li-Lee model
This undesirable result called for a method that identifies a common tendency within a group, but preserves individual deviations in the short or medium term.
Let µx,i(t) denote the death rate of age x at time t for individual population i.
Then, the Li-Lee model is specified by
µx,i(t) = αx,i+ BxKt+ βx,iκt,i+ εx,t,i. (15)
The parameters can be interpreted in a similar fashion as the parameters of the LC model. We can identify a common factor BxKt that describes the change over time
in mortality for all populations in the group and an individual factor βx,iκx,i that
allows for a short- or medium-term deviation of the common trend. αx,i is then an
individual population’s general shape of mortality schedule across age.
Similar to Lee and Carter (1992), the authors use SVD to find the parameter estimates of (15). First, the parameters of the common trend are estimated and in a second step the specific factor of the individual population is obtained such that it describes the individual population’s residual matrix. To ensure the uniqueness of the parameters, the restrictions P
xBx = 1, P tKt = 0 and P xβx,i = 1, P tκt,i = 0
(for all i) are imposed.
The time-varying parameters Kt and κtare again modelled using standard
time-series methods. The authors choose to model Ktas a random walk with drift (RWD).
κt,i is then modelled as a random walk without drift,
with et,i ∼ N (0, σ2i), or a first-order autoregressive process,
κt,i = ci+ aiκt−1,i+ et,i, (17)
with constant ci, |ai| < 1 and et,i ∼ N (0, σi2). The random walk (RW) and first-order
autoregressive process (AR(1)) ensure that κt,i tends towards a constant value on the
long run and therefore that the individual population’s trend does not diverge from the group’s common trend.
3.2.3 Poisson log-bilinear model
Brouhns et al. (2002) propose a Poisson regression approach on stochastic mortality modelling of a single population. The authors argue that the assumption of ho-moskedastic errors in the LC model is a serious drawback because at older ages there are a smaller absolute number of deaths, resulting in a more variable log force of mortality at higher ages than at lower ages. The number of deaths Dx,t at age x and
time t is recognised as a counting random variable and therefore assumed to have a Poisson structure. Specifically, Dx,t is modelled by
Dx,t ∼ Poisson (Ex,tµx(t)) with µx(t) = exp (αx+ βxκt) , (18)
determined by maximising the log-likelihood8
l (α, β, κ) =X
x,t
{Dx,t(αx+ βxκt) − Ex,texp (αx+ βxκt)} + constant. (19)
It is worth mentioning that the authors find similar trends in the parameter estimates of the Poisson log-bilinear model and the classical Lee-Carter model, and that the Poisson model outperforms the LC model for higher ages (especially for ages over 90). The time trend κt is again modelled and forecasted using an ARIMA time series
model.
3.2.4 AG2016 model
The AG2016 model is a stochastic model which distinguishes a common trend be-tween mortality in the Netherlands and other similar European countries. Moreover, it incorporates correlations in mortality rates between men and women. For age x = {0, 1, 2, . . . , 90}, time t = {1970, . . . , 2015} and gender g = {M, F } the force of mortality µgx(t) is modelled by:
log µgx(t) = log µg,EUx (t) + log µg,N Lx (t), with log µg,EUx (t) = Agx+ BxgKtg and
log µg,N Lx (t) = αgx+ βxgκgt.
(20)
For the sake of readability we will occasionally drop the superscripts g, EU and NL. Following Brouhns et al. (2002) it is assumed that the observed number of deaths Dx,t is a counting random variable with a Poisson distribution as in (18). First
8The log-likelihood function of a random sample X = (X
1, X2, . . . , Xn) iid
∼ Poisson(λ) is given by l(λ; x) =Pn
the parameters for the European force of mortality Ag
x, Bxg and K g
t are estimated for
each gender separately such that the Poisson (log-)likelihood of observed deaths is maximised, given the observed number of exposures:
l(A, B, K) = X
x,t
DEU
x,t (Ax+ BxKt) − Ex,tEUexp(Ax+ BxKt) + constant. (21)
For identification, the restrictions P
tKt= 0 and
P
xBx= 1 are imposed. Since
population data was not available for t = 2015, K2015could not be estimated directly
and its estimate is derived using linear extrapolating ˆK2014.9 Second, the parameters
Ax, Bx and Kt are taken as given and the log-likelihood for the Dutch data,
DN Lx,t ∼ Poiss Ex,tN Lµˆx(t)
with ˆµx(t) = exp( ˆAx+ ˆBxKˆt+ αx+ βxκt), (22)
is maximised to find estimates for αx, βx and κt. The parameters αx, βx and κt are
normalised as well. The dynamics of the stochastic processes Kt and κt are then
modelled as
Kt= Kt−1+ θ + εt and
κt= aκt−1+ δt.
(23)
That is, Kt is modelled by a random walk with drift and κt by a first-order
au-toregressive process. The AR(1) process ensures that on the long run, deviations of Dutch mortality from the European trend will diminish. Kt and κt could also be
drift or an AR process that yields a bounded short-term trend. When assuming the above specified dynamics for Ktgand κgt, the stochastic variable Zt = εtM, δMt , εFt, δtF
has a multivariate normal distribution with mean 0 and given covariance matrix C. That is Zt
i.i.d.
∼ N (0, C). The AR(1) and RWD parameters θM, θF, aM, aF and C
are again estimated using maximum likelihood. This last step introduces correlation between men and women in the stochastic mortality model.
Due to limitations in the reliability of mortality data for ages 91 and above the closure method of Kannisto (1994) is used until age 120. The method is based on a logistic regression on the last n = 11 estimated ages yk ∈ {80, 81, . . . , 90}. For ages
x ≥ 91 define log µx(t) = L n X k=1 wk(x)L−1 µyk(t) ! , with L and L−1 the logistic and inverse logistic functions10
L(x) = 1 1 + e−x and L −1 (x) = − log 1 x − 1 ,
and regression weights
wk(x) = 1 n + (yk− ¯y) (x − ¯y) Pk j=1(yj− ¯y) 2 = 1 11+ (yk− 85) (x − 85) 110 .
For ages x > 120 it is assumed that µx(t) = µ120(t).
We can simulate a sample path of the time series Kt+hg and κgt+hby drawing h inde-pendent samples from the multivariate normal distribution of Zt= εMt , δMt , εFt, δtF.
Then the forecasted hazard rate of an individual with gender g and age x at time t + h is ˆµg x(t + h) = exp ˆAgx+ ˆBxgKˆ g t+h+ ˆαgx+ ˆβxgˆκ g t+h .
Assuming the force of mortality does not change within a year, µx+s(t+s) = µx(t)
for 0 ≤ s < 1, the corresponding one-year mortality rates are given by qx(t) =
1−exp−R1
0 µx+s(t + s) ds
= 1−exp (−µx(t)). To find the mortality rates that are
published in the AG2016 projection table, one simply has to set εMt , δMt , εFt, δFt = (0, 0, 0, 0) for all t, i.e. simulation is not required. The mortality rates used by the pension fund are forecasted using an adjustment of the 2016 projection table from the AG. In line with the guidelines proposed by the Dutch Central Bank (DNB) in 2012, the adjustments are age specific and based on past mortality in the pension fund and socioeconomic status. Based on a participant’s postal code, we can pool individuals into different socioeconomic groups. Using past mortality from several pension funds we can determine an optimal correction factor for a given mortality table. These individual (age dependent) correction factors can be combined to find a pension fund-wide correction factor. The model that is used is owned by an external party and therefore will not be discussed in more detail in this thesis. It is worth mentioning that mortality rates of the fund’s participants are, on average, lower than the Dutch population’s.
3.3
The Bootstrap Approach
All uncertainty in the mortality projections of the above mentioned models come from the uncertainty from forecasting the time-varying parameter(s) κt(and Kt). In order
sample B times with replacement from a given dataset to create replicate datasets. Performing statistical analysis on the separate replicated datasets gives a variability of parameter estimates from which we can produce confidence intervals to assess their uncertainty. We will examine two bootstrap methods: a residual bootstrap and a semi-parametric Poisson bootstrap. The two methods will be explained in the upcoming two sections. A third (parametric) approach that was proposed by Brouhns et al. (2002), which samples from the multivariate normal distribution of the relevant parameters, will not be discussed as we are dealing with two sets of parameters (for the European trend and the Dutch deviation from this trend) and their multivariate distribution is not well specified and studied. We refer the interested reader to Renshaw and Haberman (2008) for a comparison of the three approaches to the Poisson Lee-Carter single population model.
3.3.1 Residual Bootstrap
In this section we will explain the residual bootstrap approach by Koissi et al. (2006). The authors propose to sample with replacement from the matrix of deviance residu-als rD of the Poisson log-bilinear model (18) to find B replications rbD. The deviance
residuals are then inverted to find the corresponding matrices of number of deaths Db
x,t (with b = 1, . . . , B).11 Applying the Poisson log-bilinear model to Dbx,t with
the original matrix of exposures Ex,t yields B sets of parameter estimates for αx, βx
and κt. Each ˆκbt is then modelled by an ARIMA(0,1,0) process and the B sets of
corresponding predicted death rates and life expectancy are derived. The percentile
11Find the matrix Db
approach (Efron and Tibshirani, 1994) is then applied to find confidence intervals for the predicted death rates and life expectancy. For the residual bootstrap to perform adequately, it is important that the residuals are independent and identically dis-tributed. It is therefore crucial to thoroughly inspect the deviance residuals before drawing random samples.
3.3.2 Poisson Bootstrap
The second bootstrap approach that we will investigate is the semi-parametric boot-strap by Brouhns et al. (2005). Since the observed number of deaths Dx,tare assumed
to come from a Poisson distribution, Brouhns et al. (2005) suggest to apply a Pois-son noise to Dx,t to produce B replicate datasets. Dbx,t is then a realisation from
the Poisson distribution with mean Dx,t. The corresponding Lee-Carter parameters
are then estimated for each replicate dataset and from there death rates and life expectancy are predicted with the corresponding confidence intervals.
Since we are working with two data sets and one of the quantities we are modelling is the Dutch deviation from a European mortality trend, we have to be cautious with replicating the datasets DEU
x,t and Dx,tN L. Since Dx,tEUalso contains the number of deaths
in the Netherlands, we first split the data from the other European countries and the Netherlands: Drestx,t = DEUx,t − DN L
x,t . Adding a Poisson noise to Drestx,t and Dx,tN L
independently could create an artificial large deviation from the observed values. We therefore choose for the following approach:
and parameters D rest x,t Db,EUx,t , Dx,tN L Db,EUx,t .
This way we ensure that the Poisson noise that is applied to DEU
x,t is allocated
pro-portionally as Db,EUx,t = Dx,tb,rest+ Dx,tb,N L. We then fit the AG2016 mortality model to find B sets of parameter estimates with which we calculate the measures of interest.
3.4
The Simulation Model
In order to investigate the reinsurance contract for the mortality risk in the variable pension agreement, we will combine the AG2016 mortality model with fund-specific corrections and a Monte Carlo simulation model for individual deaths.
3.4.1 Pension Fund Data
For the simulation of individual deaths we will use data provided by the pension fund. We observe N = 27,265 individuals at 30 September 2017 from whom we know the gender, date of birth, accrued capital, annual wage and corresponding part-time factor (ptf). Since our model will start at 1 January 2018, we assume that no-one will die between 30 September and 31 December. All workers (wage > 0) still need to pay the last quarter’s pension premium. The premium percentage (pr%) is age-dependent, increasing by age, and one only pays a premium for wages above the franchise12 (45,000 Euro) and below the fiscal maximum (in 2017: 103,317 Euro).
The remaining premium for 2017 (3 out of 12 months) that has to be paid is then:
premium = 3/12 · ptf · max {min {103,317; wage/ptf } − 45,000; 0} · pr%.
To speed up calculations, we will drop all individuals who are both not active at the pension fund (wage = 0) and have a low accrued capital. We define a capital of 8,000 Euro as the minimum, which is approximately the buy-off limit of pensions in the Netherlands. Moreover, we exclude workers with a low wage (wage/ptf < 45,000) and a low accrued capital.13 Even though we have dropped 8,543 participants we have
only lost 2.5% of total capital. It is reasonable to assume that the individuals that are dropped neither choose for a variable pension agreement nor have a significant effect on the outcomes as their accrued capital is so small. Furthermore, we assume that participants do not decease before retirement and continue working at the same employer until retirement. This assumption can be justified as employees who decease or quit jobs are typically replaced by people with similar characteristics (age and gender) and accrued capital that is transferred to the new pension fund. Moreover, for simplicity all birthdays are on 1 January and all individuals retire on their 67th birthday.
Table 2 provides summary statistics of the individuals at 1 January 2018. 20.8% of our population are females, the average age is 49.4 years, an average working par-ticipant has an annual wage of 67,505 Euro (for privacy reasons we have truncated wages at 105,075 Euro) and an accrued capital of 55,768 Euro. The sum of the re-tirement capital is 1.04 billion Euro, of which 87% is owned by the male participants.
Table 2: Summary Statistics of Pension Fund’s Participants
Mean St. Dev. Min Max
Gender 0.21 - -
-Age 49.4 8.7 24 67
Wage∗ 67,505 19,904 24,494 105,075
Capital 55,768 79,065 36 1,001,826
* Only for 0 < wage < 105,075.
3.4.2 Simulation Setup
At the start of each year t, we have an accrued capital K(t). For the calculations we will divide K(t) into three sets: K1(t) that is reserved for the retiree, K2(t) that
is reserved for his/her spouse and a negative K3(t) on the life of the first death.
K(t) = K1(t) + K2(t) + K3(t) holds. K1(t) will be used to provide a whole life
annuity due on the life of the retiree (K1/¨ax), K2(t) provides a whole life annuity
due on the life of the spouse (K2/¨ay) and K3(t) is used to provide a negative income
stream (K3/¨axy) that offsets K2/¨ay when the retiree is still alive. Note that if (x)
has no spouse (or decides at retirement that he does not want a pension for his spouse) K2 = 0 and K3 = 0. Figures 4 and 5 show how the total pension payment
is composed and what happens when the retiree (x) or the spouse (y) dies.
0 1 2 3 Time ¨ ax ¨ ay ¨ axy 1000 700 -700 1000 700 -700 1000 700 0 -700 0 1000 0 0 Figure 5: Time-line diagram: Spouse (y) dies at t ∈ (1, 2]
The dynamics of the pension capital K1(t) in the retirement phase is as follows:
K1(t + 1) = (K1(t) − C1(t)) (1 + rt)
1 + qx
px
Each year, the next year’s retirement benefit C1(t) = K1(t)/¨axis subtracted from
K1(t) and the rest will be invested in the capital market with return rt. If the retiree
is still alive at the end of the year, he receives a life bonus of (K1(t) − C1(t))(1 + rt)pqxx.
Immediately after acquiring the life bonus, the next year’s pension benefit will be subtracted and paid out. Similar dynamics hold for K2(t) and K3(t).
In the event of death of (x), K1(t)+K3(t) goes to the reinsurer and in the event of
death of (y), the reinsurer receives K2(t) + K3(t). If we have a total of N reinsurance
contracts with individuals i = 1, . . . , N we define14
and
V (t) =X
i
I[Tx,i≤ t]K1,i(t)
+X
i
I[Ty,i≤ t]K2,i(t)
+X
i
I[Txy,i ≤ t]K3,i(t),
where B(t) is the total life bonus that is paid out at time t and V (t) is the to-tal amount of capito-tal that will be received by the reinsurer. Tx,i, Ty,i and Txy,i
are the analogues of Tx, Ty and Txy for pension contract i. The same holds for
K1,i(t), K2,i(t), K3,i(t), and qx,i, qy,i, qxy,i, and px,i, py,i, pxy,i. Note that if (x) deceases
at time t ∈ (t1, t2], then K1(t3) = 0, i.e. his capital is zero after death: the reinsurer
has claimed it.
Individual deaths of (x) and (y) will be simulated by independent Bernoulli ran-dom variables with corresponding parameters qx and qy. Moreover, individuals pay a
pension premium for ages x = 21, . . . , 66 that will be added to their pension capital K(t). Furthermore we assume that during the individual’s working career he also receives a (rather small) life bonus from the pension fund. We will build on the idea that the life bonus at time t ensures that the retirement benefit at time t will exactly equal the retirement benefit at time t − 1 if the retiree would invest in a risk-free asset. That is, C(t) = C(t − 1) should hold if the retiree survives until t.15
We will consider the following three scenario’s. In scenario 1 the annuity factors and life bonuses are based on the projected mortality rates from the AG2016 model
and these rates are also the input of the Bernoulli random death generator. Profits or losses on the reinsurance contract solely come from the uncertainty of actual deaths not corresponding to projected deaths (for a finite population).
In scenario 2 we will use a worst-case scenario projection table as input for the mortality rates, but the annuity factors and life bonuses are based on the point-estimates of the AG2016 projection table (including a fund specific correction). As the annuity factors and life bonuses are actually too high, we know (almost certainly) that the reinsurance will result in a severe loss for the reinsurer. In this scenario there is no possibility for the reinsurer to adjust its rates if he realises he has made a (severe) miscalculation for the projection of death rates. Since this is an unrealistic scenario (especially if the reinsurer is the same pension fund) and the reinsurance is only for the mortality risk (changes in life expectancy are not reinsured), we also consider a third scenario that does allow for future adjustments. Scenario 2 is therefore mainly for illustrative purposes of the significant effects of uncorrected longevity.
For a long period of losses (actual death rates are lower than projected), this will correspond to a decrease in pension payments and bonuses. The opposite holds for a long period of profits.
For all three scenario’s we will assume a constant return on stocks of 6.5%. The returns on bonds are defined by the one-year forward rates derived from the yield curve provided by the DNB on 31-12-2017. The investment mix is based on a life-cycle which invests 83% in stocks until age 45, and decreases its exposure to stocks to 32% at ages 67 and higher. Moreover, the franchise and the maximum pensionable wage in 2018 are 45,000 and 105,075 Euro respectively and increase by 2% per year. Wages increase by 2.5% annually. We also assume that approximately 70% of the retirees will choose for a spouse pension. Furthermore, we will work with gender-neutral annuity factors as described in Section 2.2.4 and gender-gender-neutral life bonuses. The true probabilities of death that are the input for our simulation model are however gender specific. Lastly, we abstract away from any administration costs.
We would like to emphasise that we deem the chosen parameters appropriate in this particular setting. The return on stocks and bonds are in line with the param-eters provided by the Dutch Central Bank16 and the other parameters are typical for analyses for this particular pension fund. In, for example, an extensive asset and liability study one should undoubtedly apply a wider range of parameter values. However, as our primary subject of interest is mortality and time was restricted we consider it appropriate to limit the research on this front. Furthermore we expect that choosing other parameter values will have a limited effect on the outcomes. A
higher return on stocks would result in larger absolute losses or gains but since the total amount of capital in the fund would increase in a similar pattern, we do not expect a substantial relative impact.
4
Results
The following subsections, 4.1 and 4.2, discuss the appropriateness of the residual and Poisson bootstrap for the data at hand with regard to the construction of confidence intervals for the parameters of the AG2016 model. Subsection 2.3 provides the results of the simulation study for the reinsurance contract for the various scenarios.
4.1
Residual Bootstrap
0 20 40 60 80 −20 0 10 20 30
Deviance Residuals EU Males by Age
Age (0−90) Residual 0 20 40 60 80 −5 0 5 10
Deviance Residuals NL Males by Age
Age (0−90)
Residual
Figure 6: Deviance Residuals by Age: Males
0 10 20 30 40 −20 0 10 20 30
Deviance Residuals EU Males by Year
Year (1970−2014) Residual 0 10 20 30 40 −5 0 5 10
Deviance Residuals NL Males by Year
Year (1970−2015)
Residual
age and years, we conduct the Levene’s test for variance (Levene, 1960). The results can be found in Table 3. For all age-groups and years we reject the null-hypothesis that the deviance residuals are homoskedastic. Omitting age 0 from the analysis yields similar results.
Table 3: Levene’s test for equality of variances
W df1 df2 p-value Over ages European Males 13.82 90 4004 <0.001 Dutch Males 5.33 90 4095 <0.001 European Females 10.30 90 4004 <0.001 Dutch Females 3.05 90 4095 <0.001 Over years European Males 3.18 44 4050 <0.001 Dutch Males 2.80 45 4140 <0.001 European Females 2.28 44 4050 <0.001 Dutch Females 3.16 45 4140 <0.001
If p > 0.05, then homoscedasticity holds.
In addition, we conduct a visual inspection of contour plots of the deviance resid-uals. Figure 8 shows the deviance residuals for all ages and years for European females. We can strongly identify non-random trends in the residuals that violate independence of the residuals. Figure 9 shows a similar contour plot for Dutch fe-males. Even though trends are less clear they are still persistent. Similar but less prominent conclusions can be drawn from contour plots for males (Figure 15 and Figure 16 in Appendix B).
−15 −10 −5 0 5 10 15 20 1970 1980 1990 2000 2010 0 20 40 60 80
Deviance Residuals EU Females
Figure 8: Contour Plot of Deviance Residuals of European Females
−4 −2 0 2 4 1970 1980 1990 2000 2010 0 20 40 60 80
Deviance Residuals NL Females
the outlying residuals, as there seems to be a structural deviation for this particular cohort in both European male and female data. Mortality rates are remarkably lower for this cohort compared to its neighbouring cohorts because of the larger number of exposures-to-risk. In the Dutch dataset we find no such pattern. Whether this means that the latter is atypical or that the European dataset is flawed, we refrain from using the residual bootstrap, since the observed residuals do not appear to be the outcomes of i.i.d. random variables.
4.2
Poisson Bootstrap
For the Poisson bootstrap we set the number of bootstrap replications to be B = 10,000. Plots of the parameters, including approximate 90% confidence intervals, are found in Figure 10 (males) and Figure 11 (females). The Poisson noise that we add to the observed number of deaths induces some natural variability in the parame-ters. For the interested reader we provide tables with the parameter estimates, and corresponding 5th and 95th percentiles in Appendix C.
Due to the large absolute number of deaths for ages 60 to 90,17 the Poisson
noise has a rather minor effect on the male age-dependent parameters Ax+ αx and
Bx at high ages. For Bx with x < 60 we do find some variation, but especially for
pension providers and our reinsurance contract, low age mortality is not of particular importance. For Kt with t < 2015 (the estimated years) we cannot identify any
significant variability, indicating that the time trend is rather robust.
and κt. It is especially remarkable that the estimate of the AR(1) coefficient a in
(3.2.4) is, on average, lower than the point estimate of the AG2016 model (0.962 versus 0.980). The slope of the forecasted values of κtis therefore steeper, indicating
that the Dutch mortality trend will converge more rapidly to the European trend than the point estimates suggest. Similar results are obtained for the female parameters. For the remaining cohort life expectancy of a 65-year-old man and woman in 2016 we simulate B mortality tables with which we calculate eb
65. We will compare these
values with simulations based only on the uncertainty of future values of the time-varying parameters (KM
t , KtF, κMt , κFt ). Results are presented in Table 4 and Figure
17 (in Appendix B). We find minimal differences. The medians in the bootstrap setting of both males and females are slightly higher than in the non-bootstrap (regular) setting. Moreover, the distribution of the bootstrapped e65 for females
seem to have slightly fatter tails than its regular counterpart: the difference in the 1st percentile is approximately 1 month (0.08 years). For males we cannot identify significant differences.
As eb
65also includes uncertainty in the forecasting of the time-varying parameters,
we conclude that most randomness originates from forecasting these parameters. Similar results hold for the period life expectancy (not shown here), suggesting that even on the short term (forecasting the time-varying parameters only 1 period ahead) the large observed number of deaths for European mortality data and therefore the relatively small additional Poisson noise dominate the results.
0 20 40 60 80 −8 −6 −4 −2 Index ALPHA 0 20 40 60 80 0.005 0.010 0.015 0.020 Index BET A.EU 0 20 40 60 80 100 −200 −150 −100 −50 0 Index KAPP A.EU 0 20 40 60 80 −0.06 −0.02 0.00 0.02 0.04 0.06 Index BET A.NL 0 20 40 60 80 100 −4 −2 0 2 Index KAPP A.NL
0 20 40 60 80 −8 −6 −4 −2 Index ALPHA.F 0 20 40 60 80 0.010 0.015 0.020 Index BET A.EU.F 0 20 40 60 80 100 −150 −100 −50 0 50 Index KAPP A.EU.F 0 20 40 60 80 −0.02 −0.01 0.00 0.01 0.02 0.03 Index BET A.NL.F 0 20 40 60 80 100 −10 −5 0 5 10 15 Index KAPP A.NL.F
Figure 11: Bootstrapped Parameter Estimates Females with 90% confidence interval (dashed lines). Top Ax+ αx, mid-left Bx, mid-right Kt, bottom-left βx and
Table 4: Cohort Life Expectancy of a 65-year-old in 2016 1% 2.5% 5% 50% 95% 97.5% 99% Males eb65 18.79 19.02 19.19 20.07 20.89 21.07 21.26 e65 18.75 18.98 19.14 20.01 20.90 21.07 21.25 Females eb 65 22.15 22.32 22.45 23.18 23.85 23.97 24.11 e65 22.23 22.37 22.50 23.14 23.76 23.90 24.02
years ahead to obtain the corresponding mortality rates18) we identify no differences
between the bootstrap and regular forecast.
For the worst-case scenario projection table that will be used in the Monte Carlo simulation of the reinsurance contract, we will only incorporate the variability in the time-varying parameters Kt and κt of the AG2016 model. We simulate the
parameters 53 years ahead to produce 10,000 projection tables and we select the scenario that results in the 97.5% highest cohort life expectancy of a 67-year-old man and woman in 2018. In this scenario, the male’s remaining life expectancy is 19.62 years and for the female it is 22.27 years. As a reference, the point estimates of the life expectancy of a 67-year-old male and female are 18.47 years and 21.48 years respectively.
4.3
Reinsurance Contract
The three scenario’s are simulated 50,000 times. For every simulation and every time t = 2018, . . . , 2062 we observe the profit or loss on the reinsurance contract and the total amount of capital of the retirees. To evaluate the impact of a future profit or loss, we need to evaluate it relative to the total amount of capital in the whole fund at that point in time. A one million Euro loss relative to a total of 10 million or 1 billion evidently differ significantly. Moreover, as the participants still pay a pension premium until retirement, the total amount of capital will grow in the first years of the simulation. For the total DC pension fund we start with a total amount of capital Ktot
t of around 1 billion Euro in 2018 and it reaches its maximum in 2041
amount of capital that is available at the moment of retirement for each cohort and discount this using the DNB yield curve at 31-12-2017, its value Kt∗ is approximately 2,650 million Euro. Tables 5, 6 and 7 in Appendix A show summary statistics of the present values of the future cash flows up to a selected number of years in scenario 1, 2 and 3 respectively. Graphical comparisons are presented in Figures 18 to 20 in Appendix B.
For scenario 1 we find that the mean annual cash flow is approximately zero for all t. For the years 2018 and 2019 we find that around 75% of the simulations have resulted in a loss. Since the number of participants who retire in 2018 and 2019 is extremely small (respectively 2 and 28) there is a high probability that no-one deceases in these two years and therefore more bonus is paid out than the reinsurer receives from the deceased capital. As more participants retire (262, 349 and 445 in 2020, 2021 and 2022 respectively) the distribution of the annual cash flows becomes more symmetric. When compared to the sum of capital that the retirees have at the moment of retirement (Kt∗), we observe that the 2.5th percentile fluctuates around the value 0.005. This indicates that, for this given finite population, the probability that the reinsurance contract results in a loss larger than 0.5% of total capital is approximately 0.025. For the 0.5th percentile we find a relative loss fluctuating around 0.65%. The 25th and 75th percentile are approximately -0.15% and 0.15% respectively for all t.
the randomness of the survival of a person it is possible to experience some profits (especially in the early stage). We find that for the first 4 years this scenario does not significantly differ from scenario 1. The rise in life expectancy that follows from the simulated projection table presents itself in the fifth year as the average cash flow relative to total capital diverts from 0. After 8 years, at the end of 2025, we observe that the 75th percentile drops below 0. The mean loss rises to approximately 3.75% of total capital Kt∗ in the year 2062 (-99.1 million relative to 2.650 million).
In scenario 3, we adjust our mortality rates when subsequent losses on the rein-surance contract are observed. In the first two years we only have 30 retirees and we cannot draw any well founded conclusions from the results. For the period from 2020 to 2032 we find that the average cash flow diverts from zero, however not as rapidly as in scenario 2 since in some simulations we have already adjusted mortality rates. In 2032 we observe that the average simulation has resulted in a total loss of 3 million (approximately 0.3% of K2032∗ ). The absolute average loss increases in the subsequent years, but since Kt∗ increases at a larger rate, the loss relative to Kt∗ decreases. In the year 2062 we find that the average present value of the loss is 4.6 million, which is about 0.18% of K2062∗ . The 2.5th and 0.5th percentile are -19.5 (0.74%) and -25.1 (0.95%) million respectively.
simulation results of scenario 3 (approximately 95% of the simulations lie inside this confidence interval). This suggests that for this particular finite population the ef-fect of a worst-case scenario mortality table is minimal compared to the randomness manifested in individual deaths.
The previous results are the sum of the discounted future cash flows compared to the total amount of capital that the participants have at the moment of retirement. One could also be interested in the one-year profit or losses on the contract. In scenario 1 we find the 1st percentile of the annual cash flow to be less than 0.1% of total capital for the years 2018 until 2032. Thereafter, the 1st percentile of the one-year loss increases to 0.2% until 2049. For the years after 2049 we have a larger retired population with life bonuses that increase each year (qx increases and thus
qx/px increases) and total capital Kttot that decreases. Hence, the spread increases
while its expected value stays constant. In scenario 3 we find similar results as for the years 2020 to 2028 less than 1% of the simulations result in a one-year loss of 0.1% or more of total capital. For the years 2029 until 2049 we find that the 1st percentile does not exceed 0.2% of total capital.
If we translate this to the results of our reinsurance simulation, we find that the 2.5th percentile is a loss of 13.1 million Euro. The average loss (also known as Expected Shortfall, or Conditional Value-at-Risk) of the scenario’s in the 2.5th percentile is 15.6 million Euro. This would translate in a capital buffer of 15.6 million euro times the probability of the event 2.5%: 390,000 Euro. It is customary to charge this amount to the customer (here: participant in the variable pension agreement) and in order to relate this to a single risk premium that needs to be paid at the moment of retirement we find that the participant should pay a premium of 0.0147% of total capital to cover this buffer: an almost negligible amount.
5
Discussion and Conclusion
We have studied the AG2016 mortality model and incorporated a bootstrap approach to examine the effect of parameter uncertainty on forecasted life expectancies. In the majority of the model’s parameters we have found no significant effect of incor-porating all variability in the parameters. Moreover, the randomness that originates from the time-varying parameters dominates the randomness from the other param-eters. As a consequence, quantities that are derived from the model, such as life expectancy, do not significantly differ between the bootstrap and regular approach. One result that stands out is that the bootstrapped parameter estimates for aM and
aF differ remarkably from the regular point estimate. This indicates that the Dutch
existent.
For the reinsurance contract of mortality risk in the variable pension agreement we find that an actuarially fair price is nil. For this particular pension fund’s pop-ulation we find that a single risk premium that covers the capital reserve that is necessary to cover the possibility of an adverse scenario is minimal (0.0147% of one’s accrued capital at the moment of retirement). In addition, as the effect of an adverse future mortality scenario is minimal, this buffer would also capture a large amount of possible outcomes in our third scenario.
We would like to emphasise that we have used a particular pension fund’s pop-ulation with a given size. Performing a similar analysis with a significantly larger population would have resulted in smaller capital buffers, relative to the total amount of capital, as deviations from the actual mortality rate would be less common. But then again, a reinsurance contract would not be necessary as mortality risk could be shared collectively. Alternatively, if we had just a small number of retirees in the variable pension agreement, the inherent uncertainty of a statistical analysis would render it less than useful.
residual bootstrap. In order to investigate whether this particular cohort influences the outcomes significantly, one could for instance exclude it from the analysis or include a specific cohort dummy variable in the model.
Bibliography
Blake, David and William Burrows (2001). Survivor bonds: Helping to hedge mor-tality risk. Journal of Rik and Insurance 68 (2), 339–348.
Box, George EP and Gwilym M Jenkins (1976). Time series analysis. Forecasting and control. In Holden-Day Series in Time Series Analysis, Revised ed., San Francisco: Holden-Day, 1976.
Brouhns, Natacha, Michel Denuit, and Ingrid Van Keilegom (2005). Bootstrapping the Poisson log-bilinear model for mortality forecasting. Scandinavian Actuarial Journal 2005 (3), 212–224.
Brouhns, Natacha, Michel Denuit, and Jeroen K Vermunt (2002). A Poisson log-bilinear regression approach to the construction of projected lifetables. Insurance: Mathematics and economics 31 (3), 373–393.
Brouhns, Natacha, Michel Denuit, Jeroen K Vermunt, et al. (2002). Measuring the longevity risk in mortality projections. Bulletin of the Swiss Association of Actuaries 2, 105–130.
Cairns, Andrew JG, David Blake, and Kevin Dowd (2006). A two-factor model for stochastic mortality with parameter uncertainty: Theory and calibration. Journal of Risk and Insurance 73 (4), 687–718.
Cairns, Andrew JG, David Blake, Kevin Dowd, Guy D Coughlan, David Epstein, Alen Ong, and Igor Balevich (2009). A quantitative comparison of stochastic mortality models using data from England and Wales and the United States. North American Actuarial Journal 13 (1), 1–35.
Cairns, Andrew JG and Ghali El Boukfaoui (2017). Basis risk in index based longevity hedges: A guide for longevity hedgers. Technical report.
Currie, ID (2006). Smoothing and forecasting mortality rates with P-splines. Talk given at the Institute of Actuaries.
De Nederlansche Bank (2012). Good Practice Gebruik Fondsspecifieke Ervaring-sterfte. http://www.toezicht.dnb.nl/binaries/50-226664.pdf.
math-Efron, Bradley (1992). Bootstrap methods: another look at the jackknife. In Break-throughs in statistics, pp. 569–593. Springer.
Efron, Bradley and Robert J Tibshirani (1994). An introduction to the bootstrap. CRC press.
Kannisto, V¨ain¨o (1994). Development of oldest-old mortality 1950-1990: Evidence from 28 developed countries.
Koissi, Marie-Claire, Arnold F Shapiro, and G¨oran H¨ogn¨as (2006). Evaluating and extending the Lee–Carter model for mortality forecasting: Bootstrap confidence interval. Insurance: Mathematics and Economics 38 (1), 1–20.
Lee, Ronald D and Lawrence R Carter (1992). Modeling and forecasting US mortal-ity. Journal of the American statistical association 87 (419), 659–671.
Levene, Howard (1960). Robust tests for equality of variances. Contributions to probability and statistics 1, 278–292.
Li, Nan and Ronald Lee (2005). Coherent mortality forecasts for a group of popula-tions: An extension of the Lee-Carter method. Demography 42 (3), 575–594. Renshaw, AE and Steve Haberman (2008). On simulation-based approaches to risk
measurement in mortality with specific reference to Poisson Lee–Carter modelling. Insurance: Mathematics and Economics 42 (2), 797–816.
Renshaw, Arthur E and Steven Haberman (2006). A cohort-based extension to the Lee–Carter model for mortality reduction factors. Insurance: Mathematics and Economics 38 (3), 556–570.
Royal Actuarial Association (2016). Projection table AG2016. https://www.ag-ai. nl/view.php?Pagina_Id=731.
Sociaal Economische Raad (2016). Persoonlijk pensioenvermogen met col-lectieve risicodeling. https://www.ser.nl/~/media/db_adviezen/2010_2019/ 2016/persoonlijk-pensioenvermogen.ashx.
