# Random cohort effects and age groups dependency structure for mortality
modelling and forecasting: Mixed-effects time-series model approach

Ka Kin Lam111Corresponding author. E-mail<EMAIL_ADDRESS>(K.K.Lam)., Bo
Wang

School of Mathematics and Actuarial Science, University of Leicester,
Leicester LE1 7RH, UK

> Abstract
>
> There have been significant efforts devoted to solving the longevity risk
> given that a continuous growth in population ageing has become a severe
> issue for many developed countries over the past few decades. The Cairns-
> Blake-Dowd (CBD) model, which incorporates cohort effects parameters in its
> parsimonious design, is one of the most well-known approaches for mortality
> modelling at higher ages and longevity risk. This article proposes a novel
> mixed-effects time-series approach for mortality modelling and forecasting
> with considerations of age groups dependence and random cohort effects
> parameters. The proposed model can disclose more mortality data information
> and provide a natural quantification of the model parameters uncertainties
> with no pre-specified constraint required for estimating the cohort effects
> parameters. The abilities of the proposed approach are demonstrated through
> two applications with empirical male and female mortality data. The proposed
> approach shows remarkable improvements in terms of forecast accuracy
> compared to the CBD model in the short-, mid-and long-term forecasting using
> mortality data of several developed countries in the numerical examples.

Keywords: Demographic modelling; Mortality forecasting; Cohort effects; Time
series analysis; Mixed-effects model; CBD model

##  1 Introduction

Continuous growth in life expectancy has been lasting over the last few
decades in many developed countries. This unanticipated increase in lifespan
has posed significant challenges to many government sectors and life insurance
sectors due to the extra burdens on the health care services and the rapid
increase in pension and insurance expenditure. The problem caused by the
ageing population is also well-known as the ‘longevity risk’. Many demographic
researchers have recognised this problem and devoted significant efforts to
developing new statistical models for mortality modelling and forecasting. A
good statistical model for modelling mortality is always of great importance
to incorporate many different factors that help explain the historical
mortality patterns and capture their future trends. For example, cohort
effects in a statistical mortality model can be considered one of the most
important factors but is still underexplored and is missing from the majority
of mortality modelling literature, which we aim to discuss and address in this
article.

Cohort effects, also known as the year-of-birth effects, are arguably one of
the most important and discussed factors that influence mortality experience
among groups of individuals born in the same year or generation. Some evidence
for the existence of the cohort effects has been successfully spotted, and the
developments of the cohort studies can mainly be divided into two different
streams. One side relies on the empirical data analysis and the qualitative
analysis, such as descriptive and graphical representations, to conclude the
existence of cohort effects on mortality patterns in a population. For
example, Willets (2004) discovers that the ‘golden generations’ cohort who
were born between the year 1925 and the year 1945, have experienced
exceptionally rapid improvements in mortality rates through examining the
average annual mortality improvement rates for males and females in the
population of England and Wales. Murphy (2009) and Murphy (2010) further
explain this ‘golden generations’ cohort phenomenon that may be due to those
born during and after the Second World War with low fertility and thus facing
less competition for resources, better environmental conditions with changing
smoking habits and diets between generations. Li . (2016) demonstrate the
cohort surface across five different countries using graphical presentations
with kernel smoothing techniques and conclude that Japan and the UK are the
countries with relatively stronger cohort effects than the Netherlands, Spain
and the US.

Another side of the cohort studies depends on establishing a new statistical
model that offers quantitative approaches to identify the cohort patterns of
the observed mortality data. Early attempts can be traced back to the
classical age-period-cohort (APC) model in mortality modelling (Holford,
1983). Renshaw Haberman (2003) extend the well-known Lee-Carter model (Lee
Carter, 1992) by incorporating an age-modulated cohort effects parameter.
Cairns . (2009) extend their original two-factor Cairns-Blake-Dowd (CBD) model
(Cairns ., 2006) as a three-factor CBD model with cohort parameters. However,
the aforementioned statistical mortality models suffer the ‘age-period-cohort
identification problem’, given that the age, period and cohort effects are
intertwined222The age $x$ and period $t$ take values in
$\\{x_{1},\cdots,x_{m}\\}$ and $\\{t_{1},\cdots,t_{n}\\}$, the cohort index
$t-x$ takes values in $\\{t_{1}-x_{m},\cdots,t_{n}-x_{1}\\}$.. Including
indicator variables of each level of age, period and cohort respectively in a
fixed-effects estimation approach can lead to the perfect collinearity
problem, which means that there can be an infinite number of solutions that
fit the observed data equally well in the fixed-effects estimation approach
(O’Brien, 2014). The common solution to this age-period-cohort identification
problem relies on placing pre-specified constraints on the model coefficients
prior to the fitting procedure. However, such pre-specified constraints
meanwhile lead to the fitting results of the estimated model parameters
substantially depending on the initial constraints chosen (O’Brien, 2017).
Readers can refer to Haberman Renshaw (2011) and Currie (2016) for more
comprehensive summaries on these convergence and robustness problems among all
these fixed-effects estimation approaches applied to the morality models with
the cohort effects parameters.

In this article, we focus on the side of establishing a new statistical
mortality model that can extract and provide cohort effects from the observed
mortality data for modelling and forecasting. We propose a novel approach with
random cohort effects parameters and age groups dependency structure, which
can be seen as a mixed-effects time-series version of the three-factor CBD
model with the cohort effects parameters proposed by Cairns . (2009). The main
contributions of this article are to demonstrate the formulation and
parameters estimation of the proposed model with considerations of random
cohort effects parameters for mortality modelling and forecasting under a
united mixed-effects time-series framework with several desirable advantages
over the CBD model.

More will be discussed in detail and the rest of this article is organised as
follows. In Section 2, we firstly give a review of the CBD model. In Section
3, we introduce the formulation and the framework of the proposed model for
mortality modelling and forecasting. We then illustrate the proposed model
with empirical studies of mortality modelling and forecasting with comparisons
to the CBD model in terms of the systematic differences and forecasting
ability using mortality data of ten different developed countries in Section
4. We lastly conclude this article with discussions and remarks in Section 5.

##  2 Review of the Cairns-Blake-Dowd (CBD) model

In this section, we define some common actuarial notations used in the CBD
model first, and then go through a brief review of the CBD model.

###  2.1 Actuarial notations and their relationships

In the literature of mortality modelling, there are generally two types of
stochastic methods. One type of method models the central mortality rate. The
CBD model belongs to the other type of method which models the initial
mortality rate. For the sake of notation clarity, we start by defining the
notations using a subscript $i$ to denote the number of calendar year $t_{i}$
for $i=1,\cdots,n$, and a subscript $j$ as the number of age-specific group
$x_{j}$ for $j=1,\cdots,m$ in the CBD model. These notations include as
follows:

  * •

$D_{x_{j},t_{i}}$ is the observed number of deaths aged $x_{j}$ last birthday
in a calendar year $t_{i}$.

  * •

$E^{c}_{x_{j},t_{i}}$ is the central exposure to risk aged $x_{j}$ last
birthday in the middle of a calendar year $t_{i}$.

  * •

$m_{x_{j},t_{i}}$ is the central mortality rate which reflects the death
probability within the range of 0 to 1 for a population died in age $x_{j}$
last birthday in the middle of a calender year $t_{i}$. It is estimated by
$m_{x_{j},t_{i}}=D_{x_{j},t_{i}}/E^{c}_{x_{j},t_{i}}$.

  * •

$q_{x_{j},t_{i}}$ is the initial mortality rate as the one-year death
probability between $0$ and $1$ for a population died in age $x_{j}$ exactly
at time $t_{i}$.

The initial mortality rate $q_{x_{j},t_{i}}$ and the central mortality rate
$m_{x_{j},t_{i}}$ are typically very close to one another in practice (Dickson
., 2013). The link between the initial mortality rate and the central
mortality rate can be stated, i.e.

$q_{x_{j},t_{i}}\thickapprox 1-\text{exp}(-m_{x_{j},t_{i}}),$ (2.1)

where $0\leq q_{x_{j},t_{i}}\leq 1$.

###  2.2 Cairns-Blake-Dowd (CBD) model

The Cairns-Blake-Dowd (CBD) model is one of the most well-known approaches in
the mortality modelling literature. Cairns . (2009) extend their original two-
factor CBD model and propose a three-factor CBD model which incorporates the
cohort parameters based on the parsimonious design in their previous version.

We begin by using a subscript $i$ to denote the number of a calendar year
$t_{i}$ for $i=1,\cdots,n$, and a subscript $j$ for the number of an age-
specific group $x_{j}$ for $j=1,\cdots,m$. The CBD model has the form

$y_{x_{j},t_{i}}=\text{logit}(q_{x_{j},t_{i}})=\text{log}\bigg{(}\frac{q_{x_{j},t_{i}}}{1-q_{x_{j},t_{i}}}\bigg{)}=\kappa^{(1)}_{t_{i}}+{\color[rgb]{0,0,0}\kappa^{(2)}_{t_{i}}\times(x_{j}-\bar{x})}+\gamma^{(3)}_{t_{i}-x_{j}},$
(2.2)

where $y_{x_{j},t_{i}}$ is the logit of the initial mortality rate at age
$x_{j}$ in a calendar year $t_{i}$, $\kappa_{t_{i}}^{(1)}$ is the general time
effects intercept parameter, $\kappa_{t_{i}}^{(2)}$ is the time effects slope
parameter, $\bar{x}$ is the average of the observed mortality data age range,
i.e. $\bar{x}=\frac{1}{m}\sum_{j=1}^{m}x_{j}$, and
$\gamma^{(3)}_{t_{i}-x_{j}}$ is the cohort effects parameter that provides
mortality information on a group of people born in the same calendar year.

The CBD model relies on the linearity of the logit of the initial mortality
rate at older ages (age 60 to age 89) by assuming that $y_{x_{j},t_{i}}$ is a
response variable of a linear function. It treats $\kappa^{(1)}_{t_{i}}$ as an
intercept term, $\kappa^{(2)}_{t_{i}}$ as a slope term, and
$\gamma^{(3)}_{t_{i}-x_{j}}$ as a cohort effects term and assumes all of the
terms following stochastic processes for any given calendar year $t_{i}$ for
mortality modelling (Cairns ., 2009).

####  2.2.1 Model parameters estimation via the generalised linear modelling
approach

Instead of assuming the model error structure is normally distributed and
homoscedastic as in many ordinary least squares linear regression approaches,
the parameters in the CBD model are estimated based on an alternative approach
known as the generalised linear modelling (GLM) approach to calibrate the
historical mortality data to the model parameters.

We start with the explanation of the CBD model parameters estimation procedure
using the relationship between the initial mortality rate $q_{x_{j},t_{i}}$
and the central mortality rate $m_{x_{j},t_{i}}$ defined in Equation (2.1),
and the representation of the CBD model in Equation (2.2). It can be denoted
as

$\text{log}\bigg{(}\frac{1-\text{exp}(-m_{x_{j},t_{i}})}{1-(1-\text{exp}(-m_{x_{j},t_{i}}))}\bigg{)}=\kappa^{(1)}_{t_{i}}+{\color[rgb]{0,0,0}\kappa^{(2)}_{t_{i}}\times(x_{j}-\bar{x})}+\gamma^{(3)}_{t_{i}-x_{j}}.$
(2.3)

Rearranging the order of Equation (2.3) with respect to $m_{x_{j},t_{i}}$, we
have

$m_{x_{j},t_{i}}=\text{log}(1+\text{exp}{\color[rgb]{0,0,0}(\kappa^{(1)}_{t_{i}}+\kappa^{(2)}_{t_{i}}\times(x_{j}-\bar{x})+\gamma^{(3)}_{t_{i}-x_{j}})}).$

We denote above $m_{x_{j},t_{i}}$ as $m_{x_{j},t_{i}}(\theta_{i})$, where
$\theta_{i}=(\kappa_{t_{i}}^{(1)},\kappa_{t_{i}}^{(2)},{\color[rgb]{0,0,0}\gamma^{(3)}_{t_{i}-x_{j}}})$
represents the set of the model parameters depending on $m_{x_{j},t_{i}}$. The
CBD model assumes that the number of deaths $D_{x_{j},t_{i}}$ is a random
variable from count data following the Poisson distribution, i.e.

$D_{x_{j},t_{i}}\sim\text{Poisson}(E^{c}_{x_{j},t_{i}}\times
m_{x_{j},t_{i}}(\theta_{i})).$

We can thus state the log-likelihood function of $\boldsymbol{\theta}$
conditional on $\boldsymbol{D}$ and $\boldsymbol{E^{c}}$ across all age-
specific groups and calender years as

