Credible Regression Approaches to Forecast Mortality for Populations with Limited Data

Abstract

In this paper, we propose a credible regression approach with random coefficients to model and forecast the mortality dynamics of a given population with limited data. Age-specific mortality rates are modelled and extrapolation methods are utilized to estimate future mortality rates. The results on Greek mortality data indicate that credibility regression contributed to more accurate forecasts than those produced from the Lee–Carter and Cairns–Blake–Dowd models. An application on pricing insurance-related products is also provided.

1. Introduction

During the last decades, mortality has significantly declined in most developed countries around the world, mainly due to the continuous improvement of living conditions and the evolution of medical science and technology. This decline in mortality results in a steady increase of life expectancy, which creates higher financial responsibilities for governments and annuity providers. Consequently, finding ways to manage the mortality dynamics of a population is a very important step in building a sustainable health and pension system. In this spirit, actuaries and demographers have focused on the development of novel methods to model and forecast the mortality rates of a population.

Lee and Carter (1992) proposed a pioneer modelling method to forecast the mortality of the total population of the United States, by decomposing the mortality rates into two age and one period parameters. A remarkable variant of the Lee–Carter method, particularly designed for higher ages, was proposed by Cairns et al. (2006), who incorporated two period parameters. In the literature, we can find many extensions to these methods. Renshaw and Haberman (2006) extended the Lee–Carter model by including a cohort effect, while Cairns et al. (2009) added a cohort parameter to the Cairns et al. (2006) model. In addition, Plat (2009) proposed a model which combines preferable characteristics of the Lee and Carter (1992) and Cairns et al. (2006) models.

An issue that sometimes appears in mortality modelling is that, for some countries, there are too few data to fit. This issue affects the existing modelling methods, which inevitably base their forecasts on population datasets of a limited historical period of observations. In the literature, there are extensions of the Lee–Carter method that can be utilized when dealing with limited datasets. For instance, Li et al. (2004) extended the Lee–Carter model to be applied for Chinese and South Korean mortality data, which are available at only a few points in time and at unevenly spaced intervals. Zhao (2012) modified the Lee–Carter model by incorporating linearized cubic splines and other additive functions to approximate the model parameters and forecast mortality for short-base-period Chinese data and Huang and Browne (2017) presented a stochastic modification of the CMI (Continuous Mortality Investigation) model to project mortality improvement rates for limited Chinese data using clustering analysis techniques.

Recently, some alternative modelling approaches have also been proposed as a tool in mortality forecasting. Differently from the above Lee–Carter variants and extensions, these approaches are based on credibility theory, aiming to model the period patterns of limited mortality data for a specific age, using useful information from a wider age span. Bühlmann (1967) established the theoretical foundation of modern credibility theory (also known as greatest accuracy credibility theory), which is widely used in non-life insurance and Hachemeister (1975) introduced a credibility regression model to estimate auto-mobile bodily injury claims for various states in the USA.

Credibility regression has a long history in credibility literature, with applications mainly in non-life insurance: De Vylder (1978) proposed estimators for the structural parameters in a more general regression model; Norberg (1980) proposed empirical credibility estimators under various model assumptions and established asymptotic optimality; Ledolter et al. (1991) derived a credibility method that allows for time-varying parameters in the process; and Pitselis (2004) presented the relationship between claim amounts and a set of explanatory variables into a credibility regression model with cross-section and time effects, with applications for general insurance data. For an extensive review on credibility theory for non-life insurance, we also refer to the works of Goovaerts et al. (1990); Bühlmann and Gisler (2005) and Klugman et al. (2012).

Regarding some life insurance applications of credibility theory, Hardy and Panjer (1998) used empirical Bayes credibility theory to provide a theoretical basis for the calculation of risk measures associated with mortality risk for insurance companies. Salhi et al. (2016) proposed a credibility approach, which consists on reviewing the fitting parameters of a Makeham mortality curve, as new observations arrive. Schinzinger et al. (2016) presented a multivariate evolutionary credibility model for mortality improvement rates to describe the joint dynamics of mortality through time in several populations. Moreover, Li and Lu (2018) proposed a Bayesian non-parametric model for the mortality of a small population, when a benchmark mortality table of a larger population is also available and serves as part of the prior information. By using an adaptive smoothing procedure based on the local likelihood, Salhi and Thérond (2018) proposed a methodology to adjust the graduated mortality table based on credibility techniques. In addition, Gong et al. (2018) highlighted the importance of using credibility procedures in individual life and annuity business.

Two recent contributions to modelling mortality under a credibility framework were also made by Tsai and Lin (2017a, 2017b). In the first paper, they applied Bühlmann credibility to mortality data of Japan, the United Kingdom and the United States, while, in the second one, they incorporated Bühlmann credibility into the Lee and Carter (1992) model, the Cairns–Blake–Dowd model (Cairns et al. 2006) and the linear relational model of Tsai and Yang (2015) to improve forecasting performance for the United Kingdom dataset. However, it has been observed that the age-specific mortality rates show a clear downward trend over time. Moreover, when we have limited mortality data experience for a specific age, but extensive data experience for the entire age range, the use of credibility regression techniques should be preferred to capture mortality trends. Our work aims to exploit the advantages of credibility regression compared with the most widely used mortality models, as an alternative to Bühlmann credibility, to forecast the mortality rates, especially for populations with limited data.

The rest of this paper is organized as follows. Section 2 briefly reviews the Lee–Carter, the Cairns–Blake–Dowd and the random coefficients regression models. Section 3 proposes a credibility regression approach with randomly varying coefficients and a special case with fixed coefficients to model mortality rates. Section 4 presents the extrapolation methods used to estimate future mortality rates under the credibility regression approaches. An empirical illustration using Greek male and female data is presented in Section 5.1, in which forecasting performances of credibility regression, and the Lee–Carter and Cairns–Blake–Dowd methods are evaluated with the MAFE and RMSFE measures. A comparison between Bühlmann credibility and credibility regression forecasting methods is also presented in Section 5.2 and an application on pricing insurance-related products follows in Section 5.3. Finally, concluding remarks are discussed in Section 6.

2. Mortality Modelling: A Review of Methods

In this section, we briefly review the Lee–Carter model, the Cairns–Blake–Dowd model and the random coefficients regression models that will be utilized in next sections.

2.1. The Lee–Carter Model

In its original form, the Lee–Carter (LC) model links the natural logarithm of the observed mortality rates ${\displaystyle Y_{t,x}=logm(t,x)\mathrm{for}\mathrm{age}x=x_0,\cdots ,x_{k-1}\mathrm{and}\mathrm{year}t=t_0,\cdots ,t_{n-1}}$ with the following model predictor:

$$Y_{t,x}={\alpha }_x^{(1)}+{\alpha }_x^{(2)}{\kappa }_t+{ϵ}_{t,x},$$

(1)

where ${\displaystyle {\alpha }_x^{(1)}}$ is an age parameter that reflects the average mortality at age ${\displaystyle x}$, ${\displaystyle {\kappa }_t}$ is a period parameter which indicates the general level of mortality in year ${\displaystyle t}$ and ${\displaystyle {\alpha }_x^{(2)}}$ is an age parameter that indicates the deviation from the average mortality at age ${\displaystyle x}$, as the general level of mortality changes. The errors ${\displaystyle {ϵ}_{t,x}}$ are expected to be normally distributed, with zero mean and constant variance, reflecting specific period and age effects not captured by the model. Thus, after assuming that errors are independent and homoscedastic with zero mean, Lee and Carter (1992) suggested a close approximation to the SVD (Singular Value Decomposition) method, under the constraints ${\sum }_{x=x_0}^{x_{k-1}}{\alpha }_x^{(2)}=1$ and ${\sum }_{t=t_0}^{t_{n-1}}{\kappa }_t=0$, to obtain the following parameter estimates:

$${\widehat{\alpha }}_x^{(1)}=\frac{1}{t_{n-1}-t_0+1}\sum _{t=t_0}^{t_{n-1}}logm(t,x),$$

$${\widehat{\kappa }}_t=\sum _{x=x_0}^{x_{k-1}}[logm(t,x)-{\widehat{\alpha }}_x^{(1)}],$$

$${\widehat{\alpha }}_x^{(2)}=\frac{{\sum }_{t=t_0}^{t_{n-1}}[logm(t,x)-{\widehat{\alpha }}_x^{(1)}]{\widehat{\kappa }}_t}{{\sum }_{t=t_0}^{t_{n-1}}{\widehat{\kappa }}_t^2}.$$

Later on, to allow for heteroscedasticity in error variance, Brouhns et al. (2002) assumed that ${\displaystyle D(t,x)}$ follows a Poisson distribution with mean ${\displaystyle m(t,x)·E(t,x)}$. Under this approach, age and period parameters are estimated by maximising the log-likelihood function of (1).

After choosing one of the above estimation approaches, period estimates are extrapolated using time series methods. Lee and Carter (1992) suggested a random walk with a drift parameter ${\displaystyle \widehat{\theta }}$ to project period parameter for ${\displaystyle h=1,2,\cdots ,H}$ years ahead, according to ${\displaystyle {\widehat{\kappa }}_{t_{n-1}+h}={\widehat{\kappa }}_{t_{n-1}}+\widehat{\theta }h}$. Then, projected ${\displaystyle {\kappa }_t}$s are utilized along with the estimates of age parameters ${\displaystyle {\alpha }_x^{(1)}}$ and ${\displaystyle {\alpha }_x^{(2)}}$ to obtain the following mortality forecasts:

$${\widehat{Y}}_{t_{n-1}+h,x}={\widehat{\alpha }}_x^{(1)}+{\widehat{\alpha }}_x^{(2)}{\widehat{\kappa }}_{t_{n-1}+h}={\widehat{Y}}_{t_{n-1},x}+({\widehat{\alpha }}_x^{(2)}\widehat{\theta })h,\mathrm{for}h=1,2,\cdots ,H.$$

(2)

2.2. The Cairns–Blake–Dowd Model

The Cairns–Blake–Dowd (CBD) model links the logit transformation of one-year probabilities of death ${\displaystyle Y_{t,x}=logitq(t,x)}$ with the following model predictor:

$$Y_{t,x}=logitq(t,x)={\kappa }_t^{(1)}+(x-\overline{x}){\kappa }_t^{(2)}+{ϵ}_{t,x},$$

(3)

where ${\displaystyle {\kappa }_t^{(1)}}$ is a period parameter which indicates the general level of mortality in year ${\displaystyle t}$ and ${\displaystyle {\kappa }_t^{(2)}}$ is a period parameter that shows how mortality affects each age, while ${\displaystyle \overline{x}}$ is the mean age of the considered fitting age interval and ${\displaystyle {ϵ}_{t,x}}$ reflects specific effects not captured by the model and is expected to be normally distributed, with zero mean and constant variance. Again, we briefly present the estimates of the model parameters, which can be obtained by regressing ${\displaystyle logitq(t,x)}$ on ${\displaystyle (x-\overline{x})}$ for each ${\displaystyle t}$:

$${\widehat{\kappa }}_t^{(1)}=\frac{1}{x_{k-1}-x_0+1}\sum _{x=x_0}^{x_{k-1}}logitq(t,x)\mathrm{and}{\widehat{\kappa }}_t^{(2)}=\frac{{\sum }_{x=x_0}^{x_{k-1}}[logitq(t,x)(x-\overline{x})]}{{\sum }_{x=x_0}^{x_{k-1}}{(x-\overline{x})}^2}.$$

Alternatively, Cairns et al. (2009) assumed that deaths follow a Poisson distribution with mean ${\displaystyle m(t,x)·E(t,x)}$, where ${\displaystyle m(t,x)=-log[1-q(t,x)]}$. Then, the CBD model parameters are obtained by maximizing the log-likelihood function of (3). Assuming that period estimates are independent, each one of them is extrapolated using a random walk with a drift parameter (${\displaystyle {\widehat{\theta }}_i,i=1,2}$) and then mortality forecasts for ${\displaystyle h=1,2,\cdots ,H}$ are obtained by

