The Lee-Carter quantile mortality model

The Lee-Carter (LC) stochastic mortality model has been widely used for making future projections of mortality rates. In the framework of the LC model, the response function is non-linear in parameters. Here, we adapt this LC framework to compute conditional quantiles. The LC quantile model can be deﬁned as quantile non-linear regression conditioned to age and the calendar year. Two strategies for estimating coeﬃcients based on interior-point methods are described. We show that the LC quantile model provides additional information to that furnished by the traditional LC conditional mean. An application to Spanish mortality data is reported.


Introduction
Social progress has led to a steady increase in human longevity over the last 150 years.These gains in life expectancy in the late 19th and early 20th centuries reflected improvements in infant mortality rates, but in the last 50 years the increase has been mainly driven by a fall in death rates among the elderly [Wilmoth, 2000].And there are no signs that human longevity has reached its maximum limit; thus, for instance, average life expectancy increased by 5.5 years between 2000 and 2016 [WHO, 2019].This rising longevity and the associated ageing population represent major challenges to those concerned with the planning of pensions, healthcare provision and long-term care systems.As a result, the analysis and modeling of uncertainty in future mortality has become a major topic of interest in multiple fields of research.
A number of stochastic mortality models have been developed to analyze mortality rates.One of the most influential approaches to the stochastic modeling of such rates is the parametric non-linear regression model proposed by Lee and Carter [1992].The Lee-Carter (LC) model proposes estimating the expected mortality rate as the non-linear combination of age and period parameters.Subsequently, many attempts at developing mortality models have drawn inspiration from the LC model, including, but not limited to, Brouhns et al. [2002], Currie et al. [2004], Renshaw andHaberman [2003, 2006], Cairns et al. [2006] and Plat [2009].Two recent articles have outlined the common framework employed by most of these models [Hunt and Blake, 2015;Currie, 2016], the aim of which has been to estimate the expected mortality rate.Although important, the mean is just one of the measures of interest of the conditional response distribution.Indeed, analytic techniques that examine conditional distribution locations other than the mean are required to provide a more detailed description of the conditional distribution function of the mortality rate.
In this article, we define the Lee-Carter quantile model for estimating the conditional quantiles of the mortality rate.The LC quantile model is defined as a quantile regression that models conditional quantiles as a non-linear function of age and period parameters.Quantile regression (QR) was first introduced by Koenker and Bassett [1978] to model conditional quantiles as a function of predictors.QRs have been widely applied in mortality analysis, but most studies have tended to analyze the factors that account for the conditional quantiles.For instance, Yang et al. [2012] examine the effects of inequality across the mortality distribution in the US, and Kizhakethalackal et al. [2013] and Swain et al. [2017] apply quantile regressions to study the factors affecting the child mortality.A different approach has been adopted recently by Peters [2018] who applies QRs to the modelling of mortality time-series data.Similarly, Uribe et al. [2018] apply time-series techniques to model trends in the quantiles of the survival function.Here, we apply QRs in the framework of stochastic mortality predictive models, making it, to the best of our knowledge, a novel approach.
Modeling the conditional quantile mortality rates have appealing advantages.First, the framework employed by the LC quantile model means we can capture all the distributional features of the conditional distribution.For example, the observed pattern of the mortality rate frequently shows a different variability at different ages [Brouhns et al., 2002].This feature would be overlooked by the regression approaches that focus on the mean.In the LC quantile model the effect of the covariates on the response may vary at different percentiles of the distribution, and, thus, complex relationships of predictors and the response may be captured.Second, when the center of the conditional distribution of the mortality rate is of interest, the LC median regression model provides more robust estimators than the LC mean regression in presence of outliers.Third, the LC quantile regression model does not assume any particular form on the error distribution, and so LC quantile model is appropriate when errors are not normal distributed.Fourth, the LC quantile model is stable under monotonic transformations, so conditional quantile log mortality rate estimates and log conditional quantile mortality rate estimates may be considered indistinctly.Finally, the tails of the conditional distribution of the mortality rate can be modeled and used to construct prediction intervals.
The use of quantiles is a major topic in the actuarial science, particularly in the field of risk measures in risk management.Actuaries have been traditionally interested on the behavior of loss random variables in the tails of their distributions, evaluated by means of appropriate high quantiles [Embrechts et al., 2009;Belles-Sampera et al., 2014].Here, we argue that the conditional quantiles of the mortality rate distribution are also of great interest by actuaries.Median mortality rate forecasts may be used by actuaries to compute annuity factors for reserving and pensions, due to the appealing statistical properties of their estimators.On the other hand, low quantiles of mortality rates have a direct application in calculating stressed annuity factors in adverse scenarios.
The LC quantile model is a non-linear parametric regression model.Yet, estimating the coefficients of a non-linear quantile regression is far from straightforward and here, to address that challenge, we adopt the interior point method described in Koenker and Park [1996] to estimate our LC quantile model parameters.Two alternative calibration strategies are described: a) an iterative process of a sequence of linear optimization problems, and b) the linearization of the objective function.Subsequently, we apply the model to Spanish mortality data.Macias and Santolino [2018] previously compared the impact of two different mortality projection models on life insurance products in Spain, but here we fit only the LC quantile model to the Spanish male and female populations.
The rest of this article is structured as follows.The Lee-Carter model is described in Section (2).The quantile regression is described in Section (3).The LC quantile model is defined and the procedure to estimate the parameters outlined in Section (4).An application to Spanish mortality data is presented in Section (5).Finally, the main conclusions are reported in Section (6).