$\mathcal{L}(\boldsymbol{\theta}|{\color[rgb]{0,0,0}\boldsymbol{D},\boldsymbol{E^{c}}})=\sum_{i=1}^{n}\sum_{j=1}^{m}D_{x_{j},t_{i}}\times\text{log}{\color[rgb]{0,0,0}(}E^{c}_{x_{j},t_{i}}\times
m_{x_{j},t_{i}}(\theta_{i}){\color[rgb]{0,0,0})}-E^{c}_{x_{j},t_{i}}\times
m_{x_{j},t_{i}}(\theta_{i})-\text{log}(D_{x_{j},t_{i}}!),$ (2.4)

where $\boldsymbol{\theta}$ is the full set of the model parameters matrix,
$\boldsymbol{D}$ is the observed number of deaths matrix across all age-
specific groups and calender years, and $\boldsymbol{E^{c}}$ is the central
exposure matrix across all age-specific groups and calender years.

However, the full set of the model parameters cannot simply be calibrated by
maximising the above log-likelihood function to have a closed form solution as
the log-likelihood function follows the Poisson distribution. The model itself
also suffers model identification problems, in which the model parameters are
over-parametrised in Equation (2.2) since there can be two or more possible
model structures. For example, we can see that they are invariant to the
following transformation

$(\kappa^{(1)}_{t_{i}},\kappa^{(2)}_{t_{i}},\gamma^{(3)}_{t_{i}-x_{j}})\mapsto\Big{(}\kappa^{(1)}_{t_{i}}+\phi_{1}+{\color[rgb]{0,0,0}\phi_{2}\times(t_{i}-\bar{x})},\kappa^{(2)}_{t_{i}}-\phi_{2},\gamma^{(3)}_{t_{i}-x_{j}}-\phi_{1}-{\color[rgb]{0,0,0}\phi_{2}\times(t_{i}-x_{j})}\Big{)},$

where $\phi_{1}$ and $\phi_{2}$ are any non-zero constants.

Cairns . (2009) amend this identification problem by switching the cohort
effects parameter from $\gamma^{(3)}_{t_{i}-x_{j}}$ to
$\tilde{\gamma}^{(3)}_{t_{i}-x_{j}}=\phi_{1}+{\color[rgb]{0,0,0}\phi_{2}\times(t_{i}-x_{j})}+\gamma^{(3)}_{t_{i}-x_{j}}$
and imposing the following constraints on the cohort parameter
$\gamma^{(3)}_{t_{i}-x_{j}}$ throughout the iterative optimisation process for
the log-likelihood function in Equation (2.4), such that

$\sum_{(t_{i}-x_{j})\in C}\gamma^{(3)}_{t_{i}-x_{j}}=0,$
$\sum_{(t_{i}-x_{j})\in C}(t_{i}-x_{j})\times\gamma^{(3)}_{t_{i}-x_{j}}=0,$

where $C$ is the set of cohort years of birth $(t_{i}-x_{j})$ included in the
analysis. Cairns . (2009) explain the choice of the above constraints is
analogous to fitting a latent linear function of
$\phi_{1}+{\color[rgb]{0,0,0}\phi_{2}\times(t_{i}-x_{j})}$ to
$\gamma^{(3)}_{t_{i}-x_{j}}$ implicitly, so that the constraints ensure
$\hat{\phi}_{1}=0$ and $\hat{\phi}_{2}=0$, and the fitted
$\hat{\gamma}^{(3)}_{t_{i}-x_{j}}$ will fluctuate around zero and has no
discernible linear trend as a standard normal random variable. Cairns .
(2009), Haberman Renshaw (2011) and Currie (2016) provide further discussions
on these identifiability and parameter fitting issues for the three-factor CBD
model.

####  2.2.2 Out-of-sample forecasts of the mortality odds in the CBD model

For the out-of-sample forecasts of the initial mortality rate, Cairns . (2006)
assume that the fitted $\hat{\kappa}_{t_{i}}^{(1)}$ and
$\hat{\kappa}_{t_{i}}^{(2)}$ follow a bivariate random walk with drift which
can be described by defining
$\hat{\boldsymbol{\beta}}_{t_{i}}=\begin{pmatrix}\hat{\kappa}_{t_{i}}^{(1)}\\\
\hat{\kappa}_{t_{i}}^{(2)}\end{pmatrix}$, then we have

$\hat{\boldsymbol{\beta}}_{t_{i}}=\hat{\boldsymbol{\beta}}_{t_{i-1}}+\boldsymbol{d}+\boldsymbol{K}\boldsymbol{Z}_{{\color[rgb]{0,0,0}t_{i}}},$

where ${\boldsymbol{d}}$ is a $2\times 1$ vector of the drift factors and
$\boldsymbol{Z}_{{\color[rgb]{0,0,0}t_{i}}}$ is a $2\times 1$ vector of
bivariate normally distributed random variables, i.e.
$\boldsymbol{Z}_{{\color[rgb]{0,0,0}t_{i}}}\sim\mathcal{N}_{2\times
1}(\boldsymbol{0},\boldsymbol{I})$. Let $\boldsymbol{V}$ be a variance-
covariance matrix of $\hat{\kappa}_{t_{i}}^{(1)}$ and
$\hat{\kappa}_{t_{i}}^{(2)}$. Without loss of generality, $\boldsymbol{K}$ is
restricted to be an upper triangular matrix by Cholesky decomposition where
$\boldsymbol{V}=\boldsymbol{K}^{T}\boldsymbol{K}$. Denoting $\Delta$ as a
first-difference operator, such that
$\Delta\hat{\kappa}_{t_{i}}^{(l)}=\hat{\kappa}_{t_{i}}^{(l)}-\hat{\kappa}_{t_{i-1}}^{(l)}$
for $l=1,2$, we can estimate $\boldsymbol{d}$ and $\boldsymbol{V}$ as

$\boldsymbol{\hat{d}}=\begin{pmatrix}\mathbb{E}(\Delta\hat{\kappa}_{t_{i}}^{(1)})\\\
\mathbb{E}(\Delta\hat{\kappa}_{t_{i}}^{(2)})\end{pmatrix},$

and let
$\hat{\sigma}_{lk}=\mathbb{E}[(\Delta\hat{\kappa}_{t_{i}}^{(l)}-\mathbb{E}(\Delta\hat{\kappa}_{t_{i}}^{(l)}))^{T}(\Delta\hat{\kappa}_{t_{i}}^{(k)}-\mathbb{E}(\Delta\hat{\kappa}_{t_{i}}^{(k)}))]$
for $l,k=1,2$,

$\boldsymbol{\hat{V}}=\begin{pmatrix}\hat{\sigma}_{11}&\hat{\sigma}_{12}\\\
\hat{\sigma}_{21}&\hat{\sigma}_{22}\end{pmatrix}.$

By Cholesky decomposition of $\hat{\boldsymbol{V}}$, we have

$\boldsymbol{\hat{V}}=\boldsymbol{\hat{K}}^{T}\boldsymbol{\hat{K}}.$

The $h$-step ahead extrapolation of $\hat{\boldsymbol{\beta}}_{t_{i}}$ is then

${\color[rgb]{0,0,0}\hat{\boldsymbol{\beta}}_{t_{i+h}}=\hat{\boldsymbol{\beta}}_{t_{i+h-1}}+\boldsymbol{\hat{d}},}$

and the variance of the $h$-step ahead extrapolation of
$\hat{\boldsymbol{\beta}}_{t_{i}}$ is

${\text{Var}({\boldsymbol{\beta}}_{t_{i+h}})=h\times\boldsymbol{\hat{V}}.}$

The cohort effects parameter
$\hat{\gamma}^{(3)}_{{\color[rgb]{0,0,0}t_{i}-x_{j}}}$ is assumed to follow a
univariate random walk with drift (Haberman Renshaw, 2011), hence the $h$-step
ahead extrapolation of $\hat{\gamma}_{{\color[rgb]{0,0,0}t_{i}-x_{j}}}^{(3)}$
is

${\color[rgb]{0,0,0}\hat{\gamma}^{(3)}_{t_{i+h}-x_{j}}=\hat{\gamma}^{(3)}_{t_{i+h-1}-x_{j}}+\hat{\mu},}$

where $\hat{\mu}=\mathbb{E}(\Delta\hat{\gamma}^{(3)}_{t_{i}-x_{j}})$.

Combining all the extrapolated model parameters in the CBD model, the $h$-step
ahead out-of-sample forecast of $y_{x_{j},t_{i}}$ is

$\hat{y}_{x_{j},t_{i+h}}=\hat{\kappa}^{(1)}_{t_{i+h}}+{\color[rgb]{0,0,0}\hat{\kappa}^{(2)}_{t_{i+h}}\times(x_{j}-\bar{x})}+\hat{\gamma}^{(3)}_{t_{i+h}-x_{j}}.$

##  3 Methodology

###  3.1 Mixed-effects time-series model with age groups dependency and random
cohort effects parameters

In this section, we introduce the proposed mixed-effects time-series model
with age groups dependency and random cohort effects parameters.

Let our observed mortality dataset be $\\{t_{i},x_{j},y_{x_{j},t_{i}}\\}$ for
$i=1,\ldots,n$ and $j=1,\ldots,m$, where $n$ is the total number of observed
calendar years, $m$ is the maximum number of age-specific groups of interest,
$y_{x_{j},t_{i}}$ is the logit of the initial mortality rate for age $x_{j}$
in a calendar year $t_{i}$. The proposed mixed-effects time-series model is

$y_{x_{j},t_{i}}=\text{logit}(q_{x_{j},t_{i}})=\text{log}\bigg{(}\frac{q_{x_{j},t_{i}}}{1-q_{x_{j},t_{i}}}\bigg{)}=f_{x_{j}}(t_{i}-\bar{t})+\epsilon_{x_{j},t_{i}},\hskip
14.22636pt\epsilon_{x_{j},t_{i}}\sim\mathcal{N}(0,\sigma^{2}),$

where $f_{x_{j}}(t_{i}-\bar{t})$ is an underlying function which needs to be
estimated for a given age $x_{j}$, $\bar{t}$ is the average of the observed
calendar years, i.e. $\bar{t}=\frac{1}{n}\sum_{i=1}^{n}t_{i}$, and
$\\{\epsilon_{x_{j},t_{i}}\\}_{i,j=1}^{n,m}$ allow the observation error
varying among the different age groups and the calendar years based on the
assumption of i.i.d. normally distributed random variables with mean 0 and
constant variance $\sigma^{2}$.

The underlying function $f_{x_{j}}(t_{i}-\bar{t})$ in the proposed model under
the mixed-effects structure is further specified as

$f_{x_{j}}(t_{i}-\bar{t})=\mu(t_{i}-\bar{t})+\eta_{x_{j}}(t_{i}-\bar{t}),$

where $\mu(t_{i}-\bar{t})$ is a general linear mean function across all age
groups $\\{x_{j}\\}_{j=1}^{m}$, i.e.

$\mu(t_{i}-\bar{t})=\beta^{(1)}+\beta^{(2)}\times(t_{i}-\bar{t}),$

with $\beta^{(1)}$ is the fixed-effects intercept parameter and $\beta^{(2)}$
is the fixed-effects slope parameter multiplying the input $(t_{i}-\bar{t})$.

Meanwhile, $\eta_{x_{j}}(t_{i}-\bar{t})$ is the age-specific deviation
function from the general linear mean function across all age groups, which is
designed to be a function with the following components

$\eta_{x_{j}}(t_{i}-\bar{t})=f_{x_{j}}(t_{i}-\bar{t})-\mu(t_{i}-\bar{t})=\gamma_{x_{j}}^{(1)}+\gamma_{x_{j}}^{(2)}\times(t_{i}-\bar{t})+\gamma^{(3)}_{t_{i}-x_{j}},$
(3.1)

where we assume that $\\{\gamma_{x_{j}}^{(1)}\\}_{j=1}^{m}$ are the correlated
age-specific random effects intercept parameters following a prior
distribution
${(\gamma_{x_{1}}^{(1)},\ldots,\gamma_{x_{m}}^{(1)})}^{T}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{K}^{(1)})$,
$\\{\gamma_{x_{j}}^{(2)}\\}_{j=1}^{m}$ are the correlated age-specific random
effects slope parameters following a prior distribution
${(\gamma_{x_{1}}^{(2)},\ldots,\gamma_{x_{m}}^{(2)})}^{T}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{K}^{(2)})$,
and $\\{\gamma_{t_{i}-x_{j}}^{(3)}\\}_{t_{1}-x_{m}}^{t_{n}-x_{1}}$ are the
correlated age-specific random effects cohort parameters following a prior
distribution
$(\gamma_{t_{1}-x_{m}}^{(3)},\ldots,\gamma_{t_{n}-x_{1}}^{(3)})^{T}\sim\mathcal{N}(\boldsymbol{0},\boldsymbol{K}^{(3)})$
333Without any loss of completeness, the covariance structures of all the
random effects components $\gamma_{x_{j}}^{(1)}$, ${\gamma_{x_{j}}^{(2)}}$ and
$\gamma^{(3)}_{t_{i}-x_{j}}$ will be discussed in Section 3.2..

