Forecasting short-term mortality trends using Bernstein polynomials

Abstract Mortality data play an important role on the fields such as health, epidemiology and national planning. Most mortality models mainly focus on providing a perfect fitting, to the detriment of an exact forecasting result. In this paper, we fit the Bernstein polynomial to mortality data based on maximum likelihood based inference through the simulated annealing method. The proposed method utilizes the derivative of Bernstein polynomials to describe the pattern of mortality rates. The asymptotic behavior of the proposed model is also given on general results. The performance of the proposed method is evaluated by simulated examples and illustrated through applications to datasets provided from the Human Mortality Database (www.mortality.org). The simulated examples and real data analysis show that the proposed approach is quite satisfactory in forecasting the short-term mortality trends.


Introduction
Owing to the continuous longevity improvement during the past few decades, modeling the pattern of human mortality data has become an important issue confronted by governments, insurance companies and superannuation funds (Bohk-Ewald and Rau 2017; Doukhan et al. 2017;Li, O'hare, and Vahid 2017;Plant 2009). Therefore, a variety of methods for mortality modeling and forecasting have been introduced in the literature and have made great advances during the last decades (Bamidele and Adejumo 2014;Doukhan et al. 2017).
The representation of mortality data through a parametric model has aroused much attention for over a century (Debon, Monte, and Sala 2005;Dellaportas, Smith, and Stavropoulos 2001). Predicting mortality rates and estimating their sampling variance also have attracted much attention in the past few years (Lin, Wang, and Tsai 2015). Lee and Carter (1992) introduced the stochastic mortality model to model and forecast the age and period effects. The Lee-Carter model provided the fundamental set up that are the most commonly used in the current mortality models for mortality modeling and forecasting. Various extensions of the Lee-Carter model have been introduced. One of them was introduced for considering individual heterogeneity in each age-period cell (Li, Hardy, and Tan 2009). A modification of the Lee-Carter method was proposed for modeling the rotation of age patterns for long-term projections (Li, Lee, and Gerland 2013). A comprehensive overview of the literature of the extension of the Lee-Carter model can be found, e.g., in Bohk-Ewald and Rau (2017).
On the other hand, some methods for mortality modeling and forecasting are more flexible than the Lee-Carter model (Bohk-Ewald and Rau 2017). For example, Eilers (2004, 2006) developed the spline-based approach based on the twodimensional mortality data to estimate and predict the mortality surfaces. Camarda (2015) carried out the spline-based approach Eilers 2004, 2006) along with a two-dimensional version in the R package MortalitySmooth. Plant (2009) proposed a stochastic mortality model that combined the advantages of the existing methods, the Lee-Carter model and the spline-based approach of Currie, Durb an, and Eilers (2004. Li, O'hare, and Vahid (2017) proposed a two-dimensional Legendre orthogonal polynomials model to improve the overall accuracy of the future mortality projection. A proposition of generalized stochastic Milevsky-Promislov mortality models was suggested for estimating mortality rates (Sliwka and Socha 2018). In the recent years, Bayesian methods have gradually become an important method in mortality modeling and forecasting (Bohk-Ewald and Rau 2017). These Bayesian methods used the Markov chain Monte Carlo algorithm to estimate and predict smooth mortality curves (Alexander, Zagheni, and Barbieri 2017;Dellaportas, Smith, and Stavropoulos 2001;Goicoa et al. 2019;Li and Lu 2018).
However, most existing mortality models mainly focus on providing a perfect fitting, to the detriment of an exact forecasting result (Li, O'hare, and Vahid 2017). We propose a flexible maximum likelihood approach based on the shape-restricted regression model with the two-dimensional Bernstein polynomials for forecasting the short-term mortality trends. An easy-to-implement algorithm is proposed for specifying Bernstein coefficients for constructing the shape-restricted regression model. A simulated annealing method (Ross 2013) based on a Markov chain Monte Carlo (MCMC) sampling algorithm is introduced and connected with the Akaike information criterion (Akaike 1974) (AIC) to perform maximum likelihood based inference. Some fundamental theorems on the asymptotic behavior of the proposed model are also given on general results. A significant ground for using a polynomial in the form of Bernstein in studying the mortality trend is that the linear combinations of the Bernstein polynomials have much of their geometry dictated by their Bernstein coefficients (Chang et al. 2007;Chang et al. 2005). Their derivatives have simpler form than the B-splines. These derivatives of these polynomials describe the pattern of mortality rates. The appropriate approximation theories (Propositions 1-3) based on the Bernstein polynomials are provided for addressing computational problems. Therefore, Bernstein polynomials are more suitable than the polynomials in a different form (e.g., the B-splines) in studying the mortality trend.
The remainder of the paper is organized as follows. Section 2 introduces the proposed mortality model and studies its asymptotic properties. Section 3 uses simulated annealing method to perform maximum likelihood based inference. Section 4 examines the proposed model and its alternatives by simulations. Section 5 illustrates the method on data extracted from the Human Mortality Database (www.mortality.org). Section 6 presents a brief summary of the results and discussion.

The model specification
Let Mðx, yÞ denote the unobserved death rate within a group of population aged x in calendar year y. Let Dðx, yÞ and Eðx, yÞ stand for the number of deaths and the number of people living in calendar year y at age x, respectively. We start with a simple model and we assume that the number of deaths in calendar year y of people having age x is distributed as a Poisson random variable with mean lambda, Eðx, yÞMðx, yÞ, namely, Dðx, yÞ $ PoissonðEðx, yÞMðx, yÞÞ, where the natural logarithm of the unobserved death rate, ln ðMðx, yÞÞ, is modeled using the two-dimensional Bernstein polynomials in Equation (1) in the form ln ðM a ðx, yÞÞ ¼ x À s 11 s 12 À s 11 u j, m y À s 21 s 22 À s 21 , i.e., M a ðx, yÞ ¼ e lnðM a ðx, yÞÞ f g : Precisely, we assume that Dðx 1 , y 1 Þ, Á Á Á , Dðx K, y K Þ are independent Poisson random variables. The density of Dðx k , y k Þ, k ¼ 1, Á Á Á , K, at the corre- where the index k represents the number of the observed data and the number K represents the total number of the observed data. Given the observed data ððx k , y k Þ, z k , Eðx k , y k ÞÞjk ¼ 1, Á Á Á , K È É , the likelihood function of the Bernstein parameters ðn, m, aÞ is given by where the polynomial orders ðn, mÞ and the Bernstein coefficients a ¼ ða 0, 0 , a 0, 1 , Á Á Á , a n, m Þ can be estimated by the likelihood function in Equation (3). The mortality rate of Mðx, yÞ at a fixed calendar year of y is increasing with the age of x increasing, which in turn implies that the natural logarithm of the mortality rate of ln ðMðx, yÞÞ given a fixed value of y has an increasing shape restriction for x on the interval ½s 11 , s 12 : An increasing shape restriction of the Bernstein polynomials is thus used to fit ln ðMðx, yÞÞ by Equation (2) where @M a ðx, yÞ @a ¼ u i, n x À s 11 s 12 À s 11 u j, m y À s 21 s 22 À s 21 e X n i¼0 X m j¼0 a i, j u i, n x À s 11 s 12 À s 11 u j, m y À s 21 s 22 À s 21 ( ) , I À1 ðaÞ is the inverse function of the Fisher information matrix IðaÞ ¼ EðÀ@ 2 l=@a i 0 , j 0 @a i, j Þ for i ¼ 0, 1, Á Á Á , n, j ¼ 0, 1, Á Á Á , m, i 0 ¼ 0, 1, Á Á Á , n, j 0 ¼ 0, 1, Á Á Á , m, and l ¼ ln ðLðn, m, aÞÞ represents the log-likelihood function. A detail of the derivation of the variance estimator of the mortality at age x in year y is presented in Section 1 in Supplementary Material.

The Bernstein density
ðfng Â fmg Â B nm Þ: Let p denote a probability measure on B termed as the joint density of the polynomial orders ðn, mÞ and the Bernstein coefficients a ¼ ða 0, 0 , Á Á Á , a n, m Þ: The p is regarded as a Bernstein density and is defined by pðfng Â fmg Â B nm Þ ¼ pðajfng Â fmg Â B nm ÞpðnÞpðmÞ where pðajfng Â fmg Â B nm Þ is the conditional density of p conditionally on B nm and P 1 n¼1 pðnÞ ¼ P 1 m¼1 pðmÞ ¼ 1: Proposition 1 provides sufficient conditions under which the Bernstein polynomial ln ðM a ðx, yÞÞ in Equation (2) given a fixed calendar year of y is a monotone increasing function of the age of x on ½s 11 , s 12 , i.e., ln ðM a ðx, yÞÞ given a fixed calendar year of y has an increasing shape restriction for age x on ½s 11 , s 12 : Proposition 1. If À1 < R 1 a 0, j < a 1, j Á Á Á a n, j R 2 < 0, for every j ¼ 0, Á Á Á , m, then R 1 ln ðM a ðx, yÞÞ R 2 and @ ln ðM a ðx, yÞÞ=@x ! 0 for every x 2 ½s 11 , s 12 and y 2 ½s 21 , s 22 , where @ ln ðM a ðx, yÞÞ=@x ¼ ðn=ðs 12 À s 11 ÞÞ P nÀ1 i¼0 P m j¼0 ða iþ1, j À a i, j Þu i, nÀ1 ððx À s 11 Þ=ðs 12 À s 11 ÞÞ u j, m ððy À s 21 Þ=ðs 22 À s 21 ÞÞ: Let I be the set of all two-dimensional continuous functions defined on ½s 11 , s 12 Â ½s 21 , s 22 with the values between R 1 and R 2 and with an increasing shape restriction for x on ½s 11 , s 12 : Then, using Proposition 1 combined with the Bernstein-Weierstrass type approximations, we have Proposition 2 that complements Proposition 1.

Maximum likelihood inference
We first explain how to establish the Bernstein density pðn, m, aÞ described in Section 2.3. The likelihood function in Equation (3) is a non-linear form with the multi-dimensional parameter space considered. The simulated annealing is a stochastic algorithm for approaching the solution of a global optimization problem in a high-dimensional space ( Cern y 1985; Kirkpatrick, Gelatt, and Vecchi 1983). We then illustrate how to use the simulated annealing method to obtain the MLEs of the Bernstein parameters ðn, m, aÞ by a MCMC sampling algorithm.

Specifying the Bernstein density p
To generate a sample from the Bernstein density p, the values of the polynomial orders n and m are respectively chosen from the positive integers f1, 2, 3, Á Á Ág: Given a set of the polynomial orders n and m, the Bernstein coefficients a ¼ ða 0, 0 , Á Á Á , a n, m Þ are generated through the following processes.
The conditional density pðajfng Â fmg Â B nm Þ is defined by assuming that fa i, j ji ¼ 0, Á Á Á , n, j ¼ 0, m; i ¼ 0, n, j ¼ 1, Á Á Á , m À 1g have respective densities, which do not change the polynomial orders, n and m. Elements of fa i, j ji ¼ 0, Á Á Á , n, j ¼ 0, mg have densities of the form UniformðL i, j , U i, j Þ, i ¼ 0, Á Á Á , n, j ¼ 0, m, respectively, and these elements satisfy the restriction a 0, j a 1, j Á Á Á a n, j , for , j be the order statistics of a random sample from Uniformða 0, j , a n, j Þ, for j ¼ 1, Á Á Á , m À 1: Then, ðn, m, a 0, 0 , Á Á Á , a n, m Þ 2 B is a sample generated from the shape-restricted density p by the above processes. The selection of the parameters of these distributions depends on the datasets used for analyzing. A detail of the choice of the parameters of these distributions (e.g., L i, j and U i, j ) will be illustrated in the simulation studies and the practical data analysis example.

Estimating the Bernstein parameters a for fixed n and m
With the shape-based density p illustrated in Section 3.1, we use the simulated annealing method described in Ross (2013, pp. 290-292) to obtain the MLEs of the Bernstein parameters ðn, m, aÞ by a MCMC sampling algorithm. Let B ðnmÞ ¼ fðn, m, a 0, 0 , Á Á Á , a n, m Þjða 0, 0 , Á Á Á , a n, m Þ 2 B nm g: Let h ðrÞ ¼ ðn, m, a 0, 0 , Á Á Á , a n, m Þ 2 B ðnmÞ be the current state. A new state h ðrþ1Þ 2 B ðnmÞ is upgraded from the current state h ðrÞ 2 B ðnmÞ according to the following steps.

Simulation study
We now conduct the simulation study to explore the performance of the proposed method. To make the simulation study more relevant to the practical examples, the data in this simulation study are generated based on the Swedish mortality data provided from Human Mortality Database. Here we focus on the mortality rates among the middle-and old-aged adults. We consider the mortality rates of Sweden in the period from 1990 to 2010 at age from 50 to 79. Let Mðx, yÞ denote the mortality rate of Sweden at age x in year y and let Eðx, yÞ denote the population number of Sweden at age x in year y.
Using the Bernstein polynomials ln ðM a ðx, yÞÞ in Equation (2) based on the observed samples, the mortality rates at ages 50-79 in years 1990-2010, we intend to predict the mortality trend in the period 1990-2011 at the age of 80 and to forecast the mortality trend in the year of 2011 at ages 50-80. Here our proposed method focuses mainly on forecasting the short-term mortality rates.
The estimate of the mortality rates,M, is computed by Equation (4). The performance ofM is evaluated by the L 1 -norm, sup-norm and mean square error (MSE) of jM À Mj, where j Á j stands for the absolute value and M represents the mortality rate of Sweden. Here the L 1 -norm, sup-norm and mean square error (MSE) of jM À Mj are defined by Ð Ð jMðx, yÞ À Mðx, yÞjdxdy, sup ðx, yÞ jMðx, yÞ À Mðx, yÞj and Ð Ð ðMðx, yÞ À Mðx, yÞÞ 2 dxdy, respectively. In order to examine the effectiveness of the proposed method, the simulation studies include the comparison between the proposed method and the existing methods, the spline-based approach Eilers 2004, 2006;Eilers, Currie, and Durb an 2006) and the Lee-Carter model (Lee and Carter 1992), which are executed by the R packages MortalitySmooth (Camarda 2015) and demography (Hyndman 2017), respectively. The comparison is conducted based on the 200 simulation replicates. Table 1 shows the averages of the resulting 200 L 1 -norms, sup-norms and MSEs of jM À Mj from the proposed method, the spline-based approach and the Lee-Carter model, respectively. We note that the Lee-Carter model carried out by the R package demography (Hyndman 2017) cannot forecast the mortality rate at the age of 80 in the year of 2011, when we only consider the observed samples including the mortality rates at ages 50-79 in years 1990-2010. To present the consistent results obtained from the Lee-Carter model carried out by the R package demography (Hyndman 2017), Table 1 doesn't present the forecasting result of the mortality rate at the age of 80 in the year of 2011 based on the Lee-Carter model by the R package demography (Hyndman 2017). Table 1 shows that, in forecasting the mortality rates at age 80 in the year of 2011, the averages of the resulting 200 L 1 -norms and MSEs of jM À Mj based on the proposed method are 6.675 Â 10 À4 and 4.456 Â 10 À7 , respectively, which are larger than that based on the spline-based approach, 2.354 Â 10 À4 and 5.542 Â 10 À8 , respectively. However, in predicting the mortality rates at ages 50-79 in the year of 2011, the average of the resulting 200 L 1 -norms and MSEs of jM À Mj based on the proposed method, 2.650 Â 10 À4 and 1.435 Â 10 À7 , respectively, are less than that based on the spline-based approach, 2.830 Â 10 À4 and 1.786 Â 10 À7 , respectively. A similar result is observed with predicting the mortality rates at the age of 80 in years 1990-2010. On the other hand, in predicting the mortality rates at ages 50-79 in the year of 2011, the performance of the Lee-Carter is the best, among the three competing methods, when the measure of the sup-norm is considered.   The averages of the 200 confidence intervals (95%) for forecasting the mortality rates at the age of 80 in years 1990-2010 and the mortality rates at ages 50-79 in the year of 2011 are presented in Table S1 and Figures S1 and S2 in Section 5 in Supplementary Material. The comparison of the proportion of the 200 confidence intervals (95%) containing the true values of the mortality rates from the proposed method, the splinebased approach and the Lee-Carter model based on the simulation study is shown in Table S2 and Figure S3 in Section 5 in Supplementary Material.
From Tables S1 and S2 and Figures S1-S3, we observe that in forecasting the mortality rates at age 80 in years 1990-2010, the Lee-Carter model generally has wider confidence intervals (95%) and higher proportion of the confidence intervals (95%) containing the true values of the mortality rates based on the 200 simulation runs, in comparison with the proposed method and the spline-based approach. However, the Lee-Carter doesn't have better performance than the proposed method and the splinebased approach in forecasting the mortality rates at age 80 in years 1990-2010 Figure 2. Forecasted mortality rates using the proposed method based on the simulation study. Upper panel: age-specific mortality rates by year of death (left); year-specific mortality rates by age of death (right). Lower panel: derivative (rate of change) of age-specific mortality rates by year of death (left); derivative (rate of change) of year-specific mortality rates by age of death (right). (Table 1). Moreover, we observe that in forecasting the mortality rates at ages 50-79 in the year of 2011, the three method have similar confidence intervals.
On the other hand, through the estimate of the mortality rates,M, computed by Equation (4), Figure 2a plots the forecasted mortality rates at ages 50-80 in the calendar year 2011 and the forecasted mortality rates in the period from 1990 to 2011 at the age of 80. To further explore the trends of the mortality rates, included in Figure 2b are its partial derivatives with respect to x (age) and with respect to y (year), namely,M 0 x and M 0 y , respectively. More precisely, ðâ iþ1, j Àâ i, j Þu i, nÀ1 ðx=ðs 12 À s 11 ÞÞu j, m ðy=ðs 22 À s 21 ÞÞM a ðx, yÞ ðâ i, jþ1 Àâ i, j Þu i, n ðx=ðs 12 À s 11 ÞÞu j, mÀ1 ðy=ðs 22 À s 21 ÞÞM a ðx, yÞ: Figure 2a shows that, in the calendar year 2011, the forecasted mortality rates increase with the age group from 50 to 80. Its derivative shown in Figure 2b indicates that the slopes of the mortality rates with respect to x (age) at ages 50 and 80 in year 2011 are 6.2Â10 À3 and 1.8Â10 À1 , respectively, and it increases from age 50 to age 80. On the other hand, Figure 2a shows that the forecasted mortality rates decrease during the period from 1990 to 2011 at the age of 80. Figure 2b shows that the slopes of the mortality rates with respect to y (year) at age 80 in years 1990 and 2011 are À3.0Â10 À2 and À1.9Â10 À2 , respectively, and it is negative and increases from year 1990 to year 2011.
In brief, the simulation results show that the three methods have their respective merits in forecasting the short-term mortality trends. In addition, the proposed method is more effective than the existing methods, the spline-based approach and the Lee-Carter model, in facilely providing information about the slope (rate of change) of mortality rates, in order to farther explore more valid information on the mortality rates that we are interested in.
To furthermore examine the performance of the proposed method, additional simulation studies are presented in Section 6 in Supplementary Material with forecasting the mortality trends in the period 1990-2014 at ages 80-83 and in years 2011-2014 at ages 50-83. These additional simulation studies similarly show that the proposed method, the spline-based approach and the Lee-Carter method have their respective advantages in predicting the short-term mortality trends.

Data analysis
We now analyze the Swedish mortality data provided from the Human Mortality Database to illustrate the applicability and practicality of our proposed method. As in the simulation study, by considering the mortality rates of Sweden during the period from 1990 to 2010 at the ages from 50 to 79, we intend to predict the mortality rates of Sweden in the period 1990-2011 at age 80 and to forecast the mortality rates of Sweden in the year 2011 at ages 50-80.
Using the notations of Section 4, the design points fðx k , y k Þjk ¼ 1, Á Á Á , 630g form a lattice in ½0, 1 Â ½0, 1: The numbers of deaths and population of Sweden at the age of 50 þ 30x k in the year of 1990 þ 21y k are defined by Dðx k , y k Þ and Eðx k , y k Þ, respectively. The empirical mortality rates are denoted by m k ¼ Dðx k , y k Þ=Eðx k , y k Þ for k ¼ 1, Á Á Á , 630: We use the same way and setting described in Section 4 of the simulation study to analyze the real dataset. Table 2 shows that, in predicting the Swedish mortality trend at age 80 during the period 1990-2010 and in forecasting the mortality trend at ages 50-79 in the year of 2011, our proposed method is more effective than the spline-based approach and the Lee-Carter model, when the measures of the L 1 -norm and the MSE are considered. Moreover, in predicting the mortality rate of Swedish at age 80 in the year of 2011, the spline-based approach is more powerful than our proposed method. On the other hand, in forecasting the mortality trend at ages 50-79 in the year of 2011, the Lee-Carter model has better performance than the proposed method and the spline-based approach, when the measure of the sup-norm is considered. Figure 3 shows the AIC values with the different combinations of the MLEs of the Bernstein parameters ðn, m, aÞ from Swedish mortality data. Figure 3 indicates that when the MLEs of the polynomial orders n and m are 4 and 1, respectively, the MLEs of the Bernstein parameters ðn, m, aÞ, ðn,m,âÞ, have the lowest AIC value.
From Table S4 and Figures S9 and S10, we observe that, in forecasting the Swedish mortality rates at age 80 in years 1990-2010, the confidence intervals (95%) of the proposed method are wider than that of the spline-based approach. Thus, the confidence intervals (95%) of the proposed method have higher proportions containing the true values of the Swedish mortality rates at age 80 in years 1990-2010. The proportion of the confidence intervals (95%) based on our proposed method containing the true values of the Swedish mortality rates at age 80 in years 1990-2010 is 0.52. The proportion of the confidence intervals (95%) based on the spline-based method (Lee-Carter model) containing the true values of Swedish mortality rates at age 80 in years 1990-2010 is 0.19 (0.43).
On the other hand, in forecasting the Swedish mortality rates at ages 50-79 in the year of 2011, the three method have similar confidence intervals. The proportion of the confidence intervals (95%) based on the Lee-Carter model containing the true values of Swedish mortality rates at ages 50-79 in the year of 2011 is 0.6, while the proportion of the confidence intervals (95%) based on the proposed method (the spline-based approach) containing the true values of Swedish mortality rates at ages 50-79 in the year of 2011 is 0.4 (0.23). Figure 4a plots the forecasted mortality rates using the proposed method based on the Swedish mortality data at ages 50-80 in the calendar year of 2011 and at the age of 80 during the period from l990 to 2011. Their derivatives are presented in Figure 4b. Figure 4a shows that, in the calendar year 2011, the forecasted mortality rates of Swedish have an increasing tend at ages 50-80. Their derivatives in Figure 4b indicate that the slopes of the Swedish mortality rates with respect to x (age) at ages 50 and 80 in year 2011 are 5.9Â10 À3 and 1.8Â10 À1 , respectively, and it increases from age 50 to age 80. On the other hand, Figure 4a shows that, at the age of 80, the forecasted mortality rates have a trend toward decreasing during the period between 1990 and 2011. Their derivatives in Figure 4b indicate that the slopes of the Swedish mortality rates with respect to y (year) at age 80 in years 1990 and 2011 are À2.8Â10 À2 and À1.9Â10 À2 , respectively, and it is negative and increases from year 1990 to year 2011.
In addition, we illustrate our proposed method with forecasting the Swedish mortality trends in the period 1990-2014 at ages 80-83 and in years 2011-2014 at ages 50-83. year-specific mortality rates by age of death (right). Lower panel: derivative (rate of change) of age-specific mortality rates by year of death (left); derivative (rate of change) of year-specific mortality rates by age of death (right).
These forecasting results are similar to that based on the simulation data and hence they aren't presented here.
In summary, the results of the Swedish mortality data analysis show that the three competing methods have their respective advantages in predicting the short-term mortality trends. Moreover, the proposed method is more effective than the existing methods, the spline-based approach and the Lee-Carter model, in easily providing information about the slope (rate of change) of mortality rates, in order to further provide more valid information on the mortality trends that we are interested in.

Discussion
In this paper, a flexible maximum likelihood approach for forecasting the short-term mortality rates is introduced by using the shape-restricted regression model with the two-dimensional Bernstein polynomials. We use the property of the Bernstein coefficients to construct a regression model with a shape restriction for estimating the pattern of the mortality rates. An easy-to-implement algorithm is introduced to specify the Bernstein coefficients satisfying such a shape restriction. A simulated annealing method based on a MCMC sampling algorithm is introduced and connected with the AIC to effectively perform maximum likelihood based inference. Some general features of the asymptotic behavior of the proposed model are also discussed. Our numerical performances and real data analysis results show that the proposed method, the spline-based approach and the Lee-Carter model have their respective advantages in predicting the short-term mortality trends.
One of the important characters of the proposed method is that the shape-restricted regression model can be easily conducted by data information. Another important character is that the derivatives of Bernstein polynomials have simpler form than the polynomials in a different form, e.g., the B-splines. Their derivatives can be further explore the pattern of the mortality trends.
However, the proposed method has its limitations. In our results, we observe that the confidence interval estimates (95%) for forecasting the short-term mortality trends based on our proposed method and the other two competing methods, the spline-based approach and the Lee-Carter model, don't always include the true values of the mortality rates. In the real data analysis, we observe that the variance of the observed samples, the Swedish mortality rates at ages 50-79 in years 1990-2010, is 0.00019, while the variance of the samples of the forecasted mortality rates in the period 1990-2011 at the age of 80 and in the year of 2011 at ages 50-79 is 0.00060. We use these observed samples with smaller variance to predict the samples of the unknown mortality rates with larger variance. This may easily cause biased estimates or inappropriate confidence intervals for predicting the samples of the unknown mortality rates that have larger variance. Thus, it is necessary to propose a more effective algorithm for ameliorating the poor coverages of the confidence intervals from our proposed method in the future research. Moreover, our proposed method focuses mainly on forecasting the short-term mortality. An extension of the proposed method in predicting the medium-and long-term mortality rates will have significant benefits in practical data analysis.