Lee-Carter stochastic mortality model
The Lee-Carter mortality model, introduced in 1992 [Lee andCarter, 1992], can be defined as follows: where m x,t is the central rate of mortality at age x and year t, a x and b x are the specific age parameters and k t the time-varying index, x = 0, ..., ω and t = 1, ..., T .Finally, the error ε x,t has mean 0 and variance σ ε .Infinite solutions exist.For any scalars c and d, fitted values.To overcome the lack of identifiability, Lee and Carter [1992] proposed two constraints x b x = 1 and t k t = 0.
The conditional expectation is equal to The expectation is the value that minimizes the sum of squared errors.One strategy for estimating the parameters is to minimize the squared errors.However, this model cannot be directly estimated by ordinary least squares because the right-hand side of equation ( 1) is not linear with the parameters.
To estimate the coefficients, Lee and Carter [1992] proposed the application of singular value decomposition (SVD), that is, decomposing the matrix of log mortality rates once the average over time of log age-specific rates have been subtracted.Doing it, a vector of coefficient estimates θ is obtained, θ = (ã 0 , ..., ãω , b0 , ..., bω , k1 , ..., kT ) .In a second step, the authors proposed recalibrating kt by iterative processes to match the estimated number of deaths with the observed number of deaths in period t, where N x,t and E x,t are the observed number of deaths and the exposure-to-risk in period t and at age x.The motivation for this second-stage estimate is to avoid sizable discrepancies between the numbers of predicted and actual deaths (which are likely to occur because the first step is based on logarithms of death rates) [Lee and Carter, 1992;Brouhns et al., 2002].Brouhns et al. [2002] proposed estimating Lee-Carter model so as to take into greater variability in older than in younger ages.The logarithm of the observed mortality rates is much more variable at older ages than at younger ages because of the much smaller absolute number of deaths at older ages.They suggested modelling the number of deaths by means of a non-linear Poisson regression model with exposure-to-risk.By so doing, the variance of the mortality rate is inversely related to the exposure-to-risk.Various authors have proposed estimating parameters by maximum likelihood.Renshaw and Haberman [2006], for example, showed how maximum likelihood estimates may be obtained under the original LC Gaussian error structure or the Poisson error structure using an iterative process.Villegas et al. [2018] defined a unified framework of stochastic mortality projection models that include LC and other mortality models and the authors introduced the R package StMoMo to fit these models via maximum log-likelihood.