With all the designed components mentioned above, we can state the full form
of the proposed model as

$y_{x_{j},t_{i}}=\beta^{(1)}+{\color[rgb]{0,0,0}\beta^{(2)}\times(t_{i}-\bar{t})}+\gamma_{x_{j}}^{(1)}+{\color[rgb]{0,0,0}\gamma_{x_{j}}^{(2)}\times(t_{i}-\bar{t})}+\gamma^{(3)}_{t_{i}-x_{j}}+\epsilon_{x_{j},t_{i}}.$
(3.2)

Alternatively, we can group
$\kappa_{x_{j}}^{(1)}={\color[rgb]{0,0,0}\beta^{(1)}+\gamma_{x_{j}}^{(1)}}$
and
$\kappa_{x_{j}}^{(2)}={\color[rgb]{0,0,0}\beta^{(2)}+\gamma_{x_{j}}^{(2)}}$
for the simplicity of notation, and the proposed model is thus

$y_{x_{j},t_{i}}=\kappa_{x_{j}}^{(1)}+{\color[rgb]{0,0,0}\kappa_{x_{j}}^{(2)}\times(t_{i}-\bar{t})}+\gamma^{(3)}_{t_{i}-x_{j}}+\epsilon_{x_{j},t_{i}}.$

###  3.2 Formulation and parameters estimation of the proposed model

In this section, we begin by describing the formulation of the proposed model
and explain the procedure of the model parameters estimation.

####  3.2.1 Model formulation

Letting
$\boldsymbol{Y}_{{x}_{j}}={\color[rgb]{0,0,0}[{y}_{x_{j},t_{1}},{y}_{x_{j},t_{2}},\ldots,{y}_{x_{j},t_{n}}]^{T}}$,
$\boldsymbol{T}_{{x}_{j}}=[(t_{1}-\bar{t}),(t_{2}-\bar{t}),\ldots,(t_{n}-\bar{t})]^{T}$
and
$\boldsymbol{\varepsilon}_{x_{j}}={\color[rgb]{0,0,0}[\varepsilon_{x_{j},t_{1}},\varepsilon_{x_{j},t_{2}},\ldots,\varepsilon_{x_{j},t_{n}}]^{T}}$,
we can denote Equation (3.2) in matrix form as

$\boldsymbol{Y}=\boldsymbol{T}\boldsymbol{\beta}+\boldsymbol{Z}^{(1)}\boldsymbol{\gamma}^{(1)}_{\boldsymbol{x}}+\boldsymbol{Z}^{(2)}\boldsymbol{\gamma}^{(2)}_{\boldsymbol{x}}+\boldsymbol{Z}^{(3)}\boldsymbol{\gamma}^{(3)}_{\boldsymbol{t-x}}+\boldsymbol{\epsilon},\hskip
14.22636pt\boldsymbol{\epsilon}\sim
N(\boldsymbol{0},\sigma^{2}\boldsymbol{I}),$ (3.3)

where

$\boldsymbol{Y}=\begin{bmatrix}\boldsymbol{Y}_{{x}_{1}}\\\
\boldsymbol{Y}_{{x}_{2}}\\\ \vdots\\\
\boldsymbol{Y}_{{x}_{m}}\end{bmatrix}\in\mathbb{R}^{nm\times
1},\boldsymbol{T}=\begin{bmatrix}\boldsymbol{1}_{n\times
1}&\boldsymbol{T}_{{x}_{1}}\\\ \boldsymbol{1}_{n\times
1}&\boldsymbol{T}_{{x}_{2}}\\\ \vdots&\vdots\\\ \boldsymbol{1}_{n\times
1}&\boldsymbol{T}_{{x}_{m}}\end{bmatrix}\in\mathbb{R}^{nm\times
2},\boldsymbol{\beta}=\begin{bmatrix}\beta^{(1)}\\\
\beta^{(2)}\end{bmatrix}\in\mathbb{R}^{2\times 1},$
$\boldsymbol{Z}^{(1)}=\begin{bmatrix}\boldsymbol{1}_{n\times
1}&0&\ldots&0&0\\\ 0&\boldsymbol{1}_{n\times 1}&\ldots&0&0\\\
\vdots&\vdots&\ddots&\vdots&\vdots\\\ 0&0&\ldots&0&\boldsymbol{1}_{n\times
1}\\\ \end{bmatrix}\in\mathbb{R}^{nm\times
m},\boldsymbol{\gamma}^{(1)}_{\boldsymbol{x}}=\begin{bmatrix}\gamma^{(1)}_{x_{1}}\\\
\gamma^{(1)}_{x_{2}}\\\ \vdots\\\ \gamma^{(1)}_{x_{m}}\\\
\end{bmatrix}\in\mathbb{R}^{m\times 1},$
$\boldsymbol{Z}^{(2)}=\begin{bmatrix}\boldsymbol{T}_{x_{1}}&0&\ldots&0&0\\\
0&\boldsymbol{T}_{x_{2}}&\ldots&0&0\\\ \vdots&\vdots&\ddots&\vdots&\vdots\\\
0&0&\ldots&0&\boldsymbol{T}_{x_{m}}\\\ \end{bmatrix}\in\mathbb{R}^{nm\times
m},\boldsymbol{\gamma}^{(2)}_{\boldsymbol{x}}=\begin{bmatrix}\gamma^{(2)}_{x_{1}}\\\
\gamma^{(2)}_{x_{2}}\\\ \vdots\\\ \gamma^{(2)}_{x_{m}}\\\
\end{bmatrix}\in\mathbb{R}^{m\times 1},$
$\boldsymbol{Z}^{(3)}=\hbox{}\vbox{\kern 0.86108pt\hbox{$\kern 0.0pt\kern
2.5pt\kern-5.0pt\left[\kern
0.0pt\kern-2.5pt\kern-5.55557pt\vbox{\kern-0.86108pt\vbox{\vbox{
\halign{\kern\arraycolsep\hfil\@arstrut$\kbcolstyle#$\hfil\kern\arraycolsep&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep&&
\kern\arraycolsep\hfil$\@kbrowstyle#$\ifkbalignright\relax\else\hfil\fi\kern\arraycolsep\cr
5.0pt\hfil\@arstrut$\scriptstyle$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle
t_{1}-x_{m}$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle t_{2}-x_{m}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\ldots$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle
t_{n}-x_{m}$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle t_{n}-x_{m-1}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\ldots$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle
t_{n}-x_{1}\\\t_{1}-x_{1}$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle
0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\ldots$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle
0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\ldots$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle
0\\\t_{2}-x_{1}$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\ldots$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle
0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\ldots$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle
0\\\\\vdots$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\vdots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\vdots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\ddots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\vdots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\vdots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\ddots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\vdots\\\t_{n}-x_{1}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle
0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\ldots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle
0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\ldots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle 1\\\t_{1}-x_{2}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle
0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\ldots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle
0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\ldots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle 0\\\t_{2}-x_{2}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle
0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\ldots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle
0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\ldots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle 0\\\\\vdots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\vdots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\vdots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\ddots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\vdots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\vdots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\ddots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\vdots\\\t_{n}-x_{2}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle
0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\ldots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle
0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\ldots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle 0\\\\\vdots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\vdots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\vdots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\ddots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\vdots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\vdots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\ddots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\vdots\\\t_{1}-x_{m-1}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle
0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\ldots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle
0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\ldots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle 0\\\t_{2}-x_{m-1}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle
0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\ldots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle
0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\ldots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle 0\\\\\vdots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\vdots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\vdots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\ddots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\vdots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\vdots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\ddots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\vdots\\\t_{n}-x_{m-1}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle
0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\ldots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle
1$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\ldots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle 0\\\t_{1}-x_{m}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle 1$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle
0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\ldots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle
0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\ldots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle 0\\\t_{2}-x_{m}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle
1$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\ldots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle
0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\ldots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle 0\\\\\vdots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\vdots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\vdots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\ddots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\vdots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\vdots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\ddots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle\vdots\\\t_{n}-x_{m}$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle 0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle
0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\ldots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle 1$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle
0$\hfil\kern 5.0pt&5.0pt\hfil$\scriptstyle\ldots$\hfil\kern
5.0pt&5.0pt\hfil$\scriptstyle 0\\\$\hfil\kern
5.0pt\crcr}}}}\right]$}}\in\mathbb{R}^{nm\times(n+m-1)},$
$\boldsymbol{\gamma}^{(3)}_{\boldsymbol{t-x}}=\begin{bmatrix}\gamma^{(3)}_{t_{1}-x_{m}}\\\
\gamma^{(3)}_{t_{2}-x_{m}}\\\ \vdots\\\ \gamma^{(3)}_{t_{n}-x_{m}}\\\
\gamma^{(3)}_{t_{n}-x_{m-1}}\\\ \vdots\\\ \gamma^{(3)}_{t_{n}-x_{1}}\\\
\end{bmatrix}\in\mathbb{R}^{(n+m-1)\times 1},\text{and}\hskip
5.69046pt\boldsymbol{\epsilon}=\begin{bmatrix}\boldsymbol{\epsilon}_{{x}_{1}}\\\
\boldsymbol{\epsilon}_{{x}_{2}}\\\ \vdots\\\
\boldsymbol{\epsilon}_{{x}_{m}}\end{bmatrix}\in\mathbb{R}^{nm\times 1}.$

In aid of explaining the above full matrix representation of the proposed
model, here we consider a smaller matrix of cells where the row and column of
the matrix directly correspond to year $t$ and age $x$ without the use of
script for illustration purposes. The cohort effects parameter
$\gamma_{t-x}^{(3)}$ is indexed by the year-of-birth and its corresponding
value on each cell is displayed in Table 3.1 for the cohort effects. We can
see that the index values of $\gamma_{t-x}^{(3)}$ are the same in the ‘cohort
direction’, which are on the cells $(x,t)$, $(x+1,t+1)$ and so on. The total
number of cohort effects parameters is equal to the number of rows plus the
number of columns minus one, which is four in our illustration here. The
proposed model with the smaller size of matrix cells can be expressed as

$\begin{bmatrix}{y}_{{\color[rgb]{0,0,0}x_{1},t_{1}}}\\\
{y}_{{\color[rgb]{0,0,0}x_{1},t_{2}}}\\\
{y}_{{\color[rgb]{0,0,0}x_{2},t_{1}}}\\\
{y}_{{\color[rgb]{0,0,0}x_{2},t_{2}}}\\\
{y}_{{\color[rgb]{0,0,0}x_{3},t_{1}}}\\\
{y}_{{\color[rgb]{0,0,0}x_{3},t_{2}}}\end{bmatrix}=\begin{bmatrix}1&t_{1}\\\
1&t_{2}\\\ 1&t_{1}\\\ 1&t_{2}\\\ 1&t_{1}\\\
1&t_{2}\end{bmatrix}\begin{bmatrix}\beta^{(1)}\\\ \beta^{(2)}\\\
\end{bmatrix}+\begin{bmatrix}1&0&0\\\ 1&0&0\\\ 0&1&0\\\ 0&1&0\\\ 0&0&1\\\
0&0&1\end{bmatrix}\begin{bmatrix}\gamma^{(1)}_{x_{1}}\\\
\gamma^{(1)}_{x_{2}}\\\ \gamma^{(1)}_{x_{3}}\\\
\end{bmatrix}+\begin{bmatrix}t_{1}&0&0\\\ t_{2}&0&0\\\ 0&t_{1}&0\\\
0&t_{2}&0\\\ 0&0&t_{1}\\\
0&0&t_{2}\end{bmatrix}\begin{bmatrix}\gamma^{(2)}_{x_{1}}\\\
\gamma^{(2)}_{x_{2}}\\\ \gamma^{(2)}_{x_{3}}\\\
\end{bmatrix}+\begin{bmatrix}0&0&1&0\\\ 0&0&0&1\\\ 0&1&0&0\\\ 0&0&1&0\\\
1&0&0&0\\\ 0&1&0&0\end{bmatrix}\begin{bmatrix}\gamma^{(3)}_{-2}\\\
\gamma^{(3)}_{-1}\\\ \gamma^{(3)}_{0}\\\ \gamma^{(3)}_{1}\\\
\end{bmatrix}+\begin{bmatrix}{\epsilon}_{{\color[rgb]{0,0,0}x_{1},t_{1}}}\\\
{\epsilon}_{{\color[rgb]{0,0,0}x_{1},t_{2}}}\\\
{\epsilon}_{{\color[rgb]{0,0,0}x_{2},t_{1}}}\\\
{\epsilon}_{{\color[rgb]{0,0,0}x_{2},t_{2}}}\\\
{\epsilon}_{{\color[rgb]{0,0,0}x_{3},t_{1}}}\\\
{\epsilon}_{{\color[rgb]{0,0,0}x_{3},t_{2}}}\end{bmatrix}.$

