Investigating Mortality Uncertainty Using the Block Bootstrap

This paper proposes a block bootstrap method for measuring mortality risk under the Lee-Carter model framework. In order to take account of all sources of risk the process risk, the parameter risk, and the model risk properly, a block bootstrap is needed to cope with the spatial dependence found in the residuals. As a result, the prediction intervals we obtain for life expectancy are more accurate than the ones obtained from other similar methods.


Introduction
For actuarial pricing and reserving purposes, the mortality table needs to be projected to allow for improvement in mortality to be taken into account.It is now a well-accepted fact that mortality development is difficult to predict; therefore, a stochastic mortality modelling approach has been advocated 1, 2 .For a review of recent developments in stochastic mortality modelling, interested readers are referred to Pitacco 3 and Cairns et al. 4 .Under the framework of a stochastic mortality model, future mortality rates are random as are the other quantities derived from the mortality table.In order to manage the mortality risk properly, we need to assess the uncertainty coming from the mortality dynamics carefully.
There are three types of risk embedded in adopting a stochastic mortality model see 5 : a the "process risk", that is, uncertainty due to the stochastic nature of a given model; b the "parameter risk", that is, uncertainty in estimating the values of the parameters;

Definitions and Notation
The following terminology will be used at various points in this paper.For an individual at integer age x and calendar time t, we have the following: i D xt is the number of deaths at integer age x in calendar time t.
ii ETR xt is the exposure-to-risk at integer age x in calendar time t, that is, the personyears lived by people aged x in year t.
iii μ x s, t , 0 ≤ s ≤ 1 is the force of mortality.This is assumed to be constant within each year of age, that is, μ x s, t μ xt , 0 ≤ s ≤ 1.
iv m xt is the central mortality rate and is defined as D xt /ETR xt .Under the constant force of mortality assumption, m xt μ xt .
v q x t is the probability that an individual aged x in year t dies before reaching age x 1. p x t 1−q x t .Under the constant force of mortality assumption, p x t e −μ xt .vi e 0 is the life expectancy at birth and is calculated as In this paper, we investigate how e 0 is influenced by the uncertainty coming from projections under the Lee-Carter model.

Model Description
Lee and Carter 6 proposed a two-stage dynamic procedure to describe the secular change in the logarithm of age-specific central rates:

3.2
In the first stage, 3.1 is fit to the historical mortality data from year 1 to t 0 , i.e., t 1, . . ., t 0 to estimate the parameters a x , b x , and k t .The interpretation of the parameters is as follows.a x describes the general age shape of the log central rates, while the actual log rates in a specific year change according to a time-varying mortality index k t modulated by an age response variable b x .In other words, the profile of b x indicates which rates decline rapidly and which decline slowly over time in response to changes in k t .xt contains the error that is not captured by the model at age x in year t.
In the second stage, the estimated values of k t from the first stage are fit to a random walk 3.2 .According to 3.2 , the dynamics of k t follow a random walk with drift parameter c and normal noise ξ t with mean 0 and variance σ 2 .

Estimation of the Parameters
Instead of using least-squares estimation via a singular value decomposition SVD to estimate a x , b x , and k t as in the original Lee and Carter paper, we adopt the methodology proposed by Brouhns et al. 12 .The main drawback of the SVD method is that the errors xt in 3.1 are assumed to be homoscedastic, which is not realistic.In order to circumvent the problems associated with the SVD method, Brouhns et al. 12 proposed the Poisson assumption to be applied to the death counts D xt .That is,

3.4
We now can write the log-likelihood function for the parameters a x , b x , and k t : Maximizing 3.5 iteratively provides the MLE estimates a x , b x , and k t 12 .

