Model fitting for Malaysian mortality rate: Comparison of Heligman-Pollard and P-splines smoothing

Malaysia has been experiencing longevity risk since the last decade due to improvements of mortality rates. Longevity risk refers to the probability of a person living longer than expected. According to the Department of Statistics Malaysia, a baby born in the year 2018 is predicted to live an average life of 75 years. Since the minimum retirement age policy of 60 years had come into force in 2012, the 2018 baby would live approximately 15 more years after retirement. Therefore, this study aims to compare the Heligman-Pollard and P-splines smoothing for fitting the Malaysian mortality rate. This model fitting will give a clear picture of the mortality pattern in predicting the mortality rate accurately, especially for the baby boomer generation. The data obtained from the Department of Statistics Malaysia are split into groups of five years, from 0 to 75 years old, and time ranges from 1995 to 2018. The data set from 1995 to 2010, known as the training set is used to fit the mortality rate. After fitting the mortality rate for both methods, this study will measure the performance in the testing set from 2011 until 2018. This study uses the mean absolute percentage error (MAPE) to identify the better method to fit the Malaysian mortality rate. Based on the MAPE values, P-splines smoothing gives a relatively smaller value compared to the Heligman-Pollard. For overall performance from 1995 to 2018, P-spline smoothing has proven to fit the Malaysian data well.


Introduction
Many countries have recorded improved mortality rates, both developed and developing nations. As the mortality rate has improved, people have experienced higher life expectancy, which is associated with longevity risk [1]. The risk of longevity refers to the probability of a person living longer than expected. Note that life expectancy is associated with a mortality rate that is significant to demographers, actuaries and the insurance industry. In order to address the longevity risk, many instruments have been developed, such as an index-based longevity hedge and longevity swap [2]. However, the instruments cannot wholly overcome the longevity risk but can be reduced by adopting natural hedge resulting in downsized risk [3]. Besides that, another parsimonious way to reduce the impact of longevity risk is by selecting the best model to fit the mortality rate before forecasting the mortality rates more accurately.
In Malaysia, the National Registration Department (JPN) is responsible for updating the number of populations for all age groups. The Department of Statistics Malaysia (DOSM) will then finalize the data from JPN before providing a statistic on the number of predicted deaths, as well as the number of predicted births for specific age groups. Using the statistical data of mortality rate for all age groups, insurance providers can come up with life insurance and retirement products, especially for people in certain age groups. In addition, the government, actuaries and insurance providers are also concerned with death rates as this will have an impact on estimating life insurance premiums and pension liabilities. From the perspective of insurance companies, policyholders live longer due to improved mortality rate [4]. If the expected death rate is underestimated or overestimated, life insurance premiums would be affected. Besides that, the mortality rate fitting will also benefit other regulatory bodies in Malaysia in determining and establishing their mortality rate decisions.
The improvement of mortality rate impacts pension liabilities, life insurance and retirement products significantly, especially in the senescence group. A higher life expectancy implies a good health status that contributes to an increment of the senescence group population. Due to the increase in healthcare status, the minimum compulsory retirement age in Malaysia has been set to 60 years from 2012 onwards [5]. By implementing 60 years old as the standard retirement age in Malaysia, the government should drive attention to this group in terms of social securities and pension liabilities.
According to the DOSM, a newborn baby in 2018 is expected to live an average of 75 years, that is, until 2093. On average, the baby boomer generation is estimated to live for around 15 years longer after their retirement. Therefore, this study will fit the Malaysian mortality rate with two models, the Heligman-Pollard (HP) and P-splines smoothing (PS). This study would then suggest a better model that fits the mortality rate according to the mean absolute percentage error (MAPE). Finally, this study will provide a clear picture of the patterns of mortality that benefit demographers, actuaries and the insurance industry.