year $t$ $/$ age $x$ |  $x=1$ |  $x=2$ |  $x=3$
---|---|---|---
$t=1$ |  $\gamma_{0}^{(3)}$ |  $\gamma_{-1}^{(3)}$ |  $\gamma_{-2}^{(3)}$
$t=2$ |  $\gamma_{1}^{(3)}$ |  $\gamma_{0}^{(3)}$ |  $\gamma_{-1}^{(3)}$
Table 3.1: Values of the cohort effects parameters $\gamma_{t-x}^{(3)}$ in a
matrix of cells $(t,x)$.

For the stochastic components in the proposed model, we assume that the prior
distributions of
$\boldsymbol{\gamma}^{(1)}_{\boldsymbol{x}}\sim{\color[rgb]{0,0,0}\mathcal{N}(\boldsymbol{0},\boldsymbol{K}^{(1)})}$,
$\boldsymbol{\gamma}^{(2)}_{\boldsymbol{x}}\sim{\color[rgb]{0,0,0}\mathcal{N}(\boldsymbol{0},\boldsymbol{K}^{(2)})}$
and
$\boldsymbol{\gamma}^{(3)}_{\boldsymbol{t-x}}\sim{\color[rgb]{0,0,0}\mathcal{N}(\boldsymbol{0},\boldsymbol{K}^{(3)})}$.
$\boldsymbol{K}^{(1)}$ and $\boldsymbol{K}^{(2)}$ are both $m\times m$ squared
exponential covariance matrices which are used to model the dependence
structure among different age groups. They have the following form, for
$i=1,2$, i.e.

${\color[rgb]{0,0,0}\boldsymbol{K}^{(i)}}={\color[rgb]{0,0,0}h_{i}^{2}}\text{exp}\Bigg{(}-\begin{bmatrix}\frac{(x_{1}-x_{1})^{2}}{2{\color[rgb]{0,0,0}l_{i}}}&\frac{(x_{1}-x_{2})^{2}}{2{\color[rgb]{0,0,0}l_{i}}}&\ldots&\frac{(x_{1}-x_{m})^{2}}{2{\color[rgb]{0,0,0}l_{i}}}\\\
\frac{(x_{2}-x_{1})^{2}}{2{\color[rgb]{0,0,0}l_{i}}}&\frac{(x_{2}-x_{2})^{2}}{2{\color[rgb]{0,0,0}l_{i}}}&\ldots&\frac{(x_{2}-x_{m})^{2}}{2{\color[rgb]{0,0,0}l_{i}}}\\\
\vdots&\vdots&\ddots&\vdots\\\
\frac{(x_{m}-x_{1})^{2}}{2{\color[rgb]{0,0,0}l_{i}}}&\frac{(x_{m}-x_{2})^{2}}{2{\color[rgb]{0,0,0}l_{i}}}&\ldots&\frac{(x_{m}-x_{m})^{2}}{2{\color[rgb]{0,0,0}l_{i}}}\end{bmatrix}\Bigg{)},$

where $h_{i}$ is the response-scale amplitude parameter determining the
variation of function values , and $l_{i}$ is the characteristic length-scale
parameter that gives smooth variations in a covariate-scale and controls how
close the varying age groups should correlate.

$\boldsymbol{K}^{(3)}$ is another $(n+m-1)\times(n+m-1)$ squared exponential
covariance matrix for $\boldsymbol{\gamma}^{(3)}_{\boldsymbol{t-x}}$ modelling
the dependence structure among different groups of people born in the same
year with the following form

${\color[rgb]{0,0,0}\boldsymbol{K}^{(3)}}=c^{2}\text{exp}\Bigg{(}-\begin{bmatrix}\frac{((t_{1}-x_{m})-(t_{1}-x_{m}))^{2}}{2s}&\frac{((t_{1}-x_{m})-(t_{2}-x_{m}))^{2}}{2s}&\ldots&\frac{((t_{1}-x_{m})-(t_{n}-x_{1}))^{2}}{2s}\\\
\frac{((t_{2}-x_{m})-(t_{1}-x_{m}))^{2}}{2s}&\frac{((t_{2}-x_{m})-(t_{2}-x_{m}))^{2}}{2s}&\ldots&\frac{((t_{2}-x_{m})-(t_{n}-x_{1}))^{2}}{2s}\\\
\vdots&\vdots&\ddots&\vdots\\\
\frac{((t_{n}-x_{1})-(t_{1}-x_{m}))^{2}}{2s}&\frac{((t_{n}-x_{1})-(t_{2}-x_{m}))^{2}}{2s}&\ldots&\frac{((t_{n}-x_{1})-(t_{n}-x_{1}))^{2}}{2s}\end{bmatrix}\Bigg{)},$

where $c$ is the response-scale amplitude parameter that determines the
variation of function values, and $s$ is the characteristic length-scale
parameter that gives smooth variations in a covariate-scale and controls how
close the different groups of people born in the same year should correlate.

The main justification of choosing the squared exponential covariance matrices
here for all the stochastic components is to reflect the dependencies among
different age groups properly as we assume that people who are closer in age
groups or birth year would tend to have more similarities in their mortality
patterns than those who are in distant age ranges or were born in a different
generation.

With the prior distribution assumptions among all the random effects
components, Equation (3.3) can be restated as a multivariate normal
distribution among all its random effects components, i.e.

$\boldsymbol{Y}=\boldsymbol{T}\boldsymbol{\beta}+\boldsymbol{Z}^{(1)}\boldsymbol{\gamma}^{(1)}_{\boldsymbol{x}}+\boldsymbol{Z}^{(2)}\boldsymbol{\gamma}^{(2)}_{\boldsymbol{x}}+\boldsymbol{Z}^{(3)}\boldsymbol{\gamma}^{(3)}_{\boldsymbol{t-x}}+\boldsymbol{\epsilon},$
(3.4)

where

$\begin{bmatrix}\boldsymbol{\gamma}^{(1)}_{\boldsymbol{x}}\\\
\boldsymbol{\gamma}^{(2)}_{\boldsymbol{x}}\\\
\boldsymbol{\gamma}^{(3)}_{\boldsymbol{t-x}}\\\
\boldsymbol{\varepsilon}\end{bmatrix}\sim\mathcal{N}\Bigg{(}\begin{bmatrix}\boldsymbol{0}\\\
\boldsymbol{0}\\\ \boldsymbol{0}\\\
\boldsymbol{0}\end{bmatrix},\begin{bmatrix}{\color[rgb]{0,0,0}\boldsymbol{K}^{(1)}}&\boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0}\\\
\boldsymbol{0}&{\color[rgb]{0,0,0}\boldsymbol{K}^{(2)}}&\boldsymbol{0}&\boldsymbol{0}\\\
\boldsymbol{0}&\boldsymbol{0}&{\color[rgb]{0,0,0}\boldsymbol{K}^{(3)}}&\boldsymbol{0}\\\
\boldsymbol{0}&\boldsymbol{0}&\boldsymbol{0}&\sigma^{2}\boldsymbol{I}\\\
\end{bmatrix}\bigg{)}.$

####  3.2.2 Model parameters inference

With the model formulation specified in the previous section, we can infer the
model parameters by the Bayesian paradigm. The probability distribution of
$\boldsymbol{Y}$ in Equation (3.4) is given by
$\mathcal{N}(\boldsymbol{T}\boldsymbol{\beta},\boldsymbol{V(\theta)})$ where
$\boldsymbol{V(\theta)}=\boldsymbol{Z}^{(1)}{\color[rgb]{0,0,0}\boldsymbol{K}^{(1)}}{\boldsymbol{Z}^{(1)}}^{T}+\boldsymbol{Z}^{(2)}{\color[rgb]{0,0,0}\boldsymbol{K}^{(2)}}{\boldsymbol{Z}^{(2)}}^{T}+\boldsymbol{Z}^{(3)}{\color[rgb]{0,0,0}\boldsymbol{K}^{(3)}}{\boldsymbol{Z}^{(3)}}^{T}+\sigma^{2}\boldsymbol{I}$.
$\boldsymbol{V(\theta)}$ is the $nm\times nm$ variance-covariance matrix
parametrised by the (hyper-) parameters
$\boldsymbol{\theta}=[h^{(1)},l^{(1)},h^{(2)},l^{(2)},c,s]^{T}$ involved in
its components covariance matrices ${\color[rgb]{0,0,0}\boldsymbol{K}^{(1)}}$,
${\color[rgb]{0,0,0}\boldsymbol{K}^{(2)}}$,
${\color[rgb]{0,0,0}\boldsymbol{K}^{(3)}}$ and the random error variance
$\sigma^{2}$. Hence the log-likelihood function of $\boldsymbol{\beta}$, the
parameters of the variance-covariance matrix $\boldsymbol{\theta}$ and the
random error variance $\sigma^{2}$ is

$\mathcal{L}(\boldsymbol{\beta},\boldsymbol{\theta},\sigma^{2})=-\frac{1}{2}\text{log}|\boldsymbol{V(\theta)}|-\frac{1}{2}(\boldsymbol{Y}-\boldsymbol{T}\boldsymbol{\beta})^{T}\boldsymbol{V(\theta)}^{-1}(\boldsymbol{Y}-\boldsymbol{T}\boldsymbol{\beta})-\frac{nm}{2}\text{log}(2\pi),$
(3.5)

where $|\cdot|$ is the determinant of a matrix.

Given that there is no analytical closed-form solution to this maximisation
problem, we need to resort to some standard gradient-based numerical
optimisation techniques, such as the Newton-Raphson method or the Conjugate
Gradient method, to maximise the log-likelihood function
$\mathcal{L}(\boldsymbol{\beta},\boldsymbol{\theta},\sigma^{2})$ in Equation
(3.5) for obtaining the optimal estimates of the model parameters.

To set the form of the model parameters appropriate for the iterative process
of the gradient-based optimisation method, we seek to differentiate the log-
likelihood function
$\mathcal{L}(\boldsymbol{\beta},\boldsymbol{\theta},\sigma^{2})$ with respect
to its elements $\boldsymbol{\beta}$, $\boldsymbol{\theta}$, and $\sigma^{2}$.
This yields the following partial derivatives

$\frac{\partial\mathcal{L}(\boldsymbol{\beta},\boldsymbol{\theta},\sigma^{2})}{\partial\boldsymbol{\beta}}=-(\boldsymbol{T}^{T}\boldsymbol{V(\theta)}^{-1}\boldsymbol{T}\boldsymbol{\beta}-\boldsymbol{T}^{T}\boldsymbol{V(\theta)}^{-1}\boldsymbol{Y}),$
(3.6)
$\displaystyle\frac{\partial\mathcal{L}(\boldsymbol{\beta},\boldsymbol{\theta},\sigma^{2})}{\partial\theta_{i}}$
$\displaystyle=-\dfrac{1}{2}\bigg{\\{}\text{tr}\bigg{(}\boldsymbol{V(\theta)}^{-1}\frac{\partial\boldsymbol{V(\theta)}}{\partial\theta_{i}}\bigg{)}+(\boldsymbol{Y}-\boldsymbol{T}\boldsymbol{\beta})^{T}\bigg{(}\boldsymbol{V(\theta)}^{-1}\frac{\partial\boldsymbol{V(\theta)}}{\partial\theta_{i}}\boldsymbol{V(\theta)}^{-1}\bigg{)}(\boldsymbol{Y}-\boldsymbol{T}\boldsymbol{\beta})\bigg{\\}},$
$\frac{\partial\mathcal{L}(\boldsymbol{\beta},\boldsymbol{\theta},\sigma^{2})}{\partial\sigma^{2}}=-\dfrac{1}{2}\bigg{\\{}\text{tr}(\boldsymbol{V(\theta)}^{-1})+(\boldsymbol{Y}-\boldsymbol{T}\boldsymbol{\beta})^{T}(\boldsymbol{V(\theta)}^{-1}\boldsymbol{V(\theta)}^{-1})(\boldsymbol{Y}-\boldsymbol{T}\boldsymbol{\beta})\bigg{\\}},$

