An extension to a prominent Dutch mortality model
Geert Johannes Zittersteyn
Abstract
This thesis provides a complete explanation of an extension to the prominent and fre-quently applied mortality model from Actuarieel Genootschap (2018). We recreated the model in a cause-specific environment using Archimedean copulae to model dependence between various groups of causes of death. For this, Dutch data on cause-of-death mortal-ity and cause-specific mortalmortal-ity data from 14 comparable European countries is used. We estimate the cause-specific model and create future mortality forecasts using the model. Using the forecasts, we investigate the effects of an unexpected event in respiratory disease mortality. Moreover, we compare our model properties to other (cause-specific) mortality models. We find that the inclusion of a common factor to a cause-specific mortality context increases robustness of the forecast and we underline that cause-specific mortality forecasts foresee a more pessimistic mortality future than general mortality models.
For further inquiries, contact g.zittersteyn@gmail.com. We strongly recommend to print this paper in color.
Contents
1 Introduction 5 1.1 Preface. . . 6 1.2 Literature review . . . 7 1.3 Research objectives . . . 10 2 Methodology 12 2.1 Model Characterisation . . . 13 2.1.1 Lee-Carter Model . . . 13 2.1.2 Li-Lee Model . . . 14 2.2 Cause-of-death mortality . . . 17 2.2.1 Crude mortality . . . 18 2.2.2 Net mortality . . . 18 2.3 Competing Risk. . . 19 2.3.1 General framework . . . 192.3.2 Archimedean survivor copula . . . 20
2.4 Model Estimation. . . 22
2.4.1 Maximum Likelihood Estimation . . . 22
2.4.2 Newton-Raphson Method . . . 24
2.5 Time Dependent Component . . . 25
2.6 Old Ages . . . 26
2.7 Procedure Overview . . . 26
2.8 Population Dynamics . . . 27
3 Data 30 3.1 Data Description . . . 31
3.1.1 WHO Mortality Database Code Description . . . 31
3.1.2 Data difficulties . . . 33
3.1.4 Data completion . . . 37
3.2 Population Dynamics Data Description . . . 39
4 Results 42 4.1 Initial mortality calculations . . . 43
4.1.1 Crude mortality . . . 43
4.1.2 Marginal survival functions . . . 49
4.1.3 Net mortality . . . 49
4.2 Parameter estimates . . . 50
4.2.1 Maximum likelihood estimates . . . 50
Chapter 1
Introduction
1.1
Preface
The German phenomenologist Martin Heidegger describes in his book Sein und Zeit1 how
humans are driven by their fear of death. Behind his concept of human sentiment towards death lies the perception of randomness of human mortality. Although opinions may differ on the truthfulness of the philosopher’s first statement, the second point of view we share with him. In this research we will approach this randomness of mortality from an Econo-metric perspective and try to predict death using stochastic modelling.
Knowledge on how to predict future mortality is valuable for many reasons. Obviously, from a human standpoint we are interested in how long we are going to live, but maybe an even more important reason is that the whole social financial security system revolves around the prediction of future mortality. With the knowledge of future mortality rates and subsequently future life expectancy, financial institutions can calculate the value of mortality based annuities such as pensions and other life insurances. Such institutions adapt their products in order to accommodate longevity risk. Evaluation of future mor-tality is needed for governmental policy and health care as well.
Since 1950, the period life expectancy of Dutch males at age 65 has increased from 14.39 years to 18.97 in 20172. Female life expectancy at 65 years old experienced a change from
14.94 to 21.49. We observe that significant progress has been made in terms of mortal-ity reduction. Even more remarkable, for males at their retirement age (assumed 65) life expectancy has increased with approximately 2 years over the past decade. According to
Actuarieel Genootschap (2018), life expectancy is only expected to rise further to a pre-dicted level of around 22 and 25 years in 2050 for males and females, respectively.
The most prominent methods used for forecasting mortality in the Netherlands, such as the one performed in the research highlighted in the previous paragraph, are based on overall mortality. The models attempt estimation of the general mortality trend over time. The information provided earlier gave us a first insight into this trend for the Netherlands. The prediction of the continuation of this trend per age is based on long term movement and or short term shocks of mortality. Every additional year, new information on mortality becomes available which can be used to update mortality forecasts.
1
Heidegger and von Herrmann(1977)
1.2 Literature review 7
In the past fifty years, the mortality shocks have not been extremely drastic. Never-theless, one could imagine that for several hypothetical cases, mortality could be heavily altered. Such eventualities are difficult to account for beforehand, but could potentially severely affect our current predictions of future mortality rates. Among the hypothetical scenarios mentioned above are those that regard increases or decreases for certain diseases or, more general, causes of death. Possible examples could be a deadly viral epidemic or vice versa a solution to a major deadly disease. The effects are likely to be sizable for such large scale events for particular causes of death. Simultaneously, these scenarios are highly unlikely. In contrast, smaller shocks for a specific cause of death are more probable to occur and also affect mortality predictions, although to a lesser extent. A recent example of that is the observation made byActuarieel Genootschap(2018). For the first time their model experienced a decline of life expectancy in the future with respect to their previous model. This deviation was partially ascribed to a large increase of deaths caused by the influenza virus in recent years for people of older ages.
Such sudden big and small developments for cause-of-death mortality are not unimag-inable and their effect on future mortality could potentially be enormous. Furthermore, we have addressed earlier that many governmental social policies and products offered by financial institutions depend on longevity. The accumulation of these established facts leads us to believe that deeper research into the effects of cause-specific mortality is im-portant and could contribute to a more complete understanding of mortality. Moreover, a model that incorporates the added effects of cause of death mortality is capable of a better prediction of future mortality in case of unexpected events regarding specific causes of death.
1.2
Literature review
In this section we will provide a review of past literature on the subject of mortality fore-casting. We will address a selection of research papers describing related mortality models, some of which we have based our own model on.
The first model we will discuss is the one first introduced in Li and Lee (2005). This model is an extension of the first stochastic mortality model devised by Lee and Carter
factor which is an additional variable translating the shared effects of other comparable populations. Through this method, data is increased and parameter estimates become less sensitive to incidental deviation of individual mortality. In addition, the undesirable property of forecast divergence for individual populations is countered. The model is fully explained in section 2.1.2. For its forecast of Dutch mortality, Actuarieel Genootschap
(2018, 2016) uses the Li-Lee model using European mortality data and Dutch mortality data from 1970 to present. The European information will serve as the common factor in this case.
Enchev et al. (2017) provide a full comprehensive analysis of several mortality models. All four models approach mortality from a multi-population perspective. The considered models are the ordinary Li-Lee model, two simplified versions of this model and the common age effect model ofKleinow (2015). First, the authors introduce the model formulations and the fitting method. The authors implement a maximum likelihood estimation method in order to estimate the models. They have chosen this method in order to be able to use the BIC value. By evaluating the BIC they can appropriately compare the fit of the mod-els. Different variable constraints are implemented in order to circumvent the identifiability problems of the models. Furthermore, compared to the Li-Lee model different reestimation methods are used in order to capture the time dependent variables. They have chosen for a method that allows for coherent converging mortality forecasts. The conclusion drawn from their results is that the common age effect model is the highest-ranking model for their data set in terms of fit. The regular Li-Lee model was the second-best performing model. In addition, the authors addressed potential robustness issues of the Li-Lee model.
1.2 Literature review 9
We have seen that the common factor approach is not the only method to incorportate in-terdependence of mortality in related populations into mortality forecasting models. Where the Li-Lee model is considered an appropriate method for modeling the long-run risk ef-fects for mortality, the assumption of mortality convergence may be too restricting in the short term mortality dependence between co-integrated populations. By using a factor copula,Chen et al.(2015) provide a different approach to model the dependence between population. This straightforward copula is used to model the dependence of mortality improvements for six developed countries.
Utilizing a copula to model the dependence of related risk is not new in the context of mortality modeling. In their paper, Zheng and Klein (1995) go over the dependence of competing risks. Competing risks is the presence of censoring of the time of death from one cause in the event of death from another cause. They prove that under an assumed copula function, competing risk mortality data is sufficient to derive the marginal survival distributions of the separate risks.
Other authors also have devoted papers on cause-of-death mortality and often different methods are used to capture their dependence. Boumezoued et al. (2018) assumes in-dependence of cause-specific mortality in his research into the effects of a reduction of cause-specific mortality for the demographic composition in France. Arnold and Sherris
of competing risk effects for cause-removal analysis, leading to more rational results. The last model we address is that developed by Li and Lu (2019). In this paper the authors create a cause-specific mortality model using hierarchical and non-hierarchical Archimedean copulae. The underlying model formulation is the Lee-Carter single popu-lation model. Through the copula, the dependence between different cause-specific mor-talities is captured. In other words, the underlying marginal survival functions per cause of death is derived from the joint survival function. Hereafter, in the research the crude cause-specific mortality is forecasted using a Lee-Carter single population model on net mortality. Moreover, they evaluate the effect on life expectancy caused by impulses on cause-specific mortality. Their research underlines the importance of evaluating model un-certainty, since they show that differences between at first glance similar models can be significant under certain circumstances.
1.3
Research objectives
In this thesis we will extend the Dutch mortality model from Actuarieel Genootschap
1.3 Research objectives 11
The main research question is:
What are the benefits of the common factor cause-specific mortality model?
For this purpose we need to answer the following subquestions:
• Can we construct a sufficient data set to accommodate the model?
• How do the mortality forecasts from the common factor cause-specific mortality model differ from other cause-specific mortality models and general mortality models? • What implications do these differences have?
• What applicational advantages result from the model?
Methodology
This chapter offers a detailed description of the model, formulas and definitions we have used for the quantitative research of this thesis. In section 1.2 we provide the relevant background of our model. We want to replicate the model used byActuarieel Genootschap
2.1 Model Characterisation 13
2.1
Model Characterisation
2.1.1 Lee-Carter Model
As mentioned before, there exist many ways of modeling mortality. Lee and Carter(1992) first introduced a fully stochastic extrapolative method, applied on U.S. mortality data. Without the use of economic, behavioral and medical factors to predict future mortality rates, this method relies solely on persistent patterns in historical mortality levels. The Lee-Carter model (henceforth the LC model) is considered the "leading statistical model for mortality" (Deaton and Paxson,2004).
Before we proceed, we define the observed mortality for population 𝑖: 𝑚𝑖(𝑥, 𝑡) =
𝐷𝑖(𝑥, 𝑡) 𝐸𝑖(𝑥, 𝑡)
. (2.1.1)
𝐷𝑖(𝑥, 𝑡)denotes the amount of deaths observed in year t for people age 𝑥. 𝐸𝑖(𝑥, 𝑡)describes the corresponding age group exposure, i.e. population size per age. For the observed matrix containing the historical mortality rates 𝑚𝑖(𝑥, 𝑡) for age 𝑥 and period 𝑡 = [1, 𝑇 ], the LC model is described as follows:
𝑙𝑛(𝑚𝑖(𝑥, 𝑡)) = 𝑎𝑖(𝑥) + 𝑏𝑖(𝑥)𝑘𝑖(𝑡) + 𝑒𝑖(𝑥, 𝑡). (2.1.2) In the equation, the coefficient 𝑎𝑖(𝑥) denotes the general age-specific mortality constant. The responsiveness of the mortality rate of a certain age group 𝑥 to changes of time de-pendent mortality development is given by the coefficient 𝑏𝑖(𝑥) or, in other words, the coefficient 𝑏𝑖(𝑥) denotes the sensitivity across age groups to variation in the time-varying parameter 𝑘𝑖(𝑡). On the basis that 𝑒𝑖(𝑥, 𝑡) represents the errors of the logarithmic trans-formation of the death rates, the variance of the age-specific element 𝑒𝑖(𝑥, 𝑡) will not vary much over time (Lee and Carter,1992), although in a specific period the death rates may differ severely per age. Therefore the error term is assumed to be independent and iden-tically distributed with mean 0 and variance 𝜎2
𝑒𝑖(𝑥,𝑡). This is also a required assumption
needed for their applied singular value decomposition (SVD) estimation method, which we will address later on in this section. This assumption has been contested in other literature.
Due to the fact that the model is underdetermined (the number of unknown parameters exceed the number of equations), the authors impose two constraints:
∑︁ 𝑥 𝑏𝑖(𝑥) = 1and ∑︁ 𝑡 𝑘𝑖(𝑡) = 0, (2.1.3)
which in turn directly implies that 𝑎𝑖(𝑥)equals the average of log mortality rates over time. They then use SVD procedure on [𝑙𝑛(𝑚𝑖(𝑥, 𝑡)) − 𝑎𝑖(𝑥)] for estimation of the parameters 𝑏𝑖(𝑥)and 𝑘𝑖(𝑡). After this procedure, ˆ𝑘𝑖(𝑡)is determined such that the errors of the loga-rithm of the death rates are minimized rather than the actual deaths. The authors suggest the reestimation of ˆ𝑘𝑖(𝑡) with a time series model while assuming the estimates ˆ𝑏𝑖(𝑥) and ˆ
𝑎𝑖(𝑥) given.
In their research, Lee and Carter (1992) proved the autoregressive model of order 1 (AR(1) model) to be the best performing model for the reestimation of ˆ𝑘𝑖(𝑡). Nevertheless, they proposed the autoregressive integrated moving average (0,1,0) model (ARIMA(0,1,0) model), because of its simplicity. In addition, this model only underperformed the AR(1) model slightly. The ARIMA(0,1,0), or more commonly referred to as the random walk with drift (RWD), is defined as follows:
ˆ
𝑘𝑖(𝑡) = ˆ𝑘𝑖(𝑡 − 1) + 𝜌 + 𝑒(𝑡). (2.1.4)
In this equation we have 𝑒(𝑡) ∼ 𝒩 (0, 𝜎2
𝑒(𝑡)) and 𝜌 denotes the constant drift coefficient which is straightforwardly estimated by ˆ𝜌 = (ˆ𝑘𝑖(𝑇 ) − ˆ𝑘𝑖(1))/(𝑇 − 1), as 𝑡 = [1, 𝑇 ] is the range of years used for the estimation.
2.1.2 Li-Lee Model
2.1 Model Characterisation 15
(a) male (b) female
Figure 2.1: Period life expectancy at birth of selected EU countries (Actuarieel Genootschap,2016)
Figure 2.2: Development of life expectancy (LE) from birth on worldwide scale from 1950-2015 (United Nations, Department of Economic and Social Affairs, Population Division,
2017)
Moreover, research shows a more rapid decrease of mortality levels in populations with a lower life expectancy in comparison to the on average longer living populations (Oeppen and Vaupel, 2002). In other words, in countries with a relatively low expected lifespan, life expectancy rises more quickly. This has been pointed out in research conducted by
(2002) have shown that this need not be true. This leads to the general assumption that mortality trends between populations do not diverge over time. This can be rationalized by exchange of knowledge improving life expectancy among populations and by the pre-sumption that mortality reduction is more easily accomplished by countries with low life expectancy.
From a perspective of ever-increasing interconnection between countries and regions caused by globalization,Li and Lee(2005) have extended the original model by incorporating the data of closely connected and socioeconomically similar populations. This addition aug-ments the model by countering the inherent forecast divergence of the LC model discussed in the previous paragraph. The Li-Lee model (LL model) first determines the collective main trend of the group of populations as a whole by applying the so-called common fac-tor method. Next, they analyse the unique dissimilarities of an individual population on a short or medium horizon. This is done in order to model the historical deviation of the single population with respect to the common tendencies of the entire group. The com-bination of both stages, results in a mortality model that asserts a non-diverging relation between populations in the long run while retaining explanatory power of short/medium term variation of mortality for individual populations. The LL model is currently used by
Actuarieel Genootschapfor the calculation of the Dutch life tables.
For the first stage, we refer back to equation 2.1.2 of the LC expression. It is required that all populations within the group to have identical 𝑏𝑖(𝑥) and 𝑘𝑖(𝑡), denoted as 𝐵(𝑥) and 𝐾(𝑡) respectively. Furthermore, each population 𝑖 has an age specific effect 𝐴𝑖(𝑥) which is allowed to vary among the group members and are thus estimated individually. This is done because this variable is a constant and therefore will not lead to divergence. The ensuing common factor equation is equivalent to the LC model with the replaced variables:
𝑙𝑛(𝑚𝐶𝐹𝑖 (𝑥, 𝑡)) = 𝐴𝑖(𝑥) + 𝐵(𝑥)𝐾(𝑡) + 𝛾𝑖(𝑥, 𝑡), 𝛾𝑖(𝑥, 𝑡) i.i.d.
∼ (0, 𝜎2
𝛾𝑖(𝑥,𝑡)). (2.1.5)
Which, similarly to the LC model, can be estimated using SVD after having imposed constraints such as those described in equation 2.1.3. Enchev et al. (2017) describe the above mentioned constraints in a LC context, applied without exception for every variation on the model. Again, 𝐾(𝑡) is modeled as an RWD (2.1.4):
2.2 Cause-of-death mortality 17
In addition to the common factor,Li and Leepinpoint the individual mortality effects for population 𝑖. To put it differently, we need to model the residuals from the common factor equation (2.1.5) for population 𝑖 using the age dependent variable 𝑏𝑖(𝑥) and the time-varying term 𝑘𝑖(𝑡). The resulting extended common factor model is specified as follows:
𝑙𝑛(𝑚𝑖(𝑥, 𝑡)) = 𝑎𝑖(𝑥) + 𝐵(𝑥)𝐾(𝑡) + 𝑏𝑖(𝑥)𝑘𝑖(𝑡) + 𝛿𝑖(𝑥, 𝑡). (2.1.7) The authors suggest a stationary AR(1) model or an AR(1) process with unit root (random walk):
𝑘𝑖(𝑡) = 𝑐𝑖+ 𝜔𝑖𝑘𝑖(𝑡 − 1) + 𝜖𝑖(𝑡), (2.1.8) where 𝑐𝑖 is a constant, 𝜔𝑖 ≤ 1and 𝜖𝑖(𝑡) ∼ 𝒩 (0, 𝜎𝜖2𝑖(𝑡)). Enchev et al. (2017) highlight the potential shortcomings in the reestimation process of the time-dependent coefficient 𝑘𝑖(𝑡). They suggest a vector autoregressive model of order 1 (VAR(1)) instead of the AR(1) pro-posed byLi and Lee(2005). This is due to the sometimes diverging properties of the AR(1) model in a multi-population context. By the use of the VAR(1) model and its underlying covariance matrix, coherent forecasting can be retained. Moreover, this forecasting method does not deviate much from individual AR(1) forecasts (Enchev et al.,2017).
The LL model will be the model formulation used in our research. We apply this model to cause-specific net mortality. More on this will follow in the next sections. For the common factor estimation EU data will be used. The second part will be performed on mortality information for the Netherlands. The composition of this data is explained in chapter 3. In the subsequent section, we will discuss how we approach the cause-specific data and we elaborate on the involved formulas and definitions.
2.2
Cause-of-death mortality
equal to the real underlying mortality of the specific disease. This definition of mortality by taking into account other potential causes of death is called the net mortality. In the next part of this section we will provide mathematically how we derive these mortality variables.
2.2.1 Crude mortality
We regard both EU and and NL data as homogeneous populations (population is denoted as 𝑖) with independent and identical distributed lifetimes 𝑋𝑖 =min(𝑋𝑖1, . . . , 𝑋𝑖𝑧), where 𝑋𝑖𝑐 is the age of death due to cause 𝑐. Then let (𝑋𝑖, 𝐶) represent the time of death and the corresponding cause of death, 𝐶 ∈ {1, . . . , 𝑧}. Under the assumption of independence between causes of death, we find in literature (Li and Lu,2019;Boumezoued et al.,2018) the crude cause-of-death mortality for age 𝑥 defined as:
𝜇𝑐(𝑥) = lim 𝑢→0
P[𝑋𝑖 ≤ 𝑥 + 𝑢, 𝐶 = 𝑐 | 𝑋𝑖 > 𝑥]
𝑢 . (2.2.1)
In our case, we are dealing with a discrete data set. We therefore simply define the observed crude cause-of-death mortality the following way. Age group of age 𝑥 of population 𝑖 in year 𝑡 has size 𝐸𝑖(𝑥, 𝑡). The same year this group experiences 𝐷𝑐𝑖(𝑥, 𝑡) deaths from cause 𝑐. Then mortality is defined as:
𝑚𝑐𝑖(𝑥, 𝑡) = 𝐷 𝑐 𝑖(𝑥, 𝑡) 𝐸𝑖(𝑥, 𝑡)
. (2.2.2)
All independent mortality rates cumulatively describes the general mortality of the popu-lation: 𝑚𝑖(𝑥, 𝑡) = 𝑧 ∑︁ 𝑐=1 𝑚𝑐𝑖(𝑥, 𝑡), (2.2.3)
with the corresponding overall survival function: 𝑆𝑖(𝑥, 𝑡) =exp(− 𝑥 ∑︁ 𝑘=0 𝑧 ∑︁ 𝑐=1 𝑚𝑐𝑖(𝑘, 𝑡)). (2.2.4) 2.2.2 Net mortality
For the determination of the net mortality, we turn to the joint distribution of 𝑋𝑖. We follow the specification as given by Li and Lu (2019). The joint survival probability is given by:
2.3 Competing Risk 19
Then, for cause of death 𝑐 we define the net cause-specific mortality of 𝑋𝑐
𝑖 in year 𝑡 and age 𝑥 as: 𝜆𝑐(𝑥, 𝑡) = lim 𝑢→0 P[𝑋𝑐≤ 𝑥 + 𝑢, 𝐶 = 𝑐 | 𝑋𝑐> 𝑥] 𝑢 (2.2.6)
and accordingly we subsequently obtain the net survival function of 𝑋𝑐
𝑖 in a discrete con-text. The net survival function describes the probability of having survived a specific cause of death at age 𝑥 in period 𝑡 as follows:
𝑆𝑐(𝑥, 𝑡) =exp(− 𝑥 ∑︁ 𝑘=0 𝜆𝑐(𝑘, 𝑡)) =P[𝑋𝑐> 𝑥, 𝑋{𝐶}−𝑐> 0]. (2.2.7)
2.3
Competing Risk
As we determined in the literature review, many preceiding researchers have studied depen-dence of cause-specific mortality. The results convince us that it is relevant to incorporate the possible dependence of varying mortality categories. Hence, we will expand on the im-plemented method to assess the dependence between the separate cause-specific mortality distributions. In our explanation of the procedure, we follow the steps made by Li and Lu (2019). Tsiatis (1975) shows that observed data is insufficient for the identification of the latent dependence structure of cause-of-death mortality and therefore additional con-straints are needed. The use of survivor copula is a prevalent method applied to address the dependence. The idea of the copula, a multivariate distribution function with uniform marginal probability, was first presented inSklar(1959). The method offers a way to ex-amine the dependence of a joint distribution and its underlying marginal distributions of latent random variables.
2.3.1 General framework
For the study of cause-related mortality, we turn to the survivor copula. The survivor copula reveals the dependence of the implicit cause-specific survival functions, denoted as (𝑆1(𝑋1), . . . , 𝑆𝑧*(𝑋𝑧*)). We define the cumulative distribution function as follows:
determination of cause ‘other causes’ will be explained in the next chapter. Reversely, from equation (2.3.1) we obtain:
P[𝑋1 > 𝑥1, . . . , 𝑋𝑧*
> 𝑥𝑧*] = S(𝑆1(𝑥1), . . . , 𝑆𝑧*(𝑥𝑧*)), (2.3.2) with positive 𝑥𝑐 and 𝑐 = 𝑐*. We have determined an expression for the joint distribu-tion of the lifespan which depends on the specificadistribu-tion of joint distribudistribu-tion funcdistribu-tion S and the net survival function 𝑆𝑐(𝑋𝑐). Since the cause-specific survival functions have uniform marginals, S qualifies by definition as a copula and is henceforth referred to as the survivor copula.
Li and Lu (2019) consider the approach performed in Carriere (1994) to solve the sys-tem straight from data using the nonlinear partial differential equations. In their paper they indicate a major drawback of this approach, which is the lack of an analytical solution. In conformance with this argumentation we also opt for a closed-form expression. We will go into this process more deeply in the next section.
2.3.2 Archimedean survivor copula
First, we assume a joint Archimedean survivor copula applies for the random survival variables (𝑋1, . . . , 𝑋𝑧*
). We can now rewrite equation (2.3.2) using the formulation of the Archimedean copula introduced byLing(1965):
P[𝑋1 > 𝑥1, . . . , 𝑋𝑧* > 𝑥𝑧*] = S(𝑆1(𝑋 1), . . . , 𝑆𝑧 * (𝑋𝑧*)) = 𝜓(𝜓−1(𝑆1(𝑥1)), . . . , 𝜓−1(𝑆𝑧*(𝑥𝑧*))), (2.3.3) with positive 𝑥𝑐, 𝑐 = 𝑐*. The expression 𝜓(·) implies the copula generator function which satisfies the necessary and sufficient condition of d-monotonicity (McNeil et al.,2009). We will elaborate more on the copula generator function shortly.
2.4 Competing Risk 21
In these equations, 𝜃 denotes the dependence coefficient. For 𝜃 = 0 the copulas above equate to the independence copula. When an independence copula is assumed, net mor-tality value equals crude mormor-tality. The full definition of the survival functions as well as the first derivative and inverse of both corresponding generator functions are given in Appendix A, equations (1) - (6). These functions will be implemented later on in this section.
Under the assumption of a known pre-specified Archimedean copula we can proceed with the direct calculation of the net survival functions from the observed data. First, we derive the crude cause-specific mortality from the data. We plug these into the discrete formu-lation of the net survival function established by Wang (2012); Rivest and Wells (2001). We obtain the following expression of the marginal survival distribution for population 𝑖, cause 𝑐, age 𝑥 and year 𝑡:
𝑆𝑖𝑐(𝑥, 𝑡) = 𝜓 [︃ − 𝑥 ∑︁ 𝑠=0 exp(︁−∑︀𝑠 𝑢=0 ∑︀𝑧* 𝑘=1𝑚𝑘𝑖(𝑢, 𝑡) )︁ 𝜓′∘ 𝜓−1∘exp(︁−∑︀𝑠 𝑢=0 ∑︀𝑧* 𝑘=1𝑚𝑘𝑖(𝑢, 𝑡) )︁ · 𝑚 𝑐 𝑖(𝑠, 𝑡) ]︃ , (2.3.6)
where ‘∘’ is the notation used for expressing the composition of functions. Conversely, for age 𝑥 and for known 𝑆𝑖(𝑥, 𝑡) (through equation (2.2.4)) this equation becomes:
𝑚𝑐𝑖(𝑥, 𝑡) = −𝜓−1(𝑆𝑐 𝑖(𝑥, 𝑡)) − ∑︀𝑥−1 𝑠=1 [︁ 𝑆 𝑖(𝑠,𝑡) 𝜓′∘𝜓−1∘𝑆 𝑖(𝑠,𝑡)· 𝑚 𝑐 𝑖(𝑠, 𝑡) ]︁ [︁ 𝑆 𝑖(𝑥,𝑡) 𝜓′∘𝜓−1∘𝑆 𝑖(𝑥,𝑡) ]︁ . (2.3.7)
The cause-specific mortality at age 0 can be derived as: 𝑚𝑐𝑖(0, 𝑡) = −𝜓 −1(𝑆𝑐 𝑖(0, 𝑡)) [︁ 𝑆 𝑖(0,𝑡) 𝜓′∘𝜓−1∘𝑆 𝑖(0,𝑡) ]︁ . (2.3.8)
A full comprehensive overview of the steps taken to derive equation (2.3.7) in combination with2.3.8 is given in AppendixA, equation (7). Now that the marginal survival functions are known, we can readily determine the net mortality by solving:
𝜆𝑐𝑖(𝑥, 𝑡) = −ln(︁ 𝑆 𝑐 𝑖(𝑥, 𝑡) 𝑆𝑐 𝑖(𝑥 − 1, 𝑡) )︁ . (2.3.9)
2.4
Model Estimation
In addition to the rudimental singular value decomposition (SVD) fitting method sug-gested byLee and Carter(1992), other fitting methods have been proposed in the past. In particular cases, the procedure suggested byLee and Carterhas some undesirable features which could possibly complicate estimation. Since the model was designed for the mortal-ity forecasting of the total U.S. population, which has sufficient information on death rates for all ages, the fitting method proposed rarely yields sub-optimal coefficient estimates. In case of less rich information sets, such as the case when analyzing cause-of-death data, modification of the fitting method is needed to counter the complications that arise for death rates equal to zero (Wilmoth,1993). Similarly, the SVD approach shows disadvan-tages in case of data sets containing age groups with low death rates and other age groups with high death rates (Girosi and King,2007). In addition, the use of useful measures as AIC and BIC is not accomodated by SVD estimation. In his technical report, Wilmoth
describes two alternatives for fitting in the first stage (i.e. the stage before reestimation) of the LC model: the weighted least square (WLS) and maximum likelihood estimation (MLE) technique. Due to the fact that the WLS technique is proven to be not statistically robust (Girosi and King,2007), we will only discuss the MLE method here. On the same grounds will we use MLE fitting for our model.
2.4.1 Maximum Likelihood Estimation
The approach we use in this paper for fitting our model is that applied by Actuarieel Genootschap(2016) and by Enchev et al.(2017) for their mortality model, which follows an estimation method described byBrouhns et al.(2002).
Rather than a linear specification of the log-mortality, the distribution of deaths can be defined in a Poisson framework. We denote 𝑑𝑐
𝑖(𝑥, 𝑡) as the actually observed number of deaths for population 𝑖 by cause of death 𝑐 and let 𝐷𝑐
𝑖(𝑥, 𝑡)be the random variable for the amount of deaths at time 𝑡 for age 𝑥. Similar toBrillinger et al. (1986), we assume
2.4 Model Estimation 23
FollowingBrillinger et al.(1986), we derive the log-likelihood function from (2.4.2): ℓ𝑐𝑖(𝑥, 𝑡) =∑︁ 𝑥 ∑︁ 𝑡 [︀𝐷𝑐 𝑖(𝑥, 𝑡) ·ln{𝜆𝑐𝑖(𝑥, 𝑡)𝐸𝑖(𝑥, 𝑡)} − 𝜆𝑖𝑐(𝑥, 𝑡)𝐸𝑖(𝑥, 𝑡) −ln{𝐷𝑖𝑐(𝑥, 𝑡)!}. (2.4.3)
The maximum likelihood estimation for our model consist of 2 stages, on which we will elaborate hereafter.
Step 1:
We return to the Li-Lee formulation. First we focus on the common factor expression described in equation (2.1.5) for which MLE is to be performed. In our case this is repre-sented by the data for EU. The restriction imposed by this equation and the substitution of 𝑚𝐶𝐹
𝑖 (𝑥, 𝑡) by 𝜆𝑐𝑒𝑢(𝑥, 𝑡)leads to the following equation:
𝜆𝑐𝑒𝑢(𝑥, 𝑡) = 𝑒𝐴𝑐𝑒𝑢(𝑥)+𝐵𝑐𝑒𝑢(𝑥)𝐾𝑒𝑢𝑐 (𝑡). (2.4.4)
After inserting (2.4.4) into the log-likelihood equation from (2.4.3) yields: ℓ𝑐𝑒𝑢(𝑥, 𝑡) =∑︁ 𝑥 ∑︁ 𝑡 [︀𝐷𝑐 𝑒𝑢(𝑥, 𝑡) ·ln{𝐸𝑒𝑢(𝑥, 𝑡) · 𝑒𝐴 𝑐 𝑒𝑢(𝑥)+𝐵𝑒𝑢𝑐 (𝑥)𝐾𝑒𝑢𝑐 (𝑡)} − 𝐸𝑒𝑢(𝑥, 𝑡) · 𝑒𝐴𝑐𝑒𝑢(𝑥)+𝐵𝑒𝑢𝑐 (𝑥)𝐾𝑒𝑢𝑐 (𝑡)−ln{𝐷𝑐 𝑒𝑢(𝑥, 𝑡)!} ]︀ =∑︁ 𝑥 ∑︁ 𝑡 [︀𝐷𝑐 𝑒𝑢(𝑥, 𝑡) · (ln{𝐸𝑒𝑢(𝑥, 𝑡)} + 𝐴𝑐𝑒𝑢(𝑥) + 𝐵𝑒𝑢𝑐 (𝑥)𝐾𝑒𝑢𝑐 (𝑡)) − 𝐸𝑒𝑢(𝑥, 𝑡) · 𝑒𝐴 𝑐 𝑒𝑢(𝑥)+𝐵𝑒𝑢𝑐 (𝑥)𝐾𝑒𝑢𝑐 (𝑡)−ln{𝐷𝑐 𝑒𝑢(𝑥, 𝑡)!}]︀. (2.4.5)
The parameter estimates ˆ𝐴𝑐𝑒𝑢(𝑥), ˆ𝐵𝑒𝑢𝑐 (𝑥) and ˆ𝐾𝑒𝑢𝑐 (𝑡) are found by maximizing the log-likelihood function in equation (2.4.5). For this we use the Newton-Raphson method (Raphson,1702;Newton,1833), which we will describe in more detail in the next section. For MLE, similarly to the SVD method applied in both the LC and LL model, the restric-tions from equation (2.1.3) are applied, forcing the sum of ˆ𝐾𝑒𝑢𝑐 (𝑡)and ˆ𝐵𝑒𝑢𝑐 (𝑥) to be equal to 1 and 0, respectively.
Step 2:
Next, we turn to the unique tendencies of the individual population with respect to the common factor effect of the group. Then, we state the restriction as in equation (2.1.7) by inserting the optimal values for 𝐴𝑐
following manner:
𝜆𝑐𝑛𝑙(𝑥, 𝑡) = 𝑒𝐴^𝑐𝑒𝑢(𝑥)+ ^𝐵𝑒𝑢(𝑥) ^𝐾𝑒𝑢𝑐 (𝑡)+𝑎𝑐𝑛𝑙(𝑥)+𝑏𝑐𝑛𝑙(𝑥)𝑘𝑐𝑛𝑙(𝑡)
= ˆ𝜆𝑐𝑒𝑢(𝑥, 𝑡) · 𝑒𝑎𝑐𝑛𝑙(𝑥)+𝑏𝑐𝑛𝑙(𝑥)𝑘𝑐𝑛𝑙(𝑡).
(2.4.6) Comparable to step 1, we now plug in (2.4.6) into equation (2.4.3) and obtain the subse-quent statement: ℓ𝑐𝑛𝑙(𝑥, 𝑡) =∑︁ 𝑥 ∑︁ 𝑡 [︀𝐷𝑐 𝑛𝑙(𝑥, 𝑡) ·ln{ˆ𝜆𝑐𝑒𝑢(𝑥, 𝑡)𝐸𝑛𝑙(𝑥, 𝑡) · 𝑒𝑎 𝑐 𝑛𝑙(𝑥)+𝑏 𝑐 𝑛𝑙(𝑥)𝑘 𝑐 𝑛𝑙(𝑡)} − ˆ𝜆𝑐𝑒𝑢(𝑥, 𝑡)𝐸𝑛𝑙(𝑥, 𝑡) · 𝑒𝑎 𝑐 𝑛𝑙(𝑥)+𝑏𝑐𝑛𝑙(𝑥)𝑘𝑐𝑛𝑙(𝑡)−ln{𝐷𝑐 𝑛𝑙(𝑥, 𝑡)!} ]︀ =∑︁ 𝑥 ∑︁ 𝑡 [︀𝐷𝑐 𝑛𝑙(𝑥, 𝑡)(ln{ˆ𝜆𝑐𝑒𝑢(𝑥, 𝑡)𝐸𝑛𝑙(𝑥, 𝑡)} + 𝑎𝑐𝑛𝑙(𝑥) + 𝑏𝑐𝑛𝑙(𝑥)𝑘𝑛𝑙𝑐(𝑡))) − ˆ𝜆𝑐𝑒𝑢(𝑥, 𝑡)𝐸𝑛𝑙(𝑥, 𝑡) · 𝑒𝑎 𝑐 𝑛𝑙(𝑥)+𝑏 𝑐 𝑛𝑙(𝑥)𝑘 𝑐 𝑛𝑙(𝑡)−ln{𝐷𝑐 𝑛𝑙(𝑥, 𝑡)!}]︀. (2.4.7)
The Newton-Raphson method is used in order to attain the maximum of the log-likelihood function from (2.4.7). Furthermore, following the guidelines set from (2.1.3), we set the sum of 𝑏𝑛𝑙(𝑥) and 𝑘𝑛𝑙(𝑡) equal to 1 and 0, respectively. As a result, we obtain the fitted estimates of the involved parameters, denoted as ˆ𝑎𝑐
𝑛𝑙(𝑥), ˆ𝑏𝑐𝑛𝑙(𝑥) and ˆ𝑘𝑐𝑛𝑙(𝑡).
The parameter estimations are performed for both genders separately. We denote the resulting parameter estimates ˆ𝜃𝑐,𝑔 ∈ { ˆ𝐴𝑐,𝑔
𝑒𝑢(𝑥), ˆ𝐵𝑒𝑢𝑐,𝑔(𝑥), ˆ𝐾𝑒𝑢𝑐,𝑔(𝑡), ˆ𝑎𝑐,𝑔𝑛𝑙(𝑥), ˆ𝑏𝑐,𝑔𝑛𝑙 (𝑥), ˆ𝑘𝑛𝑙𝑐,𝑔(𝑥)}, for 𝑔 = {male, female} and all causes 𝑐.
2.4.2 Newton-Raphson Method
As mentioned earlier, we use the Newton-Raphson method (Newton,1833;Raphson,1702) as the procedure to obtain the optimal estimates for the parameters of our mortality model. The Newton-Raphson method is a numerical method, which through an iterative process approximates the roots of a function. Then for some stopping criterion and conditional on its previous estimated parameters for every iteration 𝜈, we estimate 𝜃𝑐,𝑔 recursively by estimating each of its elements successively, using the formula:
ˆ 𝜃(𝜈+1)𝑐,𝑔 = ˆ𝜃𝑐,𝑔(𝜈)− (︃ 𝜕(ℓ𝑐,𝑔(𝜈)) 𝜕(𝜃𝑐,𝑔) ⧸︃ 𝜕2(ℓ𝑐,𝑔 (𝜈)) 𝜕(𝜃𝑐,𝑔)2 )︃ . (2.4.8)
2.6 Time Dependent Component 25
and the parameters fitted. For the first step described in the previous section, the partial first and second order derivatives of function (2.4.5) can be found in AppendixA, equations (8) to (13). The first and second order derivatives, found in equations (14) to (19), are used for the second step of the estimation of the fitted parameters.
2.5
Time Dependent Component
The goal of the model estimation presented is forcasting future mortality, using the meth-ods discussed previously. In this section we will focus on the time dependent variables of our model, denoted by 𝐾𝑐,𝑔
𝑖 (𝑡)and 𝑘 𝑐,𝑔
𝑖 (𝑡). Particularized in the literature discussed prior, reestimation of these parameters is used in order to more accurately characterize the time dependent aspect of the model. Another benefit of this second stage estimation is the ability to simulate and forecast the model using the time series formalization.
In the previous section, we have fitted the parameters using MLE, including those de-scribing the two time dependent variables. The resulting subset of estimates ˆ𝜃𝑐,𝑔, 𝐾𝑐 = { ˆ𝐾𝑒𝑢𝑐,𝑀(𝑡), ˆ𝑘𝑛𝑙𝑐,𝑀(𝑡), ˆ𝐾𝑒𝑢𝑐,𝐹(𝑡), ˆ𝑘𝑛𝑙𝑐,𝐹(𝑡)}, we reestimate using time series analysis on the follow-ing equations usfollow-ing an ARIMA(p,d,q) process. Li and Lee(2005) suggest an ARIMA(0,1,0) and an ARIMA(1,0,0) model (see equations2.1.6and 2.1.8) for 𝐾𝑐,𝑔
𝑖 (𝑡)and 𝑘 𝑐,𝑔
𝑖 (𝑡) respec-tively. We will use the same time series specification for this research. The resulting system of time series equations for 𝐾𝑐 are estimated simultaneously using the method of Seemingly Unrelated Regressions (SUR) (Zellner,1962). This method allows information acquired from one equation to be used by the estimation of another, by assuming the error terms to be correlated. This possibly increases efficiency and can be supported by the idea that shocks in mortality for a specific age in a specification year for male and female will probably be strongly correlated. Lastly, we let the set of errors {𝜂𝑐,𝑀
𝑒𝑢 (𝑡), 𝜖𝑐,𝑀𝑛𝑙 (𝑡), 𝜂𝑒𝑢𝑐,𝐹(𝑡), 𝜖𝑐,𝐹𝑛𝑙 (𝑡)}follow a four-dimensional normal distribution with mean 0 and covariance matrix 𝐻𝑐.
From the estimation of the system we obtain the set of estimates {ˆ𝜌𝑐,𝑀, ˆ𝜔𝑐,𝑀, ˆ𝜌𝑐,𝐹, ˆ𝜔𝑐,𝐹} and the corresponding covariance matrix ˆ𝐻𝑐. With the use of these estimates, we can now forecast the time dependent variable set 𝐾𝑐. For this we assume time series error to be none. Combining this with the earlier obtained estimates from ˆ𝜃𝑐,𝑔 makes forecasting 𝜆𝑐𝑛𝑙(𝑥, 𝑡)possible. From these forecasts, we in turn can derive the crude mortality forecasts for the Netherlands 𝑚𝑐
2.6
Old Ages
The data set we use includes maximum age 90, whereas the eventual data set we opt to use includes ages up to 120. Therefore, similar to Actuarieel Genootschap (2018, 2016), we use the theory of Kannisto (1994) as the extrapolation method for the ages over the maximum age available in the data. We want to estimate the mortality for the ages 𝑥𝑜 = {90, . . . , 120}. This estimation is based on mortality of ages 𝑥𝑦 = {80, . . . , 90}. This is executed by performing a logistic regression on the 𝑛 = 11 observations 𝑥𝑦:
𝑚𝑐(𝑥𝑜, 𝑡) = 𝑓(︁ 𝑛 ∑︁
𝑘=1
𝑤(𝑥𝑜, 𝑘) · 𝑓−1(︀𝑚𝑐(𝑥𝑦𝑘, 𝑡))︀)︁, (2.6.1)
where 𝑓(·) and 𝑓−1(·) are defined as:
𝑓 (𝛼) = (1 + 𝑒−𝛼)−1, (2.6.2)
𝑓−1(𝛼) = −ln(︂ 1 𝛼
)︂
− 1. (2.6.3)
The weights 𝑤(𝑥, 𝑖) are defined in the following fashion: 𝑤(𝑥, 𝑖) = 1 𝑛+ (𝑥𝑦𝑖 − ¯𝑥𝑦)(𝑥 − ¯𝑥𝑦) ∑︀𝑖 𝑗=1(𝑥 𝑦 𝑗 − ¯𝑥𝑦)2 (2.6.4) = 1 11+ (𝑥𝑦𝑖 − 85)(𝑥 − 85) 110 . (2.6.5)
2.7
Procedure Overview
In this section we will provide a summary of the steps taken in the research. For a more detailed description we used reference numbers to the sections in this chapter. Note that all steps (except the reestimation step) are done for both male and female separately.
1. Crude mortality: For cause 𝑐 and population 𝑖, we obtain the crude cause-of-death mortality 𝑚𝑐
𝑖(𝑥, 𝑡)for age 𝑥 in year 𝑡 using equation (2.2.2).
2. Marginal survival: We obtain the marginal cause-specific survival function 𝑆𝑐 𝑖(𝑥, 𝑡) through equation (2.3.6). Henceforth, all steps are preformed for both Frank as well as Clayton copula specification.
3. Net mortality: The net mortality intensities 𝜆𝑐
2.8 Population Dynamics 27
by performing step 1 from section 2.4.1 for the obtained marginal mortality values of EU.
5. Population effect estimation: Using the results of the previous step and the marginal mortality from NL data, we estimate the parameters explaining the indi-vidual population mortality effects by means of step 2 (section 2.4.1).
6. Reestimation: By means of the steps described in section 2.5, we reestimate the time dependent set of variables for male and female simultaneously.
7. Forecast: We forecast the net mortality 45 years using the time series analysis from the previous step.
8. Inverse transformation: We implement equation (2.3.9) reversely and afterwards apply equation (2.3.7) to acquire crude mortality forecasts. For this we use forecast general mortality (see equation (2.2.3)) using the LL model.
9. Logistic extrapolation: Extrapolate the crude mortality forecasts using the theory of Kannisto (1994) for old age mortality.
The resulting set of information are crude mortality forecasts from 2016 to 2060, for ages 0 to 120 and all causes of death included in the data. The results are for male and female separately and we differentiate between results under the assumption of either a Frank or Clayon survivor copula.
2.8
Population Dynamics
• 𝑀𝑔(𝑥, 𝑡): net migration for males and females (𝑔) for age 𝑥 (from 0 to 99 years old) in year 𝑡. We assume migration for ages 100 to 120 to be zero.
• 𝐵𝑔(𝑡): total life births for males and females (𝑔) in year 𝑡. • 𝑚𝑔
𝑛𝑙(𝑥, 𝑡): the total crude mortality intensity for Dutch males and females (𝑔) age 𝑥 in year 𝑡. This variable is obtain through equation (2.2.3). We have assumed mortality for individuals older than 120 in year 𝑡 to equal that of a 120 year old.
• 𝐸𝑔(𝑥, 2016): the exposure for males and females on the first day of 2016, age 𝑥 (from 0 to 99 years old). We assume the first period exposure for ages above 99 to equal zero.
Before we proceed, we transform the crude mortality intensities into the probability of death 𝑞𝑔
𝑛𝑙(𝑥, 𝑡) for age 𝑥 in year 𝑡:
𝑞𝑛𝑙𝑔(𝑥, 𝑡) = 1 − 𝑒−𝑚𝑔𝑛𝑙(𝑥,𝑡). (2.8.1)
We will give a mathematical formulation of the iterative modeling procedure below for 𝑡 = (2016, . . . , 2059): (1) 𝐷𝑔(𝑥, 𝑡) =(︂ 𝑞 𝑔 𝑛𝑙(𝑥, 𝑡) + 𝑞 𝑔 𝑛𝑙(𝑥 + 1, 𝑡) 2 )︂ · 𝐸𝑔(𝑥, 𝑡), (2) 𝐷˜𝑔(𝑥, 𝑡) = ⌈︃ ⌊2𝐷𝑔(𝑥, 𝑡)⌋ 2 ⌉︃ , (3) 𝐸𝑔(𝑥, 𝑡 + 1) =max{0, 𝐸𝑔(𝑥 − 1, 𝑡) − ˜𝐷𝑔(𝑥, 𝑡) + 𝑀𝑔(𝑥, 𝑡)}, for 𝑥 ̸= 0, (4) 𝐸𝑔(0, 𝑡 + 1) = 𝐵𝑔(𝑡). (2.8.2)
After having simulated the Dutch population and its corresponding composition, we will compare the model outcomes using three different variables. Similar toBoumezoued et al.
(2018) we denote the ratio of individuals aged 65 and older and individuals aged between 15 and 65, commonly referred to as the dependency ratio:
𝑑𝑡=
𝐸([65, 120], 𝑡)
𝐸([15, 64], 𝑡) . (2.8.3)
2.8 Population Dynamics 29
longevity risk into its calculation of life expectancy for time 𝑡. This method incorporates the mortality risks of forecasted years past 𝑡 in order to calculate the corresponding ex-pected lifespan of an individual. For an extensive explanation and analysis of the two different measures we refer toAyuso et al.(2018).
Data
3.1 Data Description 31
3.1
Data Description
The data set containing cause-of-death mortality data for Europe used for the common factor part of the competing risk LL model is obtained from the World Health Organization (WHO) Mortality Database1. This data set contains mortality data in the form of yearly
number of deaths and mid-year population levels, ranging from the years 1950 to 2017. Most countries in the world have joined WHO and for these countries data on mortality is bundled together in the raw data files of the data base. In order to accommodate further research of the common factor model incorporating competing risk, we report in this section the steps taken to obtain the data set we used for the model estimation in a detailed manner. The data set for the Dutch cause-of-death mortality and mid-year population, from now on referred to as NL, is acquired from StatLine of Centraal Bureau voor de Statistiek (CBS)2.
The same cause-of-death specifications as for the WHO Mortality Database (see section
3.1.1) hold. We have chosen to use the CBS data instead of the WHO data since, unlike the WHO data set, this data set does not involve all shortcomings described in section
3.1.2.
3.1.1 WHO Mortality Database Code Description
Following the example set by the Actuarieel Genootschap, we have chosen 14 European countries showing similar mortality trends (Figure2.1). In line withActuarieel Genootschap, we have included data from 1970 to the latest year available for all countries, namely 2015. The countries have been chosen based on their GDP. The countries are, in alphabetical order: Austria, Belgium, Denmark, Germany, Finland, France, Iceland, Ireland, Luxem-bourg, The Netherlands, Norway, Sweden, Switzerland and the United Kingdom. Keeping in mind the in many cases slim volume of cause-of-death data, we have chosen to include mortality numbers from Germany between 1970 and 1989, when the current country of Germany was split into the Federal Republic of Germany (West Germany/FRG) and the German Democratic Republic (East Germany/GDR). This is in contrast to the general mortality data used by Actuarieel Genootschap. Moreover, in preceding cause-of-death research byArnold and Sherris (2015) mortality data of split Germany was not included, since no data is available for the German Democratic Republic before the year of 1969. This does not pose a problem in our research, because the first year of our data set is 1970 and therefore we accumulated the mortality knowledge of the German subsections to represent total German mortality for the years preceding the fall of the Berlin Wall. In
1
https://www.who.int/healthinfo/statistics/mortality_rawdata/en/
table3.1we have presented the corresponding country codes as given in the documentation handbook provided by WHO.
Austria 4010 Iceland 4160 Belgium 4020 Ireland 4170 Denmark 4050 Luxembourg 4190 Finland 4070 Netherlands 4210 France 4080 Norway 4220 Germany 4085 Sweden 4290 GDR 4090 Switzerland 4300 FRG 4100 United Kingdom 4308
Table 3.1: Countries included and the corresponding codes.
For the purpose of consistency of the cause-of-death indication between countries, the mor-tality database provided by WHO employs cause classification according to the guidelines of International Classification of Diseases (ICD), provided by WHO since 1948. The most recent revision of the classification guidelines, ICD-10, was endorsed by the Forty-third World Health Assembly in 1990 and is now used by all the European countries we have included for our research. For these countries, four different classifications have been used over the period of 1970 to 2017. These are classifications: ICD-7, ICD-8, ICD-9 and ICD-10. Since for all the different countries these classification guidelines were adopted in different years, we must for now process and refine each country individually. In table3.2
we have presented the last year before the next classification method was adopted for each country within the time period 1970-2015.
3.1 Data Description 33
Country ICD-7 ICD-8 ICD-9 ICD-10
Austria - 1979 2001 2015 Belgium - 1978 1997 2015 Denmark - 1993 - 2015 Finland - 1986 1995 2015 France - 1978 1999 2015 Germany - - 1990-1997 2015 GDR - 1978 1989 -FRG - 1978 1989 -Iceland 1970 1980 1995 2015 Ireland - 1978 2006 2015 Luxembourg 1970 1978 1997 2015 Netherlands - 1978 1995 2015 Norway - 1985 1995 2015 Sweden - 1986 1996 2015 Switzerland - 1994 - 2015 United Kingdom - 1978 1999 2015
Table 3.2: Final year of implementation of the ICD codes per country for the total cause-specific data set 1970-2015.
Code Causes of death ICD-7 ICD-8 ICD-9 ICD-10
1 Circulatory system A079-A086 A080-A088 B25-B30 I00-I99
2 Cancer A044-A060 A045-A061 B08-B17 C00-D48
3 Respiratory system A087-A097 A089-A096 B31-B32 J00-J99
4 External causes A138-A150 A138-A150 B47-B56 V00-Y89
5 Infectious and parasitic diseases A001-A043 A080-A088 B01-B07 A00-B99
Total A000 A000 B00 AAA
Table 3.3: Cause-of-death coding of the different classification systems.
3.1.2 Data difficulties
Comparability
In table 3.2 we pinpointed the transition moments from one classification to another. As a consequence of adopting new methods of recording mortality data, the data of the successive classifications are analogous. In figuresB.1a toB.4d3 we find that at the points
of transition for some countries and some causes of death the data is not comparable, which appears as a steep drop or a steep rise at the according year. These disturbances are countered via the same method applied byGaille and Sherris (2011), namely with the application of comparability ratios. The comparability ratio is acquired by the condition that the average cause-of-death specific mortality in the last two years before the transition should equal the mean cause-of-death specific mortality of the two years following the transition. Suppose the first after the transition is 𝑡 = 𝐴 and the last year before the transition is 𝑡 = 𝐵, then mathematically:
𝜇𝐵+ 𝜇𝐵−1
2 = 𝑐 ·
𝜇𝐴+ 𝜇𝐴+1
2 , (3.1.1)
where 𝑐 denotes the comparability ratio. After having calculated the comparability ratio, we multiply cause-of-death specific mortality for all years that fall under the new system. We apply this for every country and every cause of death to maintain continuity over dif-ferent classifications. The results are found in figuresB.1a to B.4d3.
The simplicity of this method is advantageous, but the transformation described in equa-tion ()3.1.1) includes a considerable downside. Figures B.1a to B.4d3 show that in most
cases, in particular for the transition to ICD-10, reported mortality lies significantly lower than is the case for the preceding phase. This, in turn, leads to a high corresponding comparability ratio. As a direct result, any deviation in the phase that is less rich in infor-mation is enhanced, increasing the volatility of mortality in this time period. Next section we discuss the cases in which this aspect becomes problematic and we bring forward a solution.
Incomplete and insufficient data
This section we will discuss the defects inherent with the raw data set as provide by WHO. For the countries GDR, Switzerland and United Kingdom the data for years 1979, 1994 and 2000 are missing, respectively. The missing data is most likely caused by the tran-sition from one classification system to another which has occurred in all three specific
3.1 Data Description 35
cases. In order to cope with the absence of information for these specific countries in these specific years, we have assumed the number of deaths to be the same as the year before. The sometimes non-conform death levels between two classification periods, as discussed in the section on comparability, is the argumentation behind not choosing to use linear interpolation.
Furthermore, we previously observed that within our data set for some countries there exists a lack of sufficient information on some cause-of-death specific mortality within the last classification scheme, ICD-10. Generally, this issue is overcome with the use of the earlier discussed comparability ratio. Nevertheless, in some cases, there still exists too little data on mortality and in others no information at all. The former could result in a too large comparability ratio leading to an unrealistically high mortality volatility and the latter even leads to the absence of a comparability ratio. For the instances involving no information, we opted for either linear interpolation between the two enclosing years that contained sufficient information or linear extrapolation using the gradient as the preceding years. For the occurrences with too little information, we added the last year with suffi-cient information as a basis to the years with the shortage in data. In table3.4we created an overview of these instances.
Table 3.4: Details of manipulated instances for insufficient data.
Country Cause of death Period Method*
Austria 4 2002 - 2015 Basis
Iceland 4 1996 - 2009 Interpolation
Ireland 4 2007 - 2015 Extrapolation
Luxembourg 4 1998 - 2015 Basis
*The three methods used to substitude the lack of data are: linear interpolation between surrounding years,
linear extrapolation based on preceiding years or add the level of deaths from the last year with sufficient information to the insufficient years.
Figure 3.1: Cause-specific mortality intensities of the 14 selected countries in the data set after coping with the inherent data difficulties (1970-2015).
3.1.3 Age groups
The WHO data set provides an extensive amount of information on cause-of-death mortal-ity. The number of deaths and population volumes per age are reported in varying ways. For certain years and certain countries the used age specific formatting distinguishes be-tween the age-groups in a more detailed fashion than another country in the same year or the same country in a different year. For the whole WHO data file there are ten different age formatting methods used, however for the 14 selected European countries only three of these apply (see tableB.1), namely format 00, format 01 and format 02.
Format 02 has population and death levels divided in five-year age-groups with the ex-ception of age group 0-4, which is split up in the age 0 and the age-group 1-4. The last group, constituting the highest ages, contains ages 85 and above. This amounts to a total of 19 separate age groups for format 02. Format 01 offers more detailed information about the age group 1-4 by splitting this group into four groups with a one year interval. Lastly, format 00 is an expansion of format 01 due to the addition of age groups 85-89 and 90-94. The final group now becomes ages 95 and above, resulting in a total of 24 separate age groups.
3.1 Data Description 37
the randomness of the formatting within countries and years and the ensuing complexity because of it we have opted for an overall format 02 by bundling age groups 1, 2, 3 and 4 for format 00 and format 01 data points and by bundling age groups 85-89, 90-94 and 95+ for format 00 data points.
For our research, we aim to use mortality on a age interval of one year. Since the data sets at hand, both NL and EU, contain less detailed information due to the 5-year age groups, we use interpolation in order to expand the number of data points. We interpolate following a similar method as the one provided by the HMD Method Protocol byWilmoth et al.(2007). We take a lower boundary 𝑙 = 0 and an upper boundary 𝑢 = 91. This means that we assume that all death in age group 85+ occur before age 91. Next, we create a variable representing the cumulative distribution of deaths per year up to age 𝑥:
𝐶(𝑥, 𝑡) = 𝑥−1 ∑︁
𝑘=0
𝐷(𝑘, 𝑡), where 𝑥 ∈ [𝑙, 𝑢]. (3.1.2)
We apply the cubic spline interpolation method as described byMcNeil et al.(2011), which employs the ‘Hyman filter’ (Dougherty et al.,1989;Hyman,1983) to ensure monotonicity and therefore implicitly non-negative deaths. This yields a cumulative distribution of deaths on a one-year interval from age 0 to 91. We proceed by applying the following formula to obtain the number of deaths per age 𝑥:
ˆ
𝐷(𝑥, 𝑡) = ˆ𝐶(𝑥 + 1, 𝑡) − ˆ𝐶(𝑥, 𝑡), where 𝑥 ∈ [𝑙, 𝑢 − 1]. (3.1.3) The resulting data set ˆ𝐷(𝑥, 𝑡) ranges from ages 0 to 90 for the years 1970-2015.
3.1.4 Data completion
negative number of deaths under cause Other. Suppose in year 𝑡 this is the case for age 𝑖. We construct the new assumed mortality as follows:
𝐷6(𝑥, 𝑡) = ⎧ ⎪ ⎨ ⎪ ⎩ 𝐷6(0,𝑡) 𝐷6(0,𝑡−1)· 𝐷6(𝑥, 𝑡 − 1), if 𝐷𝑇 𝑜𝑡𝑎𝑙(𝑥, 𝑡) − ( ∑︀5 𝑐=1𝐷𝑐(𝑥, 𝑡)) ≤ 0 𝐷𝑇 𝑜𝑡𝑎𝑙(𝑥, 𝑡) − (∑︀5𝑐=1𝐷𝑐(𝑥, 𝑡)), otherwise. (3.1.4) (a) EU (b) NL
Figure 3.2: The levels of population/deaths for the Netherlands and the 14 selected Eu-ropean countries, for the data sets from this research (green) andActuarieel Genootschap
(2018) (red).
(a) EU (b) NL
3.2 Population Dynamics Data Description 39
the NL case we see that over the years the level of mortality in our data set increasingly exceeds that of Actuarieel Genootschap (2016). Some of this we can account for by the transformation introduced in equation (3.1.4), but most of the difference in slope is inherent to the cause-of-death data set we have used. In addition, we see that also the population size of our NL data set also exceeds its reference level increasingly, which compensates for the mortality surplus in terms of deaths per exposure ratio. In figure4.10later on in this research the correspondence of our data set in confirmed once more. Therefore, due to the equivalence with a trustworthy and commonly used data set used for socioeconomic studies and by prominent financial institutions, we proceed with the research with confidence in the accurateness of the underlying data set. Figure 3.3illustrates per cause of death the level of deaths over the time period 1970-2015 for the ultimate data sets. For the time period 1970 to 2015, the contribution of each cause to total mortality is presented in table
3.5.
Table 3.5: Contribution to overall mortality per cause-of-death, for males and females in EU and NL (1970-2015).
EU NL
Code Causes of death Male (%) Female* (%) Male (%) Female (%)
1 Circulatory system 32.1 37.0 29.6 31.1
2 Cancer 20.8 17.1 24.6 20.9
3 Respiratory system 4.9 4.2 7.4 6.2
4 External causes 5.8 3.8 4.2 3.2
5 Infectious and parasitic diseases 0.8 0.8 0.8 0.9
6 Other 35.6 37.2 33.4 37.7
*Due to rounding, the summation of all percentages does not equal 100%.
3.2
Population Dynamics Data Description
In this section we will elaborate on the data used for the simulation of future Dutch population up to the year 2060. We want to investigate the composition of the population in terms of age groups. More specifically, we want to research the weight of the older generations with respect to younger generations within Dutch society. We measure this dependency using equation (2.8.3), which renders the ratio of people older than 65 and people between 15 and 65 years old. Since the emphasis of this research lies on the mortality model, we have relied on a straightforward data set to model the population dynamics. We have acquired all variables from StatLine of CBS4, the most prominent statistical
institution in the Netherlands. From this source, we have used the most detailed data
available describing the variables we included in the model. Moreover, we will compare our forecasting results with the prognoses of CBS in a later stage. Therefore, the use of CBS data is an obvious choice as a foundation of the population model.
(a) Forecast birth. (b) Forecast net migration.
(c) Historical distribution net migration per age (1977-2015).
Figure 3.4: Graphical representation of used demographical data for the forecast of Dutch population, for male (blue) and female (pink).
Multiple demographic fluxes have an effect on the dynamics of a population. Together they constitute a yearly renewed balance of exposure, i.e. population size per age. We have opted to incorporate several demographic flows in order to predict future population composition. The first are the two variables discussed in the previous sections. One is the starting exposure. The exposure is taken from January the first 2016 and equates to nearly 17 million individuals of which 8.4 million are male. Another influence is the number of death each year per age group. We will calculated the number of deaths using the obtained crude mortality forecasts from previous sections. The mortality rates are from age 0 to 120 through application of equation (2.6.1).
3.2 Population Dynamics Data Description 41
males and the other 48.73% are females. The birth rate forecast is shown in figure 3.4a. In addition, we have chosen to include net migration per age to the model. For this we have again used the forecast performed by CBS. We then apply the historical distribution of net migration per age based on data from 1977 to 2015 to obtain net migration per age, between the ages of 0 and 99 years old. The distribution of the contribution per age to the net migration is illustrated in figure3.4c.
Results
4.1 Initial mortality calculations 43
4.1
Initial mortality calculations
The first part of this section we will present and discuss the crude mortality. After the steps performed in the previous chapter, we obtain the crude cause-of-death (henceforth referred to as ‘cod’) mortality directly from the ultimate data set. Next, the crude mortality rates are transformed for each cause of death to the marginal survival functions. Based on Li and Lu(2019), this is done under the assumption of both a Clayton and Frank survivor copula with parameter 𝜃 = 2. The subsequent net survival functions are transformed into net mortality intensities using equation (2.3.6) and again compared.
4.1.1 Crude mortality
We apply equation 2.2.2 and receive the crude mortality for EU data. In figure 4.1 the development is illustrated. The year 1970 is represented by the red line and the last year, i.e. 2015, is drawn in violet. The years in between are represented by a rainbow color ramp between these two colors. Moreover, a comparative overview for male and female is given in figure4.3.
Firstly, in general males mortality is higher than their female counterparts. Male mor-tality exceeds that of females for circulatory, cancerous and respiratory deseases (cod 1, cod 2 and cod 3). Especially for cancer, mortality levels are double and for some ages/years even triple. In total, it can be concluded that most deaths fall under the categories circula-tory diseases and other causes. Next, we see that cancer has the second highest mortality rate. For newborns, mortality is mostly comprised from death caused by infectious diseases or external/other causes (cod 4, cod 5 and cod 6). We also observe a much higher mortality rate for males from external causes (cod 4) around the adolescent ages, which is a known phenomenon commonly referred to as the ‘accident hump’.
the same age. Deaths from external causes (cod 4) have developed similarly for male and female. Overall trends for this cause imply mortality to be increasing after a period of decrease. For females we observe a mortality increase from infectious or parasitic diseases (cod 5). Only for newborns the mortality intensity has decreased. The same has happened for males, with the exception of the age group above 75, which has experienced a large mortality increase. Lastly, we see a general decrease of deaths under the category other causes (cod 6).
We perform the same procedure for NL data and receive the data presented in figures
4.2 and 4.4. Male mortality exceeds that of female, in line with the conclusions drawn from EU data. We again observe the biggest portion of mortality comprised by circulatory causes (cod 1). Instead of other causes for EU data (cod 6), the second biggest mortality contributor for NL data is cancer (cod 2). Again, we see that most newborn deaths are caused by infectious diseases or external/other causes (cod 4, cod 5 and cod 6) and that adolescent males are more prone to die from external causes then their female counterparts, as noted earlier for EU data.
(a) Crude mortality, male (1970-2015) (b) Crude mortality, female (1970-2015) Figure 4.1
EU: crude mortality per age from historical data, for causes of death 1-6.
A rainbow color ramp is used to represent the course of years from red (1970) to violet (2015).
(a) Crude mortality, male (1970-2015) (b) Crude mortality, female (1970-2015)
Figure 4.2
NL: crude mortality per age from historical data, for causes of death 1-6.
(b) Crude mortality, causes of death 4-6 (1970-2015) Figure 4.3
4.1 Initial mortality calculations 47
(a) Crude mortality, causes of death 1-3 (1970-2015)
(b) Crude mortality, causes of death 4-6 (1970-2015) Figure 4.4
(c) Male (NL) (d) Female (NL) Figure 4.5
4.1 Initial mortality calculations 49
4.1.2 Marginal survival functions
We will only shortly address the observed marginal survival functions in this section since next section we will discuss the net mortality intensities. We have explained in section
2.3.2 that the net mortality intensity is derived directly from the marginal survival func-tion. Hence, the inference of the net mortality is identical to that of the marginal survival function. Therefore, we omitted a lengthly explanation of the outcomes of the marginal survival functions in this thesis. After the application of equation (2.3.6) we obtain the graphs shown in figures C.1 to C.2 in Appendix C. Again, a rainbow color ramp is used to indicate the course of mortality over time with red representing 1970 and violet repre-senting 2015. We see that for both EU and NL data, the survival functions corresponding to a Clayton copula decrease much further with a significantly steeper slope than those corresponding to the Frank copula. This would indicate higher net mortality rates in the Clayton case compared to the Frank copula model.
4.1.3 Net mortality
Following the transformation of the marginal survival function described in equation (2.3.9), we obtain the marginal mortality intensities. For both EU and NL as well as for male and female, we have presented the results in figure4.5. We have excluded the results for other causes (cod 6), as this category is assumed independent and therefore its net mortality matches its crude mortality.
When we compare the marginal mortality rates to the crude mortality rates, we observe that the marginal mortality exceeds the crude mortality in all cases. This is in line with ex-pectations. The empirical data consisting of exposure and cause specific deaths is directly translated in the crude mortality. When death occurs from some cause, the individual im-plicitly can not die from another. We can regard this mutual exclusion as censoring of the real underlying stochastic process. In literature this is referred to as frailty (Oakes,1989;
net mortality intensities. Since we have eliminated the ‘distortion’, we obtain mortality rates that exceed the crude mortality rates. Only in the case of cod 1 net mortality for males and females we observe net mortality rates that approach the corresponding crude mortality (both EU and NL). This can be explained by the fact that this cause (circulatory diseases) is the main cause of mortality and is therefore less sensitive to distortion of other causes.
For older ages we see that the relative differences between crude and net mortality in-creases. This can be contributed to the fact that mortality is higher for these ages. The potential censoring probabilities become bigger for this reason and therefore net and crude diverge increasingly.
The trends of the net mortality under the Frank copula assumption matches the Clay-ton model mortality trends almost one to one. We observe that for males in the earliest observations of both NL and EU the Clayton mortality surpasses the Frank mortality for old ages. For females we see the same property in lesser degree. As time proceeds the curves tend to run parallel. The Frank copula model yields higher mortality in general than the Clayton copula. Moreover, we observe that the increase of mortality starts with a bit of lag under the Clayton assumption compared to the Frank assumption.
4.2
Parameter estimates
4.2.1 Maximum likelihood estimates
Through the formulation of the LL model from section 2.1.2 using maximum likelihood estimation method described in section2.4.1 we model the marginal mortality intensities. The resulting set of parameters is ˆ𝜃𝑐,𝑔 ∈ { ˆ𝐴𝑐,𝑔
𝑒𝑢(𝑥), ˆ𝐵𝑒𝑢𝑐,𝑔(𝑥), ˆ𝐾𝑒𝑢𝑐,𝑔(𝑡), ˆ𝑎𝑐,𝑔𝑛𝑙(𝑥), ˆ𝑏𝑐,𝑔𝑛𝑙(𝑥), ˆ𝑘𝑐,𝑔𝑛𝑙 (𝑥)}. The estimates are presented in figures C.3 to C.8. First we will discuss the differences in estimates between an assumed Clayton or Frank survivor copula. After that we will highlight some estimate properties that stand out.