$${\widehat{Y}}_{t_{n-1}+h,x}=({\widehat{\kappa }}_{t_{n-1}}^{(1)}+{\widehat{\theta }}_{1}h)+(x-\overline{x})({\widehat{\kappa }}_{t_{n-1}}^{(2)}+{\widehat{\theta }}_{2}h)={\widehat{Y}}_{t_{n-1},x}+[{\widehat{\theta }}_1+(x-\overline{x}){\widehat{\theta }}_2]h.$$

(4)

Remark 1.

We can easily observe that expressions in Equations (2) and (4) are both linear functions of the forecasting horizon ${\displaystyle h}$, where their intercepts are equal to the fitted rates of the last observed year and their slopes are the products of the estimated age parameters with the drift terms.

2.3. The Random Coefficients Regression Model

Empirical data indicate that mortality in each age ${\displaystyle x=x_0,\cdots ,x_{k-1}}$ decreases linearly over time. Especially in higher ages, mortality rates have been significantly improving over the last few years. We are interested in a model structure able to capture the improvement trends and describe the mortality evolution through time. For this reason, we consider a regression structure with random coefficients, aiming to capture the underlying mortality effects that are not included in the explanatory variables.

For each age ${\displaystyle x}$, the regression model with random coefficients is defined by $Y_{t,x}={\mathbf{\beta }}_{1t,x}+{\sum }_{k=2}^p{\mathbf{\beta }}_{kt,x}{Z}_{kt,x}+{ϵ}_{0t,x},\mathrm{for}t=t_0,t_1,\cdots ,t_{n-1}$, where ${\displaystyle Y_{t,x}}$ is the response variable, ${\displaystyle {\mathbf{\beta }}_{kt,x},k=1,2,\cdots ,p}$ are the randomly varying coefficients and ${\displaystyle Z_{kt,x}}$ are the explanatory variables. Then, each coefficient element can be decomposed in ${\displaystyle {\mathbf{\beta }}_{kt,x}={\mathbf{\beta }}_{k,x}+{ϵ}_{kt,x}}$, for all ${\displaystyle t}$ and ${\displaystyle k}$, with ${\displaystyle {\mathbf{\beta }}_k,x}$ and ${\displaystyle {ϵ}_{kt,x}}$ being the fixed and random parts, respectively, assuming that ${\displaystyle E({ϵ}_{kt,x})=0}$, ${\displaystyle \mathrm{Var}({ϵ}_{kt,x})={\sigma }_{k,x}^2}$ for all ${\displaystyle t}$ and ${\displaystyle \mathrm{Cov}({ϵ}_{kt,x},{ϵ}_{k^{{}^{\prime }}{t}^{{}^{\prime }},x})=0}$ for ${\displaystyle k\ne k^{{}^{\prime }}}$ and ${\displaystyle t\ne t^{{}^{\prime }}}$. For more details on regression models with random coefficients, we refer to the works of Hildreth and Houck (1968); Hsiao (1986) and Greene (2012).

The above formulation means that the unknown regression coefficients can take different values over an observed period. Actually, mortality dynamics for a specific age can vary over time, due to unknown or exogenous1 factors.

Nevertheless, the random coefficients regression model may be reduced to a fixed coefficients model with heteroscedastic variances, defined as

$$Y_{t,x}={\mathbf{\beta }}_{1,x}+{\displaystyle \sum _{k=2}^p}{\mathbf{\beta }}_{k,x}{Z}_{kt,x}+v_{t,x},\mathrm{with}{v}_{t,x}=({ϵ}_{0t,x}+{ϵ}_{1t,x})+{\sum }_{k=2}^p{Z}_{kt,x}{ϵ}_{kt,x},$$

(5)

where

$$E(v_{t,x})=0,\mathrm{Var}(v_{t,x})=E(v_{t,x}^2)=({\sigma }_{0,x}^2+{\sigma }_{1,x}^2)+{\sum }_{k=2}^p{\sigma }_{k,x}^2{Z}_{kt,x}^2\mathrm{and}\mathrm{Cov}(v_{t,x},v_{t^{{}^{\prime }},x})=0,$$

(6)

for all ${\displaystyle x}$ and ${\displaystyle t}$, with ${\displaystyle t\ne t^{{}^{\prime }}}$.

We have to point out that error variances ${\displaystyle {\sigma }_{0,x}^2}$ and ${\displaystyle {\sigma }_{1,x}^2}$ cannot be identified separately, while the sum (${\displaystyle {\sigma }_{0,x}^2+{\sigma }_{1,x}^2}$) can. Therefore, without loss of generality, ${\displaystyle {\sigma }_{0,x}^2}$ is dropped and the above variance is simplified to $\mathrm{Var}(v_{t,x})={\sigma }_{1,x}^2+{\sum }_{k=2}^p{\sigma }_{k,x}^2{Z}_{kt,x}^2$. Note that variance heteroscedasticity is still present even if ${\displaystyle {\sigma }_{k,x}^2={\sigma }_x^2}$ for ${\displaystyle k=1,2,\cdots ,p}$, due to the existence of squared explanatory variables ${\displaystyle Z_{kt,x}^2}$.

3. Credible Regression Mortality Models

In this section, we propose a mortality modelling approach embedded, for the first time, in a credibility regression framework with varying coefficients. The parameter estimation procedure is described and a special case with fixed coefficients is also provided.

3.1. A Credibility Regression Approach with Randomly Varying Coefficients

Denote ${\displaystyle D(t,x)}$ as the observed number of deaths at age ${\displaystyle x}$ in year ${\displaystyle t}$ and ${\displaystyle E(t,x)}$ as the average population aged ${\displaystyle x}$ during year ${\displaystyle t}$ (also called as population exposure to risk). Then, age-specific mortality rates ${\displaystyle m(t,x)}$ are obtained by the ratio ${\displaystyle D(t,x)/E(t,x)}$ and one-year probabilities of death can be derived from the identity ${\displaystyle q(t,x)=1-exp[-m(t,x)]}$, which is implied by the assumption of a constant force of mortality over each year of integer age and over each calendar year.

We assume that response variable ${\displaystyle Y_{t,x}}$ refers to an appropriate transform (${\displaystyle log}$ or ${\displaystyle logit}$) of a mortality measure [${\displaystyle m(t,x)}$ or ${\displaystyle q(t,x)}$] for age ${\displaystyle x=x_0,\cdots ,x_{k-1}}$ of year ${\displaystyle t=t_0,\cdots ,t_{n-1}}$, where variable ${\displaystyle x}$ corresponds to consecutive integer ages (${\displaystyle k}$ in total) and ${\displaystyle t}$ corresponds to consecutive calendar years (${\displaystyle n}$ in total). We also consider ${\displaystyle A_x}$ as an age-related random risk parameter, ${\displaystyle {\mathit{Y}}_x={(Y_{t_0,x},Y_{t_1,x},\cdots ,Y_{t_{n-1},x})}^{{}^{\prime }}}$ as a mortality vector and ${\displaystyle {\mathit{Z}}_x}$ as the design matrix of explanatory variables. We note that, in general, the design matrix could consist of various explanatory variables that reflect mortality characteristics. For instance, in a medical study, mortality may depend on various factors, such as the genetic background of an individual aged ${\displaystyle x}$, the life style, the nutrition, the toxicity of the environment, a possible infectious cause (bacteria, parasites, or fungi) or other socio-demographic factors that should affect mortality dynamics. Therefore, the pair that describes mortality evolution in age ${\displaystyle x}$ is ${\displaystyle (A_x,{\mathit{Y}}_x)}$, under the following assumptions:

(i)

The pairs ${\displaystyle (A_{x_0},{\mathit{Y}}_{x_0})}$, ${\displaystyle (A_{x_1},{\mathit{Y}}_{x_1})}$, ${\displaystyle \cdots }$, ${\displaystyle (A_{x_{k-1}},{\mathit{Y}}_{x_{k-1}})}$ are independent and ${\displaystyle A_{x_0},\cdots ,A_{x_{k-1}}}$ are independent and identically distributed.

(ii)

${\displaystyle E({\mathit{Y}}_x|A_x)={\mathit{Z}}_x\mathbf{\beta }(A_x)}$, where ${\displaystyle {\mathit{Z}}_x}$ is a fixed ${\displaystyle n\times p}$ design matrix of full rank ${\displaystyle p(<n)}$ and ${\displaystyle \mathbf{\beta }(A_x)}$ is an unknown regression vector of length ${\displaystyle p}$.

(iii)

${\displaystyle \mathrm{Cov}({\mathit{Y}}_x|A_x)=\mathrm{diag}\left[d_{t_0{t}_0}(A_x),\cdots ,d_{t_{n-1}{t}_{n-1}}(A_x)\right]}$,

${\displaystyle \mathrm{where}{d}_{tt}(A_x)={\sigma }_1^2(A_x)+\sum _{k=2}^p{\sigma }_k^2(A_x)Z_{kt,x}^2,\mathrm{with}{\sigma }_1^2(A_x)={\sigma }_{01}^2(A_x)+{\sigma }_{11}^2(A_x),}$

or in matrix formulation:

${\displaystyle \mathrm{Cov}({\mathit{Y}}_x|A_x)}$ = ${\displaystyle \left(\begin{array}{cc}{\sigma }_1^2(A_x)+{\displaystyle \sum _{k=2}^p}{\sigma }_k^2(A_x)Z_{k{t}_0,x}^2& 0\\ \hfill \ddots \\ 0& {\sigma }_1^2(A_x)+{\displaystyle \sum _{k=2}^p}{\sigma }_k^2(A_x)Z_{k{t}_{n-1},x}^2\end{array}\right)}$.

The structural parameters are defined as follows:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill \mathit{b}=E(\mathbf{\beta }(A_x)),\mathbf{\Phi }=\mathrm{Cov}\left[\mathbf{\beta }(A_x)\right],& {\mathit{s}}^2=E\left[{\sigma }^2(A_x)\right]=E\left[{\left({\sigma }_1^2(A_x),\cdots ,{\sigma }_p^2(A_x)\right)}^{{}^{\prime }}\right]\mathrm{and}\hfill \\ {\Delta }_x=E\left[\mathrm{Cov}({\mathit{Y}}_x|A_x)\right].\hfill \end{array}\end{array}$$

(7)

In such a regression setting, ${\displaystyle {\Delta }_x}$ has to be estimated. Consequently, instead of the ordinary least squares method, regression coefficients are estimated with the generalised least squares method (GLS). Then, an individual estimator of ${\displaystyle \mathbf{\beta }(A_x)}$ can be obtained by

$${\widehat{\mathbf{\beta }}}_x={({\mathit{Z}}_x^{{}^{\prime }}{\Delta }_x^{-1}{\mathit{Z}}_x)}^{-1}{\mathit{Z}}_x^{{}^{\prime }}{\Delta }_x^{-1}{\mathit{Y}}_x\mathrm{and}\mathrm{Cov}({\widehat{\mathbf{\beta }}}_x|A_x)={({\mathit{Z}}_x^{{}^{\prime }}{\Delta }_x^{-1}{\mathit{Z}}_x)}^{-1}.$$

(8)

Proposition 1.

Under the above assumptions, the credibility estimator of ${\displaystyle \mathbf{\beta }(A_x)}$ is given by

$${\mathit{B}}_x^{RC}={\mathit{C}}_x{\widehat{\mathbf{\beta }}}_x+(\mathit{I}-{\mathit{C}}_x)\mathit{b},$$

(9)

with

$${\mathit{C}}_x=\mathbf{\Phi }{({\mathbf{\Xi }}_x+\mathbf{\Phi })}^{-1},$$

(10)

where ${\displaystyle {\widehat{\mathbf{\beta }}}_x}$ is given in (8), ${\displaystyle \mathit{b}}$ and ${\displaystyle \mathbf{\Phi }}$ are defined in (7), ${\displaystyle {\mathbf{\Xi }}_x=E\left[\mathrm{Cov}({\widehat{\mathbf{\beta }}}_x|A_x)\right]}$ and ${\displaystyle \mathit{I}}$ is the ${\displaystyle p\times p}$ identity matrix.

Proof. 