Uncertainties in the Lee-Carter Mortality Projections
Let a x , b x , k t , c, and σ represent all the estimated parameters obtained from fitting the Lee-Carter model to the data from year 1 to t 0 .In order to forecast the time-varying index k t n years ahead given all the data up to time t 0 , formula 3.2 is extrapolated as Here, ξ j represents the uncertainty coming from forecasting k t .
Accordingly, the future central mortality rates m xt in the calendar year t n t 0 n are given by log m xt n a x b x k t n , 3.7 where a x and b x remain constant.From these rates, we can construct projected life tables and calculate the associated life expectancies using the relations given in Section 2.2.Under the Lee-Carter model, k t follows a stochastic process.Further, according to model 3.7 , all future mortality rates in the same year are affected by the same timevarying random variable k t , where ξ t is the sole error source of the model.Under this circumstance, the projected life expectancy given in formula 2.1 can be viewed as a sum of a sequence of comonotonic variables, since all survival probabilities can be written as monotone transformations of the same underlying random variable which is the noise ξ t in this context so that they are perfectly positively dependent.Comonotonicity allows us to express the quantiles for projected life expectancies in terms of the quantiles of ξ t analytically.According to Denuit 13 , the p-percentile for e 0 in the forecast year t n taking account of only the uncertainty from k t can be computed as The Denuit quantile formula 3.8 will be used in this paper as a benchmark to compare with our simulated prediction intervals.
Of course, uncertainty may also come from errors in parameter estimation and/or from model discrepancy.If the model is properly specified, it is reasonable to expect that these other sources of error are minor.Under the assumption that different sources of error at a given age are uncorrelated and that the error sources other than k t are uncorrelated across age, Lee and Carter 6, in Appendix B estimated that these other types of error contribute less than 2% of the forecast error for life expectancy e 0 for forecast horizons greater than 10 years.However, this assertion is not supported by the data; inspection of the residuals reveals substantial correlations across age and year 9, 14 .In other words, spatial dependence has been found in the residuals.It has also been found that the actual coverage of the prediction intervals based solely on the process risk of k t is lower than their nominal level see 14 .This indicates the need for appropriate methods to take into account different types of uncertainty.

The Deviance Residuals
As mentioned above, we need to check the residuals to assess model adequacy.Under the Poisson Lee-Carter model, deviance residuals can be used for this purpose in a similar manner to how ordinary residuals are used in regression analysis.For reference, see Maindonald 15 .Deviance residuals are calculated as 3.9

Ordinary Bootstrap Method
It is not possible to analytically obtain prediction intervals for e 0 taking account of all sources of variability.One numerical method that has been employed for this purpose is the bootstrap.The basic idea of the bootstrap is to artificially generate resamples of the data which can be used as a basis for approximating the sampling distribution of model parameter estimates.Koissi et al. 9 give a brief overview of the bootstrap and also some references in the context of applying the bootstrap method to the Lee-Carter model.They use an ordinary deviance residual-based bootstrap for this problem.
Briefly, their approach is to fit the Lee-Carter model to the data and to obtain deviance residuals.They then sample with replacement from these residuals to obtain resampled residuals.These resampled residuals are transformed back to get resampled death counts.The Lee-Carter model parameters are estimated from this artificially generated data set, and the corresponding k t process is also fit and resimulated.As a result, we obtain a projected mortality table see Section 3.3 for projection method .The process is repeated a large number of times say, N 5000 giving a collection of projected mortality tables.Approximate 90% prediction intervals for e 0 are given by the 5% percentile and 95% percentile of the e 0 's derived from these tables.Better bootstrap prediction intervals have been proposed see, e.g., 16 , but any improvement in accuracy will depend heavily on the validity of the Lee-Carter model assumption.See also D'Amato et al. 17 for a stratified bootstrap sampling method which reduces simulation error.