Quantile regression
Let Y be a continuous random variable (r.v.) with finite expectation and cumulative distribution function F Y defined by F Y (y) = P (Y ≤ y).The inverse function of F Y is known as quantile function, Q.The quantile of order α is defined as ≥ α} where α ∈ (0, 1).The quantile is a left-continuous increasing function.If F Y is continuous and strictly increasing, the mathematical expectation can be represented as Quantile regressions were introduced by Koenker and Bassett [1978].Let us define the conditional quantile where X is a row vector of exogenous variables and F Y |X the conditional cumulative distribution.Suppose that the conditional quantile can be defined as for a certain function f , where θ α is a vector of coefficients which depends on α.In the particular case of the quantile linear regression, we have that f (X, θ α ) = X θ α .
For a given covariate set, the conditional quantile function fully characterizes the entire conditional distribution function, unlike the mean which is just one of its characteristic quantities.Parameters of quantile regression are estimated by least absolute techniques.Quantile regression has many of the appealing properties of the ordinary sample quantiles [Koenker, 2005].Thus, least absolute regression estimates are less sensitive to the presence of outliers than ordinary least square regression estimates.Another appealing property is that the QR is stable under monotonic transformations.For any monotone function g, it holds that ).These two characteristics are specially appealing in the context of mortality rates since there can be outliers (wars, pandemics, etc.) and often the log transformation is the function of interest.

Quantile regression in presence of heteroscedasticity
The constant error variance (homoscedasticity) assumption of the original Lee-Carter model in (1) has been criticized as unrealistic because of the observed pattern of the mortality rates showing a different variability at different ages [Brouhns et al., 2002].Quantile regression models are particularly interesting in the presence of heteroscedasticity.
Let us assume that Y can be represented as for a certain function f , where f (X, θ) can be thought of as the conditional mean, which depends on regressors X and the vector of coefficients θ, and ε is the error.Assume that errors are independent and identically distributed, then the conditional quantile function (2) can be written as where the vector of coefficients θ does not depend on α.In that case, the quantile functions are simply a location displacement of one another [Koenker, 2005].Departing from the assumption of identically distributed errors, we suppose now that Y is represented as Y = f (X, θ) + σ(X)ε, where σ(X) is a conditional scale that depends on X.The conditional quantile function would be, in that case, The more general expression of the quantile function (2) includes both of these situations, since θ α is not restricted to be equal for all α.
Summarizing, quantile regressions can be used in those contexts that the conditional distribution of Y depends on the regressors in complex ways in which the error is not independent of regressors.Asymptotic theory for quantile linear regression models under independent but not identically distributed errors is well developed [Koenker and Bassett, 1982;White and Kim, 2002;Parente and Santos Silva, 2016].In fact, a number of tests for heterocedasticity in linear regression models are based on the analysis of variations of estimated coefficients of quantile linear regressions [Koenker and Bassett, 1982;Wilcox and Keselman, 2006;Machado and Silva, 2011].In the context of quantile non-linear regression models, the pairs-bootstrap estimator of the covariance matrix of coefficient estimates is widely used [Buchinsky, 1995].
where superscript α indicates the α-quantile associated with the parameters.As in the case of the Lee-Carter mean regression model, to overcome the lack of identifiability, two constraints are established, namely x b α x = 1 and t k α t = 0.
Under these constraints, the interpretation of the coefficients is equivalent to that in the traditional Lee-Carter regression model but here adapted to quantiles.Specifically, • a α x indicates the expected α quantile of log(m x,t ) for age x over time.

• b α
x indicates variations to the expected α quantile for age x in relation to the time trend.
• k α t is the time index that captures the trend over time of the α-quantile for all ages.