The mean square error of (9) can be defined in terms of the norm ${\displaystyle {\parallel .\parallel }_E^2}$ as

$$\begin{array}{ccc}\hfill Q& =& \parallel \mathbf{\beta }(A_x)-{\mathit{B}}_x^{RC}{\parallel }_E^2\hfill \\ & =& E\left\{{[\mathbf{\beta }(A_x)-{\mathit{B}}_x^{RC}]}^{{}^{\prime }}[\mathbf{\beta }(A_x)-{\mathit{B}}_x^{RC}]\right\}\hfill \\ & =& E\left[{\left[{\mathbf{\beta }}^0(A_x)\right]}^{{}^{\prime }}{\mathbf{\beta }}^0(A_x)+{({\mathbf{\beta }}_x^0)}^{{}^{\prime }}{\mathit{C}}_x^{{}^{\prime }}{\mathit{C}}_x{\mathbf{\beta }}_x^0-{\left[{\mathbf{\beta }}^0(A_x)\right]}^{{}^{\prime }}{\mathit{C}}_x{\mathbf{\beta }}_x^0-{({\mathit{C}}_x{\mathbf{\beta }}_x^0)}^{{}^{\prime }}{\mathbf{\beta }}^0(A_x)\right],\hfill \end{array}$$

(11)

where ${\displaystyle {\mathbf{\beta }}^0(A_x)=\mathbf{\beta }(A_x)-\mathit{b}}$ and ${\displaystyle {\mathbf{\beta }}_x^0={\mathbf{\beta }}_x-\mathit{b}.}$ Using the product rule and differentiating with respect to matrix ${\displaystyle {\mathit{C}}_x}$, we have

$$\begin{array}{ccc}\hfill \frac{\partial Q}{\partial ({\mathit{C}}_x)}& =& -2E\left[{\mathbf{\beta }}^0(A_x){({\mathbf{\beta }}_x^0)}^{{}^{\prime }}-{\mathit{C}}_x{\mathbf{\beta }}_x^0{({\mathbf{\beta }}_x^0)}^{{}^{\prime }}\right].\hfill \end{array}$$

(12)

By substituting the values of ${\displaystyle {\mathbf{\beta }}^0(A_x)}$ and ${\displaystyle {\mathbf{\beta }}_x^0}$ and setting (12) equal to zero, we obtain