where $\theta_{i}$ is a generic notation of the $i^{\text{th}}$ element of
parameters in the vector $\boldsymbol{\theta}$.

Setting Equation (3.6) equal to zero and solving it give us the closed-form
maximum likelihood estimator $\hat{\boldsymbol{\beta}}$, i.e.

$\hat{\boldsymbol{\beta}}=(\boldsymbol{T}^{T}\boldsymbol{V(\theta)}^{-1}\boldsymbol{T})^{-1}\boldsymbol{T}^{T}\boldsymbol{V(\theta)}^{-1}\boldsymbol{Y}.$
(3.7)

The typical iterative procedure in the gradient-based optimiser for estimating
${{\color[rgb]{0,0,0}\boldsymbol{\beta}}}$, $\boldsymbol{\theta}$ and
$\sigma^{2}$ with the above partial derivatives involves the following four
steps:

  1. 1.

Assign initial values to
${\boldsymbol{\hat{\theta}}}_{(1)}=[\hat{h}_{(1)}^{(1)},\hat{l}_{(1)}^{(1)},\hat{h}_{(1)}^{(2)},\hat{l}_{(1)}^{(2)},\hat{c}_{(1)},\hat{s}_{(1)}]^{T}$
and $\hat{\sigma}^{2}_{(1)}$.

  2. 2.

Substitute $\boldsymbol{\hat{\theta}}_{(k)}$ and $\hat{\sigma}^{2}_{(k)}$ to
solve for $\hat{\boldsymbol{\beta}}_{(k)}$ in Equation (3.7) and calculate the
value of the following log-likelihood function in each $k^{\text{th}}$ step
iteration:

$\displaystyle\mathcal{L}(\boldsymbol{\hat{\beta}}_{(k)},\boldsymbol{\hat{\theta}}_{(k)},\hat{\sigma}^{2}_{(k)})=-\frac{1}{2}\text{log}|\boldsymbol{V}(\boldsymbol{\hat{\theta}}_{(k)})|-\frac{1}{2}(\boldsymbol{Y}-\boldsymbol{T}\boldsymbol{\hat{\beta}}_{(k)})^{T}\boldsymbol{V}(\boldsymbol{\hat{\theta}}_{(k)})^{-1}(\boldsymbol{Y}-\boldsymbol{T}\boldsymbol{\hat{\beta}}_{(k)})-\frac{nm}{2}\text{log}(2\pi).$

  3. 3.

Update
${\hat{\theta}}_{i,(k+1)}={\hat{\theta}}_{i,(k)}+\frac{{\partial\mathcal{L}(\boldsymbol{\hat{\beta}}_{(k)},\boldsymbol{\hat{\theta}}_{(k)},\hat{\sigma}^{2}_{(k)})/\partial\hat{\theta}_{i,(k)}}}{{\partial^{2}\mathcal{L}(\boldsymbol{\hat{\beta}}_{(k)},\boldsymbol{\hat{\theta}}_{(k)},\hat{\sigma}^{2}_{(k)})/\partial^{2}\hat{\theta}_{i,(k)}}}$
and
$\hat{\sigma}^{2}_{(k+1)}=\hat{\sigma}^{2}_{(k)}+\frac{{\partial\mathcal{L}(\boldsymbol{\hat{\beta}}_{(k)},\boldsymbol{\hat{\theta}}_{(k)},\hat{\sigma}^{2}_{(k)})/\partial\hat{\sigma}^{2}_{(k)}}}{\partial^{2}\mathcal{L}(\boldsymbol{\hat{\beta}}_{(k)},\boldsymbol{\hat{\theta}}_{(k)},\hat{\sigma}^{2}_{(k)})/\partial^{2}\hat{\sigma}^{2}_{(k)}}$
based on the gradient ascent criteria.

  4. 4.

Repeat Step 2 and Step 3 until the above log-likelihood function converges
444The performances of the gradient-based optimisation approaches may depend
on the initial values supplied. Several optimisations running can be applied
with different initial values may be necessary if no prior knowledge about the
function to optimise..

Once the estimated $\boldsymbol{\hat{\theta}}_{(p)}$, $\hat{\sigma}^{2}_{(p)}$
and $\boldsymbol{V}(\boldsymbol{\hat{\theta}}_{(p)})$ are obtained in the
$p^{\text{th}}$ step iteration where the log-likelihood function converges, we
have

$\begin{bmatrix}\boldsymbol{Y}\\\
\boldsymbol{\gamma}^{(1)}_{\boldsymbol{x}}\\\
\boldsymbol{\gamma}^{(2)}_{\boldsymbol{x}}\\\
\boldsymbol{\gamma}^{(3)}_{\boldsymbol{t-x}}\\\
\end{bmatrix}\sim\mathcal{N}\Bigg{(}\begin{bmatrix}\boldsymbol{T}\hat{\boldsymbol{\beta}}_{(p)}\\\
\boldsymbol{0}\\\ \boldsymbol{0}\\\
\boldsymbol{0}\end{bmatrix},\begin{bmatrix}\boldsymbol{V}(\boldsymbol{\hat{\theta}}_{(p)})&\boldsymbol{Z}^{(1)}{\color[rgb]{0,0,0}\hat{\boldsymbol{K}}^{(1)}}&\boldsymbol{Z}^{(2)}{\color[rgb]{0,0,0}\hat{\boldsymbol{K}}^{(2)}}&\boldsymbol{Z}^{(3)}{\color[rgb]{0,0,0}\hat{\boldsymbol{K}}^{(3)}}\\\
{\color[rgb]{0,0,0}\hat{\boldsymbol{K}}^{(1)}}{\boldsymbol{Z}^{(1)}}^{T}&{\color[rgb]{0,0,0}\hat{\boldsymbol{K}}^{(1)}}&\boldsymbol{0}&\boldsymbol{0}\\\
{\color[rgb]{0,0,0}\hat{\boldsymbol{K}}^{(2)}}{\boldsymbol{Z}^{(2)}}^{T}&\boldsymbol{0}&{\color[rgb]{0,0,0}\hat{\boldsymbol{K}}^{(2)}}&\boldsymbol{0}\\\
{\color[rgb]{0,0,0}\hat{\boldsymbol{K}}^{(3)}}{\boldsymbol{Z}^{(3)}}^{T}&\boldsymbol{0}&\boldsymbol{0}&{\color[rgb]{0,0,0}\hat{\boldsymbol{K}}^{(3)}}\\\
\end{bmatrix}\Bigg{)},$

and the conditional distributions of
$\boldsymbol{\gamma}^{(1)}_{\boldsymbol{x}}$,
$\boldsymbol{\gamma}^{(2)}_{\boldsymbol{x}}$,
$\boldsymbol{\gamma}^{(3)}_{\boldsymbol{t-x}}$ given $\boldsymbol{Y}$ follow
the multivariate Gaussian distribution, which are
$\boldsymbol{\gamma}^{(1)}_{\boldsymbol{x}}|\boldsymbol{Y}\sim\mathcal{N}(\boldsymbol{\hat{\gamma}}^{(1)}_{\boldsymbol{x}},\text{Var}(\boldsymbol{{\gamma}}^{(1)}_{\boldsymbol{x}}))$,
$\boldsymbol{\gamma}^{(2)}_{\boldsymbol{x}}|\boldsymbol{Y}\sim\mathcal{N}(\boldsymbol{\hat{\gamma}}^{(2)}_{\boldsymbol{x}},\text{Var}(\boldsymbol{{\gamma}}^{(2)}_{\boldsymbol{x}}))$
and
$\boldsymbol{\gamma}^{(3)}_{\boldsymbol{t-x}}|\boldsymbol{Y}\sim\mathcal{N}(\boldsymbol{\hat{\gamma}}^{(3)}_{\boldsymbol{t-x}},\text{Var}(\boldsymbol{{\gamma}}^{(3)}_{\boldsymbol{t-x}}))$,
where

$\boldsymbol{\hat{\gamma}}^{(1)}_{\boldsymbol{x}}={\color[rgb]{0,0,0}\hat{\boldsymbol{K}}^{(1)}}{\boldsymbol{Z}^{(1)}}^{T}\boldsymbol{V}(\boldsymbol{\hat{\theta}}_{(p)})^{-1}(\boldsymbol{Y}-\boldsymbol{T}\hat{\boldsymbol{\beta}}_{(p)}),$
(3.8)
$\text{Var}(\boldsymbol{{\gamma}}^{(1)}_{\boldsymbol{x}})={\color[rgb]{0,0,0}\hat{\boldsymbol{K}}^{(1)}}-[{\color[rgb]{0,0,0}\hat{\boldsymbol{K}}^{(1)}}{\boldsymbol{Z}^{(1)}}^{T}\boldsymbol{V}(\boldsymbol{\hat{\theta}}_{(p)})^{-1}{\boldsymbol{Z}^{(1)}}{\color[rgb]{0,0,0}\hat{\boldsymbol{K}}^{(1)}}],$
$\boldsymbol{{\hat{\gamma}}}^{(2)}_{\boldsymbol{x}}={\color[rgb]{0,0,0}\hat{\boldsymbol{K}}^{(2)}}{\boldsymbol{Z}^{(2)}}^{T}\boldsymbol{V}(\boldsymbol{\hat{\theta}}_{(p)})^{-1}(\boldsymbol{Y}-\boldsymbol{T}\hat{\boldsymbol{\beta}}_{(p)}),$
(3.9)
$\text{Var}(\boldsymbol{{\gamma}}^{(2)}_{\boldsymbol{x}})={\color[rgb]{0,0,0}\hat{\boldsymbol{K}}^{(2)}}-[{\color[rgb]{0,0,0}\hat{\boldsymbol{K}}^{(2)}}{\boldsymbol{Z}^{(2)}}^{T}\boldsymbol{V}(\boldsymbol{\hat{\theta}}_{(p)})^{-1}{\boldsymbol{Z}^{(2)}}{\color[rgb]{0,0,0}\hat{\boldsymbol{K}}^{(2)}}],$
$\boldsymbol{\hat{\gamma}}^{(3)}_{\boldsymbol{t-x}}={\color[rgb]{0,0,0}\hat{\boldsymbol{K}}^{(3)}}{\boldsymbol{Z}^{(3)}}^{T}\boldsymbol{V}(\boldsymbol{\hat{\theta}}_{(p)})^{-1}(\boldsymbol{Y}-\boldsymbol{T}\hat{\boldsymbol{\beta}}_{(p)}),$
(3.10)
$\text{Var}(\boldsymbol{{\gamma}}^{(3)}_{\boldsymbol{t-x}})={\color[rgb]{0,0,0}\hat{\boldsymbol{K}}^{(3)}}-[{\color[rgb]{0,0,0}\hat{\boldsymbol{K}}^{(3)}}{\boldsymbol{Z}^{(3)}}^{T}\boldsymbol{V}(\boldsymbol{\hat{\theta}}_{(p)})^{-1}{\boldsymbol{Z}^{(3)}}{\color[rgb]{0,0,0}\hat{\boldsymbol{K}}^{(3)}}].$

The linear fixed effects estimator $\hat{\boldsymbol{\beta}}$ and its variance
$\text{Var}({\boldsymbol{\beta}})$ in the proposed model is

$\hat{\boldsymbol{\beta}}=(\boldsymbol{T}^{T}\boldsymbol{V}(\boldsymbol{\hat{\theta}}_{(p)})^{-1}\boldsymbol{T})^{-1}\boldsymbol{T}^{T}\boldsymbol{V}(\boldsymbol{\hat{\theta}}_{(p)})^{-1}\boldsymbol{Y},$
$\text{Var}({\boldsymbol{\beta}})=\hat{\sigma}^{2}_{(p)}\boldsymbol{I}(\boldsymbol{T}^{T}\boldsymbol{V}(\boldsymbol{\hat{\theta}}_{(p)})^{-1}\boldsymbol{T})^{-1}.$

We lastly combine all the estimated components to obtain the complete matrix
form of the estimated mean of the proposed model, i.e.

$\hat{\boldsymbol{Y}}=\boldsymbol{T}\hat{\boldsymbol{\beta}}+{\boldsymbol{Z}^{(1)}}\boldsymbol{\hat{\gamma}}^{(1)}_{\boldsymbol{x}}+{\boldsymbol{Z}^{(2)}}\boldsymbol{\hat{\gamma}}^{(2)}_{\boldsymbol{x}}+{\boldsymbol{Z}^{(3)}}\boldsymbol{\hat{\gamma}}^{(3)}_{\boldsymbol{t-x}},$

and the variance of the estimated mean in the proposed model is