Model fitting
In the same way as the mean is the value that minimizes the sum of squared deviations, the median minimizes the sum of absolute deviations.A widely used strategy to estimate the parameters of the mean regression involves minimizing the squared errors by means of least squares techniques.Similarly, parameters of quantile regression are estimated by least absolute techniques.
The problem of finding the α sample quantile of an empirical distribution τ α may be written as a problem of minimization of a loss function ρ α as follows, where n is the sample size.The loss function is defined as ρ α (r i ) = r i (α − I r i <0 ) for α ∈ (0, 1), where I r i <0 is an indicator function such that I r i <0 = 1 if r i < 0 and zero otherwise.Note that in the median case (α = 0.5), this optimization problem is equivalent to min n i=1 |y i − τ 0.5 |.Let us consider the LC conditional α quantile model defined in expression (3).The set of parameters is estimated by solving the following optimization problem min where The set of parameters θ α is a vector of dimension p = 2 + 2 • ω + T and the sample size n is equal to n = (1 + ω) • T .The quantile regression optimization problem can be represented as min where u i = g i (θ α ) if g i (θ α ) > 0 or 0 otherwise, and Two alternative strategies can be adopted to estimate parameters based on the interior point method optimization algorithm explained in Koenker and Park [1996].Under the first strategy (Method A), coefficients are estimated by means of an iterative process of a sequence of linear optimization problems.The interior point method algorithm for linear problems is used in each sequence of the optimization process.Under the second strategy (Method B ), the interior point method algorithm for non-linear problems is applied directly.The two optimization methods are described in Appendix A. The algorithm used in both strategies is based on Maketon method [Koenker and Park, 1996].Alternative algorithms may be used to resolve quantile regression problems by interior point methods, some of them implemented in standard software.For instance, the algorithms implemented in the R package quantreg for quantile linear regressions are based on Mehrotra's Predictor-Corrector method [Portnoy and Koenker, 1997;Koenker, 2019].In the case of non-linear quantile regressions, the algorithm implemented in the function nlrq is based on the Meketon algorithm, but this function can only be applied to full-rank design matrices (which is not the case of the LC model framework).
To conclude, a number of studies has recommended to smooth coefficient estimates of LC conditional mean regression models [Currie et al., 2004;Renshaw and Haberman, 2003].Smoothing techniques may be applied in the context of quantile regression models.For example, Koenker et al. [1994] and Yu and Jones [1998] used local lineal fitting and splines, respectively.Wang et al. [2009] analyzed varying coefficient models for dependent data.The optimization methods described in Appendix A can be adapted to smooth coefficients estimates.Smoothing methods can be applied after each linear minimization problem is optimized (Method A) or after each iterate (Method B).In this article coefficient estimates are not smoothed in order to evaluate the performance of the two optimization methods without the impact of exogenous factors, as the application of smoothing techniques.Estimation techniques of quantile non-linear models are not as well studied as compared to methods of estimating conditional mean non-linear models.So, the analysis of the parameter estimation methods of the LC quantile model is of great interest itself.In addition, LC quantile coefficient estimates can be easily compared with LC mean coefficient estimates.

Forecasting the quantile mortality index
This section follows the same approach as that taken in the case of LC conditional mean models [Lee and Carter, 1992;Brouhns et al., 2002;Renshaw and Haberman, 2006].The dynamics of the mortality quantiles are captured by the set of estimated mortality indexes (k α 1 , ..., k α T ).If the goal of the LC quantile model is to project mortality quantile rates into the future, then time-series techniques are required to model quantile mortality indexes.
Our objective in this section is to find the appropriate univariate ARIMA model for the mortality index k α t , t = 1, .., T .In this approach, k α t follows an ARIMA(p,d,q) with drift, so that where ∆ is the difference operator, d the number of non-seasonal differences needed for stationarity and η t is a Gaussian noise process with variance σ η .The set of parameters is (δ 0 , φ 1 , .., φ p , δ 1 , ..., δ q ), with δ 0 the drift parameter, p the number of autoregressive terms and q the number of lagged forecast errors.

Application to the Spanish data
In what follows the LC quantile model is fitted to Spanish mortality data to illustrate its application.The number of deaths observed, exposures and central mortality rates for the Spanish population by gender were obtained directly from the Human Mortality Database HMD [2019] for the observation period 1908-2016 and for ages 0-100.

