Multiple mortality modeling in Poisson Lee–Carter framework

ABSTRACT The academic literature in longevity field has recently focused on models for detecting multiple population trends (D'Amato et al., 2012b; Njenga and Sherris, 2011; Russolillo et al., 2011, etc.). In particular, increasing interest has been shown about “related” population dynamics or “parent” populations characterized by similar socioeconomic conditions and eventually also by geographical proximity. These studies suggest dependence across multiple populations and common long-run relationships between countries (for instance, see Lazar et al., 2009). In order to investigate cross-country longevity common trends, we adopt a multiple population approach. The algorithm we propose retains the parametric structure of the Lee–Carter model, extending the basic framework to include some cross-dependence in the error term. As far as time dependence is concerned, we allow for all idiosyncratic components (both in the common stochastic trend and in the error term) to follow a linear process, thus considering a highly flexible specification for the serial dependence structure of our data. We also relax the assumption of normality, which is typical of early studies on mortality (Lee and Carter, 1992) and on factor models (see e.g., the textbook by Anderson, 1984). The empirical results show that the multiple Lee–Carter approach works well in the presence of dependence.


Introduction
The whole financial system is dramatically threatened by the systematic improvements in longevity phenomenon, especially regarding the welfare and public pensions. The empirical evidence shows the presence of common factors that would affect survival rates across multiple populations in a similar way. The reason behind is that populations of the world are becoming more closely linked by communication, transportation, trade, technology, and disease.
Modeling mortality co-movements for multiple populations would have significant implications for longevity risk management.
In a multi-population mortality model, the analysis is focused on more than one population by taking into account a joint framework. The main idea behind these models lies on the convergence indefinitely over time observed between populations having similar socioeconomic features. The development of multiple setting is related to the several applications.
First of all, securitization of longevity risk, i.e., the transfer of longevity risk in the capital markets, which typically occurs through the creation of derivatives or securities whose cash flows are linked to the survival of a reference population. A more effective risk management can be obtained throughout the study of the mortality correlations between populations. The relationship between populations is involved in the estimation of the actual basis risk which arises from the difference in mortality improvements between the insured population and the population to which the standardized longevity hedging instruments are linked. Hedging instruments have been developed in order to help pension funds to protect themselves against longevity risk, in particular by reinsurance and longevity derivatives.
Second, the multiple population mortality model allows for enhancing the estimate reliability by modeling the smaller population jointly with a larger population (Li et al., 2010).
Furthermore, the multiple approach can also be used for enforcing greater consistency of the sex differentials when an analysis on both genders has to be performed.
Recently, significant developments in multi-population modeling have been recorded. The seminal interest in studying cross-country longevity common trends focused on modeling the interdependence of the mortality rates of two populations. Li and Lee (2005) developed an augmented common factor model for modeling convergent mortality dynamics. Li and Hardy (2011) propose a model that measures the population basis risk involved in a longevity hedge. Cairns et al. (2011) represents the short trends by using a mean-reverting spread, where the long-run improvements are parallel according to biological principle. Jarner and Kryger (2011) develop a similar two-step approach modeling the mortality of larger reference population and the mortality spread between the two populations. The current studies investigate on long-run equilibrium relationships across countries, by considering more than two populations in the mortality framework, for capturing valuable information about the factors driving changes in mortality (Njienga et al., 2011;Russolillo et al., 2011). Other authors model mortality dependence (structure) across countries using a dynamic copula approach (Yang et al., 2013). In particular, they employ time-varying copula introducing mortality dependence and demonstrate symmetric dependence. In the same context, MacMinn et al. (2013) are the first to apply the factor copula model to mortality fitting for multiple populations. They focus on the residual risk (tail dependence) by setting out an efficient approach for high-dimensional data. Kleinow (2013) shows that the period effects of different populations in a CAE model are better comparable with each other since the impact of different age effects is eliminated. This aspect can be relevant in effectiveness of hedge positions.
In order to study cross-country longevity common trends, it is essential to consider tools to quantify, compare, and model the strength of dependence. Therefore, it is necessary to take into account either the dependence for adjacent age groups, or the dependence structure across time in a single population setting: a sort of intra-dependence structure (D' Amato et al., 2012c). At the same time, it is important to consider the dependence across multiple populations, what we call inter-dependence, for capturing common long-run relationships between countries.
In a previous contribution (D' Amato et al., 2013), we kept into account the issue of crosssectional and time dependence, by proposing an algorithm based on the Lee-Carter framework (1992), but relaxing the assumption of normality, which is typical of early studies on mortality (Lee and Carter, 1992) and on factor models (see e.g., the textbook by Anderson, 1984). In this article, we extend that formulation by taking into consideration a different fitting model for mortality data, the Poisson log-bilinear mortality model, which allows to overcome the problems associated with the ordinary least square (OLS) method in the fitting procedure.
We also consider a bootstrap procedure for dependent data, thus preserving both the historical parametric structure and the intra-group error correlation structure. In particular, we apply a Sieve bootstrap algorithm (Bühlmann, 1997) to the vector autoregression (VAR henceforth) model containing the estimated common factors (both stationary and non stationary). However, in a previous paper (D' Amato et al., 2012c), we stated that in our context we cannot apply a standard Sieve bootstrap algorithm, since, when resampling the estimated common factors, a generated regressor problem arises. In this order of ideas, our article is based on Trapani (2012), which develops the full-blown theory to apply Sieve bootstrap to the context of non stationary panel factor series, developing selection rules for the order of the VAR and showing the superior performance of Sieve bootstrap compared to first-order asymptotic. In particular, the article is structured as follows: Section 2 faces the issue of multiple alignment. In Section 3, we present the Poisson log-bilinear mortality model. Section 4 discusses the methodology to generate the bootstrap sample by proposing the multiple Poisson Lee-Carter Panel Sieve Algorithm. Section 5 shows the results of the Numerical Application. Section 6 contains the Concluding Remarks.

Multiple setting
The representation of multiple populations bases on similar mortality behaviors typical for people which share analogous living conditions. Studying mortality experience for a group of populations with similar mortality behaviors might improve the stability of mortality modeling (Yang et al., 2011). Moreover, it could allow for solving the problem of small population. Indeed, some authors propose the replication of the data by mixing appropriately the mortality data from neighboring countries (Olivieri, 2011). The crucial point working with different populations is that the dependence structure analyzed in previous works for a single dataset (D' Amato et al., 2012a) becomes very complex and has to be taken into account under a multidimensional approach. In fact, three kinds of dependence have to be captured: the crosssectional dependence for adjacent age groups, across countries, and serial/time dependence. In this case, the classical Sieve bootstrap cannot be applied to the three-dimensional dataset, due to a too large number of cross-sectional units. D' Amato et al. (2013) originally solve this problem looking for a rational reduction of the dataset. To this aim, they consider that the common trends between countries are captured by the time-varying parameters k t of the LC model and propose the following alignment of the data (we refer to Lee and Carter, 1992 for a deepened description of the parameters involved in the Lee-Carter method). They fit separately the LC on some mortality dataset of M different populations, composed by the same ages x = a, a + 1, . . . , a + N and years t = b, b + 1, . . . , b + T , where a represents the first age, fixed equal to zero and b the first time, respectively. N and T represent the last age and the length of the period considered, respectively. Once they obtain the k t 's for each country, they arrange the M time series of k t in a matrix, generating a panel data in which the single units are represented by the different populations and are collected in rows. In this way, they obtain a reduced dataset on which it is possible to design a Sieve bootstrap. The framework they propose is very flexible and lends itself to interesting extensions and more accurate formulations. One of this is the possibility to take into account a different fitting model for mortality data. An example could be the Poisson log-bilinear mortality model, described in the next section.

The Poisson log-bilinear mortality model
D' Amato et al. (2012c) develop the idea of first fitting Lee-Carter parametric model, and then re-sampling a particular class of the residuals, the so-called centered residuals, according to the Sieve scheme, through an autoregressive approximation for generating bootstrap replications of the data. They first consider the Lee-Carter parametric framework because of its well-known properties (Deaton and Paxson, 2004); however, this model making use of the least-square tool for mortality estimation is based on the hypothesis of homoscedastic errors, but this assumption is not confirmed by the empirical evidence. Instead, the Poisson log-bilinear model allows to overcome the problems associated with the OLS method in the fitting procedure. Because the number of deaths is a counting random variable, according to Brillinger (1986), the Poisson assumption appears to be plausible. It has been argued that the number of deaths, when the central exposed-to-risk is given, may be assumed to follow a Poisson distribution. At the same time, the promising estimates may be obtained by fitting the Poisson regression (see Renshaw and Haberman, 1997;Renshaw and Haberman, 1996;Sithole et al., 2000): where the parameters are still subjected to the constraints t k t = 0 and x b x = 1 and the force of mortality is thus assumed to have the log-bilinear form:

Algorithm: Multiple Poisson Lee-Carter panel sieve
This section discusses the methodology to generate the bootstrap sample. The preliminary step is the construction of matrix K through the fitting of the LC model on different populations separately.
For each population i, we fit the Poisson version of the LC model according to eq. (1), determiningk ti . The parameters k t,i are then arranged in a matrix K MXT . A VAR is fitted to this matrix. The VAR is the statistical tool employed to represent how different populations are related each other.
Let K t = [k 1t , …,k Mt ] denote an (Mx1) vector of time series variables. The basic q-lag VAR(q) has the form: where A is the matrix of coefficient of the selected VAR(q) model. Hence the bootstrapping algorithm is as follows: Step 1. (2.2) Compute the residualsê t,q = k t − qK j=1Â q, j k t− j and center them around their mean, defining themē t,q .
In this way, we obtain b matrix K b MXT from which we can generate the pseudo-sample {μ b xt,i } T t=1 according to eq. (1).

Numerical application
In the present section, we apply the methodology described in the Section 4 to the historical mortality data for five countries expected to have experienced common longevity improvements, on the basis of similar socioeconomic features: United Kingdom (henceforth UK), France, Italy, Spain, and Belgium. The study is performed for each country on total population (composed by male and female) ranging from 1950 to 2006, for ages from 0 up to 110 years, considered by single calendar year and by single year of age, where the class of age above 100 years is collected in an open age group 100+ (the data were downloaded from the Human Mortality Database). The first step is the fitting of the Poisson version of the LC on the datasets of the five selected countries. In Fig. 1, we show the estimates of the trend parameters for each country: Observing Fig. 1, we note that the estimated trends of the parameter k t are quite similar for each country, showing a decreasing slope. We refer to D' Amato et al. (2013) for a detailed description of similarities and measurement of dependence in mortality between the selected countries. In the second step of this numerical application, we fit the VAR model to the k t of each country and calculate the residuals. The number of lags of the VAR model is selected through the application of the Akaike criterion, useful to avoid the overparameterization of the model. The result is that the selected lags are two so that each lag term Kt-q, m i for country i affect  the response K t , m j with i = j and i ࣔ j, q = 1,2. In other words, not only own country lags are relevant, but each country influences the response of the others for two legs. Figures 2-6 display, for each country, a diagram of fit, a residual plot, the auto-correlation, and partial auto-correlation function of the residuals.
Finally, we have implemented the bootstrap simulation according to the algorithm proposed in this article with a number of replications equal to 1000 and we have obtained, for each period, the sample average k t . Then, we have projected them by using ARIMA models and calculated confidence intervals. Figure 7 displays the mean of the simulated k t and the projections with confidence intervals for UK, Belgium, Spain, France, and Italy.
In Fig. 7, the solid lines represent the k t trend for UK, Belgium, Spain, France, and Italy, with the confidence intervals highlighted in yellow. As it is clear by the diagrams, the trend is quite similar between countries, which have experienced similar trends in mortality reduction due  to common social-economic factors and improvements in medical research. But if we look at the width of the confidence intervals, we can see they are quite unequal between countries. In particular, as result of our application, it is relevant to notice the wider CI's for Belgium, which is the smallest of the considered countries. In the Introduction, we stated that by considering a multiple model we allow enhancing the reliability for small-population mortality estimates. In other words, the knowledge about countries with small population can be enriched by the consideration of common features between different populations. Further research will aim to verify the forecast goodness of the multiple Lee-Carter approach throughout different measures of forecast accuracy and also to test the improvement of the proposed model with respect to the classical Lee-Carter model.

Concluding remarks
In the last few years, there has been an increasing interest in the actuarial literature on the issue of modeling multiple populations characterized by similar socioeconomic and living conditions. It has also been shown that the problems arising when considering simultaneously different countries are strictly related to the dependence analysis. The existence of dependence in mortality data involves the interactions for adjacent age groups between age and time and also across the different populations.
From a methodological point of view, the research presented in this contribution is oriented to consider the flexibility of the methodology proposed in a previous article (D' Amato et al., 2013). In particular, we build on the Lee-Carter method in its Poisson version to develop an approach within this general framework. In particular, we extend the basic structure to include some cross-dependence in the error term on the basis of a tailor-made bootstrap algorithm, explained in detail in Section 4.
Another attractive methodology for our application could have been Bayesian (empirical, full) with the k t as random effects. Bayesian simulation differs from classical simulation analysis in that probability distributions are used to represent the uncertainty about model parameters, rather than point estimates and confidence intervals. In the Bayesian approach, the prior information about the distribution of parameters is modified into posterior one after the observation of sample. In the language of uncertainty, classical simulation models only aleatory uncertainty (the randomness of the system itself), while Bayesian simulation models both the aleatory and epistemic uncertainty (the lack of knowledge about the system). However, as it has been observed (Hastie et al., 2009), "bootstrap distribution represents an approximate non-parametric posterior based on a particular choice of noninformative prior … and is typically much simpler to carry out. " In other words, referring to our case the prior from which we start is the fitted kt from the model and the bootstrap algorithm approximate the posterior information about the distribution of the parameter offering point estimates and confidence intervals; in this way, we reach a quicker results than that we should obtain with Bayesian simulation. In actuarial application, it is important to reduce the time of simulation to make really applicable the algorithm for practical uses; moreover, for practitioners the quantities of interest are just the confidence intervals to estimate the impact of longevity risk on insurance liabilities.