Literature Review
To express the law of mortality that fits the post-war Australian mortality data, Heligman and Pollard developed a mathematical expression known as the Heligman-Pollard model, which consists of eight parameters [6]. The Heligman-Pollard model fits well with a range of populations, especially in capturing three age mortality patterns from three different age groups: childhood, young adulthood and senescence [7]. The first term consists of three parameters during the early years of life represented by an exponential function with a significant decrease to describe mortality. While the next term, which also includes three parameters, is presented with lognormal distribution expressing the young adult mortality as 'accident hump' (natural curve mortality). The final term is comprised of two parameters, distributed by the Gomperts exponential function, implying the degradation in the senescence group.
Since the Heligman-Pollard model's establishment, this technique has been used to fit and forecast mortality rate by many researchers [8][9][10][11]. Besides that, the Heligman-Pollard model also performs better than the Lee-Carter model for Dutch and Nigerian populations [12,13]. In Malaysia, the Heligman-Pollard model has been used to predict and model life tables in a series of studies [14][15][16][17]. It was seen that over the years, Malaysians are getting access to better healthcare, leading to the improvement of the mortality rate over time [18]. However, the Heligman-Pollard model seems to underestimate the probability of death at higher ages, where it is inappropriate once extrapolated beyond its range of applicability [19]. Besides that, the Heligman-Pollard model is often encountered with overparameterization issues [20,21]. Due to this overparameterization issues, this method is compared to another approach, namely P-splines smoothing.
P-splines smoothing was established in 1996 and used to smooth the mortality rate in projecting the mortality rate [22]. The authors mentioned integrated B-splines and penalties on the basis of adjacent coefficient finite order differences, then named as smoothing of P-splines. One of the benefits of this approach is that it can precisely fit polynomial data. Over the past few years, P-splines smoothing was 3 a popular method because this approach is effective in practical applications and simulations [23]. Since the development of P-spline smoothing, this technique has also been used in data fitting and forecasting in various studies [24][25][26].
Moreover, P-splines smoothing, as one of the smoothers, can be adapted to smoothing and forecasting two-dimensional mortality tables. Due to its flexibility, this smoothing has also been built in the R-package software, which is applied to model aggregate mortality data in both one and twodimensional contexts [27]. This smoothing proved to be efficient for forecasting broader range models [28]. Furthermore, P-spline model has been extended in terms of ANOVA-type and flexible Bayesian to provide great reliability by considering smooth space and time parameters [29,30]. P-splines smoothing method gives researchers the ability to smooth an unequally distributed data and selects the number of knots as well as the degree of polynomials to fit the actual data depending on the penalty order [31].
To the best of the author's knowledge, there is still a scarcity of research involving P-splines smoothing in fitting the mortality rate in Malaysia. Therefore, this study compares the Heligman-Pollard model and P-splines smoothing to fit the mortality rate in Malaysia. This comparison will determine which approach for both genders is better equipped to produce a more accurate result.

Methodology and Data Sources
There are two types of data required to measure the death rates, which are the number of exposures and the number of deaths for age groups. In determining the number of exposures, the average population at the beginning and end is estimated for a given year. The mortality rates are calculated from the year 1995 until 2018. The data is extracted from the DOSM where the population age is between 0 and 84 years old. As a whole, for both genders, this analysis fits the mortality rate for 17 age groups. For cross-validation purposes, the data is divided into two sets which are the training and testing set. The training set begins from 1995 to 2010, while the testing set commences from 2011 to 2018. Based on the training and testing set, the MAPE is computed to determine either the Heligman-Pollard or P-splines smoothing is a better model in fitting the Malaysian mortality rate. For each gender, the MAPE is calculated on a yearly basis for the training and testing set correspondingly. The analysis is carried out by using R software for fitting the mortality rate.

The Heligman-Pollard
The formulation of the Heligman-Pollard model is given by: where is the mortality rate for a specific age group of . Equation (1) is the Heligman-Pollard model formula that was established by [6] which consists of eight parameters and represented by A, B, C, D, E, F, G and H. The first term of equation (1) reflects the mortality rate from 0 to 9 years of age during early childhood. While the second term of equation (1) defines the rate of adult mortality from 10 to 40 years of age. The senescence mortality rate for 40 years of age and above is described in the last term of equation (1).
The available data is fitted using the Heligman-Pollard model for the training set to obtain eight parameters. Due to the complexity of the model, the built-in Levenberg-Marquardt algorithm and Poisson loss function are selected in the 'MortalityLaws' package to optimize the model. After that, using the parameters obtained from the training set, the testing set parameters are estimated. The testing set parameters are estimated using Autoregressive Moving Average (ARIMA), which can be expressed in the form of (p,d,q). The first term, p, is known as the lag order that indicates the number of lag observations included in the model. The second term, d, is called as the degree of differencing, representing the number of times the raw data are differenced. Finally, the last term, q, which is known as the moving average order, stipulates the size of the moving average. The ARIMA model is evaluated separately for each parameter. Then, the testing set's mortality rates are computed using the predicted parameters using equation (1). The computed value represents the fitted mortality rate and is used in computing the MAPE value.