$$\begin{array}{ccc}\hfill {\mathit{C}}_x& =& E\left\{\left[\mathbf{\beta }(A_x)-\mathit{b}\right]{\left[{\mathbf{\beta }}_x-\mathit{b}\right]}^{{}^{\prime }}\right\}{\left\{E\left[({\mathbf{\beta }}_x-\mathit{b}){({\mathbf{\beta }}_x-\mathit{b})}^{{}^{\prime }}\right]\right\}}^{-1}\hfill \\ & =& \mathrm{Cov}[\mathbf{\beta }(A_x),{\mathbf{\beta }}_x]{\left[\mathrm{Cov}({\mathbf{\beta }}_x)\right]}^{-1}\hfill \\ & =& \left\{E\left\{\mathrm{Cov}[\mathbf{\beta }(A_x),{\mathbf{\beta }}_x|A_x]\right\}+\mathrm{Cov}\{E\left[\mathbf{\beta }(A_x)\right|A_x],E\left[{\mathbf{\beta }}_x\right|A_x]\}\right\}{\left[\mathrm{Cov}({\mathbf{\beta }}_x)\right]}^{-1}\hfill \\ & =& \left\{0+\mathrm{Cov}\left[\mathbf{\beta }(A_x)\right]\right\}\{E{[\mathrm{Cov}({\mathbf{\beta }}_x|A_x)+\mathrm{Cov}\left[E({\mathbf{\beta }}_x|A_x)\right]\}}^{-1},\hfill \end{array}$$

(13)

which yields (10). □

Then, the credibility estimator of future mortality rates ${\displaystyle Y_{t_{n-1}+h,x},h=1,2,\cdots ,H}$ may be compactly written as

$${\mathit{Y}}_x^{n+h}={\mathit{Z}}_x^{n+h}{\mathit{B}}_x^{RC},$$

where ${\displaystyle {\mathit{Z}}_x^{n+h}}$ denotes the design matrix of future periods.

3.2. Estimation of Structural Parameters

To estimate the structural parameters of the random coefficients credibility regression model, we can proceed similarly as in Hildreth and Houck (1968). Let ${\displaystyle {\mathit{r}}_x={(r_{t_0,x},\cdots ,r_{t_{n-1},x})}^{{}^{\prime }}}$ be the vector of the least squares residuals from the regression of ${\displaystyle {\mathit{Y}}_x}$ on ${\displaystyle {\mathit{Z}}_x}$ given ${\displaystyle A_x}$, which is obtained by

$$\begin{array}{c}\hfill {\mathit{r}}_x={\mathit{Y}}_x-{\mathit{Z}}_x{\widehat{\mathbf{\beta }}}_x={\mathit{M}}_x{\mathit{v}}_x,\end{array}$$

(14)

where ${\displaystyle {\widehat{\mathbf{\beta }}}_x={({\mathit{Z}}_x^{{}^{\prime }}{\mathit{Z}}_x)}^{-1}{\mathit{Z}}_x^{{}^{\prime }}{\mathit{Y}}_x}$ is the least squares estimator of coefficients in ordinary regression, ${\displaystyle {\mathit{M}}_x=\mathit{I}-{\mathit{Z}}_x{({\mathit{Z}}_x^{{}^{\prime }}{\mathit{Z}}_x)}^{-1}{\mathit{Z}}_x^{{}^{\prime }}}$ is a symmetric and idempotent matrix of order ${\displaystyle n\times n}$ and ${\displaystyle {\mathit{v}}_x={\mathit{Y}}_x-{\mathit{Z}}_x\mathbf{\beta }(A_x)}$ is the error term. Then, given ${\displaystyle A_x}$, the variance matrix of ${\displaystyle {\mathit{r}}_x}$, via (6), becomes

$$\begin{array}{c}\hfill E({\mathit{r}}_x{\mathit{r}}_x^{{}^{\prime }}|A_x)=E({\mathit{M}}_x{\mathit{v}}_x{\mathit{v}}_x^{{}^{\prime }}{\mathit{M}}_x|A_x),\end{array}$$

(15)

from which we can get

$$E({\dot{\mathit{r}}}_x|A_x)={\dot{\mathit{M}}}_x{\dot{\mathit{Z}}}_x{\sigma }^2(A_x),$$

(16)

where ${\displaystyle {\dot{\mathit{r}}}_x={(r_{t_0,x}^2,\cdots ,r_{t_{n-1},x}^2)}^{{}^{\prime }}}$, ${\displaystyle {\dot{\mathit{M}}}_x={\left\{m_{ts,x}^2\right\}}_{t,s=t_0,\cdots ,t_{n-1}}}$ and ${\displaystyle {\dot{\mathit{Z}}}_x={\left\{Z_{kt,x}^2\right\}}_{k=1,\cdots ,p,t=t_0,\cdots ,t_{n-1}}}$ are the Hadamard products of matrices ${\displaystyle {\mathit{r}}_x}$, ${\displaystyle {\mathit{M}}_x}$ and ${\displaystyle {\mathit{Z}}_x}$, respectively, while ${\displaystyle {\sigma }^2(A_x)}$ is as defined in (7). In addition, (16) implies that, for given ${\displaystyle A_x}$, least squares residuals ${\displaystyle {\dot{\mathit{r}}}_x}$ are regressed on ${\displaystyle {\sigma }^2(A_x)}$, which yields

$${\dot{\mathit{r}}}_x={\dot{\mathit{M}}}_x{\dot{\mathit{Z}}}_x{\sigma }^2(A_x)+{\mathit{e}}_x={\mathit{G}}_x{\sigma }^2(A_x)+{\mathit{e}}_x,$$

(17)

where ${\displaystyle {\mathit{G}}_x={\dot{\mathit{M}}}_x{\dot{\mathit{Z}}}_x}$ and ${\displaystyle {\mathit{e}}_x}$ is a ${\displaystyle n\times 1}$ disturbance vector, such that ${\displaystyle E({\mathit{e}}_x|A_x)=0}$. Hence, its variance-covariance matrix is given by

$$\begin{array}{ccc}\hfill \mathrm{Cov}({\mathit{e}}_x|A_x)& =& E\left\{[{\dot{\mathit{r}}}_x-E({\dot{\mathit{r}}}_x|A_x)]{[{\dot{\mathit{r}}}_x-E({\dot{\mathit{r}}}_x|A_x)]}^{{}^{\prime }}\right|A_x\}\hfill \\ & =& E({\dot{\mathit{r}}}_x|A_x){\left[E({\dot{\mathit{r}}}_x|A_x)\right]}^{{}^{\prime }}+2E({\mathit{r}}_x{\mathit{r}}_x^{{}^{\prime }}|A_x)\ast E({\mathit{r}}_x{\mathit{r}}_x^{{}^{\prime }}|A_x)-E({\dot{\mathit{r}}}_x|A_x){\left[E({\dot{\mathit{r}}}_x|A_x)\right]}^{{}^{\prime }}\hfill \\ & =& 2{\dot{\mathbf{\Psi }}}_x,\hfill \end{array}$$

(18)

where ${\displaystyle {\dot{\mathbf{\Psi }}}_x}$ represents the Hadamard product of matrix ${\displaystyle {\mathbf{\Psi }}_x}$ by itself, with

$${\mathbf{\Psi }}_x=E({\mathit{r}}_x{\mathit{r}}_x^{{}^{\prime }}|A_x)=E({\mathit{M}}_x{\mathit{v}}_x{({\mathit{M}}_x{\mathit{v}}_x)}^{{}^{\prime }}|A_x)={\mathit{M}}_{x}E({\mathit{v}}_x{\mathit{v}}_x^{{}^{\prime }}|A_x){\mathit{M}}_x={\mathit{M}}_x{\Delta }_x{\mathit{M}}_x.$$

Then, if ${\displaystyle {\sigma }_k^2}$s are known, the GLS estimator of ${\displaystyle {\sigma }^2(A_x)}$ in (17) is obtained by minimising the criterion function ${\displaystyle {[{\dot{\mathit{r}}}_x-{\mathit{G}}_x{\sigma }^2(A_x)]}^{{}^{\prime }}{(2{\dot{\mathbf{\Psi }}}_x)}^{-1}[{\dot{\mathit{r}}}_x-{\mathit{G}}_x{\sigma }^2(A_x)]}$, which gives

$${\widehat{\sigma }}_x^2={({\mathit{G}}_x^{{}^{\prime }}{\dot{\mathbf{\Psi }}}_x^{-1}{\mathit{G}}_x)}^{-1}{\mathit{G}}_x^{{}^{\prime }}{\dot{\mathbf{\Psi }}}_x^{-1}{\dot{\mathit{r}}}_x.$$

(19)

However, estimators of ${\displaystyle \mathbf{\beta }(A_x)}$ in (8) and ${\displaystyle {\sigma }^2(A_x)}$ in (19) are non-operational, because the variance-covariance matrices ${\displaystyle {\Delta }_x}$ and ${\displaystyle 2{\dot{\mathbf{\Psi }}}_x}$ are functions of unknown variances. Therefore, operational estimators of ${\displaystyle \mathbf{\beta }(A_x)}$ and ${\displaystyle {\sigma }^2(A_x)}$ can be obtained by replacing unknown matrices with estimators ${\displaystyle {\widehat{\Delta }}_x}$ and ${\displaystyle {\widehat{2\dot{\mathbf{\Psi }}}}_x}$, respectively. A least squares estimator of the unknown variances ${\displaystyle {\sigma }^2(A_x)}$ is directly obtained from (17) as follows:

$$\begin{array}{ccc}\hfill {\widehat{\sigma }}_x^{\mathbf{2}}& =& {({\mathit{G}}_x^{{}^{\prime }}{\mathit{G}}_x)}^{-1}{\mathit{G}}_x^{{}^{\prime }}{\dot{\mathit{r}}}_x\hfill \\ & =& {\left[{({\dot{\mathit{M}}}_x{\dot{\mathit{Z}}}_x)}^{{}^{\prime }}({\dot{\mathit{M}}}_x{\dot{\mathit{Z}}}_x)\right]}^{-1}{({\dot{\mathit{M}}}_x{\dot{\mathit{Z}}}_x)}^{{}^{\prime }}{\dot{\mathit{r}}}_x\hfill \\ & =& {({\dot{\mathit{Z}}}_x^{{}^{\prime }}{\dot{\mathit{M}}}_x^2{\dot{\mathit{Z}}}_x)}^{-1}{\dot{\mathit{Z}}}_x^{{}^{\prime }}{\dot{\mathit{M}}}_x{\dot{\mathit{r}}}_x,\hfill \end{array}$$

(20)

where equality ${\displaystyle {\dot{\mathit{M}}}_x^{{}^{\prime }}={\dot{\mathit{M}}}_x}$ holds true, since ${\displaystyle {({\mathit{M}}_x\ast {\mathit{M}}_x)}^{{}^{\prime }}={\mathit{M}}_x\ast {\mathit{M}}_x}$ for a symmetric matrix ${\displaystyle {\mathit{M}}_x}$.

Remark 2.

In the actuarial literature, there are many other types of estimators for variance in (17). For instance, Hildreth and Houck (1968) suggested the unbiased estimator ${\displaystyle {\widehat{\sigma }}_x^{\mathbf{2}(\mathrm{alt}1)}={({\dot{\mathit{Z}}}_x^{{}^{\prime }}{\dot{\mathit{M}}}_x{\dot{\mathit{Z}}}_x)}^{-1}{\dot{\mathit{Z}}}_x^{{}^{\prime }}{\dot{\mathit{r}}}_x}$, while Rao (1973) proposed the so-called “Minimum Norm Quadratic Unbiased Estimator” (MINQUE), given by ${\displaystyle {\widehat{\sigma }}_x^{\mathbf{2}(\mathrm{alt}2)}={({\dot{\mathit{Z}}}_x^{{}^{\prime }}{\dot{\mathit{M}}}_x{\dot{\mathit{Z}}}_x)}^{-1}{\dot{\mathit{Z}}}_x^{{}^{\prime }}{\dot{\mathit{M}}}_x{\dot{\mathit{r}}}_x}$.

The random coefficients (RC) credibility estimator of ${\displaystyle \mathbf{\beta }(A_x)}$, denoted as ${\displaystyle {\widehat{\mathit{B}}}_x^{RC}={({\widehat{B}}_{1x}^{RC},\cdots ,{\widehat{B}}_{px}^{RC})}^{{}^{\prime }}}$, is given by

$${\widehat{\mathit{B}}}_x^{RC}={\widehat{\mathit{C}}}_x{\widehat{\widehat{\mathbf{\beta }}}}_x+(\mathit{I}-{\widehat{\mathit{C}}}_x)\widehat{\mathit{b}},$$

(21)

where

$${\widehat{\widehat{\mathbf{\beta }}}}_x={({\mathit{Z}}_x^{{}^{\prime }}{\widehat{\Delta }}_x^{-1}{\mathit{Z}}_x)}^{-1}{\mathit{Z}}_x^{{}^{\prime }}{\widehat{\Delta }}_x^{-1}{\mathit{Y}}_x$$

and ${\displaystyle {\widehat{\Delta }}_x=diag\left({\widehat{\delta }}_{t_0{t}_0}^x,...,{\widehat{\delta }}_{t_{n-1},t_{n-1}}^x\right),}$ with ${\displaystyle {\widehat{\delta }}_{tt}^x={\widehat{s}}_1^2+\sum _{k=2}^p{\widehat{s}}_k^2{Z}_{kt,x}^2,t=t_0,\dots ,t_{n-1},}$ obtained according to (7), by using the mean of the estimated variances in (20). Future mortality estimates follow from

$${\widehat{\mathit{Y}}}_x^{n+h}={\mathit{Z}}_x^{n+h}{\widehat{\mathit{B}}}_x^{RC}={\mathit{Z}}_x^{n+h}{\widehat{\mathit{C}}}_x{\widehat{\widehat{\mathbf{\beta }}}}_x+{\mathit{Z}}_x^{n+h}(\mathit{I}-{\widehat{\mathit{C}}}_x)\widehat{\mathit{b}},h=1,2,\cdots ,H,$$

(22)

where ${\displaystyle {\widehat{\mathit{C}}}_x=\widehat{\mathbf{\Phi }}{({\widehat{\mathbf{\Xi }}}_x+\widehat{\mathbf{\Phi }})}^{-1},x=x_0,\cdots ,x_{k-1},}$ is the corresponding credibility factor. We suggest the following estimators for parameters ${\displaystyle \mathit{b}}$, ${\displaystyle {\mathbf{\Xi }}_x}$ and ${\displaystyle \mathbf{\Phi }}$ to obtain De Vylder’s (1978) optimality (minimum variance within the class of unbiased estimators):

$$\begin{array}{c}\hfill \widehat{\mathit{b}}={(\sum _{x=x_0}^{x_{k-1}}{\widehat{\mathit{C}}}_x)}^{-1}\sum _{x=x_0}^{x_{k-1}}{\widehat{\mathit{C}}}_x{\widehat{\widehat{\mathbf{\beta }}}}_x,\end{array}$$

(23)

$$\begin{array}{c}\hfill {\widehat{\mathbf{\Xi }}}_x=\frac{1}{x_{k-1}-x_0+1}\sum _{x^{\prime }=x_0}^{x_{k-1}}{({\mathit{Z}}_x^{{}^{\prime }}{\widehat{\Delta }}_{x^{\prime }}^{-1}{\mathit{Z}}_x)}^{-1}\end{array}$$

(24)

and

$$\begin{array}{c}\hfill \widehat{\mathbf{\Phi }}=\frac{1}{x_{k-1}-x_0}\sum _{x=x_0}^{x_{k-1}}{\widehat{\mathit{C}}}_x({\widehat{\widehat{\mathbf{\beta }}}}_x-\widehat{\mathit{b}}){({\widehat{\widehat{\mathbf{\beta }}}}_x-\widehat{\mathit{b}})}^{{}^{\prime }}.\end{array}$$

(25)

Note that the estimators of ${\displaystyle \widehat{\mathbf{\Phi }}}$ and ${\displaystyle \widehat{\mathit{b}}}$ are implicit functions of the parameter to be estimated and should be calculated iteratively, by imposing ${\displaystyle (\widehat{\mathbf{\Phi }}+{\widehat{\mathbf{\Phi }}}^{{}^{\prime }})/2=0}$ to retain symmetry after each iteration.

3.3. Credibility Regression with Fixed Coefficients and Weights: A Special Case

In the case of fixed regression’s coefficients, the previous model reduces to a special case of Hachemeister’s (1975) model with no weighs, i.e., ${\displaystyle {\mathit{W}}_x=\mathit{I}}$. In particular, some weights may appear in each regression line of ${\displaystyle A_x}$. For instance, population exposures ${\displaystyle E(t,x)}$, for ${\displaystyle t=t_0,\cdots ,t_{n-1}}$ can be used as weights. In this case, we have the standard regression case of Hachemeister’s model. To proceed, we follow the same Assumptions (i) and (ii) as in the random coefficients case, but covariance matrix in Assumption (iii) is simplified to ${\displaystyle \mathrm{Cov}({\mathit{Y}}_x|A_x)={\sigma }^2(A_x){\mathit{W}}_x}$, where ${\displaystyle {\mathit{W}}_x}$ is a fixed ${\displaystyle n\times n}$ positive definite diagonal matrix, with weights ${\displaystyle {\mathit{W}}_x=diag\left[E(t_0,x),\cdots ,E(t_{n-1},x)\right]}$. The structural parameters are now defined as

$$\mathit{b}=E\left[\mathbf{\beta }(A_x)\right],\mathit{U}=\mathrm{Cov}\left[\mathbf{\beta }(A_x)\right]\mathrm{and}{s}^2=E\left[{\sigma }^2(A_x)\right]$$

(26)

and the ordinary least squares estimator of the coefficients vector ${\displaystyle \mathbf{\beta }(A_x)}$ is given by

$${\widehat{\mathbf{\beta }}}_x={({\mathit{Z}}_x^{{}^{\prime }}{\mathit{W}}_x^{-1}{\mathit{Z}}_x)}^{-1}{\mathit{Z}}_x^{{}^{\prime }}{\mathit{W}}_x^{-1}{\mathit{Y}}_x$$

(27)

and the variance-covariance matrix is obtained by ${\displaystyle \mathrm{Cov}({\widehat{\mathbf{\beta }}}_x|A_x)={\sigma }^2(A_x){({\mathit{Z}}_x^{{}^{\prime }}{\mathit{W}}_x^{-1}{\mathit{Z}}_x)}^{-1}}$, while its expected value is given by ${\displaystyle E\left[\mathrm{Cov}({\widehat{\mathbf{\beta }}}_x|A_x)\right]=E\left[{\sigma }^2(A_x){({\mathit{Z}}_x^{{}^{\prime }}{\mathit{W}}_x^{-1}{\mathit{Z}}_x)}^{-1}\right]=s^2{({\mathit{Z}}_x^{{}^{\prime }}{\mathit{W}}_x^{-1}{\mathit{Z}}_x)}^{-1}.}$

Based on the above assumptions, the credibility estimator ${\displaystyle {\widehat{\mathit{B}}}_x^{FC}={({\widehat{B}}_{1x}^{FC},\cdots ,{\widehat{B}}_{px}^{FC})}^{{}^{\prime }}}$ of ${\displaystyle \mathbf{\beta }(A_x)}$ for the fixed coefficients (FC) model is given by

$${\widehat{\mathit{B}}}_x^{FC}={\widehat{\mathit{K}}}_x{\widehat{\mathbf{\beta }}}_x+(\mathit{I}-{\widehat{\mathit{K}}}_x)\widehat{\mathit{b}},$$

(28)

where ${\displaystyle {\widehat{\mathit{K}}}_x=\widehat{\mathit{U}}{[{\widehat{s}}^2{({\mathit{Z}}_x^{{}^{\prime }}{\mathit{W}}_x^{-1}{\mathit{Z}}_x)}^{-1}+\widehat{\mathit{U}}]}^{-1}}$ is the estimated credibility factor. Similarly, for the derivation of (28), we refer to Bühlmann and Gisler (2005). To recapture De Vylder’s (1978) optimality, we use the following estimators:

$${\widehat{s}}^2=\frac{1}{(x_{k-1}-x_0+1)(t_{n-1}-t_0+1-p)}\sum _{x=x_0}^{x_{k-1}}{({\mathit{Y}}_x-{\mathit{Z}}_x{\widehat{\mathbf{\beta }}}_x)}^{{}^{\prime }}{\mathit{W}}_x^{-1}({\mathit{Y}}_x-{\mathit{Z}}_x{\widehat{\mathbf{\beta }}}_x),$$

(29)

$$\widehat{\mathit{U}}=\frac{1}{x_{k-1}-x_0}\sum _{x=x_0}^{x_{k-1}}{\widehat{\mathit{K}}}_x({\widehat{\mathbf{\beta }}}_x-\widehat{\mathit{b}}){({\widehat{\mathbf{\beta }}}_x-\widehat{\mathit{b}})}^{{}^{\prime }},$$

(30)

$$\widehat{\mathit{b}}={(\sum _{x=x_0}^{x_{k-1}}{\widehat{\mathit{K}}}_x)}^{-1}\sum _{x=x_0}^{x_{k-1}}{\widehat{\mathit{K}}}_x{\widehat{\mathbf{\beta }}}_x.$$

(31)

Again, the estimators of ${\displaystyle \widehat{\mathit{U}}}$ and ${\displaystyle \widehat{\mathit{b}}}$ should be calculated iteratively, imposing ${\displaystyle (\widehat{\mathit{U}}+{\widehat{\mathit{U}}}^{{}^{\prime }})/2=0}$ after each iteration.

4. Extrapolation Methods for Estimating Future Mortality Rates

In this section, we fit the random coefficients (RC) and the fixed coefficients (FC) credibility regression models to mortality rates for age ${\displaystyle x=x_0,\cdots ,x_{k-1}}$ of year ${\displaystyle t=t_0,\cdots ,t_{n-1}}$. For both models, the fitted rates up to year ${\displaystyle t_{n-1}}$ can be compactly written as ${\displaystyle {\widehat{\mathit{Y}}}_x={\mathit{Z}}_x{\widehat{\mathbf{\beta }}}_x}$. As we noted before, design matrix ${\displaystyle {\mathit{Z}}_x}$ could consist of various independent variables that reflect risk factors for any given age ${\displaystyle x}$, but due to lack of specific data, we assume that ${\displaystyle {\mathit{Y}}_x}$s for each given age ${\displaystyle x}$, depend only on the period effects of each calendar year, i.e., ${\displaystyle {\mathit{Z}}_x=\mathit{Z}}$. However, if specific data are available, for instance in case of life insurance datasets, then more explanatory variables can be incorporated in the regression model. Henceforth, we consider the same design matrix ${\displaystyle \mathit{Z}={\displaystyle {\left(\begin{array}{cccc}1& 1& \cdots & 1\\ 1& 2& \cdots & n\end{array}\right)}^{{}^{\prime }}}}$ for all ${\displaystyle {\mathit{Y}}_x}$s.

4.1. Standard Extrapolation Method (SEM)

Based on current fitting data of the response variable ${\displaystyle \widehat{{\mathit{Y}}_x}={(Y_{t_0,x},Y_{t_1,x},\cdots ,Y_{t_{n-1},x})}^{{}^{\prime }}}$, mortality rates for one-year ahead are estimated by

$${\widehat{Y}}_{t_{n-1}+1,x}={\widehat{B}}_{1x}^c+{\widehat{B}}_{2x}^c(t_{n-1}-t_0+2),\mathrm{where}c=\mathrm{RC}\mathrm{or}\mathrm{FC}.$$

(32)

Similarly, estimates of future mortality rates for age ${\displaystyle x=x_0,\cdots ,x_{k-1}}$ are given by extrapolating one-year ahead estimates in (32) to ${\displaystyle {\widehat{Y}}_{t_{n-1}+h,x}={\widehat{B}}_{1x}^c+{\widehat{B}}_{2x}^c(t_{n-1}-t_0+1+h)}$, for ${\displaystyle h=2,3,\cdots ,H}$, where the credibility estimators ${\displaystyle {\widehat{\mathit{B}}}_x^c={({\widehat{B}}_{1x}^c,{\widehat{B}}_{2x}^c)}^{{}^{\prime }}}$ are obtained by (21) for the RC or (28) for the FC model. Hence, under this method, future estimates are based on the mortality data of the initial fitting span ${\displaystyle [t_0,t_{n-1}]}$.

4.2. Other Extrapolation Methods

In practice, two additional methods can also be used to extrapolate mortality rates over a given forecasting horizon ${\displaystyle h=1,2,\cdots ,H}$. Thus, for each one of the RC and FC models, one-year ahead estimates ${\displaystyle {\widehat{Y}}_{t_{n-1}+1,x}}$ can be embedded to the existing fitting span, with ${\displaystyle Y_{t_0,x}}$ simultaneously excluded from it, so that the fitting year span is moved forward by one year each time to ${\displaystyle [t_1,t_{n-1}+1]}$, ${\displaystyle [t_2,t_{n-1}+2]}$, ${\displaystyle [t_3,t_{n-1}+3]}$, ${\displaystyle \cdots }$. Then, after repeating the estimation procedure, we can consecutively obtain ${\displaystyle {\widehat{Y}}_{t_{n-1}+2,x}}$, ${\displaystyle {\widehat{Y}}_{t_{n-1}+3,x},{\widehat{Y}}_{t_{n-1}+4,x},\cdots ,{\widehat{Y}}_{t_{n-1}+H,x}}$. Under this “moving extrapolation method (MEM)”, future estimates are based on more recent mortality trends.

Alternatively, one-year ahead estimates ${\displaystyle {\widehat{Y}}_{t_{n-1}+1,x}}$ can be embedded to the existing fitting span, without removing ${\displaystyle Y_{t_0,x}}$, so that the fitting year span is extended by one year each time to ${\displaystyle [t_0,t_{n-1}+1]}$, ${\displaystyle [t_0,t_{n-1}+2]}$, ${\displaystyle [t_0,t_{n-1}+3]}$, ${\displaystyle \cdots }$. Hence, in each estimation step, credibility regression models are fitted on a continuously extended response variable, to obtain ${\displaystyle {\widehat{Y}}_{t_{n-1}+2,x}}$, ${\displaystyle {\widehat{Y}}_{t_{n-1}+3,x}}$, ${\displaystyle {\widehat{Y}}_{t_{n-1}+4,x},\cdots ,{\widehat{Y}}_{t_{n-1}+H,x}}$. Under this “extended extrapolation method (EEM)”, future mortality trends are based on both the initial mortality rates and the recent ones that have been obtained after each estimation step. Similar practical approaches have also been adopted by Luan (2015). The numerical results in the following section justify that all methods can be efficiently applied in actuarial practice.

Remark 3.

Similar extrapolation methods may be used in other regression or time series contexts, but here are customized to be used with the credibility regression models presented in Section 3.

5. Empirical Illustration

In this section, the Lee–Carter (LC), the Cairns–Blake–Dowd (CBD) and the credibility regression models are fitted on Greek mortality data. Then, forecasting results are evaluated using the mean absolute forecast error (MAFE) and the root mean of squared forecast error (RMSFE) measures. Greek data have a limited number of historical mortality observations (1981–2013), which are available on the Human Mortality Database (2017), structured by year, age and gender. Furthermore, in life insurance datasets similar limitations frequently exist. Credibility regression can efficiently capture the underlying data trends, especially in cases where there is limited mortality experience for a specific age, but extensive experience for the entire age range (the case of Greek data). Of course, credibility regression methods can also be used for larger datasets.

Mortality evolution for the period 1981–2010 in Greece is illustrated in Figure 1 and Figure 2 for ${\displaystyle logm(t,x)}$ and ${\displaystyle logitq(t,x)}$, respectively. Both mortality measures show a linearity for discrete ages ${\displaystyle x=40,60,80}$ of males (left panels of Figure 1 and Figure 2) and females (middle panels of Figure 1 and Figure 2). In addition, for both genders, average mortality decline shows a clear downward trend over time (right panels of Figure 1 and Figure 2).

5.1. Forecasting Results

For the numerical illustration that follows, we used the empirical age-specific mortality rates ${\displaystyle m(t,x)}$ from 1981 to 2010, for males and females at the ages of 15 to 84. This age span choice is in accordance with similar studies (Tsai and Lin 2017a, 2017b) as it corresponds to the age of a young adult up to the overall level of life expectancy in developed countries. To ensure robustness, relative to changes in the fitting range of data, we used two age and three period spans to extract forecasts for a 10-year ${\displaystyle (H=10)}$ forecasting horizon, presented in

Table 1. Selected fitting and forecasting periods.

Fitting Ages	Fitting Period	Forecasting Period
${\displaystyle [x_0,x_{k-1}]}$	${\displaystyle [t_0,t_{n-1}]}$	${\displaystyle [t_{n-1}+1,t_{n-1}+H]}$
${\displaystyle [15,84]}$	${\displaystyle [1981,2000]}$	${\displaystyle [2001,2010]}$
${\displaystyle [15,84]}$	${\displaystyle [1986,2000]}$	${\displaystyle [2001,2010]}$
${\displaystyle [15,84]}$	${\displaystyle [1991,2000]}$	${\displaystyle [2001,2010]}$
${\displaystyle [55,84]}$	${\displaystyle [1981,2000]}$	${\displaystyle [2001,2010]}$
${\displaystyle [55,84]}$	${\displaystyle [1986,2000]}$	${\displaystyle [2001,2010]}$
${\displaystyle [55,84]}$	${\displaystyle [1991,2000]}$	${\displaystyle [2001,2010]}$

. In particular, for the FC model, we used ${\displaystyle {\mathit{W}}_x=\mathit{I}}$ as weights. The credibility regression methods, as well as the LC and the CBD mortality models, were implemented in R (R Core Team 2017). In particular, for the Poisson LC and CBD fitting methods, we used the “LifeMetrics” R package2.

To retain linearity over each corresponding fitting period, the logarithmic transform ${\displaystyle Y_{t,x}=logm(t,x)}$ was used for the age-specific mortality rates and the logit transform ${\displaystyle Y_{t,x}=logitq(t,x)=log\frac{q(t,x)}{1-q(t,x)}}$ for the one-year probabilities of death. Forecast errors were then evaluated over the 10-year forecasting horizon using MAFE and RMSFE measures3, where smaller values indicate a better forecasting performance. Averaged (avg) MAFE and RMSFE values are obtained by using

$$MAF{E}_{avg}=\frac{1}{H\times (x_{k-1}-x_0+1)}\sum _{h=1}^H\sum _{x=x_0}^{x_{k-1}}\left|\widehat{m}(t_{n-1}+h,x)-m(t_{n-1}+h,x)\right|\times 100$$

(33)

and

$$RMSF{E}_{avg}=\sqrt{\frac{1}{H\times (x_{k-1}-x_0+1)}\sum _{h=1}^H\sum _{x=x_0}^{x_{k-1}}{\left[\widehat{m}(t_{n-1}+h,x)-m(t_{n-1}+h,x)\right]}^2}\times 100.$$

(34)

Similarly, in the case of using ${\displaystyle Y_{t,x}=logitq(t,x)}$ as response variable, ${\displaystyle m(t,x)}$ should be replaced by ${\displaystyle q(t,x)}$ in above formulas. Forecast accuracy results at percentage (%) scales are evaluated over the period ${\displaystyle [2001,2010]}$. MAFE and RMSFE values for fitting ages ${\displaystyle [15,84]}$, using ${\displaystyle Y_{t,x}=logm(t,x)}$ are illustrated in

Table 2. MAFE and RMSFE values of forecast errors over the period ${\displaystyle [2001,2010]}$ for ages ${\displaystyle [15,84]}$.

(a) MAFE Values
${\displaystyle {\mathit{MAFE}}_{[\mathbf{15},\mathbf{84}]}}$		Lee–Carter		Random Coefficients (RC)			Fixed Coefficients (FC)
Fitting Period	Gender	LC	LC-Poisson	${\displaystyle \mathit{SEM}}$	${\displaystyle \mathit{MEM}}$	${\displaystyle \mathit{EEM}}$	${\displaystyle \mathit{SEM}}$	${\displaystyle \mathit{MEM}}$	${\displaystyle \mathit{EEM}}$
${\displaystyle [1981,2000]}$	Male	0.1513	0.1569	0.1338	0.1205	0.1322	0.1352	0.1256	0.1361
	Female	0.0831	0.0861	0.0702	0.0740	0.0711	0.0691	0.0657	0.0690
${\displaystyle [1986,2000]}$	Male	0.1684	0.1514	0.1175	0.1196	0.1158	0.1203	0.1221	0.1206
	Female	0.0625	0.0799	0.0650	0.0696	0.0758	0.0608	0.0651	0.0613
${\displaystyle [1991,2000]}$	Male	0.1468	0.1681	0.1275	0.1257	0.1280	0.1288	0.1289	0.1289
	Female	0.0763	0.0959	0.0705	0.0678	0.0750	0.0622	0.0663	0.0669
Average		0.1147(7)	0.1231(8)	0.0974(5)	0.0962(3)	0.0997(6)	0.0961(2)	0.0956(1)	0.0971(4)
(b) RMSFE Values
${\displaystyle {\mathit{RMSFE}}_{[\mathbf{15},\mathbf{84}]}}$		Lee–Carter		Random Coefficients (RC)			Fixed Coefficients (FC)
Fitting Period	Gender	LC	LC-Poisson	${\displaystyle \mathit{SEM}}$	${\displaystyle \mathit{MEM}}$	${\displaystyle \mathit{EEM}}$	${\displaystyle \mathit{SEM}}$	${\displaystyle \mathit{MEM}}$	${\displaystyle \mathit{EEM}}$
${\displaystyle [1981,2000]}$	Male	0.3165	0.3220	0.2661	0.2349	0.2629	0.2716	0.2511	0.2745
	Female	0.1791	0.1825	0.1398	0.1594	0.1457	0.1376	0.1365	0.1374
${\displaystyle [1986,2000]}$	Male	0.3543	0.3200	0.2257	0.2265	0.2204	0.2362	0.2364	0.2375
	Female	0.1307	0.1742	0.1410	0.1509	0.1700	0.1264	0.1385	0.1288
${\displaystyle [1991,2000]}$	Male	0.3180	0.4010	0.2478	0.2457	0.2470	0.2570	0.2551	0.2516
	Female	0.1694	0.2415	0.1580	0.1511	0.1707	0.1302	0.1438	0.1476
Average		0.2447(7)	0.2735(8)	0.1964(4)	0.1948(3)	0.2028(6)	0.1932(1)	0.1936(2)	0.1962(5)

(a) and (b), respectively, while the corresponding values for ages ${\displaystyle [55,84]}$ with ${\displaystyle Y_{t,x}=logitq(t,x)}$ are presented in

Table 3. MAFE and RMSFE values of forecast errors over the period ${\displaystyle [2001,2010]}$ for ages ${\displaystyle [55,84]}$.

(a) MAFE Values
${\displaystyle {\mathit{MAFE}}_{[\mathbf{55},\mathbf{84}]}}$		Mortality Models				Random Coefficients (RC)			Fixed Coefficients (FC)
Fitting Period	Gender	LC	LC-Poisson	CBD	CBD-Poisson	${\displaystyle \mathit{SEM}}$	${\displaystyle \mathit{MEM}}$	${\displaystyle \mathit{EEM}}$	${\displaystyle \mathit{SEM}}$	${\displaystyle \mathit{MEM}}$	${\displaystyle \mathit{EEM}}$
${\displaystyle [1981,2000]}$	Male	0.3191	0.3322	0.2924	0.3247	0.2885	0.2642	0.2846	0.2871	0.2673	0.2870
	Female	0.1884	0.1933	0.1694	0.1884	0.1624	0.1458	0.1611	0.1629	0.1448	0.1627
${\displaystyle [1986,2000]}$	Male	0.2928	0.3186	0.2682	0.2988	0.2506	0.2547	0.2494	0.2544	0.2581	0.2541
	Female	0.1577	0.1769	0.1618	0.1708	0.1287	0.1377	0.1344	0.1289	0.1351	0.1288
${\displaystyle [1991,2000]}$	Male	0.3091	0.3622	0.2790	0.3348	0.2483	0.2461	0.2464	0.2538	0.2493	0.2525
	Female	0.1723	0.2126	0.1659	0.1868	0.1324	0.1350	0.1363	0.1363	0.1382	0.1361
Average		0.2399(8)	0.2660(10)	0.2228(7)	0.2507(9)	0.2018(3)	0.1973(1)	0.2020(4)	0.2039(6)	0.1988(2)	0.2035(5)
(b)RMSFEValues
${\displaystyle {\mathit{RMSFE}}_{[\mathbf{55},\mathbf{84}]}}$		Mortality Models				Random Coefficients (RC)			Fixed Coefficients (FC)
Fitting Period	Gender	LC	LC-Poisson	CBD	CBD-Poisson	${\displaystyle \mathit{SEM}}$	${\displaystyle \mathit{MEM}}$	${\displaystyle \mathit{EEM}}$	${\displaystyle \mathit{SEM}}$	${\displaystyle \mathit{MEM}}$	${\displaystyle \mathit{EEM}}$
${\displaystyle [1981,2000]}$	Male	0.4616	0.4848	0.3904	0.4467	0.4041	0.3644	0.3963	0.4065	0.3786	0.4061
	Female	0.2795	0.2842	0.2221	0.2512	0.2260	0.1996	0.2213	0.2304	0.2010	0.2299
${\displaystyle [1986,2000]}$	Male	0.4320	0.4872	0.3551	0.4073	0.3506	0.3522	0.3419	0.3631	0.3653	0.3618
	Female	0.2340	0.2699	0.2165	0.2244	0.1805	0.1940	0.1895	0.1803	0.1897	0.1800
${\displaystyle [1991,2000]}$	Male	0.4671	0.6129	0.3698	0.4625	0.3484	0.3423	0.3389	0.3660	0.3501	0.3616
	Female	0.2652	0.3721	0.2202	0.2510	0.1866	0.1888	0.1912	0.1961	0.1930	0.1954
Average		0.3566(9)	0.4185(10)	0.2957(7)	0.3405(8)	0.2827(4)	0.2736(1)	0.2799(3)	0.2904(6)	0.2796(2)	0.2891(5)

(a) and (b), respectively. Note that CBD model is included only for comparisons in fitting ages ${\displaystyle [55,84]}$, as it has been particularly designed for higher ages.

For both genders, accuracy results in Table 2 (a), (b) for fitting ages ${\displaystyle [15,84]}$ and Table 3 (a), (b) for ages ${\displaystyle [55,84]}$ indicate that, for each fitting period, credibility regression models outperform LC and CBD models for both error measures. Average values over the whole period are given in the last rows of each measure’s subtable. More precisely, for ages ${\displaystyle [15,84]}$, the FC-MEM and FC-SEM produce the smallest average MAFE and RMSFE, while for ages ${\displaystyle [55,84]}$, RC-MEM performs better in average under both measures, which indicates that forecasts for higher ages are based on more recent mortality trends. Moreover, we observe that errors are getting evidently larger, when shortening the age fitting span to ${\displaystyle [55,84]}$. This is due to the fact that both ${\displaystyle \left|\widehat{m}(t_{n-1}+h,x)-m(t_{n-1}+h,x)\right|}$ in (33) and ${\displaystyle {\left[\widehat{m}(t_{n-1}+h,x)-m(t_{n-1}+h,x)\right]}^2}$ in (34) are generally increasing with age ${\displaystyle x}$. Therefore, ${\displaystyle MAF{E}_{avg}}$ and ${\displaystyle RMSF{E}_{avg}}$ for ages ${\displaystyle [55,84]}$ are larger than those for ${\displaystyle [15,84]}$.

We note that, for our comparison, we used the Lee–Carter (1992) and Cairns–Blake–Dowd (2006) models, which incorporate only age and period effects. Models with cohort parameters were intentionally excluded from our analysis to be consistent with the age-period structure of the proposed credibility regression methods that model the period dynamics of mortality across the ages. For a modelling comparison study on Greek data that allows for models with cohort effects, we refer to the work of Bozikas and Pitselis (2018).

Credibility Effects on Mortality Modelling

In the preceding section, we used the proposed credibility regression methods to estimate the actual mortality trend for a specific age, by weighting the mortality trend for this age and the mean trend over a wider group of ages that encompasses much more information. Figure 3 illustrates the linear trend of the actual (observed) ${\displaystyle logitq(t,x)}$ for Greek males (left panel) and females (right panel), aged 55, 65 and 75 over the period 1981–2010. The intuition behind using credibility regression is that the proposed methods could potentially lead us to more accurate estimates for the intercept and the slope of the mortality curve for a given age ${\displaystyle x=x_0,\cdots ,x_{k-1}}$. To assure this, we used the absolute forecast errors by age (${\displaystyle {\mathrm{AFE}}_x}$) to compare the linear trend (intercept and the slope) of the ${\displaystyle logitq(2000+h,x)}$, ${\displaystyle h=1,\cdots ,10}$ between the actual rates and the rates produced from the best performing models for both genders over years ${\displaystyle [2001,2010]}$, with and without credibility, for pension ages ${\displaystyle [65,84]}$, fitted for ${\displaystyle [1981,2000]}$. For each model, ${\displaystyle {\mathrm{AFE}}_x}$ can be obtained by

$$AF{E}_x=\left|logit\widehat{q}(2000+h,x)-logitq(2000+h,x)\right|\times 100.$$

(35)

Figure 4 displays the ${\displaystyle AF{E}_x}$ comparison results, which indicate that, almost for all ages, credibility regression methods (dot lines) perform better than the LC (solid lines) and CBD (dashed lines) models. An alternative way to see how close the credibility forecasts are to the the actual mortality trend, Figure 5 illustrates the intercept and the slope of the actual rates and the forecasted ones for some ages, under the best performing methods (based on ${\displaystyle {\mathrm{AFE}}_x}$) with credibility (FC-MEM for males and RC-MEM for females) and without credibility (LC, CBD).

The trend lines for the RC-MEM and FC-MEM forecasts can be easily extracted using the ordinary least squares method. Recall that, the intercept and the slope for the LC and CBD models is given by Equations (2) and (4), respectively (Remark 1), while for the credibility method RC by (21) and for FC by (28). The illustrated results in Figure 5 indicate that intercepts and slopes of the FC-MEM (for males) and RC-MEM (for females) lines are closer to the actual ones, which set the best starting point for the forecasts.

5.2. Applying the Bühlmann Credibility Approach

Tsai and Lin (2017a) proposed a Bühlmann credibility approach to forecast mortality rates for both genders in Japan, the United Kingdom and the United States. This model can be directly obtained from the more general regression model, presented in Section 3.3, if we set ${\displaystyle {\mathit{Z}}_x={\displaystyle {\left(\begin{array}{cccc}1& 1& \cdots & 1\end{array}\right)}^{{}^{\prime }}}}$ and ${\displaystyle {\mathit{W}}_x=\mathit{I}}$ for ${\displaystyle x=x_0,\cdots ,x_{k-1}}$. Then, from (27), ${\displaystyle {\mathbf{\beta }}_x}$ is equal to ${\displaystyle {\overline{\mathit{Y}}}_x}$ and the model parameters, which are scalars now, can be estimated by

$${\widehat{s}}^2=\frac{1}{(x_{k-1}-x_0+1)(t_{n-1}-t_0)}\sum _{x=x_0}^{x_{k-1}}\sum _{t=t_0}^{t_{n-1}}\left(Y_{t,x}-{\overline{\mathit{Y}}}_x\right),$$

(36)

$$\widehat{b}=\frac{1}{x_{k-1}-x_0+1}\sum _{x=x_0}^{x_{k-1}}{\overline{\mathit{Y}}}_x=\frac{1}{(x_{k-1}-x_0+1)(t_{n-1}-t_0+1)}\sum _{x=x_0}^{x_{k-1}}\sum _{t=t_0}^{t_{n-1}}{Y}_{t,x}=\overline{Y},$$

(37)

$$\widehat{U}=\frac{1}{x_{k-1}-x_0}\sum _{x=x_0}^{x_{k-1}}\left({\overline{\mathit{Y}}}_x-\overline{Y}\right)-\frac{{\widehat{s}}^2}{t_{n-1}-t_0+1},$$

(38)

$$\widehat{K}=\left[(t_{n-1}-t_0+1)\widehat{U}\right]{\left[{\widehat{s}}^2+(t_{n-1}-t_0+1)\widehat{U}\right]}^{-1}.$$

(39)

The Bühlmann credibility estimates for one year ahead can be obtained by

$${\widehat{Y}}_{t_{n-1}+1,x}=\widehat{K}{\overline{\mathit{Y}}}_x+(1-\widehat{K})\overline{Y},\mathrm{for}x=x_0,\cdots ,x_{k-1}.$$

(40)

In contrast to the credibility regression approaches, which aim to capture the downward trend of ${\displaystyle m(t,x)}$s over ${\displaystyle t}$, for the Bühlmann credibility approach to be applied, this downward trend must be eliminated. For this reason, Tsai and Lin (2017a) applied the Bühlmann credibility model on the time series of mortality rate changes rather than the mortality rate levels, i.e., ${\displaystyle Y_{t,x}=logm(t,x)-logm(t-1,x)}$, for ${\displaystyle t_1,\cdots ,t_{n-1}}$. Then, they proposed two strategies for estimating ${\displaystyle Y_{t+h,x}}$, ${\displaystyle h=2,\cdots ,H}$. The first strategy expands fitting window (EW) by one year, similarly with the EEM regression method, described in Section 4 and the second one moves fitting window (MW) by one year, similarly with the MEM regression method. In what follows, we compare the forecasting performance between the Bühlmann and the credibility regression methods on Greek data. To be consistent with the Bühlmann modelling framework of Tsai and Lin (2017a), age fitting spans [21, 85] and [56, 85] were selected and forecast errors were also evaluated under the averaged MAPFE values, which is defined by ${\displaystyle MAPF{E}_{avg}=\frac{1}{H\times (x_{k-1}-x_0+1)}\sum _{h=1}^H\sum _{x=x_0}^{x_{k-1}}\frac{\left|\widehat{m}(t_{n-1}+h,x)-m(t_{n-1}+h,x)\right|}{\left|m(t_{n-1}+h,x)\right|}\times 100.}$

Error values for each gender were evaluated by fitting ${\displaystyle Y_{t,x}}$s on [1982, 2000], [1986, 2000], and [1990, 2000] period spans. Comparison of averaged MAFE, RMSFE and MAPFE4 results between Bühlmann and credibility regression methods is given for both genders in

Table 4. MAFE, RMSFE and MAPFE values of forecast errors over the period ${\displaystyle [2001,2010]}$ for ages ${\displaystyle [21,85]}$.

(a) MAFE Values
${\displaystyle {\mathit{MAFE}}_{[\mathbf{21},\mathbf{85}]}}$		Bühlmann Methods		Regression Methods – RC			Regression Methods – FC
Fitting Period	Gender	EW	MW	${\displaystyle \mathit{SEM}}$	${\displaystyle \mathit{MEM}}$	${\displaystyle \mathit{EEM}}$	${\displaystyle \mathit{SEM}}$	${\displaystyle \mathit{MEM}}$	${\displaystyle \mathit{EEM}}$
${\displaystyle [1982,2000]}$	Male	0.2348	0.2334	0.1404	0.1287	0.1379	0.1444	0.1361	0.1460
	Female	0.0930	0.0931	0.0816	0.0882	0.0843	0.0791	0.0780	0.0790
${\displaystyle [1986,2000]}$	Male	0.2170	0.2294	0.1329	0.1321	0.1306	0.1364	0.1358	0.1373
	Female	0.0918	0.0919	0.0782	0.0852	0.0909	0.0741	0.0805	0.0747
${\displaystyle [1990,2000]}$	Male	0.2355	0.2258	0.1399	0.1392	0.1369	0.1434	0.1422	0.1423
	Female	0.0954	0.0933	0.0836	0.0839	0.0879	0.0798	0.0818	0.0802
Average		0.1613(8)	0.1612(7)	0.1094(2)	0.1096(4)	0.1114(6)	0.1095(3)	0.1091(1)	0.1099(5)
(b)RMSFEValues
${\displaystyle {\mathit{RMSFE}}_{[\mathbf{21},\mathbf{85}]}}$		Bühlmann Methods		Regression Methods – RC			Regression Methods – FC
Fitting Period	Gender	EW	MW	${\displaystyle \mathit{SEM}}$	${\displaystyle \mathit{MEM}}$	${\displaystyle \mathit{EEM}}$	${\displaystyle \mathit{SEM}}$	${\displaystyle \mathit{MEM}}$	${\displaystyle \mathit{EEM}}$
${\displaystyle [1982,2000]}$	Male	0.4980	0.4948	0.2633	0.2342	0.2566	0.2756	0.2564	0.2799
	Female	0.1795	0.1795	0.1613	0.1884	0.1730	0.1540	0.1581	0.1541
${\displaystyle [1986,2000]}$	Male	0.4584	0.4861	0.2447	0.2387	0.2386	0.2564	0.2532	0.2591
	Female	0.1767	0.1772	0.1633	0.1781	0.1999	0.1484	0.1643	0.1502
${\displaystyle [1990,2000]}$	Male	0.4997	0.4765	0.2578	0.2574	0.2472	0.2704	0.2666	0.2668
	Female	0.1849	0.1802	0.1640	0.1761	0.1767	0.1567	0.1674	0.1570
Average		0.3329(8)	0.3324(7)	0.2091(1)	0.2122(5)	0.2153(6)	0.2103(2)	0.2110(3)	0.2112(4)
(c)MAPFEValues
${\displaystyle {\mathit{MAPFE}}_{[\mathbf{21},\mathbf{85}]}}$		Bühlmann Methods		Regression Methods – RC			Regression Methods – FC
Fitting Period	Gender	EW	MW	${\displaystyle \mathit{SEM}}$	${\displaystyle \mathit{MEM}}$	${\displaystyle \mathit{EEM}}$	${\displaystyle \mathit{SEM}}$	${\displaystyle \mathit{MEM}}$	${\displaystyle \mathit{EEM}}$
${\displaystyle [1982,2000]}$	Male	11.90	11.86	11.97	11.24	11.80	11.95	11.43	12.00
	Female	13.75	13.76	11.66	11.83	11.69	11.66	11.54	11.66
${\displaystyle [1986,2000]}$	Male	11.30	11.71	12.05	10.76	11.73	11.86	10.72	11.90
	Female	13.71	13.72	11.52	11.82	11.89	11.56	11.73	11.55
${\displaystyle [1990,2000]}$	Male	11.93	11.60	10.81	9.81	10.60	10.57	9.71	10.52
	Female	13.83	13.77	11.84	11.79	12.08	12.05	11.83	11.99
Average		12.73(7)	12.74(8)	11.64(6)	11.21(2)	11.63(5)	11.61(4)	11.16(1)	11.60(3)

(a)–(c) for ages [21, 85] and

Table 5. MAFE, RMSFE and MAPFE values of forecast errors over the period ${\displaystyle [2001,2010]}$ for ages ${\displaystyle [56,85]}$.

(a) MAFE Values
${\displaystyle {\mathit{MAFE}}_{[\mathbf{56},\mathbf{85}]}}$		Bühlmann Methods		Regression Methods – RC			Regression Methods – FC
Fitting Period	Gender	EW	MW	${\displaystyle \mathit{SEM}}$	${\displaystyle \mathit{MEM}}$	${\displaystyle \mathit{EEM}}$	${\displaystyle \mathit{SEM}}$	${\displaystyle \mathit{MEM}}$	${\displaystyle \mathit{EEM}}$
${\displaystyle [1982,2000]}$	Male	0.3599	0.3503	0.3272	0.3012	0.3210	0.3262	0.3036	0.3255
	Female	0.1686	0.1623	0.1717	0.1633	0.1735	0.1711	0.1595	0.1709
${\displaystyle [1986,2000]}$	Male	0.3233	0.3430	0.2893	0.2958	0.2886	0.2946	0.2991	0.2937
	Female	0.1481	0.1539	0.1534	0.1601	0.1617	0.1495	0.1573	0.1511
${\displaystyle [1990,2000]}$	Male	0.3745	0.3641	0.2958	0.2934	0.2937	0.2999	0.2954	0.2973
	Female	0.1670	0.1646	0.1617	0.1616	0.1613	0.1601	0.1625	0.1615
Average		0.2569(8)	0.2564(7)	0.2332(3)	0.2293(1)	0.2333(4)	0.2336(6)	0.2296(2)	0.2334(5)
(b) MAPFE Values
${\displaystyle {\mathit{RMSFE}}_{[\mathbf{56},\mathbf{85}]}}$		Bühlmann Methods		Regression Methods – RC			Regression Methods – FC
Fitting Period	Gender	EW	MW	${\displaystyle \mathit{SEM}}$	${\displaystyle \mathit{MEM}}$	${\displaystyle \mathit{EEM}}$	${\displaystyle \mathit{SEM}}$	${\displaystyle \mathit{MEM}}$	${\displaystyle \mathit{EEM}}$
${\displaystyle [1982,2000]}$	Male	0.5411	0.5261	0.4670	0.4213	0.4524	0.4700	0.4305	0.4679
	Female	0.2358	0.2282	0.2368	0.2242	0.2366	0.2389	0.2202	0.2381
${\displaystyle [1986,2000]}$	Male	0.4852	0.5159	0.4065	0.4138	0.3987	0.4221	0.4224	0.4185
	Female	0.2120	0.2178	0.2151	0.2235	0.2271	0.2089	0.2192	0.2107
${\displaystyle [1990,2000]}$	Male	0.5636	0.5472	0.4139	0.4130	0.4072	0.4291	0.4184	0.4195
	Female	0.2338	0.2307	0.2243	0.2246	0.2236	0.2217	0.2257	0.2232
Average		0.3786(8)	0.3777(7)	0.3273(4)	0.3201(1)	0.3243(3)	0.3318(6)	0.3227(2)	0.3297(5)
(c)MAPFEValues
${\displaystyle {\mathit{MAPFE}}_{[\mathbf{56},\mathbf{85}]}}$		Bühlmann Methods		Regression Methods – RC			Regression Methods – FC
Fitting Period	Gender	EW	MW	${\displaystyle \mathit{SEM}}$	${\displaystyle \mathit{MEM}}$	${\displaystyle \mathit{EEM}}$	${\displaystyle \mathit{SEM}}$	${\displaystyle \mathit{MEM}}$	${\displaystyle \mathit{EEM}}$
${\displaystyle [1982,2000]}$	Male	9.53	9.34	9.48	9.17	9.54	9.29	8.97	9.31
	Female	9.93	9.72	9.98	9.81	10.36	9.65	9.43	9.69
${\displaystyle [1986,2000]}$	Male	8.82	9.20	8.78	8.97	8.99	8.61	8.85	8.66
	Female	9.23	9.45	9.14	9.42	9.62	8.84	9.26	8.98
${\displaystyle [1990,2000]}$	Male	9.82	9.61	8.85	8.74	9.00	8.62	8.74	8.78
	Female	9.88	9.81	9.49	9.33	9.46	9.32	9.37	9.48
Average		9.54(8)	9.52(7)	9.29(5)	9.24(4)	9.50(6)	9.06(1)	9.10(2)	9.15(3)

(a)–(c) for ages [56, 85].

The results indicate that credibility regression methods produce the smallest MAFE, RMSFE and MAPFE values for the majority of the selected fitting periods for both age spans. More precisely, the FC-MEM method has the best average performance according to MAFE and MAPFE values for ages [21, 85], while the RC-MEM method seems to be more appropriate to capture future mortality trends for older ages [56, 85]. We note that the smallest values in average are produced by different regression methods, depending on which measure is used. Such inconsistencies are expected due to the nature of MAFE, RMSFE and MAPFE formulas. That was also pointed out by (Tsai and Yang 2015, p. 9).

5.3. Application in Insurance-Related Products

In this section, we apply the mortality forecasts obtained from the Lee–Carter, the Cairns–Blake–Dowd and the credibility regression models to calculate life premiums, reflecting the appropriateness of each model in pricing applications. Denote ${\displaystyle A_{{\stackrel{1}{t}}_{n-1}+1,x:\overline{K\hspace{.2em}|}}}$ as the fully discrete life insurance premium, payable at the end of the year of death, if it occurs within a term of ${\displaystyle K}$ years and ${\displaystyle A_{t_{n-1}+1,x:\stackrel{1}{\overline{K\hspace{.2em}|}}}}$ as the pure endowment, payable at the end of ${\displaystyle K}$ years in case of being alive. Both products are issued to an insured aged ${\displaystyle x}$ in year ${\displaystyle t_{n-1}+1}$. Net premiums (NP) are obtained (see Bozikas and Pitselis 2018) by

$$A_{{\stackrel{1}{t}}_{n-1}+1,x:\overline{K\hspace{.2em}|}}=\sum _{k=0}^{K-1}{}_k{p}_{t_{n-1}+1,x}.q(t_{n-1}+1+k,x+k).{(1+i)}^{-(k+1)},$$

(41)

$$A_{t_{n-1}+1,x:\stackrel{1}{\overline{K\hspace{.2em}|}}}={}_K{p}_{t_{n-1}+1,x}.{(1+i)}^{-K},$$

(42)

where ${\displaystyle {}_k{p}_{t_{n-1}+1,x}}$ denotes the k-year survival probability for age ${\displaystyle x}$ in year ${\displaystyle t_{n-1}+1}$, while its estimate is given by ${\displaystyle {}_k{\widehat{p}}_{t_{n-1}+1,x}={\widehat{p}}_{t_{n-1}+1,x}.\cdots .{\widehat{p}}_{t_{n-1}+1+k-1,x+k-1}}$, ${\displaystyle k=1,\cdots ,K-1}$ and similarly for ${\displaystyle {}_K{\widehat{p}}_{t_{n-1}+1,x}}$, where ${\displaystyle i}$ is the interest rate and ${\displaystyle {}_0{\widehat{p}}_{t_{n-1}+1,x}=1}$. In addition, to see the performance on a life annuity product, typically used for pension applications, denote ${\displaystyle {\ddot{a}}_{t_{n-1}+1,x:\overline{K\hspace{.2em}|}}}$ as the discrete life annuity-due at age ${\displaystyle x}$ in year ${\displaystyle t_{n-1}+1}$, payable annually for up to ${\displaystyle K}$ years. Its actuarial present value (APV) can be obtained by

$${\ddot{a}}_{t_{n-1}+1,x:\overline{K\hspace{.2em}|}}=\sum _{k=0}^{K-1}{}_k{p}_{t_{n-1}+1,x}.{(1+i)}^{-k}.$$

(43)

Then, we apply the estimated mortality rates obtained from the LC, CBD and credibility methods, fitted to 1981–2000 rates, to calculate the NPs and APVs for ages 55–74 with ${\displaystyle K=10}$, assuming ${\displaystyle i=4\%}$. The errors between forecasted values and those produced from the observed mortality rates for the years 2001–2010 are evaluated using MAFE and RMSFE, which are defined by

$$MAF{E}_{avg}=\frac{1}{20}\sum _{x=55}^{74}\left|{\widehat{A}}_{\stackrel{1}{2}001,x:\overline{10\hspace{.2em}|}}-A_{\stackrel{1}{2}001,x:\overline{10\hspace{.2em}|}}\right|\times 100,$$

(44)

$$RMSF{E}_{avg}=\sqrt{\frac{1}{20}{\displaystyle \sum _{x=55}^{74}}{({\widehat{A}}_{\stackrel{1}{2}001,x:\overline{10\hspace{.2em}|}}-A_{\stackrel{1}{2}001,x:\overline{10\hspace{.2em}|}})}^2}\times 100.$$

(45)

Similarly, MAFE and RMSFE formulas are adjusted for pure endowment or annuity products by replacing ${\displaystyle A_{{\stackrel{1}{t}}_{n-1}+1,x:\overline{K\hspace{.2em}|}}}$ with ${\displaystyle A_{t_{n-1}+1,x:\stackrel{1}{\overline{K\hspace{.2em}|}}}}$ or ${\displaystyle {\ddot{a}}_{t_{n-1}+1,x:\overline{K\hspace{.2em}|}}}$ in Equations (44) and (45).

Table 6. MAFE and RMSFE values (ranking order in brackets) for a 10-year forecasted life insurance, a pure endowment and a life annuity for males and females of ages 55–74 during 2001–2010.

(a) Life Insurance
${\displaystyle MAF{E}_{avg}}$	Mortality Models				Random Coefficients (RC)			Fixed Coefficients (FC)
Gender	LC	LC-Poisson	CBD	CBD-Poisson	${\displaystyle SEM}$	${\displaystyle MEM}$	${\displaystyle EEM}$	${\displaystyle SEM}$	${\displaystyle MEM}$	${\displaystyle EEM}$
Male	1.6019(8)	1.5640(7)	1.7151(10)	1.6794(9)	1.5000(6)	1.4169(2)	1.4924(5)	1.4735(3)	1.3932(1)	1.4741(4)
Female	1.0264(6)	1.0269(7)	1.2141(10)	1.1079(9)	1.0262(5)	0.9317(2)	1.0346(8)	0.9898(3)	0.8840(1)	0.9910(4)
${\displaystyle RMSF{E}_{avg}}$	Mortality Models				Random Coefficients (RC)			Fixed Coefficients (FC)
Gender	LC	LC-Poisson	CBD	CBD-Poisson	${\displaystyle SEM}$	${\displaystyle MEM}$	${\displaystyle EEM}$	${\displaystyle SEM}$	${\displaystyle MEM}$	${\displaystyle EEM}$
Male	1.8423(8)	1.8043(7)	1.9401(10)	1.9089(9)	1.7125(6)	1.6143(2)	1.7043(5)	1.6871(3)	1.5989(1)	1.6875(4)
Female	1.2320(8)	1.2294(7)	1.4023(10)	1.2918(9)	1.2133(5)	1.0965(2)	1.2215(6)	1.1744(3)	1.0494(1)	1.1756(4)
(b) Pure Endowment
${\displaystyle MAF{E}_{avg}}$	Mortality Models				Random Coefficients (RC)			Fixed Coefficients (FC)
Gender	LC	LC-Poisson	CBD	CBD-Poisson	${\displaystyle SEM}$	${\displaystyle MEM}$	${\displaystyle EEM}$	${\displaystyle SEM}$	${\displaystyle MEM}$	${\displaystyle EEM}$
Male	1.1439(8)	1.1139(7)	1.2417(10)	1.2044(9)	1.0722(6)	1.0153(2)	1.0681(5)	1.0512(3)	0.9942(1)	1.0518(4)
Female	0.7181(7)	0.7192(8)	0.8894(10)	0.7923(9)	0.7340(5)	0.6717(2)	0.7463(6)	0.7026(3)	0.6297(1)	0.7038(4)
${\displaystyle RMSF{E}_{avg}}$	Mortality Models				Random Coefficients (RC)			Fixed Coefficients (FC)
Gender	LC	LC-Poisson	CBD	CBD-Poisson	${\displaystyle SEM}$	${\displaystyle MEM}$	${\displaystyle EEM}$	${\displaystyle SEM}$	${\displaystyle MEM}$	${\displaystyle EEM}$
Male	1.3274(8)	1.2975(7)	1.4104(10)	1.3786(9)	1.2347(6)	1.1659(2)	1.2303(5)	1.2150(3)	1.1535(1)	1.2154(4)
Female	0.8745(8)	0.8717(7)	1.0310(10)	0.9319(9)	0.8774(5)	0.7968(2)	0.8889(6)	0.8440(3)	0.7552(1)	0.8451(4)
(c) Life Annuity
${\displaystyle MAF{E}_{avg}}$	Mortality Models				Random Coefficients (RC)			Fixed Coefficients (FC)
Gender	LC	LC-Poisson	CBD	CBD-Poisson	${\displaystyle SEM}$	${\displaystyle MEM}$	${\displaystyle EEM}$	${\displaystyle SEM}$	${\displaystyle MEM}$	${\displaystyle EEM}$
Male	5.4466(8)	5.2602(7)	6.3032(10)	5.7857(9)	5.1642(6)	4.9260(2)	5.1561(5)	5.0331(3)	4.7893(1)	5.0369(4)
Female	2.9140(7)	2.9361(8)	4.4471(10)	3.5932(9)	3.1527(5)	2.9479(2)	3.2151(6)	2.9656(3)	2.7024(1)	2.9721(4)
${\displaystyle RMSF{E}_{avg}}$	Mortality Models				Random Coefficients (RC)			Fixed Coefficients (FC)
Gender	LC	LC-Poisson	CBD	CBD-Poisson	${\displaystyle SEM}$	${\displaystyle MEM}$	${\displaystyle EEM}$	${\displaystyle SEM}$	${\displaystyle MEM}$	${\displaystyle EEM}$
Male	6.6138(7)	6.4342(8)	7.2583(10)	6.8729(9)	6.1786(6)	5.8730(2)	6.1608(5)	6.0681(3)	5.7919(1)	6.0704(4)
Female	3.7013(7)	3.7218(8)	5.1510(10)	4.3300(9)	3.9187(5)	3.6344(2)	3.9846(6)	3.7084(3)	3.3878(1)	3.7155(4)

presents the averaged error values in ranking order for a 10 year forecasted life insurance, pure endowment and life annuity for both genders, aged 55–74 in 2001–2010. In addition, Figure 6 illustrates the absolute forecast error values against each corresponding age (${\displaystyle {\mathrm{AFE}}_x}$) for the top LC, CBD, RC and FC credibility regression methods for males and females, according to Table 6 values. For each model, ${\displaystyle {\mathrm{AFE}}_x}$ is obtained from ${\displaystyle AF{E}_x=\left|{\widehat{A}}_{\stackrel{1}{x},2001:\overline{K\hspace{.2em}|}}-A_{\stackrel{1}{x},2001:\overline{K\hspace{.2em}|}}\right|\times 100}$.

According to MAFE and RMSFE values for both genders and insurance products in Table 6, credibility regression models produce better insurance-related forecasts in comparison with the LC and CBD modelling methods. We can easily observe that for each gender, error measures coincide in the same ranking order for all insurance products. In particular, measures show that credibility regression methods under a moving fitting span outperform LC and CBD methods in aggregate, with FC-MEM being dominant and RC-MEM following. This fact is also evident in Figure 6, where absolute error values against age for the MEM regression models lie on the lower levels for all the insurance products. Nevertheless, the FC-SEM should also be a good modelling choice for pricing insurance-related products.

6. Concluding Remarks

Credibility regression techniques seem to be of special interest and particularly useful for mortality datasets of a relatively short historical period of observations (limited data), as they can efficiently capture the underlying mortality trend for a given age, using all the information gained from populations of other ages. This paper proposes mortality modelling approaches embedded, for the first time, in a credibility regression framework. In our illustration on Greek data, credibility regression approaches resulted in better forecasts for both genders (in terms of MAFE and RMSFE measures), compared to the Lee–Carter and Cairns–Blake–Dowd models, as well as the Bühlmann credibility approach (Tsai and Lin 2017a). Finally, their performance was also evaluated on insurance-related applications.

Specifically, in Section 3, we proposed a credibility regression mortality framework with randomly varying coefficients and a special case with fixed coefficients. To estimate future mortality rates, we presented extrapolation methods for each credibility approach in Section 4. The applicability of our modelling approaches was comparatively illustrated on Greek male and female data in Section 5, accompanied with an explanation of the credibility effects in mortality modelling and a pricing application on insurance-related products. From our analysis, we concluded that, in aggregate, credibility modelling methods performed better than the LC and CBD methods. Forecasting accuracy results indicate that, for the whole age fitting span, fixed coefficients credibility methods performed better on average, while, for higher ages, the RC-MEM should also be a good choice. In addition, the FC-MEM performed a bit better in aggregate on pricing insurance-related products.

Furthermore, we noted that FC-SEM credibility forecasts were closer to observed rates for the same periods, when we used population exposure to risk as weights, i.e., ${\displaystyle {\mathit{W}}_x=diag\left[E(t_0,x),\cdots ,E(t_{n-1},x)\right],\mathrm{for}x=x_0,\cdots ,x_{k-1}}$, but weighted regression is restricted for use only under the SEM, as ${\displaystyle E(t,x)}$s are practically unknown for the upcoming years. Additionally, during the estimation procedure for the random regression models, we observed that, if we use the MINQUE estimator (Remark 2) instead of (20), error values for all the credibility modelling methods become even smaller for both genders.

For the sake of comparability, the Bühlmann credibility approach (Tsai and Lin 2017a) was applied on our dataset in Section 5, where the credibility regression methods resulted to better forecasts based on MAFE, RMSFE and MAPFE measures. In addition, credibility regression methods had a very good forecasting performance, when we applied them to the datasets of other countries (with a relatively small population5) for a limited selected fitting period (1980–2000), such as Belgium, Finland, Norway, Ireland, Slovakia and New Zealand. A further forecasting comparison between datasets of other countries has been left for future work.

Finally, we have to mention that our numerical illustration yielded results that are fully applicable and provide encouragement that credibility modelling approaches, including those of Tsai (2017a, 2017b), could contribute to future mortality projection studies.