Mean vs. Median
Mean and median give information of the center of the distribution.While the mean minimizes the sum of squared errors, the median minimizes the sum of absolute errors.So, LC quantile regression coefficient estimates are expected to be close to the LC mean regression coefficient estimates in the case of the median.
We compare the performance of the LC conditional mean model with that of the LC conditional median (50%-quantile) model.Three alternatives LC mean models were considered: the original LC model introduced by Lee and Carter [1992](LC Mean); the LC mean model with Poisson residuals (LC Mean Poisson); and, the LC mean model with Gaussian residuals (LC Mean Gaussian).The LC Mean model was estimated by SVD and k t 's were recalibrated in a second step.The LC Mean Poisson and LC Mean Gaussian models were estimated by maximum likelihood of the generalized non-linear model [Turner and Firth, 2018].Models were separately calibrated by gender.The sum of absolute errors and the sum of squared errors of the fitted models are reported in Table (1).
Almost identical results were obtained when using the two LC median model estimation methods.A comparison with the LC mean models shows that the sum of absolute errors is lower for the two LC median models estimated.Note that the performance of the LC median models is reasonably good, even when the sum of squared errors is analyzed.More specifically, only the LC Mean Gaussian model presents a lower sum of squared errors.In contrast, the original LC Mean model and LC Mean Poisson model show the poorest performance in all cases, albeit that the latter performs slightly better than the former.
The section that follows reports the results of the LC α-quantile models for the Spanish male population.Results for the Spanish female are in Appendix B.

Coefficient estimates and α quantile mortality rates
Coefficient estimates of the LC α-quantile model are shown in Figure (1).LC quantile models are estimated with Method B for the set of α's equal to (0.10, 0.50, 0.90).Results are compared with estimates for the LC conditional mean model (with Poisson error distribution).Although this is not the LC mean model that performs best, there are two reasons that justify its selection.First, the LC model with Poisson distributed errors is widely used and, second, the performance of the LC mean model with Gaussian errors is similar to that of the LC median model.Thus, the comparison between LC conditional quantile and conditional mean models is enriched by considering the Poisson error distribution.

Note that a α
x increases with α for all ages.This result is expected as a α x collects the estimated α quantile for age x over time.The largest gap between a 0.10 x and a 0.90 x is ob-  served in the late twenties-early thirties.Similar values are observed for a 0.50 x and a mean x .In the case of the b α x and k α t , a clear pattern of variations of coefficients with α is not observed.In both cases, the estimates seem to present higher oscillations when α is 0.90.
Let us focus on time-varying index estimates.Three peaks are evident over time for k 0.90 t , i.e., in the late 1910s (Spanish flu epidemic), around 1940 (Spanish civil war), and in the late '80s and early '90s (mortality associated with HIV, [CNE, 2011]).The first two peaks are also appreciated for the other time index estimates, particularly for k 0.50 t .Note that the trend of the time-varying index provides information of the variation of the (log) mortality rate over time.A similar decreasing trend (degree of inclination) is observed in all k α t 's.Thus, in relative terms, the reductions over time of all three α quantile (log) mortality rate estimates are similar.However, in those periods that mortality rises, the relative increase is more accentuated for high quantile (log) mortality rate estimates.Comparing k α t 's with the mean time-varying index k t , the latter seems to show a slightly faster decreasing trend in the last two decades of the time series.
Quantile (log) mortality rate estimates for four different years (1908, 1975, 2000, and 2016) are reported in Figure (2).It is apparent that observed mortality rates lay between the estimated 10 and 90% quantile rates in most cases.Interestingly, in 2000 and 2016, the estimated mean rates were lower than the observed rates for ages below 40.The distance is particularly notable for ages below 20 and for calendar year 2016.

Analysis of heteroscedasticity
The 80%-confidence intervals of the estimated mean mortality rates at ages 20, 40, 60, 80 and 100 are reported in Figure (3).Confidence intervals were based on 5,000 semiparametric boostrap samples of the Lee-Carter model.The semiparametric bootstrap technique is described in Brouhns et al. [2005].We use the function bootstap.fitStMoMofrom the  package StMoMo [Villegas et al., 2018].Figure (3) shows that parameter uncertainty at 80% confidence level has a moderate impact.Confidence intervals could be estimated to a higher degree of confidence level to see the impact of parameter uncertainty.Additionally, it is important to consider all sources of uncertainty in mortality projection models.In Section (5.5) it is shown the high impact of the uncertainty arising from projections of the time-varying index.
The estimated median, 10 and 90% quantile mortality rates at ages 20, 40, 60, 80 and 100 are reported in Figure ( 4).Estimated quantile mortality rates seem to capture better the shape of the observed mortality data.Now, we can see that most of observations lay on the inter-decile area.The width of the inter-decile range varies with age and calendar year.For example, at age 100 a declining trend is observed for the the width of the range between the 10 and 90% quantile rates.This declining trend is not observed for the rest of ages.At age 20 the inter-decile range increases from the middle of 80's onwards.The peaks between the estimated 10 and 90% quantile mortality rates are observed in the late 1910s, around 1940, and in the late '80s and early '90s (flu, civil war and HIV).Another interesting result is that the distance of the median with the 10 and the 90% quantiles is not symmetric and varies over time.