Block Bootstrap Method
The ordinary bootstrap is not valid in case of spatial dependence, since it is based on simple random sampling with replacement from the original sample.In order for the bootstrap sampling distribution of the statistic of interest to be a valid approximation of the true sampling distribution, it is thus necessary for the original sample to be a random sample from the underlying population.Spatial dependence represents a violation of this assumption; it can cause either a systematic overestimation or systematic underestimation of the amount of uncertainty in a parameter estimate, depending upon whether negative or positive dependence dominates.
Confidence and prediction intervals constructed from bootstrap samples in the presence of positive spatial autocorrelation will thus be too narrow so that their actual coverage probabilities will be lower than their nominal levels.In the presence of negative autocorrelation, actual coverage levels will systematically exceed nominal levels.
One remedy for this problem is the block bootstrap, which gives improved approximate prediction intervals, taking account of the spatial dependence, at least locally.The basic idea is to resample rectangular blocks of observations, instead of individual observations one at a time.A full bootstrap resample is obtained by concatenating randomly sampled blocks of observations.Thus, apart from the "seams" between the blocks, some of the dependence between neighboring observations will be retained in a bootstrap resample.If the blocks are taken large enough, then most of the dependence between observations located near each other will be retained in the bootstrap resample.Of course, block sizes should not be taken too large; otherwise, the resamples begin to resemble the original sample too closely, and this will again cause uncertainties to be underestimated.The method provides reasonable approximations for large samples in much the way that the m-dependent central limit theorem provides a normal distribution approximation for averages of locally dependent data.It should be noted that the method will fail without adjustment for data with long memory.A full account of the block bootstrap for time series is given by Davison  To set up a new artificial set of residuals, we start with an array, which has the same dimensions as the original matrix of residuals.The empty array is then partitioned into smaller rectangular blocks.Each block is replaced by a block of the same size, which is randomly selected from the original matrix.These random blocks are constructed as follows.First, uniformly randomly select an element from the original matrix.The associated block consists of all residuals in the appropriately sized rectangle to the southeast of the chosen point.In cases where part of the rectangular block falls outside the original matrix, a periodic extension of the matrix is used to fill out the remaining cells of such a rectangle.As stated earlier, more details can be found, for example, in Nordman et al. 18 .
Nordman et al. 18 give advice on choice of square block size for situations where the parameter of interest is a smooth function of the mean.The optimal block size is based on a nonparametric plug-in principle and requires 2 initial guesses see the appendix .Our parameters of interest do not completely satisfy their assumptions, but we have experimented with their recommendations.In the absence of firm theoretical guidance, we have found it useful to plot a correlogram of the original raw residuals and compare with the resampled residuals.The correlogram is a plot of the spatial autocorrelation against distance.See Venables and Ripley 20 for more information.If the correlograms match reasonably well, this gives us confidence in our block choice.

Fitting Performance on Swedish Males 1921 to 1960
We fit the Lee-Carter model to the Swedish male data.The estimated parameters a x , b x , and k t are displayed in Figures 1, 2, and 3. We then computed the deviance residuals.Contour maps are plotted in Figure 4, while the correlogram plots are given in Figure 5.For comparison purposes, we put together the results based on the original deviance residuals and the resampled residuals of different block sizes.
Both plots are highly suggestive of spatial dependence, because of the occurrence of large patches of large positive and large negative residuals, that is, clustering in Figure 4 a .Evidence of short range dependence is also seen in the correlogram of the raw residuals.To see if such clustering could be due purely to chance, we simulated from the Poisson Lee-Carter model three times and fit the model to the simulated data and plotted the deviance residuals.The results are shown in Figures 6 and 7.These plots show no sign of spatial dependence; clusters appear to be much smaller than for the real data.We conclude that the ordinary bootstrap will not give valid prediction intervals.Thus, we employed a block bootstrap as described in the previous section.

Prediction Intervals for e 0 Based on Block Bootstrap
We considered a number of block sizes.The theoretical block sizes suggested by Nordman et al. 18 depend on the actual parameter being estimated, and they also depend  on the initial guesses.The optimal sizes estimated depend on different forecast horizons.Table 1 gives two sets of block sizes: one marked as m 1 is based on initial guesses of 5 × 5 and 10 × 10, and the other marked as m 2 on 10 × 10 and 20 × 20.Note that the results vary a lot, depending on the initial guess as well as the forecast horizon.Thus, it is not clear whether these block sizes are really the best, and we checked correlograms for some different block sizes including rectangular ones.Figure 5 shows correlograms for some of the better cases raw residuals, and 8 × 4, 10 × 5, 15 × 10, and 20 × 12 block sizes .We have chosen to work with 15 × 10 blocks, since that case matches the raw residuals most closely.As further confirmation, revisit the residual plots in Figure 4.These correspond to a raw residuals, b 1 × 1, c 8 × 4, and d 15 × 10 block sizes.Quantitatively, one looks for the same kind of "patchiness" in the bootstrap residuals as can be seen in the raw residuals.There are obvious discrepancies when the block size is too small; there are no large patches of large positive or large negative residuals at the 1 × 1 block size, for the larger block sizes, we see patterns which are more like the patterns seen in the raw residuals.
Prediction intervals are displayed in Figure 8.We have plotted the 90% pointwise prediction bands based on the Denuit quantile formula, the ordinary bootstrap and the 15×10 block bootstrap.The observed life expectancy e 0 is also shown.We comment on these results in the next section.