$\text{Var}({\boldsymbol{Y}})=\boldsymbol{T}\text{Var}({\boldsymbol{\beta}})\boldsymbol{T}^{T}+\boldsymbol{Z}^{(1)}\text{Var}(\boldsymbol{{\gamma}}^{(1)}_{\boldsymbol{x}}){\boldsymbol{Z}^{(1)}}^{T}+\boldsymbol{Z}^{(2)}\text{Var}(\boldsymbol{{\gamma}}^{(2)}_{\boldsymbol{x}}){\boldsymbol{Z}^{(2)}}^{T}+\boldsymbol{Z}^{(3)}\text{Var}(\boldsymbol{{\gamma}}^{(3)}_{\boldsymbol{t-x}}){\boldsymbol{Z}^{(3)}}^{T}+\hat{\sigma}^{2}_{(p)}\boldsymbol{I}.$

####  3.2.3 Out-of-sample forecasts and the prediction intervals

In the proposed model, we can extrapolate a $h$-step forecast of the logit of
the initial mortality rate from the last observed calendar year $t_{n}$ to the
calendar year $t_{n+h}$ by extending the row of covariate matrices
${\boldsymbol{T}}$, ${\boldsymbol{Z}^{(1)}}$, ${\boldsymbol{Z}^{(2)}}$ from
the size of $n\times m$ to $(n+h)\times m$ with
$[1_{n+1},1_{n+2},\ldots,1_{n+h}]^{T}$, new input
$[t_{n+1},t_{n+2},\ldots,t_{n+h}]^{T}$ and reconstructing
${\boldsymbol{Z}^{(3)}}$ as a new $((n+h)\times m)\times(n+h+m-1)$ cohort
pattern dummy covariate matrix in Equation (3.3) denoted by
${\boldsymbol{T}^{*}}$, ${\boldsymbol{Z}^{*(1)}}$, ${\boldsymbol{Z}^{*(2)}}$,
${\boldsymbol{Z}^{*(3)}}$, such that

$\hat{\boldsymbol{Y}}^{*}=\boldsymbol{T}^{*}\hat{\boldsymbol{\beta}}+{\boldsymbol{Z}^{*(1)}}\boldsymbol{\hat{\gamma}}^{(1)}_{\boldsymbol{x}}+{\boldsymbol{Z}^{*(2)}}\boldsymbol{\hat{\gamma}}^{(2)}_{\boldsymbol{x}}+{\boldsymbol{Z}^{*(3)}}\boldsymbol{\hat{\gamma}}^{(3)}_{\boldsymbol{t+h-x}},$

where

$\boldsymbol{\hat{\gamma}}^{(3)}_{\boldsymbol{t+h-x}}={\color[rgb]{0,0,0}\hat{\boldsymbol{K}}^{*(3)}}{\boldsymbol{Z}^{(3)}}^{T}\boldsymbol{V}(\boldsymbol{\hat{\theta}}_{(p)})^{-1}(\boldsymbol{Y}-\boldsymbol{T}\boldsymbol{\hat{\beta}}),$
(3.11)

and $\hat{\boldsymbol{K}}^{*(3)}$ is the size of $(n+h+m-1)\times(n+m-1)$
covariance matrix for extrapolating
$\boldsymbol{\hat{\gamma}}^{(3)}_{\boldsymbol{t-x}}$ to $h$-step ahead , i.e.

${\color[rgb]{0,0,0}\hat{\boldsymbol{K}}^{*(3)}}={\color[rgb]{0,0,0}\hat{c}_{(p)}^{2}}\text{exp}\Bigg{(}-\begin{bmatrix}\frac{((t_{1}-x_{m})-(t_{1}-x_{m}))^{2}}{2{\color[rgb]{0,0,0}\hat{s}_{(p)}}}&\frac{((t_{1}-x_{m})-(t_{2}-x_{m}))^{2}}{2{\color[rgb]{0,0,0}\hat{s}_{(p)}}}&\ldots&\frac{((t_{1}-x_{m})-(t_{n}-x_{1}))^{2}}{2{\color[rgb]{0,0,0}\hat{s}_{(p)}}}\\\
\frac{((t_{2}-x_{m})-(t_{1}-x_{m}))^{2}}{2{\color[rgb]{0,0,0}\hat{s}_{(p)}}}&\frac{((t_{2}-x_{m})-(t_{2}-x_{m}))^{2}}{2{\color[rgb]{0,0,0}\hat{s}_{(p)}}}&\ldots&\frac{((t_{2}-x_{m})-(t_{n}-x_{1}))^{2}}{2{\color[rgb]{0,0,0}\hat{s}_{(p)}}}\\\
\vdots&\vdots&\ddots&\vdots\\\
\frac{(({\color[rgb]{0,0,0}t_{n+h}}-x_{1})-(t_{1}-x_{m}))^{2}}{2{\color[rgb]{0,0,0}\hat{s}_{(p)}}}&\frac{(({\color[rgb]{0,0,0}t_{n+h}}-x_{1})-(t_{2}-x_{m}))^{2}}{2{\color[rgb]{0,0,0}\hat{s}_{(p)}}}&\ldots&\frac{(({\color[rgb]{0,0,0}t_{n+h}}-x_{1})-(t_{n}-x_{1}))^{2}}{2{\color[rgb]{0,0,0}\hat{s}_{(p)}}}\end{bmatrix}\Bigg{)}.$

The $h$-step out-of-sample forecast variance by the proposed model is

$\text{Var}({{\boldsymbol{Y}}^{*}})=\boldsymbol{T}^{*}\text{Var}({\boldsymbol{\beta}}){\boldsymbol{T}^{*}}^{T}+{\boldsymbol{Z}^{*(1)}}\text{Var}(\boldsymbol{{\gamma}}^{(1)}_{\boldsymbol{x}}){\boldsymbol{Z}^{*(1)}}^{T}+{\boldsymbol{Z}^{*(2)}}\text{Var}(\boldsymbol{{\gamma}}^{(2)}_{\boldsymbol{x}}){\boldsymbol{Z}^{*(2)}}^{T}+{\boldsymbol{Z}^{*(3)}}\text{Var}(\boldsymbol{{\gamma}}^{(3)}_{\boldsymbol{t+h-x}}){\boldsymbol{Z}^{*(3)}}^{T}+\hat{\sigma}^{2}_{(p)}\boldsymbol{I},$

where

$\text{Var}(\boldsymbol{{\gamma}}^{(3)}_{\boldsymbol{t+h-x}})={\color[rgb]{0,0,0}\hat{\boldsymbol{K}}^{**(3)}}-{\color[rgb]{0,0,0}\hat{\boldsymbol{K}}^{**(3)}}{\boldsymbol{Z}^{(3)}}^{T}\boldsymbol{V}(\boldsymbol{\hat{\theta}}_{(p)})^{-1}{\boldsymbol{Z}^{(3)}}{\color[rgb]{0,0,0}\hat{\boldsymbol{K}}^{**(3)}},$

and $\hat{\boldsymbol{K}}^{**(3)}$ is the size of $(n+h+m-1)\times(n+h+m-1)$
covariance matrix, i.e.

${\color[rgb]{0,0,0}\hat{\boldsymbol{K}}^{**(3)}}={\color[rgb]{0,0,0}\hat{c}_{(p)}^{2}}\text{exp}\Bigg{(}-\begin{bmatrix}\frac{((t_{1}-x_{m})-(t_{1}-x_{m}))^{2}}{2{\color[rgb]{0,0,0}\hat{s}_{(p)}}}&\frac{((t_{1}-x_{m})-(t_{2}-x_{m}))^{2}}{2{\color[rgb]{0,0,0}\hat{s}_{(p)}}}&\ldots&\frac{((t_{1}-x_{m})-({\color[rgb]{0,0,0}t_{n+h}}-x_{1}))^{2}}{2{\color[rgb]{0,0,0}\hat{s}_{(p)}}}\\\
\frac{((t_{2}-x_{m})-(t_{1}-x_{m}))^{2}}{2{\color[rgb]{0,0,0}\hat{s}_{(p)}}}&\frac{((t_{2}-x_{m})-(t_{2}-x_{m}))^{2}}{2{\color[rgb]{0,0,0}\hat{s}_{(p)}}}&\ldots&\frac{((t_{2}-x_{m})-({\color[rgb]{0,0,0}t_{n+h}}-x_{1}))^{2}}{2{\color[rgb]{0,0,0}\hat{s}_{(p)}}}\\\
\vdots&\vdots&\ddots&\vdots\\\
\frac{(({\color[rgb]{0,0,0}t_{n+h}}-x_{1})-(t_{1}-x_{m}))^{2}}{2{\color[rgb]{0,0,0}\hat{s}_{(p)}}}&\frac{(({\color[rgb]{0,0,0}t_{n+h}}-x_{1})-(t_{2}-x_{m}))^{2}}{2{\color[rgb]{0,0,0}\hat{s}_{(p)}}}&\ldots&\frac{(({\color[rgb]{0,0,0}t_{n+h}}-x_{1})-({\color[rgb]{0,0,0}t_{n+h}}-x_{1}))^{2}}{2{\color[rgb]{0,0,0}\hat{s}_{(p)}}}\end{bmatrix}\Bigg{)}.$

With the normality assumptions on the model error and the estimated
$\text{Var}({{\boldsymbol{Y}}^{*}})$, a $100(1-\alpha)\%$ prediction interval
for $\hat{\boldsymbol{Y}}^{*}$ can be calculated as
$\hat{\boldsymbol{Y}}^{*}\pm
z_{\alpha}(\text{diag}(\text{Var}({{\boldsymbol{Y}}^{*}})))^{1/2}$, where
$z_{\alpha}$ is the $(1-\alpha/2)$ quantile of the standard normal
distribution.

##  4 Applications

In this section, we demonstrate the usefulness of the proposed model by
modelling and forecasting two sets of sex-specific empirical mortality data.
We apply the proposed model to the sex-specific mortality data of Japan for
illustration purposes first, then compare and evaluate the quality of the
fitted sex-specific mortality curves by the proposed approach with the CBD
model using sex-specific mortality data of ten different developed countries.

###  4.1 Male mortality data

The male mortality data of Japan is available from the year 1947 to the year
2016 (70 years in total) from the Human Mortality Database (2020). The
database consists of the age-specific male central mortality rate by a single
calendar year of age up to 110 years old. Given that the CBD model is designed
in particular for annuities and pensions schemes for older people, we restrict
the numerical examples for the same examined age groups from age 60 to age 89
(30 age groups in total) as the CBD model (Cairns ., 2009).

We convert the central mortality rate to the initial mortality rate by
Equation (2.1) and present the observed logit initial male mortality rates as
univariate time series with 5-year age intervals in Figure 4.1 and the
observed logit initial male mortality curves from age 60 to age 89 from the
year 1947 to the year 2016 in Japan in Figure 4.1. We can see that there is a
general decrease in male mortality rates among all the selected age groups
during the examined period. The changing patterns of the time-series male
mortality rates look reasonably similar among all the selected groups over the
examined period in Figure 4.1. Figure 4.1 presents the same set of data from
the age side as a bunch of the mortality curves in a time-ordering indicated
by the colours of the rainbow from red to violet. Its vertical pattern shows
the general trends and variations of the observed initial male mortality rates
across the examined period. We can see that the historical male mortality
curves move downwards over time in general, due primarily to the advances in
medical technology.

(a)

(b)

Figure 4.1: (a) Univariate time series of the observed logit initial male
mortality rates with 5-year age intervals and (b) the logit initial male
mortality curves from age 60 to age 89 from the year 1947 to the year 2016 in
Japan.

###  4.2 Male mortality modelling and forecasting

As the demonstration, we aim to display 10-years-ahead out-of-sample forecasts
of the male mortality rates. We first split the dataset into a training
dataset with the observed mortality rates from the year 1947 to the year 2006
and a test dataset with the remaining observed mortality data from the year
2007 to the last sample data in the year 2016.