P-splines Smoothing
P-splines smoothing is a combination of B-splines and penalties. B-splines are linked from polynomial points at specific values, known as knots [32,33]. Selection of the number of knots plays an important role since too many knots may lead to overfitting data and too few knots may cause the data to be underfitted. In order to choose the best fit curve, [34] proposed using the overfitting approach and addressed the overfitting data by introducing penalty with the integral of a squared higher derivative of the fitted curve as defined in equation (2): (2) Equation (2) compromises of two terms, which are the least square function and penalty terms. Instead of using the integral penalty, [22] proposed a simple difference penalty on the coefficient of adjacent B-splines to incorporate with overfitting data which can be expressed in equation (3): where; : fitted curve for i th data, : distance between knots, ( ) : derivative of the fitted function, and : smoothing parameter.
The iteratively reweighted least squares (IRWLS) are used in order to minimize the objective function in equation (3), which can be formulated as follows: The estimation of the coefficient ̃ in equation (4) is acquired from the IRWLS approach. Since the mortality rate requires the number of deaths that can be defined in a discrete variable, the errors produced are expressed in a Poisson distribution. Therefore, the diagonal matrix of weight is simplified into ̃= (̃) [35]. While ⊺̃ is obtained by following the standard procedure for fitting a generalized linear model with the B-splines regressors via the modification by [36]. The penalty is weighted by a positive regularization parameter, and described in equation (5). where; : d th order differences of the coefficient, and = ∆ .
Equation (5) depicts the relationship between the penalty, and the coefficient of order differences. To analyze the P-splines, the data is rearranged according to the particular age group. Each group is categorized by age groups of five years, i.e., 0-4, 5-9, …, 80-84. For each specific age group, there are 24 samples represent by 24 years of the study period. Next, the data is fitted by minimizing the Bayesian information criteria (BIC) for the training set. From the fitted data, the suggested smoothing parameter is generated. After that, the mortality rates are estimated by considering the smoothing parameter's recommended value. The estimated and the actual mortality rates of the testing set are plotted on the same axis to determine the forecasted value's deviance. The smoothing parameter is adjusted accordingly to achieve the best curve, which satisfies the goodness of fit and the smoothness of curve. Finally, the data is rearranged according to a yearly basis before computing the MAPE value.

Results and Discussion
MAPE is a practical approach to assess the average model performance quantitatively [37]. Therefore, this section presents the MAPE values for both genders' training and testing set. The MAPE values are illustrated in line graphs for both methods, the Heligman-Pollard (HP) and P-splines smoothing (PS). The MAPE values that show a lower value indicate that the model performed well from the other model. The MAPE values are presented in percentages and computed by using the formula as follows: where; : number of observations, : the actual value of mortality rates for i th data, and : the fitted value of mortality rates for i th data.

Training set
In this study, the MAPE computation is conducted for the male and female mortality rate separately. Then, the MAPE values are compared to identify each of the model's performance.   Figure 1 illustrates the MAPE values for the training set starting from the year 1995 to 2010 for the two methods, the Heligman-Pollard (HP) and P-splines smoothing (PS). Figure 1(a) reveals that Psplines smoothing outperformed the Heligman-Pollard model for the male mortality rate, with the exception of several years (1995, 1998, 1999 and 2000). In 1999, the maximum MAPE value was about 8% for the P-splines smoothing model. Whereas the minimum MAPE value was around 2% for P-splines smoothing in 2004. Figure 1(b) exhibits that P-splines performed better than the Heligman-Pollard model for the female mortality rate, except for in the years 1999, 2000, 2001, and 2002. The maximum MAPE value was about 8% for the Heligman-Pollard model in 1998. In 2005, the minimum MAPE value was about 1% for P-splines smoothing. It can be deduced that both models seem to fit well with the Malaysian mortality data because the MAPE for both genders' mortality rate varies from 1% to 8%. Regardless, the P-splines smoothing model performs better than the Heligman-Pollard model overall.

Testing set
This section provides the MAPE values for the Heligman-Pollard (HP) and P-splines smoothing (PS) for the male mortality rate for the testing set. For P-splines smoothing analysis, the data is organized according to the specific age group. Then, the smoothing curve is generated for the training set by minimizing the BIC as tabulated in table 1. In order to determine the best smoothing parameter, the fitted and the actual mortality rates for the testing set are plotted on the same axis. By examining the deviance between the actual mortality rates and the fitted curve of the testing set, the smoothing parameter is adjusted accordingly. In order to select the best smoothing parameter, the smoothing parameters produced must satisfy both of the goodness of fit and the smoothness of curve. After acquiring the best smoothing curve, the fitted mortality rates are recorded. Before computing the MAPE values, the data is rearranged on a yearly basis.    Year PS