Analysis of residuals
The diagnostic of residuals is carried out by estimating raw residuals, i.e., g i ( θα ) = y i − f (x i , t i , θα ) for i = 1, . . ., n. Quantile residuals are then compared with the raw ), the performance of the the LC quantile model is good for all ages and α's.The percentage of observations below the α quantile values remains very stable around α for all ages.By contrast, some oscillations are observed for the mean of residuals for young ages in the LC mean model.Remember that the expected mean of residuals should be null in the case of the LC mean model.
In general, more oscillations are observed when residuals are analyzed by the calendar year (Figure (5) [right]), especially for the LC median model and mean models.In both cases, results show an increasing trend from the 1960s onwards.In the case of the LC 10%-quantile model, the percentage shows a slightly increasing trend from the 2000s onwards.Finally, the period of greatest instability for the 90% quantile is between the 1940s and 1980s.These results seem to indicate that the LC model framework does not fully capture the trend over time for the entire period under observation.To improve the model's performance, a potential solution could involve shortening the observation window.Another option would be to consider a more complex model framework.

Forecasting mortality indexes and α quantile mortality rates
Forecasting the time index requires time-series techniques.First, we need to select the best time series process to model k α t in the framework of the expression (6).We fit ARIMA models with drift to data and select the best model for each time series based on AIC criteria.Table (2) shows the selected ARIMA model for each time index and coefficient estimates.Poisson LC mean model fitted to the Spanish male population for ages 0-100 and the period 1908-2016.

LC
Taking the estimated coefficients as shown in Table (2), each time index is projected for 34 years in the future.Figure (6) shows the projections of each time index with a 95% confidence interval.Note that randomness -that is, the variance of the noise processincreases with α as indicated in Table (2).Figure (6) shows that the width of the confidence intervals of k α t 's increases with α.
Finally, Figure (7) shows predicted α quantile log mortality rates and mean log mortality rates for 2025 (Figure ( 7) [left]) and 2050 (Figure (7)[right]).Mortality is expected to improve, especially for the young and middle-aged (under 50 years).In contrast, improvements are more moderate among those of a more advanced age.Note that in both cases, the inter-quantile range between 10 and 90% reaches a maximum at ages below 10, the distance thereafter decreasing with age.In both cases the predicted mean mortality rate is clearly lower than the predicted median for ages under 50 years.The difference is particularly notable at young ages, where the predicted mean mortality rate is lower than the predicted 10% quantile.This result reflects the fact that the conditional mean only gives a partial portrait of the mortality rate distribution and emphasizes the need to have information for the whole conditional distribution of the mortality rate, especially that of the tails.

Illustration: annuities
The application of quantile mortality rates is illustrated in the valuation of term immediate annuities.Consider an annuity of 1 per year payable annually at the end of the year, conditional on the survival of the annuitant to the payment dates.We compute the present value (PV) of a term immediate annuity.We follow the notation given by Dickson et al. [2009].The present value of a term immediate annuity is represented as a x:n| , where x is the age of the annuitant and n the number of terms.Under this annuity payments of 1 are made at times k = 1, 2, ..., n, conditional on the survival of the annuitant.The conditional probability that annuitant at age x in the year t survives at least age x + 1 is computed as p x,t = e −mx,t .Quantile and mean mortality rates are used in the valuation of term immediate annuities.It is worth mentioning that Q 1−α (p x,t ) = e −Qα(mx,t) , while in general [mx,t] .
Figure (8) shows the present value of the term immediate annuities.Actuaries are primarily concerned with the trends in adult ages in extreme scenarios.The a x:n| is computed for annuitants at ages from 67 to 100 years and the last payment is made when the annuitant reaches the age of 101, n = (101 − x).The valuation of the PV is made at the beginning of the year 2017.Mortality rates used in the computation of the PV are quantile mortality rate forecasts at levels 5%, 50%, 95%.Results are compared with the PV when mean mortality rate forecasts are applied.Interest rates used in the valuation are spot rates provided by the European Central Bank.Negative interest rates were replaced by zero.