We fit the male mortality data in the training dataset into the proposed model
stated in Equation (3.3). The mean components and random components of the
proposed model can be estimated as discussed in Section 3.2. Thanks to the age
group dependency design in the random components of the proposed model, we can
observe how the mortality rates in different age groups correlate with each
other from the estimated random effects intercept coefficients and the
estimated random effects slope coefficients stated in Equation (3.8) and
Equation (3.9) with the 95% confidence intervals in Figure 2(a) and 2(b). We
can see that they are both increasing along with the increase in age,
indicating that age is a key factor to increase the male mortality rate.
Figure 2(c) presents the estimated random cohort effects coefficients stated
in Equation (3.10). The shape of the estimated random cohort effects
coefficients reflects the improvements or deteriorations of the male mortality
among groups of people born within the same year. It shows that there is a
significant improvement in mortality for those groups of people born in
between the years 1880 to 1900, while a downturn happened for those groups of
people born in around the years 1910 to 1920 in Japan. This may be due to the
fact that this group of people were sent to the main army force in the Second
World War. We then extrapolate the estimated random cohort effects
coefficients for a 10-years-ahead out-of-sample forecast using (3.11) as shown
in Figure 2(c). Figure 4.3 presents the predicted logit initial mate mortality
rates of some selected age groups by the proposed method. We can see that the
proposed model can capture the linear patterns of the logit initial mate
mortality rates accurately with no significant misprediction. Figure 4.4 gives
an example of the 10-years-ahead forecast results of the logit initial mate
mortality curves from age 60 to age 89 (with RMSE = 0.0590) with 95%
prediction intervals by the proposed approach for the year 2016 based on the
observations from 1947 to 2006 in Japan. The root mean square error (RMSE)
here is defined as

$\text{RMSE}=\sqrt{\frac{1}{30}\sum_{j=1}^{30}\bigg{(}y_{x_{j},t_{2016}}-\hat{y}_{x_{j},t_{2016}}\bigg{)}^{2}},$

where $y_{x_{j},t_{2016}}$ is the logit initial male mortality rate aged
$x_{j}$ in the year 2016.

(a)

(b) (c) Figure 4.2: (a) Estimated random effects intercept coefficients, (b)
the estimated random effects slope coefficients and (c) the estimated random
cohort effects coefficients using the observed logit initial male mortality
rates of Japan from the year 1947 to the year 2006. Dashed lines are the 95%
confidence intervals and the vertical line indicates the start point of the
10-years-ahead forecasts of the estimated random cohort effects coefficients.

[1ex]

Figure 4.3: Predicted logit initial male mortality rates of the selected age
groups from ages 60 to 89 with 5-year age intervals using the proposed model
from the year 2007 to the year 2016 based on the observations from the year
1947 to the year 2006 in Japan. The circles are the observed values, the solid
lines are the predictions, and the dashed lines are the 95% prediction
intervals. The vertical line indicates the start point of the predictions.
Figure 4.4: Predicted logit initial male mortality curve (with RMSE = 0.0793)
from age 60 to age 89 with the 95% prediction intervals using the proposed
model for the year 2016 based on the observations from the year 1947 to the
year 2006 in Japan. The circles are the observed logit initial male mortality
rates, the solid line is the prediction and the dashed lines are the 95%
prediction intervals.

###  4.3 Female mortality data

We move to the second case study of modelling and forecasting female mortality
data. Female mortality data of Japan are also available from the year 1947 to
the year 2016 in the Human Mortality Database (2020). We maintain the same
selection for the age groups and also convert the central female mortality
rate to the initial female mortality rate by Equation (2.1) as in the male
mortality case for demonstration in this section.

We present the historical mortality data of Japan as separate univariate time
series of the logit initial female mortality rates with 5-year age intervals
in Figure 4.5 and as the logit initial mortality curves from age 60 to age 89
from the year 1947 to the year 2016 in a time-ordering indicated by the
colours of the rainbow from red to violet in Figure 4.5. We can see that the
downward movement of the time-series female mortality looks fairly similar
among all the selected groups as well as the patterns of the male mortality
over the examined period in Figure 4.5. Figure 4.5 presents the same set of
data from the age side as a bunch of the mortality curves in a time-ordering
indicated by the colours of the rainbow from red to purple. Its vertical
pattern shows the general trends and variations of the observed initial female
mortality rates across the examined period. We can also see that the
historical female mortality curves move downwards over time in general due
mainly to the advances in medical technology and longer longevity as in the
case of male mortality discussed in the previous section.

###  4.4 Female mortality modelling and forecasting

For the task of female mortality modelling and forecasting, we again attempt
to make 10-years-ahead out-of-sample forecasts of the logit initial female
mortality rates for demonstration and maintain the same setting for the split
of dataset as for the male mortality case in this section.

The female mortality data in the training dataset is fitted into the proposed
model stated in Equation (3.3). Figure 6(a) and Figure 6(b) demonstrate the
estimated random effects intercept coefficients and the estimated random
effects slope coefficients with the 95% confidence intervals for the fitted
female mortality data. We can see that they also have increasing patterns
along with the increase in age, reflecting that age is also an essential
element to increase the female mortality rate in Japan. Figure 6(c) presents
the estimated random cohort effects coefficients stated in Equation (3.10).
The shape of the estimated random cohort effects coefficients reflects the
improvements or deteriorations of the female mortality among groups of people
born within the same year. Similar to the pattern of the estimated random
cohort effects coefficients in the male mortality case, there is also a
significant improvement in female mortality for those groups of people born in
between the years 1880 to 1900, while a downturn happened for those groups of
people born in around the years 1910 to 1920 in Japan. We then extrapolate the
estimated random cohort effects coefficients for a 10-years-ahead out-of-
sample forecast using Equation (3.11) as shown in Figure 6(c). Figure 4.7
presents the predicted logit initial female mortality rates of the selected
age groups by the proposed model. We can observe that the proposed model can
catch the linear patterns of the logit initial female mortality rates
reasonably well with no significant mistake. Figure 4.8 provides the
demonstration of the 10-years-ahead forecast results of the logit initial
female mortality curves from age 60 to age 89 (with RMSE = 0.0842) with the
95% prediction intervals by the proposed model for the year 2016 based on the
observations from 1947 to 2006 in Japan.

(a)

(b)

Figure 4.5: (a) Univariate time series of the observed logit initial female
mortality rates with 5-year age intervals and (b) the logit initial female
mortality curves from age 60 to age 89 from the year 1947 to the year 2016 in
Japan. (a)

(b) (c) Figure 4.6: (a) Estimated random effects intercept coefficients, (b)
the estimated random effects slope coefficients and (c) the estimated cohort
effects coefficients using the observed logit initial female mortality rates
of Japan from the year 1947 to the year 2006. Dashed lines are the 95%
confidence intervals and the vertical line indicates the start point of the
10-years-ahead forecasts of the estimated random cohort effects coefficients.

[1ex]

Figure 4.7: Predicted logit initial male mortality rates of the selected age
groups from age 60 to age 89 with 5-year age intervals using the proposed
model from the year 2007 to the year 2016 based on the observations from the
year 1947 to the year 2006 in Japan. The circles are the observed values, the
solid lines are the predictions and the dashed lines are the 95% prediction
intervals. The vertical line indicates the start point of the predictions.
Figure 4.8: Predicted logit initial female mortality curve (with RMSE =
0.0842) from age 60 to age 89 with the 95% prediction intervals using the
proposed model for the year 2016 based on the observations from the year 1947
to the year 2006 in Japan. The circles are the observed logit initial female
mortality rates, the solid line is the prediction and the dashed lines are the
95% prediction intervals.

###  4.5 Forecast accuracy evaluations and comparisons with the CBD model

We now compare and evaluate the forecast performance of the proposed mixed-
effects time series model with the CBD model reviewed in Section 2. In order
to have a comprehensive investigation of the forecast accuracy among these two
models, we consider ten other developed countries for which data are also
available in the Human Mortality Database (2020). We restrict data periods of
all selected countries commencing in the year 1947 up to the year 2016 (70
years in total) for a unified purpose. We examine and quantify the forecasting
performance of the models by the rolling window analysis, which is frequently
used for assessing the consistency of a model’s forecasting ability by rolling
a fixed size prediction interval (window) throughout the observed period
(Zivot Wang, 2007). We hold the sample data from the initial year up to the
year $t_{l}$, where $t_{l}<t_{n}$, as holdout samples and produce the forecast
for the $t_{l+h}$ year where $h$ is the forecast horizon. The forecasts errors
are then determined by comparing the out-of-sample forecast results with the
actual data. We increase one rolling-window (1 year ahead) in year
${\color[rgb]{0,0,0}t_{l+1}}$ to make the same procedure again for the year
$t_{l+1+h}$ until the rolling-window analysis covers all available data in
year $t_{n}$. We consider four different forecast horizons $(h=5,10,15$ and
$20)$ with ten sets of rolling-window to examine the short-term, the mid-term,
and the long-term forecast abilities of the models. In our mortality rolling-
window experiments, the RMSE is defined as

$\text{RMSE}^{(i)}_{c}(h)=\sqrt{\frac{1}{10\times
30}\sum_{w=0}^{9}\sum_{j=1}^{30}\bigg{(}y^{(i)}_{x_{j},t_{l+w+h}}-\hat{y}^{(i)}_{x_{j},t_{l+w+h}}\bigg{)}^{2}},$

where $c$ is the selected country, $w$ is the index of the rolling-window
sets, $i$ is for male $(i=M)$ and for female $(i=F)$, and $j$ is the index of
the examined age groups covering from age 60 to age 89.

Table 4.1 presents the average RMSEs of ten sets of the rolling window
analysis across the ten selected countries in four different forecast horizons
in the male and female mortality curves fitting experiments. The back-testing
results confirm that the forecasting performances and accuracy of the proposed
model are much more desirable than the CBD model555The CBD model was
implemented using R package ‘StMoMo’ in the numerical experiment (Villegas .,
2015).. Among all the ten selected countries and the four different forecast
horizons, the proposed model scored the majority of the lowest RMSEs in both
genders, different countries and forecast horizons in the numerical
experiment.

Country | CBD model |  | Proposed model
---|---|---|---
| M | F | $\frac{\text{M+F}}{2}$ |  | M | F | $\frac{\text{M+F}}{2}$
$h=5$ |  |  |  |  |  |  |
Australia | 0.1174 | 0.1603 | 0.1388 |  | 0.1398 | 0.1103 | 0.1251
Belgium | 0.1039 | 0.1547 | 0.1293 |  | 0.0783 | 0.0557 | 0.0670
Canada | 0.0938 | 0.1607 | 0.1272 |  | 0.0481 | 0.0547 | 0.0514
France | 0.1929 | 0.2614 | 0.2271 |  | 0.0780 | 0.0580 | 0.0680
Italy | 0.1076 | 0.1544 | 0.1310 |  | 0.1109 | 0.0624 | 0.0866
Japan | 0.1229 | 0.1659 | 0.1444 |  | 0.0682 | 0.0759 | 0.0721
Netherlands | 0.1219 | 0.1769 | 0.1494 |  | 0.1050 | 0.1423 | 0.1236
Spain | 0.1106 | 0.1570 | 0.1338 |  | 0.0839 | 0.0744 | 0.0792
UK | 0.1381 | 0.1280 | 0.1330 |  | 0.1237 | 0.0584 | 0.0911
USA | 0.1065 | 0.1729 | 0.1397 |  | 0.0558 | 0.0600 | 0.0579
Average | 0.1215 | 0.1692 | 0.1454 |  | 0.0892 | 0.0752 | 0.0822
$h=10$ |  |  |  |  |  |  |
Australia | 0.1644 | 0.2270 | 0.1957 |  | 0.1660 | 0.1635 | 0.1648
Belgium | 0.1196 | 0.2054 | 0.1625 |  | 0.1153 | 0.1467 | 0.1310
Canada | 0.1337 | 0.2467 | 0.1902 |  | 0.1042 | 0.1174 | 0.1108
France | 0.3247 | 0.4384 | 0.3815 |  | 0.0787 | 0.0846 | 0.0817
Italy | 0.1562 | 0.2079 | 0.1821 |  | 0.1194 | 0.0731 | 0.0963
Japan | 0.1756 | 0.2093 | 0.1925 |  | 0.1023 | 0.1134 | 0.1078
Netherlands | 0.1580 | 0.2572 | 0.2076 |  | 0.1504 | 0.1765 | 0.1635
Spain | 0.1563 | 0.2234 | 0.1899 |  | 0.1191 | 0.1039 | 0.1115
UK | 0.1915 | 0.1812 | 0.1864 |  | 0.1360 | 0.1493 | 0.1426
USA | 0.1735 | 0.2953 | 0.2344 |  | 0.1308 | 0.0742 | 0.1025
Average | 0.1754 | 0.2492 | 0.2123 |  | 0.1222 | 0.1203 | 0.1213
Table 4.1: The average RMSEs in ten sets of the rolling-windows analysis of the predicted mortality curves across the ten selected countries using the CBD model and the proposed model. The minimal forecast errors are highlighted, and the lowest averaged forecast errors among models in different forecast horizons are highlighted in bold. Country | CBD model |  | Proposed model
---|---|---|---
| M | F | $\frac{\text{M+F}}{2}$ |  | M | F | $\frac{\text{M+F}}{2}$
$h=15$ |  |  |  |  |  |  |
Australia | 0.2373 | 0.2679 | 0.2526 |  | 0.2899 | 0.2039 | 0.2469
Belgium | 0.1412 | 0.2126 | 0.1769 |  | 0.1373 | 0.1286 | 0.1330
Canada | 0.1753 | 0.2851 | 0.2302 |  | 0.1373 | 0.0631 | 0.1002
France | 0.3858 | 0.4793 | 0.4326 |  | 0.1672 | 0.1009 | 0.1340
Italy | 0.2212 | 0.2460 | 0.2336 |  | 0.2446 | 0.0899 | 0.1672
Japan | 0.2271 | 0.2268 | 0.2270 |  | 0.1504 | 0.1079 | 0.1292
Netherlands | 0.2454 | 0.2680 | 0.2567 |  | 0.2276 | 0.1969 | 0.2123
Spain | 0.2092 | 0.2848 | 0.2470 |  | 0.1536 | 0.1360 | 0.1448
UK | 0.2698 | 0.2300 | 0.2499 |  | 0.1540 | 0.2130 | 0.1835
USA | 0.2056 | 0.3218 | 0.2637 |  | 0.1153 | 0.0523 | 0.0838
Average | 0.2318 | 0.2822 | 0.2570 |  | 0.1777 | 0.1292 | 0.1535
$h=20$ |  |  |  |  |  |  |
Australia | 0.3387 | 0.2887 | 0.3137 |  | 0.3878 | 0.2783 | 0.3330
Belgium | 0.1835 | 0.2280 | 0.2057 |  | 0.2052 | 0.1462 | 0.1757
Canada | 0.2250 | 0.2680 | 0.2465 |  | 0.3091 | 0.0582 | 0.1836
France | 0.4166 | 0.4128 | 0.4147 |  | 0.1876 | 0.1287 | 0.1582
Italy | 0.3180 | 0.3037 | 0.3109 |  | 0.2875 | 0.1204 | 0.2039
Japan | 0.2945 | 0.2127 | 0.2536 |  | 0.2551 | 0.1258 | 0.1904
Netherlands | 0.3457 | 0.2922 | 0.3190 |  | 0.3111 | 0.1729 | 0.2420
Spain | 0.2628 | 0.3555 | 0.3091 |  | 0.1972 | 0.1787 | 0.1880
UK | 0.3711 | 0.2950 | 0.3331 |  | 0.2866 | 0.2538 | 0.2702
USA | 0.2173 | 0.2906 | 0.2539 |  | 0.1706 | 0.0532 | 0.1119
Average | 0.2973 | 0.2947 | 0.2960 |  | 0.2598 | 0.1516 | 0.2057
Table 4.2:  (continued) The average RMSEs in ten sets of the rolling-windows
analysis of the predicted mortality curves across the ten selected countries
using the CBD model and the proposed model. The minimal forecast errors are
highlighted, and the lowest averaged forecast errors among models in different
forecast horizons are highlighted in bold.