Discussion
In Figure 8, we note that the theoretical prediction intervals calculated from the Denuit quantile formula is the narrowest.This makes sense since those intervals have only accounted for the uncertainty in forecasting k t 's future random path.With bootstrap methods, sampling errors in the parameter estimates and partial model uncertainty have also been considered.Therefore, the prediction intervals from simulations are expected to be wider than the theoretical ones, which has been confirmed by our results as well.Secondly, while the prediction intervals from the ordinary bootstrap are only marginally wider than their theoretical counterparts, we note that the prediction intervals from the block bootstrap make a substantial difference.This is due to the fact that the deviance residuals are not pattern-free random noise.When this is the case, the residuals may carry important information about the correlations among the estimated parameters, and/or may indicate some type of model risk.Ordinary bootstrap resampling with replacement then destroys this information and thereby diminishes the variability.Therefore, the block bootstrap is needed to obtain accurate prediction intervals.By taking into account the other types of risk properly, we obtain wider prediction intervals which serve as a better representation of the uncertainty associated with mortality predictions.
Another observation concerns the importance of the k t as a component of the variability.Previously, it has been claimed that the uncertainty from k t dominated all other errors.This is again an argument that is only true under restricted conditions see Section 3.3 for more details .Figure 9 presents k t 's share in the total forecast error and shows that other sources of variability are not negligible even when the forecast horizon becomes very long.
Finally, we note that most of the observed life expectancies lie within our prediction intervals, but some of them mainly in the later years fall outside, indicating that there may still be factors influencing mortality that are not explained by the Lee-Carter model, even with dependence accounted for.This is certainly a cause for concern: model misspecification is a problem when we use the Lee-Carter model for long-term mortality prediction.Sometimes introducing more variables may improve model fitting and remove correlations in the residuals to some extent.Examples of those models include Renshaw and Haberman 21 that allows for the cohort period effect and Renshaw and Haberman 22 that adds a second or even third bilinear term to the Lee-Carter for better fit.However, as shown in Dowd et al. 23 and Cairns et al. 24 , those models may suffer from nonrobustness in parameter estimation and sometimes generate implausible or unstable predictions.It is our belief that simpler models are preferable for long-term prediction, since they run less risk of being overfitted.
However, we remark that the block bootstrap is not able to solve the fundamental problem of model misspecification, which is not the purpose of the paper nor a problem that can be easily fixed.This paper tackles the problem from a different perspective-we adopt the Lee-Carter model framework but work on deriving more accurate measures of the uncertainty degree in mortality prediction.There are other approaches using a similar idea-for example, the Bayesian method and MCMC simulation as discussed by Cairns et al. 25 .The advantage of our proposed block bootstrap method is in that it is distribution-free.Without knowing the exact form of correlation among errors, we are able to provide more accurate prediction intervals for mortality predictions.A.3

Figure 1 :
Figure 1: Parameter a x based on the Swedish male 1921-1960 mortality data.

Figure 2 :Figure 3 :
Figure 2: Parameter b x based on the Swedish male 1921-1960 mortality data.

where σ 2 b
is the sample variance of the bootstrap estimates of the parameter estimates based on a b × b block, andB 0 2b 2 σ 2 b 2 − σ 2 2b 2 .
and Hinkley 16 , and the paper by Nordman et al. 18 describes the spatial block bootstrap we will employ.See Taylor and McGuire 19 for another actuarial application of bootstrap for dependent data.

Table 1 :
The theoretical block sizes for the forecast horizons 1 to 46 years.