Remark
. Aggregation properties of quantiles have been extensively studied in the actuarial field [Embrechts et al., 2009;Emmer et al., 2013;Dhaene et al., 2006].Let us consider two mortality rates, m x,t and m x+1,t+1 .Under comonotonicity (perfect positive dependence between m x,t and m x+1,t+1 ), it holds that ).In case that the subadditivity property is satisfied, defined as ) is the maximum possible value of the quantile of the aggregation of m x,t and m x+1,t+1 .Consequently, Q 1−α ( 2 p x,t ) = e −(Qα(mx,t)+Qα(m x+1,t+1 )) is the minimum (1 − α) quantile of the probability of the annuitant at age x and time t survives at least age x + 2. Formally, the subadditive property of quantiles cannot be generalized to all contexts, but the violation of subadditiviy is infrequent in practice [de Vries et al., 2005].Some interesting cases satisfying the subadditivity property of quantiles include elliptical distributions, Archimedean survival dependence structure or independent and identically distributed and positively regularly varying random variables.

Conclusions
In this article, we have proposed a quantile non-linear regression model to predict conditional quantile mortality rates.We have adapted the Lee-Carter model framework to the context of conditional quantiles and we describe two strategies for estimating the regression coefficients.We show that the LC quantile model fully characterizes the entire conditional distribution function of the mortality rates, unlike the LC mean model which only provides the expected mean value.
The application of the methodology is illustrated by means of an example using Spanish mortality data.The performance of the LC conditional median model was compared with that of the LC conditional mean model.The sum of absolute errors and the sum of squared errors of the fitted models was analyzed, and the LC median model was found to perform well even when the sum of squared errors was compared.The residuals of the LC quantile models were evaluated conditioned to age and the calendar year.The results when conditioned to age, especially for males, were highly satisfactory.However, greater instability was observed when conditioned to the year.The performance of the LC mean model was no better and, on occasions, was markedly worse, indicating that additional time index terms may be required to model Spanish mortality data.

A Appendix. Optimization strategies A.1 Sequence of linear optimization problems
Method A draws inspiration from Wilmoth [1993] and Renshaw and Haberman [2006] who propose an iterative process to estimate sequentially the parameters of the traditional LC regression model.Similarly, in this section we convert the non-linear optimization problem (4) in a sequence of three minimization linear problems.The iterative process ends when the objective function no longer improves.The algorithm can be stated as follows:
At steps (i), (ii) and (iii), a linear minimization problem has to be optimized.The quantile linear regression optimization problem in matrix notation can be represented as min θ α ∈R p ,u≥0,v≥0 where Y = (log(m 0,1 ), ..., log(m ω,T )), θ α is the set of parameters to estimate and X is an observed matrix of binary variables (obviously, θ α and X will depend on the step in the algorithm).The dual version of this optimization problem is max To solve the linear problem (A.2) the Meketon algorithm described in Koenker and Park [1996] is applied.We start with any feasible point d 0 = {d i } i=1,..,n in the interior of the constraint set: The parameter ν is bounded in [0, 1].Koenker and Park [1996] recommend using a value close to one and propose setting ν = 0.97.We use this value in our application.

B Appendix. Analysis for Spanish female population
Results for Spanish female data are shown in the Appendix.According to AIC criteria, the selected ARIMA process to model k 0.10 t was ARIMA(2,1,0) with drift.In the rest of cases, including the LC mean regression, the best fit was observed for an ARIMA(1,1,0) model with drift.

Figure 1 :
Figure 1: Coefficient estimates for LC α-quantile models fitted to the Spanish male population for ages 0-100 and the period 1908-2016.Green dash-dotted lines represent coefficient estimates of the LC 10-quantile model, black dotted lines represent coefficient estimates of the LC median model, and red dashed lines represent coefficient estimates of the LC 90-quantile model.Black dashed lines represent coefficient estimates of the Poisson LC mean model.