##  5 Discussion and conclusion remarks

In this article, we have introduced a novel mortality model under the time-
series and mixed-effects structure. We have presented the formulation, the
model assumptions and the estimation of the parameters as well as the
differences between the CBD model and the proposed model. We lastly compared
the forecasting abilities among the two models using the rolling window
analysis across the mortality data of ten selected countries. The proposed
mixed effects model has four remarkable advantages over the CBD model as
summarised below.

  1. 1.

Circumventing the age-period-cohort identification problem The proposed mixed
effects model provides an alternative way to avoid the age-period-cohort
identification problem when estimating the model cohort effects parameters.
Unlike the CBD model which estimates the model cohort effects parameters as
the fixed-effects parameters by placing the specified constraints through the
iterative optimisation process, the proposed model under the mixed-effects
structure can produce a unique best fitting solution without the need for
introducing any pre-defined constraints on the model cohort effects
parameters. It is because the model cohort effects parameters are treated as
random variables with their own specific distributions. The structure of the
proposed model also allows us to examine how different age groups correlate
with each other and disclose underlying mortality information about the age
groups correlation patterns in a more compendious fashion than the CBD model.

  2. 2.

Unified time series framework for out-of-sample extrapolation The CBD model
obtains the estimates of the model parameters under a fixed-effects regression
based approach first, then applies the random walk with drift processes method
to do the model coefficients extrapolation in the second step. In contract,
our proposed model can extrapolate the out-of-sample mortality rates under a
unified single step framework. Indeed, De Jong Tickle (2006), Sweeting (2011),
Fung . (2017) and Li O’Hare (2017) doubt about the appropriateness of using
the random walk with drift processes method to extrapolate the coefficients of
the CBD model as a two-step approach. Leng Peng (2016) also argue that the
two-steps approach is a somewhat ad-hoc procedure, and it would be more ideal
to perform the estimation and the extrapolation under a consistent and
rigorous single universal framework from the statistical point of view. The
proposed model in this article aims to provide the unified time series
framework and avoids any potential pitfalls using the two-step estimation
procedure that the CBD model may fall.

  3. 3.

Natural confidence intervals under the normal distribution assumption The CBD
model is constructed based on the Poisson distribution assumption, and it is
difficult to quantify the uncertainties of the estimated model parameters due
to the equal Poisson mean and variance and the over-dispersion problems.
However, the proposed model can provide natural confidence intervals that
include the estimated model parameters uncertainties and the prediction
intervals for the out-of-sample extrapolation given that the proposed model is
established under the Bayesian paradigm with the normal distribution
assumption. The ability to quantifying the parameter risks can enhance the
reliability of the prediction, and this ability is especially practical for
valuing the uncertainties for life insurance portfolios and pension schemes in
the insurance industry.

  4. 4.

Improved forecasting accuracy The forecasting performance of the proposed
model is better than the CBD model. In our window-rolling experiments across
the mortality data of the ten developed countries, the accuracies of the
predicted mortality curves are significantly improved by the proposed model in
short-term, mid-term and long-term forecasting. The improvements may be mainly
due to the time-series and mixed-effects framework of the proposed model,
which can directly extrapolate the data forward without the need of placing
any extra assumptions on the model parameters.

Despite the advantages of the proposed CBD model listed above, both our model
and the CBD model only work when there is a linear pattern in the latter part
of the mortality curve above the pension age. Neither of them are capable of
handling any non-linear mortality patterns for young ages. Although some
attempts have been made to extend the existing CBD model design while keeping
its simplicity to cover the whole age range for mortality modelling, see, for
example, Plat (2009), Cairns . (2009), further research is still needed to
consider about extending the proposed model with more factors which can handle
more comprehensive age ranges or embedding some time weightings approaches in
the model estimation procedure so that it can capture more recent mortality
patterns for better forecasting performance.

## References

  * Cairns . (2006) cairns2006twoCairns, AJ., Blake, D. Dowd, K.  2006\.  A two-factor model for stochastic mortality with parameter uncertainty: Theory and calibration A two-factor model for stochastic mortality with parameter uncertainty: Theory and calibration.  Journal of Risk and Insurance734687–718.
  * Cairns . (2009) cairns2009quantitativeCairns, AJ., Blake, D., Dowd, K., Coughlan, GD., Epstein, D., Ong, A. Balevich, I.  2009\.  A quantitative comparison of stochastic mortality models using data from England and Wales and the United States A quantitative comparison of stochastic mortality models using data from England and Wales and the United States.  North American Actuarial Journal1311–35.
  * Currie (2016) currie2016fittingCurrie, ID.  2016\.  On fitting generalized linear and non-linear models of mortality On fitting generalized linear and non-linear models of mortality.  Scandinavian Actuarial Journal20164356–383.
  * De Jong Tickle (2006) de2006extendingDe Jong, P. Tickle, L.  2006\.  Extending Lee–Carter mortality forecasting Extending Lee–Carter mortality forecasting.  Mathematical Population Studies1311–18.
  * Dickson . (2013) dickson2013actuarialDickson, DC., Hardy, M., Hardy, MR. Waters, HR.  2013\.  Actuarial Mathematics for Life Contingent Risks Actuarial Mathematics for Life Contingent Risks.  Cambridge University Press.
  * Fung . (2017) fung2017unifiedFung, MC., Peters, GW. Shevchenko, PV.  2017\.  A unified approach to mortality modelling using state-space framework: Characterisation, identification, estimation and forecasting A unified approach to mortality modelling using state-space framework: Characterisation, identification, estimation and forecasting.  Annals of Actuarial Science112343–389.
  * Haberman Renshaw (2011) haberman2011comparativeHaberman, S. Renshaw, A.  2011\.  A comparative study of parametric mortality projection models A comparative study of parametric mortality projection models.  Insurance: Mathematics and Economics48135–55.
  * Holford (1983) holford1983estimationHolford, TR.  1983\.  The estimation of age, period and cohort effects for vital rates The estimation of age, period and cohort effects for vital rates.  Biometrics311–324.
  * Human Mortality Database (2020) human2017universityHuman Mortality Database.  2020\.  University of California, Berkeley (USA), and Max Planck Institute for Demographic Research (Germany). Available online: https://www.mortality.org/. Accessed on 20 April 2020 University of California, Berkeley (USA), and Max Planck Institute for Demographic Research (Germany). Available online: https://www.mortality.org/. Accessed on 20 April 2020\.
  * Lee Carter (1992) lee1992modelingLee, RD. Carter, LR.  1992\.  Modeling and forecasting US mortality Modeling and forecasting US mortality.  Journal of the American Statistical Association87419659–671.
  * Leng Peng (2016) leng2016inferenceLeng, X. Peng, L.  2016\.  Inference pitfalls in Lee–Carter model for forecasting mortality Inference pitfalls in Lee–Carter model for forecasting mortality.  Insurance: Mathematics and Economics7058–65.
  * Li O’Hare (2017) li2017semiLi, H. O’Hare, C.  2017\.  Semi-parametric extensions of the Cairns–Blake–Dowd model: A one-dimensional kernel smoothing approach Semi-parametric extensions of the Cairns–Blake–Dowd model: A one-dimensional kernel smoothing approach.  Insurance: Mathematics and Economics77166–176.
  * Li . (2016) li2016twoLi, H., O’Hare, C. Vahid, F.  2016\.  Two-Dimensional Kernel Smoothing of Mortality Surface: An Evaluation of Cohort Strength Two-dimensional kernel smoothing of mortality surface: An evaluation of cohort strength.  Journal of Forecasting356553–563.
  * Murphy (2009) murphy2009goldenMurphy, M.  2009\.  The ‘Golden Generations’ in historical context The ‘Golden Generations’ in historical context.  British Actuarial Journal151–184.
  * Murphy (2010) murphy2010reexaminingMurphy, M.  2010\.  Reexamining the dominance of birth cohort effects on mortality Reexamining the dominance of birth cohort effects on mortality.  Population and Development Review362365–390.
  * O’Brien (2014) o2014estimableO’Brien, RM.  2014\.  Estimable functions in age-period-cohort models: A unified approach Estimable functions in age-period-cohort models: A unified approach.  Quality & Quantity481457–474.
  * O’Brien (2017) o2017mixedO’Brien, RM.  2017\.  Mixed models, linear dependency, and identification in age-period-cohort models Mixed models, linear dependency, and identification in age-period-cohort models.  Statistics in Medicine36162590–2600.
  * Plat (2009) plat2009stochasticPlat, R.  2009\.  On stochastic mortality modeling On stochastic mortality modeling.  Insurance: Mathematics and Economics453393–404.
  * Renshaw Haberman (2003) renshaw2003leeRenshaw, AE. Haberman, S.  2003\.  Lee–Carter mortality forecasting with age-specific enhancement Lee–Carter mortality forecasting with age-specific enhancement.  Insurance: Mathematics and Economics332255–272.
  * Sweeting (2011) sweeting2011trendSweeting, PJ.  2011\.  A trend-change extension of the Cairns-Blake-Dowd model A trend-change extension of the Cairns-Blake-Dowd model.  Annals of Actuarial Science52143–162.
  * Villegas . (2015) villegas2015stmomoVillegas, A., Kaishev, VK. Millossovich, P.  2015\.  StMoMo: An R package for stochastic mortality modelling StMoMo: An R package for stochastic mortality modelling.  7th Australasian Actuarial Education and Research Symposium. 7th Australasian Actuarial Education and Research Symposium.
  * Willets (2004) willets2004cohortWillets, RC.  2004\.  The cohort effect: insights and explanations The cohort effect: insights and explanations.  British Actuarial Journal833–898.
  * Zivot Wang (2007) zivot2007modelingZivot, E. Wang, J.  2007\.  Modeling Financial Time Series with S-Plus® Modeling Financial Time Series with S-Plus® ( 191).  Springer Science & Business Media.