Figure 2 :Figure 3 :
Figure 2: Quantile (log) mortality rate estimates for years1908, 1975, 2000 and 2016  from LC αquantile models fitted to the Spanish male population.Green dash-dotted lines represent 10-quantile (log) mortality rate estimates, black dotted lines represent median (log) mortality rate estimates, and red dashed lines represent 90-quantile (log) mortality rate estimates.Black dashed lines represent mean (log) mortality rate estimates with the Poisson LC mean model and black thick lines represent observed (log) mortality rates.

Figure 4 :Figure 5 :
Figure4: Quantile mortality rate estimates at ages 20, 40, 60, 80 and 100 from the LC α-quantile model fitted to the Spanish male population.Dots show observed mortality rates and solid grey lines show the corresponding fitted median.Black (short) dashed lines represent 10% quantile mortality rate estimates and black (large) dashed lines represent 90% quantile mortality rate estimates.

Figure 6 :
Figure 6: Forecasting k α t indexes of LC α-quantile models and k t index of Poisson LC mean model for the period 2017-2050.Blue lines represent point forecasts and shaded area represent 95% prediction intervals.

Figure 7 :
Figure7: Forecasting α quantiles of (log) mortality rates for years 2025 and 2050 from Lee-Carter αquantile models fitted to the Spanish male population.Green dash-dotted lines represent 10-quantile predictions, black dotted lines represent median predictions, and red dashed lines represent 90-quantile predictions.Black dashed lines represent mean predictions with the Poisson LC mean model.
Figure B.1: Coefficient estimates for LC α-quantile models fitted to the Spanish female population for ages 0-100 and the period 1908-2016.Green dash-dotted lines represent coefficient estimates of the LC 10-quantile model, black dotted lines represent coefficient estimates of the LC median model, and red dashed lines represent coefficient estimates of the LC 90-quantile Black dashed lines represent coefficient estimates of the Poisson LC mean model.
Figure B.2: Quantile (log) mortality rate estimates for years1908, 1975, 2000 and 2016  from LC quantile models fitted to the Spanish female population.Green dash-dotted lines represent 10-quantile (log) mortality rate estimates, black dotted lines represent median (log) mortality rate estimates, and red dashed lines represent 90-quantile (log) mortality rate estimates.Black dashed lines represent mean (log) mortality rate estimates with the Poisson LC model and black thick lines represent observed (log) mortality rates.
Figure B.3: Mortality female rate estimates and 80% confidence interval for the Poisson LC mean model (left) and α quantile mortality female rates estimates for LC α-quantile models (right) at ages 20, 40, 60, 80 and 100 and α = (0.10, 0.50, 0.90).Dots show observed mortality rates and solid grey lines show the corresponding fitted rates.

Figure B. 5 :
Figure B.5: Forecasting k α t indexes of LC α-quantile models and k t index of Poisson LC mean model for the period 2017-2050.Blue lines represent point forecasts and shaded area represent 95% prediction intervals.

Table 1 :
Comparison of LC conditional median models and LC conditional mean models fitted to the Spanish male population for ages 0-100 and the period 1908-2016.LC Mean represents the original LC mean model introduced by Lee and Carter [1992], LC Mean Poisson represents the LC mean model with Poisson residuals and LC Mean Gaussian represents the LC mean model with Gaussian residuals.LC Median (Method A) and LC Median (Method B) represent the LC conditional median models estimated by method A and method B.

Table 2 :
Coefficient estimates of the ARIMA model for the time index of LC α-quantile models and . For any natural number h > 1, the conditional survival probability h p x,t is Figure 8: Present value (PV) of a term immediate annuity at ages from 67 to 100 years and the last payment at age of 101.Dash-dotted line represents PV's based on 5-quantile mortality rate predictions, dotted line represents PV's based on median mortality rate predictions, and dashed line represents PV's based on 95-quantile mortality rate predictions.Solid line represents PV's based on mean mortality rate predictions with the Poisson LC mean model.computed as h p x,t = e −(mx,t+m x+1,t+1