2.1. Definitions and notation

Hereafter, we use the following standard definitions from actuarial literature on mortality modelling (Dickson et al., Reference Dickson, Hardy and Waters2009; Pitacco et al., Reference Pitacco, Denuit, Haberman and Olivieri2009). Let T _x be a random variable representing the remaining lifetime of a person aged x. The cumulative distribution function and survival function of T _x are written as $$_{\tau } q_{x} \, {\equals}\, P\left( {T_{x} \leq \tau } \right)$$ and $$_{\tau } p_{x} \, {\equals}\, P\left( {T_{x} \,\gt\,\tau } \right)$$ , respectively. For a person aged x, the force of mortality at age x+τ is defined as

(1) $$\mu _{{x{\plus}\tau }} \,\colon\! {\equals}\, \mathop {\rm lim}\limits_{h\to0} {1 \over h}P\left( {T_{x} \,\lt\,\tau {\plus}h\!\mid\!T_{x} \,\gt\,\tau } \right)\, {\equals}\, \, {\minus}\, {d \over {d\tau }}{\rm ln}_{\tau } p_{x} $$

Let f _x (t) be the density function of T _x , then from (1) we have $$_{\tau } q_{x} \, {\equals}\, {\int}_0^\tau \,{f_{x} (s)\:ds} \, {\equals}\, {\int}_0^\tau \,{{}_{s}p_{x} \,\mu _{{x{\plus}s}} \:ds} $$ . The central death rate for a x-year-old, where $$x\in{\Bbb N}$$ , is defined as

(2) $$m_{x} \,\colon\! {\equals} \, {{q_{x} } \over {{\int}_0^1 \,{_{s} \:p_{x} ds} }}\, {\equals} \, {{{\int}_0^1 \,{_{s} p_{x} \mu _{{x{\plus}s}} \:ds} } \over {{\int}_0^1 \,{_{s} p_{x} \:ds} }}$$

which is a weighted average of the force of mortality (here $$q_{x} \,\colon\! {\equals} \, {}_{1}q_{x} $$ ). Under the so-called piecewise constant force of mortality assumption, that is $$\mu _{{x{\plus}s}} \, {\equals} \, \mu _{x} $$ , where 0≤s<1 and $$x\in{\Bbb N}$$ , we have, from (2), m _x =μ _x . Moreover, if a Poisson assumption is made for the actual number of deaths, then the resulting maximum likelihood estimate of the force of mortality $$\hat{\mu }_{x} $$ (and hence $$\hat{m}_{x} $$ ) is given by $$\hat{\mu }_{x} \, {\equals} \, {{D_{x} } \mathord{\left/ {\vphantom {{D_{x} } {E_{x} }}} \right. \kern-\nulldelimiterspace} {E_{x} }}\, {\equals} \, \hat{m}_{x} $$ , where D _x is the number of deaths recorded at age x last birthday and the exposure-to-risk E _x the average number of people aged x last birthday, during the observation year. Note that E _x is approximated by an estimate of the population aged x last birthday in the middle of the observation year. We refer to $$\hat{m}_{x} $$ as the crude death rate.

In the above setup it is assumed that the force of mortality μ is deterministic. The stochastic case can be handled by the intensity-based framework where death time is modelled as the first jump time of a doubly stochastic process (Biffis, Reference Biffis2005). Hereafter, we treat the force of mortality μ _x+t (t), the central death rate m _x,t and the crude death rate $$\hat{m}_{{x,t}} $$ as stochastic processes. For a detailed discussion of the background of stochastic mortality modelling in discrete-time and continuous-time, see Cairns et al. (Reference Cairns, Blake and Dowd2008).

2.2. Stochastic mortality models

One of the most widely considered examples of stochastic factor model in the context of mortality modelling is the approach first presented in Lee & Carter (Reference Lee and Carter1992) who proposed a stochastic mortality model for the age-specific crude death rate $$\hat{m}_{{x,t}} $$ , where x=x ₁, … , x _p and t=1, … , T represent age (or age group) and year (time), respectively. Under the model, the dynamics of the log crude death rates, $$y_{{x,t}} \, {\equals} \, {\rm ln}\,\hat{m}_{{x,t}} $$ , is given byFootnote ²

(3) $$y_{{x,t}} \, {\equals} \, \alpha _{x} {\plus}\beta _{x} \kappa _{t} {\plus}{\varepsilon}_{{x,t}} ,\quad {\varepsilon}_{{x,t}} \mathop{\,\sim\,}\limits^{{iid}} N\left( {0,\sigma _{{\varepsilon}}^{2} } \right)$$

where $$N\left( {0,\sigma _{{\varepsilon}}^{2} } \right)$$ denotes a Gaussian distribution with zero mean and variance $$\sigma _{{\varepsilon}}^{2} $$ . The vector $${\mib \bialpha} \, {\equals} \, \alpha _{{x_{1} \! \colon\,\! x_{p} }} \,\colon\! {\equals} \, \left[ {\alpha _{{x_{1} }} ,\,\ldots\!,\alpha _{{x_{p} }} } \right]$$ represents the age-profile of the log death rates and $${\mib \bibeta } \, {\equals} \, \beta _{{x_{1} \!\colon\! \,x_{p} }} $$ measures the sensitivity of death rates for different age group to a change of the time series κ _t . The period effect, κ _t , for forecasting purpose, is assumed to satisfy the equation

(4) $$\kappa _{t} \, {\equals} \, \kappa _{{t\, {\minus}\, 1}} {\plus}\theta {\plus}\omega _{t} ,\quad \omega _{t} \mathop{\,\sim\,}\limits^{{iid}} N(0,\sigma _{\omega }^{2} )$$

where ε _x,t and ω _t are independent.

Under this specification, it is clear that the LC model is not identifiable, since (3) is invarifollowing standard definitions from actuarial literature on mortality modellinagnt up to some linear transformations of the parameters:

(5) $${\mib{y}} _{t} \, {\equals} \, {\mib \bialpha } {\plus}{\mib \bibeta } \kappa _{t} {\plus}{\mib {\bivarepsilon}} _{t} \, {\equals} \, {\mib \bialpha } {\plus}{\mib \bibeta } c{\plus}{{\mib \bibeta } \over d}\left( {\left( {k_{t} \, {\minus}\, c} \right)d} \right){\plus}{\mib {\bivarepsilon}} _{t} \, {\equals} \, {\tilde{\bialpha }} {\plus}{\tilde{\bibeta }} \tilde{\kappa }_{t} {\plus}{\mib {\bivarepsilon}} _{t} $$

where $${\tilde{\bialpha }} \, {\equals} \, { \bialpha } {\plus}{ \bibeta } c$$ , $${\tilde{\bibeta }} \, {\equals} \, {{\bibeta } \mathord{\left/ {\vphantom {{\beta } d}} \right. \kern-\nulldelimiterspace} d}$$ and $$\tilde{\kappa }_{t} \, {\equals} \, \left( {\kappa _{t} \, {\minus}\, c} \right)d$$ .

To overcome this identification issue when estimating the LC model, one has to impose a non-unique choice of constraints to restrict the model to an identifiable class. It is standard practice in actuarial literature to consider the following two constraints:

(6) $$\mathop{\sum}\limits_{x\, {\equals} \, x_{1} }^{x_{p} } \beta _{x} \, {\equals} \, 1,\quad \mathop{\sum}\limits_{t\, {\equals} \, 1}^T \kappa _{t} \, {\equals} \, 0$$

as suggested in Lee & Carter (Reference Lee and Carter1992) to remedy the identifiability issue. This choice of constraints is equivalent to fixing $$c\, {\equals} \, (1\,/\,T)\mathop{\sum}\nolimits_{t\, {\equals} \, 1}^T \kappa _{t} $$ and $$d\, {\equals} \, \mathop{\sum}\nolimits_{x\, {\equals} \, x_{1} }^{x_{p} } \beta _{x} $$ . Consequently we have $$\mathop{\sum}\nolimits_{t\, {\equals} \, 1}^T \,\tilde{\kappa }_{t} \, {\equals} \, 0$$ and $$\mathop{\sum}\nolimits_{x\, {\equals} \, x_{1} }^{x_{p} } \,\tilde{\beta }_{x} \, {\equals} \, 1$$ . The reason for these particular form of identification constraints relates to the fact that the constraint on the path space of the stochastic factor κ ₁, … , κ _T is intended to have the effect of centring the κ _t values over the range t∈{1, … , T}, such that the structure is designed to capture age–period effects with the α _x terms incorporating the main age effects, averaged over time and the bilinear terms β _x κ _t incorporating the age-specific period trends (relative to the main age effects).

Since the introduction of the LC model it has found a widespread uptake of this class of factor model in both practice, where the LC model is now used as a benchmark methodology by the US Bureau of the Census, and in academia where a range of stochastic mortality model extensions have been proposed in the literature, see Table 1.

Table 1 Several popular stochastic mortality models.

We note here that Renshaw & Haberman (Reference Renshaw and Haberman2003) and Renshaw & Haberman (Reference Renshaw and Haberman2006) introduces multi-period $$\left( {\mathop{\sum}\nolimits_{i\, {\equals} \, 1}^k \beta _{x}^{{(i)}} \kappa _{t}^{{(i)}} } \right)$$ and cohort factor $$\left( {\zeta _{{t\, {\minus}\, x}} } \right)$$ , respectively, to the LC method. Currie (Reference Currie2009) considers a simplified version of the model in Renshaw & Haberman (Reference Renshaw and Haberman2006). Cairns et al. (Reference Cairns, Blake and Dowd2006) propose to model $${\rm logit}\left( {q_{{x,t}} } \right)\,\colon\! {\equals} \, {\rm ln}\left( {q_{{x,t}} \,/\,(1\, {\minus}\, q_{{x,t}} )} \right)$$ instead of log death rates and $$\bar{x}$$ is the average age in the sample range. An addition of cohort factor is studied in Cairns et al. (Reference Cairns, Blake, Dowd, Coughlan, Epstein, Ong and Balevich2009). Plat (Reference Plat2009) introduces a model which combines the desirable features of the previous models and include a term $$\left( {\bar{x}\, {\minus}\, x} \right)^{{\plus}} \,\colon\! {\equals} \, {\rm max}\left( {\bar{x}\, {\minus}\, x,0} \right)$$ to capture better young age mortality. The specification of identification constraints for the LC type models, that is for those where the log death rate is being modelled in Table 1, is discussed in Hunt & Villegas (Reference Hunt and Villegas2015).

2.3. Two-stage estimation approaches: frequentist view

Several “classical” approaches to LC model estimation have been proposed in the literature, though they typically involve a two-stage procedure looking first at the observation equation as a regression (ignoring the latent factor structure explicitly) and then in the second stage they fit time series models to the latent factor structures. A good overview of such methods is obtained in Pitacco et al. (Reference Pitacco, Denuit, Haberman and Olivieri2009). This two-stage procedure is at odds with modern state-space modelling procedures which have been progressively moving towards joint parameter estimation and latent state estimation in frequentist and Bayesian formulations, which will be discussed in subsequent sections. This is reflected in the first attempt to improve the calibration approaches as reflected in the comment in Cairns et al. (Reference Cairns, Blake, Dowd, Coughlan and Khalaf-Allah2011) where they highlight that the “… key element of the proposed framework is our single-stage approach to model fitting and process parameter estimation”. Such sentiments, relating to consistent single-stage joint estimation are also echoed in the work of Czado et al. (Reference Czado, Delwarde and Denuit2005).

2.3.1. Multi-factor LC SVD-based two-stage calibration

One of the most commonly adopted approaches to estimate stochastic mortality models is via SVD. We use the multi-period (k-factor) LC model (Renshaw & Haberman, Reference Renshaw and Haberman2003) with identification constraints given by

(7) $$\mathop{\sum}\limits_{t\, {\equals} \, 1}^T \kappa _{t}^{{(i)}} \, {\equals} \, 0,\quad \mathop{\sum}\limits_{x\, {\equals} \, x_{1} }^{x_{p} } \beta _{x}^{{(i)}} \, {\equals} \, 1$$

where i=1, … , k, as an example to illustrate the methodology below (Koissi et al., Reference Koissi, Shapiro and Hognas2006).

Stage 1a – Observation Equation Estimation Stage: We first notice that the constraint $$\mathop{\sum}\nolimits_{t\, {\equals} \, 1}^T \kappa _{t}^{{(i)}} \, {\equals} \, 0$$ will lead to an estimator for the level α given by

(8) $$\hat{\bialpha }_{x} \, {\equals} \, {1 \over T}\mathop{\sum}\limits_{y\, {\equals} \, 1}^T y_{{x,t}} $$

Stage 1b – Observation Equation Estimation Stage: The next stage is to de-trend the observations { y _1:T } by the level estimate $${\hat{\bialpha }} $$ and then to perform a SVD on the resulting (p×T) matrix of residual observations to obtain the decomposition:

(9) $${\rm SVD}\left[ {{\mib{y}} _{{1\: \colon\: T}} \, {\minus}\, { \hat{\bialpha }} } \right]\, {\equals} \, \mathop{\sum}\limits_{i\, {\equals} \, 1}^h \rho _{i}\; {\mib \bimu } _{i} {\mib \binu } _{i}^{ \top } $$

where ⊤ denotes transposition and ρ _i , for i∈{1, … , h}, are the descending singular values where h is the rank of the data matrix. Here μ _i and ν _i are the corresponding left and right singular vectors of the singular value ρ _i with dimension p and T, respectively. For a k-rank, where k≤h, approximation of the matrix, we have

(10) $${\mib{y}} _{{1\: \colon\: T}} \, {\minus}\, {\hat{\alpha }} \, {\equals} \, \mathop{\sum}\limits_{i\, {\equals} \, 1}^k \rho _{i} {\mib \bimu } _{i} {\mib \binu } _{i}^{ \top } {\plus}\bivarphi _{{1\: \colon\: T}} $$

where $${\mib \bivarphi } _{{1\: \colon\: T}} \, {\equals}\, \mathop{\sum}\nolimits_{i\, {\equals} \, k{\plus}1}^h \,\rho _{i} {\bimu } _{i} {\binu } _{i}^{ \top } $$ is the k-rank residuals. We then identify $${ \tilde{\beta }} ^{{(i)}} \, {\equals} \, { \bimu } _{i} $$ and $${ \tilde{\kappa }} ^{{(i)}} \, {\equals} \, \rho _{i} {\mib \binu } _{i} $$ , for i=1, … , k. One then performs the transformation

(11) $$\kappa _{t}^{{(i)}} \, {\equals} \, \tilde{\kappa }_{t}^{{(i)}} \mathop{\sum}\limits_x \tilde{\beta }_{x}^{{(i)}} ,\quad \beta _{x}^{{(i)}} \, {\equals} \, {{\tilde{\beta }_{x}^{{(i)}} } \over {\mathop{\sum}\limits_x \tilde{\beta }_{x}^{{(i)}} }}$$

to ensure the constraints $$\mathop{\sum}\nolimits_x \beta _{x}^{{(i)}} \, {\equals} \, 1$$ , for i=1, … , k, are satisfied.

Stage 2 – Latent Process Factor Estimation Stage: At this stageFootnote ³ , the estimation of the latent factors can be performed by specifying a time series model structure such as ARIMA model for each of the factors:

(12) $$\kappa _{t}^{{(j)}} \, {\equals} \, \theta ^{{(j)}} {\plus}\mathop{\sum}\limits_{r\, {\equals} \, 1}^p \kappa _{{t\,\, {\minus}\, \,r}}^{{(j)}} {\plus}\mathop{\sum}\limits_{s\, {\equals} \, 1}^q {\varepsilon}_{{t\,\, {\minus}\, \,s}}^{{(j)}} {\plus}{\varepsilon}_{t} $$

or alternatively one could fit the equivalent vector autoregressive model structure not treating each factor as independent in the time series specification. One would typically perform this stage of estimation via the Yule–Walker equations, see for instance discussions in Tsay & Tiao (Reference Tsay and Tiao1984). Under such specifications, one then obtain closed form distributions and estimators for period effect latent factor forecasts that can be substituted into the observation model for forecasts of the mortality by age in future forecast horizons and used to construct life tables.

2.3.2. Regression-based approaches

It is important to note that the SVD approach assumes homoscedasticity in the error structure. Therefore, to account for heteroscedasticity in mortality data for different ages, Brouhns et al. (Reference Brouhns, Denuit and Vermunt2002) propose to model death counts, instead of death rates, via Poisson regression where the addition error term in the LC approach is replaced by Poisson random variation. Specifically, the number of death D _x,t is modelled as

(13) $$D_{{x,t}} \:\sim\:{\rm Poisson}\left( {E_{{x,t}} \:m_{{x,t}} (\Phi )} \right)$$

where E _x,t is the death exposure, m _x,t (Φ) is a model of the central death rate and Φ is the parameter vector according to the model being used, including time dynamic factors such as period and cohort effect, see, for example, Table 1. The parameter vector is then estimated by maximising the log-likelihood function, which is given by

(14) $$l(\Phi;\,D,E)\, {\equals} \, \mathop{\sum}\limits_t \mathop{\sum}\limits_x \left( {D_{{x,t}} {\rm ln}\left( {E_{{x,t}} \:m_{{x,t}} (\Phi )} \right)\, {\minus}\, E_{{x,t}} \:m_{{x,t}} (\Phi )\, {\minus}\, {\rm ln}\left( {D_{{x,t}} \,!\,} \right)} \right)$$

where D _x,t ! indicates the factorial of D _x,t . Times series models are then used to model the time dynamic factors forming a second-stage estimation procedure for forecasting purpose. Note that the CBD type models can be estimated under this approach since we have q _x,t =1−exp{−m _x, t} (Cairns et al., Reference Cairns, Blake, Dowd, Coughlan, Epstein, Ong and Balevich2009).

Remark 2.1 In all the discussed cases above, there is the general idea that the two-stage estimation approaches (SVD and regression) treat the unobserved factors corresponding to for instance a period effect κ_t and a cohort effect ζ _t−x as parameters. For forecasting purpose, these dynamics factors are then modelled as time series, typically under the ARIMA framework. In this paper we argue that a more consistent approach involves embedding the specification of the model formally within a state-space model structure and to perform the estimation via a joint combination of filtering and static parameter estimation, which can be achieved either in Bayesian (posterior-based) or frequentist (likelihood-based) settings. We will demonstrate both in this paper.

2.4. Estimation approaches: Bayesian view

From the Bayesian modelling perspective there are few papers that study stochastic mortality models, the main papers in this area involve the works of Czado et al. (Reference Czado, Delwarde and Denuit2005), Pedroza (Reference Pedroza2006), Kogure et al. (Reference Kogure, Kitsukawa and Kurachi2009) and Cairns et al. (Reference Cairns, Blake, Dowd, Coughlan and Khalaf-Allah2011). As observed in these studies, there are many possible advantages to adopting a Bayesian approach for mortality modelling, especially in the context of small populations which may also have substantial quantities of missing data.

An important point to note is that all Bayesian model formulations to date in the mortality modelling literature, that we are aware of, have utilised what would, in modern statistical approaches be considered rudimentary sampling-based approaches to performing Bayesian estimation of the LC type models. The criticism here can be levelled in two ways:

1. 1. The first relates to the fact that in these Bayesian formulations the latent dynamic process states are still treated in the MCMC sampling procedures as if they were a set of static model parameters. The issues with doing this have been mentioned in numerous places, see for example, Carter & Kohn (Reference Carter and Kohn1994). Recently new approaches to such inference in Bayesian models have been developed to avoid having to make univariate conjugate Gibbs or Metropolis-within-Gibbs steps for the latent processes. The reason for this is that it is known in general to be very inefficient in performing inference and can be prone to misleading posterior inference results due to poor mixing performance of the Markov chain for a finite computational budget. Detailed discussions have been provided on such problems in Andrieu et al. (Reference Andrieu, Doucet and Holenstein2010) and subsequently in work such as Chopin et al. (Reference Chopin, Jacob and Papaspiliopoulos2013), and the specific case to population-based state-space models in ecology in Peters et al. (2010 Reference Peters, Hosack and Hayesa ). Pedroza (Reference Pedroza2006) also follows this more efficient sampling approach for mortality modelling based on Carter & Kohn (Reference Carter and Kohn1994).

2. 2. Second, all existing MCMC sampling-based approaches we are aware of for Bayesian inference in the mortality modelling literature tends to neglect the issue of model identification in the likelihood which can cause issues in the Bayesian formulation. In fact, some approaches implement identification constraint in the Bayesian model and develop an MCMC sampler that tries to impose the identification constraint in such a manner that the resulting Markov chain may not be consistent with preserving the correct invariant stationary distribution if one applies the constraints inappropriately. We investigate this issue in a separate paper (Peters et al., Reference Peters, Fung and Shevchenko2017).

These two considerations need to be resolved to update the approaches to more efficient sampling approaches with enhanced specifications of the model formulation to deal with such issues directly. In particular, modern approaches to such model estimations are to treat the latent unobserved process not as static parameters but as a state-space model in which filtering-based methods (Kalman filter variants, SMC) can be utilised for the latent process estimations jointly with consistent estimation of the “static” model parameters. We will detail such estimation procedures which are also consistent with imposing specific identification constraints of relevance to the LC model formulations, that are developed to ensure the correct invariant Bayesian posterior model is preserved by the Markov chain sampler and filters developed.

3.1. LC model with heteroscedasticity (LC-H model)

We present here a state-space formulation of the LC model with heteroscedasticity structure. In this context, the heteroscedasticity refers to a relaxation of the constant single degree of freedom diagonal covariance assumption typically made on the observation vector for each year t across the panel of age group stratifications $${\mib{y}} _{t} \, {\equals} \, \left( {y_{{x_{1} ,t}} ,y_{{x_{2} ,t}} ,\,\ldots\!,y_{{x_{p} ,t}} } \right)$$ . Within this state-space model structure we propose an alternative identification constraint which is tailored for the estimation under the state-space approach.

The LC model with heteroscedasticity structure can be written in state-space form by combining the processes $${\mib{y}} _{t} \, {\equals} \, \left( {y_{{x_{1} ,t}} ,\,\ldots\!,y_{{x_{p} ,t}} } \right)$$ and κ _t into one dynamical system

(17a) $${\mib{y}} _{t} \, {\equals} \, {\mib \bialpha } {\plus}{\mib \bibeta } \kappa _{t} {\plus}{\mib {\bivarepsilon}} _{t} ,\quad {\mib {\bivarepsilon}} _{t} \mathop{\,\sim\,}\limits^{{iid}} N(0,\Sigma )$$

(17b) $$\kappa _{t} \, {\equals} \, \kappa _{{t\, {\minus}\, 1}} {\plus}\theta {\plus}\omega _{t} ,\quad \omega _{t} \mathop{\,\sim\,}\limits^{{iid}} N\left( {0,\sigma _{\omega }^{2} } \right)$$

where $${\mib \bialpha } \, {\equals} \, \alpha _{{x_{1} \!\colon\,x_{p} }} $$ , $${\mib \bibeta } \, {\equals} \, \beta _{{x_{1} \!\colon\,x_{p} }} $$ and κ _t is the latent state of the resulting linear Gaussian state-space model. Here, Σ is a p by p diagonal matrix with $$\sigma _{{{\varepsilon},x_{1} \,\colon\,x_{p} }}^{2} $$ on the diagonal. We refer to this model as LC-H model, and the special case with $$\sigma _{{{\varepsilon},x_{i} }}^{2} \, {\equals} \, \sigma _{{\varepsilon}}^{2} $$ , i∈{1, … , p}, as LC model.

Instead of the identification constraint (6), we suggest an alternative constraint which is simpler and more readily applicable to Monte Carlo-based procedures such as MCMC and SMC. Our formulation of the identification constraint is given by setting

(18) $$\alpha _{{x_{1} }} \, {\equals} \, {\rm constant},\quad \beta _{{x_{1} }} \, {\equals} \, {\rm constant}$$

Such a choice is a valid identification constraint since if one of the elements of each α and β are known (here we have arbitrarily chosen $$\alpha _{{x_{1} }} $$ and $$\beta _{{x_{1} }} $$ ; in general one can choose $$\alpha _{{x_{i} }} $$ and $$\beta _{{x_{j} }} $$ , where i, j∈{1, … , p}), then a non-trivial linear transformation in (5) is not allowed anymore; that is, we must have c=0 and d=1. Note that implementing the proposed constraint is straightforward in both maximum likelihood and Bayesian setting compared to the constraint (6).

3.2. Two-factor LC model with age-based heteroscedasticity (LC2-H model)

A natural extension of the LC-H model is to include a second stochastic factor for the cohort effect. We denote this model by LC2-H model. The cohort effect (Renshaw & Haberman, Reference Renshaw and Haberman2006) can be modelled under the state-space framework as follows:

(19) $$\left[ {\matrix{ {y_{{x_{1} ,t}} } \cr {y_{{x_{2} ,t}} } \cr \vdots \cr {y_{{x_{p} ,t}} } \cr } } \right]\, {\equals} \, \left[ {\matrix{ {\alpha _{{x_{1} }} } \cr {\alpha _{{x_{2} }} } \cr \vdots \cr {\alpha _{{x_{p} }} } \cr } } \right]{\plus}\left[ {\matrix{ {\beta _{{x_{1} }}^{{(1)}} } &#x0026; {\beta _{{x_{1} }}^{{(2)}} } &#x0026; 0 &#x0026; \cdots &#x0026; 0 \cr {\beta _{{x_{2} }}^{{(1)}} } &#x0026; 0 &#x0026; {\beta _{{x_{2} }}^{{(2)}} } &#x0026; \cdots &#x0026; 0 \cr \vdots &#x0026; \vdots &#x0026; \vdots &#x0026; \ddots &#x0026; \vdots \cr {\beta _{{x_{p} }}^{{(1)}} } &#x0026; 0 &#x0026; 0 &#x0026; \cdots &#x0026; {\beta _{{x_{p} }}^{{(2)}} } \cr } } \right]\left[ {\matrix{ {\kappa _{t} } \cr {\zeta _{t}^{{x_{1} }} } \cr {\zeta _{t}^{{x_{2} }} } \cr \vdots \cr {\zeta _{t}^{{x_{p} }} } \cr } } \right]{\plus}\left[ {\matrix{ {{\varepsilon}_{{x_{1} ,t}} } \cr {{\varepsilon}_{{x_{2} ,t}} } \cr \vdots \cr {{\varepsilon}_{{x_{p} ,t}} } \cr } } \right]$$

where $$\zeta _{t}^{x} \,\colon\! {\equals} \, \zeta _{{t\, {\minus}\, x}} $$ . The state equation can be expressed as

(20) $$\left[ {\matrix{ {\kappa _{t} } \cr {\zeta _{t}^{{x_{1} }} } \cr {\zeta _{t}^{{x_{2} }} } \cr \vdots \cr {\zeta _{t}^{{x_{{p\, {\minus}\, 1}} }} } \cr {\zeta _{t}^{{x_{p} }} } \cr } } \right]\, {\equals} \, \left[ {\matrix{ 1 &#x0026; 0 &#x0026; 0 &#x0026; \cdots &#x0026; 0 &#x0026; 0 \cr 0 &#x0026; \vartheta &#x0026; 0 &#x0026; \cdots &#x0026; 0 &#x0026; 0 \cr 0 &#x0026; 1 &#x0026; 0 &#x0026; \cdots &#x0026; 0 &#x0026; 0 \cr 0 &#x0026; 0 &#x0026; 1 &#x0026; \cdots &#x0026; 0 &#x0026; 0 \cr \vdots &#x0026; \vdots &#x0026; \vdots &#x0026; \ddots &#x0026; \vdots &#x0026; \vdots \cr 0 &#x0026; 0 &#x0026; 0 &#x0026; \cdots &#x0026; 1 &#x0026; 0 \cr } } \right]\left[ {\matrix{ {\kappa _{{t\, {\minus}\, 1}} } \cr {\zeta _{{t\, {\minus}\, 1}}^{{x_{1} }} } \cr {\zeta _{{t\, {\minus}\, 1}}^{{x_{2} }} } \cr \vdots \cr {\zeta _{{t\, {\minus}\, 1}}^{{x_{{p\, {\minus}\, 1}} }} } \cr {\zeta _{{t\, {\minus}\, 1}}^{{x_{p} }} } \cr } } \right]{\plus}\left[ {\matrix{ \theta \cr 0 \cr 0 \cr \vdots \cr 0 \cr 0 \cr } } \right]{\plus}\left[ {\matrix{ {\omega _{t}^{\kappa } } \cr {\omega _{t}^{{\zeta _{1} }} } \cr 0 \cr \vdots \cr 0 \cr 0 \cr } } \right]$$

Here, we assume κ _t is a random walk with drift process and an autoregressive model of order 1 (AR(1)) process is assumed for the cohort effect, that is $$\zeta _{t}^{{x_{1} }} \, {\equals} \, \vartheta \zeta _{{t\, {\minus}\, 1}}^{{x_{1} }} {\plus}\omega _{t}^{{\zeta _{1} }} $$ , where |ϑ|<1. Note that, from (20), we have $$\zeta _{t}^{{x_{i} }} \, {\equals} \, \zeta _{{t\, {\minus}\, 1}}^{{x_{{i\, {\minus}\, 1}} }} $$ for i=2, … , p, which is the defining property of the cohort effect and consequently we are only required to model the dynamics of $$\zeta _{t}^{{x_{1} }} $$ . We can write the model (19)–(20) in the following form:

(21a) $${\mib{y}} _{t} \, {\equals} \, {\mib \bialpha } {\plus}B\left[ {\kappa _{t} ,{\mib \bizeta } _{t} } \right]^{ \top } {\plus}{\mib {\bivarepsilon}} _{t} ,\quad {\mib {\bivarepsilon}} _{t} \mathop{\,\sim\,}\limits^{{iid}} N\left( {0,\Sigma } \right)$$

(21b) $$\kappa _{t} \, {\equals} \, \kappa _{{t\, {\minus}\, 1}} {\plus}\theta {\plus}\omega _{t}^{\kappa } ,\quad \omega _{t}^{\kappa } \mathop{\,\sim\,}\limits^{{iid}} N\left( {0,\sigma _{{\omega ^{\kappa } }}^{2} } \right)$$

(21c) $${\mib \bizeta } _{t} \, {\equals} \, C\,{\mib \bizeta } _{{t\, {\minus}\, 1}} {\plus}D{\mib \biomega } _{t}^{\zeta } ,\quad \omega _{t}^{{\zeta _{1} }} \mathop{\,\sim\,}\limits^{{iid}} N\left( {0,\sigma _{{\omega ^{\zeta } }}^{2} } \right)$$

where B is the p by p+1 matrix in (19), C the corresponding p by p sub-matrix in (20) and D a zero p by p matrix except for the (1, 1) element with value 1.

For further details on the formulation, estimation and forecasting of cohort models in state-space framework, see Fung et al. (Reference Fung, Peters and Shevchenko2017).

3.3. Three-factor LC model with dynamic time-based heteroscedasticity (LC3-H2 model)

We can further extend the LC2-H model by adding a third factor for the dynamic of volatility in the observation vector over time. Specifically, a time dependent stochastic factor $$\sqrt {\gamma _{t}^{y} } $$ , which is independent to age (or age group), is incorporated into the observation noise term as follows:

(22a) $${\mib{y}} _{t} \, {\equals} \, {\mib \bialpha } {\plus}B\left[ {\kappa _{t} ,{\mib \bizeta } _{t} } \right]^{ \top } {\plus}\sqrt {\gamma _{t}^{y} } {\bf {\bivarepsilon}}_{t} ,\quad {\bf {\bivarepsilon}}_{t} \mathop{\,\sim\,}\limits^{{iid}} N\left( {0,\Sigma } \right)$$

(22b) $$\kappa _{t} \, {\equals} \, \kappa _{{t\, {\minus}\, 1}} {\plus}\theta {\plus}\omega _{t}^{\kappa } ,\quad \omega _{t}^{\kappa } \mathop{\,\sim\,}\limits^{{iid}} N\left( {0,\sigma _{{\omega ^{\kappa } }}^{2} } \right)$$

(22c) $${\mib \bizeta } _{t} \, {\equals} \, C{\mib \bizeta } _{{t\, {\minus}\, 1}} {\plus}D{\mib \biomega } _{t}^{\zeta } ,\quad \omega _{t}^{{\zeta _{1} }} \mathop{\,\sim\,}\limits^{{iid}} N\left( {0,\sigma _{{\omega ^{\zeta } }}^{2} } \right)$$

(22d) $$\gamma _{t}^{y} \, {\equals} \, a\left( {b\, {\minus}\, \gamma _{{t\, {\minus}\, 1}}^{y} } \right){\plus}\gamma _{{t\, {\minus}\, 1}}^{y} {\plus}\sigma \sqrt {\gamma _{{t\, {\minus}\, 1}}^{y} } {\varepsilon}_{t}^{{\gamma ^{y} }} ,\quad {\varepsilon}_{t}^{{\gamma ^{y} }} \mathop{\,\sim\,}\limits^{{iid}} N\left( {0,1} \right)$$

where $$\gamma _{t}^{y} $$ is a process obtained via an Euler discretisation of a square Bessel process corresponding to the Cox–Ingersoll–Ross process given by

$$d\gamma _{t}^{y} \, {\equals} \, a\left( {b\, {\minus}\, \gamma _{t}^{y} } \right)\:dt{\plus}\sigma \sqrt {\gamma _{t}^{y} } \:dW_{t} $$

where 2ab≥σ ² to ensure $$\gamma _{t}^{y} $$ is strictly positive. Such a dynamic volatility factor can be used to explain time-varying periods of heightened observation variance, which potentially occur in some populations over time. These may be attributed to disease, war, famine, environmental factors or shocks as well as changes in migration and immigration patterns that could influence the volatility of the observed death counts in different age groups. Note that this is one possible way to introduce stochastic volatility to mortality models where the stochastic volatility factor enters via the observation noise term. Another possible choice is to consider stochastic volatility through the state equation described in section 3.4 which will be the focus of this paper.

3.4. Multi-factor model with stochastic volatility in the latent process (LCSV and LCSV-H models)

A common assumption in mortality modelling is that the period effect is derived from a discretisation of a random walk with drift process. Such a process may be sufficient for modelling simple dynamics, but can be insufficient if time-varying periods of volatility are present in the time series. In fact, much of the literature focusses mainly on capturing the trend of the period effect κ _t for the past several decades where mortality time series for many countries are reasonably smooth.

Here, we extend the LC framework to incorporate stochastic volatility in the latent process. As a result, the impact of epidemics, natural disasters, medical breakthrough or wars on the evolution of mortality can be taken into account. This will produce different structural effects on the calibration and importantly on the forecasting when compared to the previously developed model of LC3-H2. We refer (23a)–(23c) as LC stochastic volatility model which we denote as LCSV model:

(23a) $${\mib{y}} _{t} \, {\equals} \, {\mib \bialpha } {\plus}{\mib \bibeta } \kappa _{t} {\plus}{\mib {\bivarepsilon}} _{t} ,\quad {\mib {\bivarepsilon}} _{t} \mathop{\,\sim\,}\limits^{{iid}} N\left( {{\bf{0}} ,\sigma _{{\varepsilon}}^{2} {\bf{1}} _{p} } \right)$$

(23b) $$\kappa _{t} \, {\equals} \, \kappa _{{t\, {\minus}\, 1}} {\plus}\theta {\plus}\omega _{t} ,\quad \omega _{t} \!\mid\!\gamma _{t} \!\sim\1N\left( {0,exp\left\{ {\gamma _{t} } \right\}} \right)$$

(23c) $$\gamma _{t} \, {\equals} \, \lambda _{1} \gamma _{{t\, {\minus}\, 1}} {\plus}\lambda _{2} {\plus}\eta _{t} ,\quad \eta _{t} \mathop{\,\sim\,}\limits^{{iid}} N\left( {0,\sigma _{\gamma }^{2} } \right)$$

The log-volatility process γ _t is introduced in the state equation for κ _t via the error term ω _t . The process γ _t is an AR(1) with |λ ₁|<1 and the mean reverting level is given by λ ₂/(1−λ ₁). A heteroscedasticity structure can be introduced in (23a) and will be referred to as LCSV-H model. Specifically, the LCSV-H model is given by

(24a) $${\mib{y}} _{t} \, {\equals} \, {\mib \bialpha } {\plus}{\mib \bibeta } \kappa _{t} {\plus}{\mib {\bivarepsilon}} _{t} ,\quad {\mib {\bivarepsilon}} _{t} \mathop{\,\sim\,}\limits^{{iid}} N\left( {{\bf{0}} ,\Sigma } \right)$$

(24b) $$\kappa _{t} \, {\equals} \, \kappa _{{t\, {\minus}\, 1}} {\plus}\theta {\plus}\omega _{t} ,\quad \omega _{t} \!\mid\!\gamma _{t} \: \sim\: N\left( {0,{\rm exp}\left\{ {\gamma _{t} } \right\}} \right)$$

(24c) $$\gamma _{t} \, {\equals} \, \lambda _{1} \gamma _{{t\, {\minus}\, 1}} {\plus}\lambda _{2} {\plus}\eta _{t} ,\quad \eta _{t} \mathop{\,\sim\,}\limits^{{iid}} N\left( {0,\sigma _{\gamma }^{2} } \right)$$

Cohort effect can also be incorporated into the LCSV model as follows:

(25a) $${\mib y}_{t} \, {\equals} \, {\mib \bialpha } {\plus}B\left[ {\kappa _{t} ,{\mib \bizeta } _{t} } \right]^{ \top } {\plus}{\mib {\bivarepsilon}} _{t} ,\quad {\mib {\bivarepsilon}} _{t} \mathop{\,\sim\,}\limits^{{iid}} N\left( {{\bf{0}} ,\sigma _{{\varepsilon}}^{2} {\bf{1}} _{p} } \right)$$

(25b) $$\kappa _{t} \, {\equals} \, \kappa _{{t\, {\minus}\, 1}} {\plus}\theta {\plus}\omega _{t} ,\quad \omega _{t} \!\mid\!\gamma _{t} \:\sim\: N\left( {0,exp\left\{ {\gamma _{t} } \right\}} \right)$$

(25c) $$\gamma _{t} \, {\equals} \, \lambda _{1} \gamma _{{t\, {\minus}\, 1}} {\plus}\lambda _{2} {\plus}\eta _{t} ,\quad \eta _{t} \mathop{\,\sim\,}\limits^{{iid}} N\left( {0,\sigma _{\gamma }^{2} } \right)$$

(25d) $${\mib \bizeta } _{t} \, {\equals} \, C{\mib \bizeta } _{{t\, {\minus}\, 1}} {\plus}D{\mib \biomega } _{t}^{\zeta } ,\quad {\mib \omega } _{t}^{{\zeta _{1} }} \mathop{\,\sim\,}\limits^{{iid}} N\left( {0,\sigma _{{\omega ^{\zeta } }}^{2} } \right)$$

Compared to the LC3-H2 model where stochastic volatility is included in the observation noise term, introducing stochastic volatility in the latent period process has the advantage, in terms of simplicity and ease of interpretation, that the variability of mortality data in the time dimension is captured purely by the latent process.

4.1. Closed form filter recursions for LC-H model

The aim of filtering is to obtain the distribution of the latest state given observations. For a general state-space model, (15)–(16), the filtering density π( ϕ _t | z _1:t ) at time t can be calculated sequentially by first assuming the filtering density π( ϕ _t−1| z _1:t−1) at time t−1 is known. Then the one-step ahead predictive density for the state is given by

(32) $$\pi \left( {{\mib \biphi } _{t} \!\mid\!{\mib{z}} _{{1\: \colon\: t\, {\minus}\, 1}} } \right)\, {\equals} \, {\int} {\pi \left( {{\mib \biphi } _{t} \!\mid\!{\mib \biphi } _{{t\, {\minus}\, 1}} } \right)\pi \left( {{\mib \biphi } _{{t\, {\minus}\, 1}} \!\mid\!{\mib{z}} _{{1\: \colon\: t\, {\minus}\, 1}} } \right)d{\mib \biphi } _{{t\, {\minus}\, 1}} } $$

From Bayes’ formula and the structure of the conditional dependency of the state-space model, one can obtain the filtering density as

(33) $$\pi \left( {{\mib \biphi } _{t} \!\mid\!{\mib{z}} _{{1\: \colon\: t}} } \right)\, {\equals} \, {{\pi \left( {{\mib \biphi } _{t} \!\mid\!{\mib{z}} _{{1\: \colon\: t\, {\minus}\, 1}} } \right)\pi \left( {{\mib{z}} _{t} \!\mid\!{\mib \biphi } _{t} } \right)} \over {{\int} {\pi \left( {{\mib \biphi } _{t} \!\mid\!{\mib{z}} _{{1\: \colon\: t\, {\minus}\, 1}} } \right)\pi \left( {{\mib{z}} _{t} \!\mid\!{\mib \biphi } _{t} } \right)d{\mib \biphi } _{t} } }}$$

For non-linear and non-Gaussian state-space models, numerical techniques such as SMC methods are required to estimate the filtering density, see Doucet et al. (Reference Doucet, De Freitas and Gordon2001), Doucet et al. (Reference Doucet, Godsill and Andrieu2000) and Liu (Reference Liu2008).

In the case of the LC-H model, since it is a linear and Gaussian state-space model, the filtering distribution can be obtained analytically via Kalman filtering. In particular, we can find the conditional distributions of the key quantities in the filtering recursions are all Gaussian distributions as follows:

(34a) $$\kappa _{{t\, {\minus}\, 1}} \!\mid\!{\mib{y}} _{{1\: \colon\: t\, {\minus}\, 1}} \: \sim\: N\left( {m_{{t\, {\minus}\, 1}} ,C_{{t\, {\minus}\, 1}} } \right)$$

(34b) $$\kappa _{t} \!\mid\!{\mib{y}} _{{1\: \colon\: t\, {\minus}\, 1}} \: \sim\: N\left( {a_{t} ,R_{t} } \right)$$

(34c) $${\mib{y}} _{t} \!\mid\!{\mib{y}} _{{1\: \colon\: t\, {\minus}\, 1}} \: \sim\: N\left( {{\mib{f}} _{t} ,{\mib{Q}} _{t} } \right)$$

(34d) $$\kappa _{t} \!\mid\!{\mib{y}} _{{1\: \colon\: t}} \: \sim\: N\left( {m_{t} ,C_{t} } \right)$$

where the recursive nature of these distributions arises from the recursions of the sufficient statistics:

(35a) $$a_{t} \, {\equals} \, m_{{t\, {\minus}\, 1}} {\plus}\theta ,\quad R_{t} \, {\equals} \, C_{{t\, {\minus}\, 1}} {\plus}\sigma _{\omega }^{2} $$

(35b) $${\mib{f}} _{t} \, {\equals} \, {\mib \bialpha } {\plus}{\mib \bibeta } a_{t} ,\quad {\mib{Q}} _{t} \, {\equals} \, {\mib \bibeta \bibeta } ^{{\rm \top }} R_{t} {\plus}\Sigma $$

(35c) $$m_{t} \, {\equals} \, a_{t} {\plus}R_{t} {\mib \bibeta } ^{{\rm \top }} {\mib{Q}} _{t}^{{\, {\minus}\, 1}} \left( {{\mib{y}} _{t} \, {\minus}\, {\mib{f}} _{t} } \right),\quad C_{t} \, {\equals} \, R_{t} \, {\minus}\, R_{t} {\mib \bibeta } ^{{\rm \top }} {\mib{Q}} _{t}^{{\, {\minus}\! 1}} {\mib \bibeta } R_{t} $$

That is, given the filtering distribution at t−1, (34a), the filtering distribution at t is given by (34d) using (35a)–(35c).

4.2. Closed form gradient-based estimation via score and Hessian recursions for LC-H model

For the LC-H model, the log-likelihood function $$\ell ({\mib \bipsi } )\,\colon\! {\equals} \, {\rm ln}\,\pi \left( {{\mib{y}} _{{1\: \colon\: T}} \!\mid\!{\mib \bipsi } } \right)$$ is given by

(36) $$\ell ({\mib \bipsi } )\, {\equals} \, {\rm ln}\,\left( {\prod\limits_{t\, {\equals} \, 1}^T \pi ({\mib{y}} _{t} \!\mid\!{\mib{y}} _{{1\: \colon\: t\, {\minus}\, 1}} ,{\mib \bipsi } )} \right)\, {\equals} \, \, {\minus}\, {{pT} \over 2}{\rm ln}\,2\pi \, {\minus}\, {1 \over 2}\mathop{\sum}\limits_{t\, {\equals} \, 1}^T \left( {{\rm ln}\,\!\mid\!{\mib{Q}} _{t} \!\mid\!{\plus}{\mib \binu } _{t}^{{\rm \top }} {\bf Q}_{t}^{{\, {\minus}\! 1}} {\mib \binu } _{t} } \right)$$

where ν _t := y _t − f _t , and $${\mib \bipsi } \, {\equals} \, \big( {\alpha _{{x_{2} \: \colon\: x_{p} }} ,\beta _{{x_{2} \: \colon\: x_{p} }} ,\theta ,\sigma _{{{\varepsilon},x_{1} \: \colon\: x_{p} }}^{2} ,\sigma _{\omega }^{2} } \big)$$ is an n-dimensional parameter vector. The log-likelihood function (54) can be derived directly from (34c).

It can be shown that (Harvey, Reference Harvey1989) the elements of the score vector and the information matrix are given in closed form for the LC-H model according to the expressions:

(37) $${{\partial \ell } \over {\partial \psi _{i} }}\, {\equals} \, {1 \over 2}\mathop{\sum}\limits_{t\, {\equals} \, 1}^T \left\{ {{\rm tr}\left[ {\left( {{\mib{Q}} _{t}^{{\, {\minus}\! 1}} {{\partial {\mib{Q}} _{t} } \over {\partial \psi _{i} }}} \right)({\bf{1}} _{p} \, {\minus}\, {\mib{Q}} _{t}^{{\, {\minus}\! 1}} {\mib \binu } _{t} {\mib \binu } _{t}^{{\rm \top }} )} \right]{\plus}2{{\partial {\mib \binu } _{t}^{{\rm \top }} } \over {\partial \psi _{i} }}{\mib{Q}} _{t}^{{\, {\minus} 1}} {\mib \binu } _{t} } \right\},\quad i\, {\equals} \, 1,\,\ldots\:,n$$

where tr[·] denotes the trace operator and

(38) $$\, {\minus}\, {\rm E}\left[ {{{\partial ^{2} \ell } \over {\partial \psi _{i} \partial \psi _{j} }}} \right]\, {\equals} \, {1 \over 2}\mathop{\sum}\limits_{t\, {\equals} \, 1}^T \left[ {{\rm tr}\left( {{\mib{Q}} _{t}^{{\, {\minus}\! 1}} {{\partial {\mib{Q}} _{t} } \over {\partial \psi _{i} }}{\mib{Q}} _{t}^{{\, {\minus}\! 1}} {{\partial {\mib{Q}} _{t} } \over {\partial \psi _{j} }}} \right)} \right]{\plus}{\rm E}\left[ {\mathop{\sum}\limits_{t\, {\equals} \, 1}^T {{\partial {\mib \binu } _{t}^{{\rm \top }} } \over {\partial \psi _{i} }}{\mib{Q}} _{t}^{{\, {\minus}\! 1}} {{\partial {\mib \binu } _{t} } \over {\partial \psi _{j} }}} \right],\quad i,j\, {\equals} \, 1,\,\ldots\:,n$$

and the expectation operator E[·] on the second term in (38) can be dropped (since the expressions are asymptotically equivalent). In order to evaluate the score vector and the information matrix, we need

(39) $${{\partial {\mib \binu } _{t} } \over {\partial \psi _{i} }}\, {\equals} \, \, {\minus}\, {{\partial {\mib \bialpha } } \over {\partial \psi _{i} }}\, {\minus}\, {{\partial {\mib \bibeta } } \over {\partial \psi _{i} }}a_{t} \, {\minus}\, {\mib \bibeta } {{\partial a_{t} } \over {\partial \psi _{i} }}$$

and

(40) $${{\partial {\mib{Q}} _{t} } \over {\partial \psi _{i} }}\, {\equals} \, {{\partial {\mib \bibeta } } \over {\psi _{i} }}R_{t} {\mib \bibeta } ^{{\rm \top }} {\plus}{\mib \bibeta } {{\partial R_{t} } \over {\partial \psi _{i} }}{\mib \bibeta } ^{{\rm \top }} {\plus}{\mib \bibeta } R_{t} {{\partial {\mib \bibeta } ^{{\rm \top }} } \over {\partial \psi _{i} }}{\plus}{{\partial \Sigma } \over {\partial \psi _{i} }}$$

The expressions (39) and (40) require, for t=1, … , T and i=1, … , n

(41) $${{\partial a_{t} } \over {\partial \psi _{i} }}\, {\equals} \, {{\partial m_{{t\, {\minus}\, 1}} } \over {\partial \psi _{i} }}{\plus}{{\partial \theta } \over {\partial \psi _{i} }}\quad {\rm and}\quad {{\partial R_{t} } \over {\partial \psi _{i} }}\, {\equals} \, {{\partial C_{{t\, {\minus}\, 1}} } \over {\partial \psi _{i} }}{\plus}{{\partial \sigma _{\omega }^{2} } \over {\partial \psi _{i} }}$$

The expressions in (41) in turn require, for t=1, … , T−1 and i=1, … , n

(42) $${{\partial m_{t} } \over {\partial \psi _{i} }}\, {\equals} \, {{\partial a_{t} } \over {\partial \psi _{i} }}{\plus}{{\partial R_{t} } \over {\partial \psi _{i} }}{\mib \bibeta } ^{{\rm \top }} {\mib{Q}} _{t}^{{\, {\minus}\! 1}} {\mib \binu } _{t} {\plus}R_{t} {{\partial {\mib \bibeta } } \over {\partial \psi _{i} }}{\mib{Q}} _{t}^{{\, {\minus}\! 1}} {\mib \binu } _{t} \, {\minus}\, R_{t} {\mib \bibeta } ^{{\rm \top }} {\mib{Q}} _{t}^{{\, {\minus}\! 1}} {{\partial {\mib Q } _{t} } \over {\partial \psi _{i} }}{\mib{Q}} _{t}^{{\, {\minus}\! 1}} {\mib \binu } _{t} {\plus}R_{t} {\mib \bibeta } ^{{\rm \top }} {\mib{Q}} _{t}^{{\, {\minus}\! 1}} {{\partial {\mib \binu } _{t} } \over {\partial \psi _{i} }}$$

and

(43) $$\eqalignno{ {{\partial C_{t} } \over {\partial \psi _{i} }}\, {\equals} \, &#x0026; {{\partial R_{t} } \over {\partial \psi _{i} }}\, {\minus}\, {{\partial R_{t} } \over {\partial \psi _{i} }}{\mib \bibeta } ^{{\rm \top }} {\mib{Q}} _{t}^{{\, {\minus}\! 1}} {\mib \bibeta } R_{t} \, {\minus}\, R_{t} {{\partial {\mib \bibeta } } \over {\partial \psi _{i} }}{\mib{Q}} _{t}^{{\, {\minus}\! 1}} {\mib \bibeta } R_{t} \cr &#x0026; {\plus}R_{t} {\mib \bibeta } ^{{\rm \top }} {\mib{Q}} _{t}^{{\, {\minus}\! 1}} {{\partial {\mib{Q}} _{t} } \over {\partial \psi _{i} }}{\mib{Q}} _{t}^{{\, {\minus}\! 1}} {\mib \bibeta } R_{t} \, {\minus}\, R_{t} {\mib \bibeta } ^{{\rm \top }} {\mib{Q}} _{t}^{{\, {\minus}\! 1}} {{\partial {\mib \bibeta } } \over {\partial \psi _{i} }}R_{t} \, {\minus}\, R_{t} {\mib \bibeta } ^{{\rm \top }} {\mib{Q}} _{t}^{{\, {\minus}\! 1}} {\mib \bibeta } {{\partial R_{t} } \over {\partial \psi _{i} }} $$

Note that $${{\partial m_{0} } \over {\partial \psi _{i} }}\, {\equals} \, {{\partial C_{0} } \over {\partial \psi _{i} }}\, {\equals} \, 0$$ for i=1, … , n and the required differentiation matrices $${{\partial {\mib \bialpha } } \over {\partial {\mib \bipsi } _{i} }}$$ , $${{\partial {\mib \bibeta } } \over {\partial {\bf \psi }_{i} }}$$ , $${{\partial {\bf \biSigma } } \over {\partial {\bf \psi }_{i} }}$$ , $${{\partial {\mib \bitheta } } \over {\partial {\mib \bipsi } _{i} }}$$ and $${{\partial {\mib \bisigma } _{{\mib \biomega} }^{{\bf 2}} } \over {\partial {\mib \bipsi } _{i} }}$$ are displayed in Appendix A. The gradient-based estimation for the LC-H model is described in Algorithm 1.

5.1. Linear Gaussian state-space inference

We develop an efficient approach involving a combined Gibbs sampling conjugate model sampler for the marginal target distributions of the static model parameters along with a forward-backward Kalman filter sampler for the latent process κ _1:T . A sample of the targeted density is obtained via Gibbs sampling where M is the number of MCMC iterations (Algorithm 2).

The general block Gibbs sampling algorithm steps require to sample from the full conditional densities π(κ _0:T | $$\bipsi $$ , y _1:T ) and π( $$\bipsi $$ |κ _0:T , y _1:T ), which are shown in the following.

5.1.1. Sampling from the full conditional density π(κ _0:T | $$\bipsi $$ , y _1:T )

Samples of the full conditional density π(κ _0:T | $$\bipsi $$ , y _1:T ) can be obtained via the so-called forward-filtering-backward sampling (FFBS) procedure (Carter & Kohn, Reference Carter and Kohn1994). This sampling approach is also utilised in Pedroza (Reference Pedroza2006). We can write

(45) $$\pi \left( {\kappa _{{0\: \colon\: T}} \!\mid\!{\mib \bipsi } ,{\mib{y}} _{{1\: \colon\: T}} } \right)\, {\equals} \, \prod\limits_{t\, {\equals} \, 0}^T \pi \left( {\kappa _{t} \!\mid\!\kappa _{{t{\plus}1\: \colon\: T}} ,{\mib \bipsi } ,{\mib{y}} _{{1\: \colon\: T}} } \right)\, {\equals} \, \prod\limits_{t\, {\equals} \, 0}^T \pi \left( {\kappa _{t} \!\mid\!\kappa _{{t{\plus}1}} ,{\mib \bipsi } ,{\mib{y}} _{{1\: \colon\: t}} } \right)$$

where the last term in the product, π(κ _T | $$\bipsi $$ , y _1:T ), is distributed as N(m _T ,C _T ) which is obtained from the last iteration of the Kalman filtering procedure.

Once we draw a sample κ _T from N(m _T ,C _T ), then (45) suggests that we can draw recursively and backwardly κ _t from π(κ _t |κ _t+1, $$\bipsi $$ , y _1:t ), where t=T−1, T−2, … , 1, 0. Moreover, we have

(46) $$\kappa _{t} \!\mid\!\kappa _{{t{\plus}1}} ,{\mib \bipsi } ,{\mib{y}} _{{1\: \colon\: t}} {\rm \, \sim\, }N\left( {h_{t} ,H_{t} } \right)$$

where

(47a) $$h_{t} \, {\equals} \, m_{t} {\plus}C_{t} R_{{t{\plus}1}}^{{\, {\minus}\! 1}} \left( {\kappa _{{t{\plus}1}} \, {\minus}\, a_{{t{\plus}1}} } \right)$$

(47b) $$H_{t} \, {\equals} \, C_{t} \, {\minus}\, C_{t} R_{{t{\plus}1}}^{{\, {\minus}\! 1}} C_{t} $$

which can be derived based on Kalman smoother (Petris et al., Reference Petris, Petrone and Campagnoli2009).

The FFBS procedure is displayed in Algorithm 3. Note that the prior distribution for κ ₀ can be set to be vague to run the Kalman filter; the output of the algorithm includes the posterior distribution of κ ₀.

5.1.2. Sampling from the full conditional density π( $$\bipsi $$ |κ _0:T , y _1:T )

The first thing to observe is that under the reparameterisation of the identification constraints (18), the following Gibbs sampling stages can be performed exactly.

We assume that the prior for $$\big( {\alpha _{{x_{2} \: \colon\: x_{p} }} ,\beta _{{x_{2} \: \colon\: x_{p} }} ,\theta ,\sigma _{{{\varepsilon},x_{1} \: \colon\: x_{p} }}^{2} ,\sigma _{\omega }^{2} } \big)$$ are given by

(48a) $$\alpha _{x} \: \sim\: N\left( {\tilde{\bimu }_{\alpha } ,\tilde{\sigma }_{\alpha }^{2} } \right),\quad \beta _{x} \: \sim\: N\big( {\tilde{\bimu }_{\beta } ,\tilde{\sigma }_{\beta }^{2} } \big),\quad \theta \: \sim\: N\left( {\tilde{\bimu }_{\theta } ,\tilde{\sigma }_{\theta }^{2} } \right)$$

(48b) $$\sigma _{{{\varepsilon},x}}^{2} \: \sim\: {\rm IG}\big( {\tilde{a}_{{\varepsilon}} ,\tilde{b}_{{\varepsilon}} } \big),\quad \sigma _{\omega }^{2} \: \sim\: {\rm IG}\big( {\tilde{a}_{\omega } ,\tilde{b}_{\omega } } \big)$$

where $${\rm IG}\big( {\tilde{a}_{\omega } ,\tilde{b}_{\omega } } \big)$$ denotes an inverse-γ distribution with mean $${{\tilde{b}_{\omega } } /{\left( {\tilde{a}_{\omega } {\minus}\, 1} \right)^{2}}$$ and variance $${{\tilde{b}_{\omega }^{2} } {/{{{\big( {\big( {\tilde{a}_{\omega } \, {\minus}\, 1} \big)^{2} \big( {\tilde{a}_{\omega } \, {\minus}\, 2} \right)} \big)}}} \big. \kern-\nulldelimiterspace} }$$ for $$\tilde{a}_{\omega } \,\gt\,2$$ . We assume that the priors for all parameters are independent. In this case the posterior densities of parameters are of the same type as the prior densities, a so-called conjugate prior. The posterior distribution for each parameter is given by (we write, for ease of notation, y = y _{1 :T} , κ =κ _0:T and family $$\bipsi $$ _−λ means “parameter vector $$\bipsi $$ without the parameter λ”):

(49) $$\alpha _{x} \!\mid\!{\mib{y}} _{{1\: \colon\: T}} ,{\bikappa } ,{\mib \bipsi } _{{\, {\minus}\, \alpha _{x} }} {\rm \: \sim\: }N\left( {{{\tilde{\mu }_{\alpha } \sigma _{{{\varepsilon},x}}^{2} {\plus}\tilde{\sigma }_{\alpha }^{2} \mathop{\sum}\nolimits_t \left( {y_{{xt}} \, {\minus}\, \beta _{x} \kappa _{t} } \right)} \over {\tilde{\sigma }_{\alpha }^{2} T{\plus}\sigma _{{{\varepsilon},x}}^{2} }},{{\tilde{\sigma }_{\alpha }^{2} \sigma _{{{\varepsilon},x}}^{2} } \over {\tilde{\sigma }_{\alpha }^{2} T{\plus}\sigma _{{{\varepsilon},x}}^{2} }}} \right)$$

(50) $$\beta _{x} \!\mid\!{\mib{y}} _{{1\: \colon\: T}} ,{\mib \bikappa } ,{\mib \bipsi } _{{\, {\minus}\, \beta _{x} }} \: \sim\: N\left( {{{\tilde{\sigma }_{\beta }^{2} \mathop{\sum}\nolimits_t \left( {y_{{xt}} \, {\minus}\, \alpha _{x} } \right)\kappa _{t} {\plus}\tilde{\mu }_{\beta } \sigma _{{{\varepsilon},x}}^{2} } \over {\tilde{\sigma }_{\beta }^{2} \mathop{\sum}\nolimits_t \kappa _{t}^{2} {\plus}\sigma _{{{\varepsilon},x}}^{2} }},{{\tilde{\sigma }_{\beta }^{2} \sigma _{{{\varepsilon},x}}^{2} } \over {\tilde{\sigma }_{\beta }^{2} \mathop{\sum}\nolimits_t \kappa _{t}^{2} {\plus}\sigma _{{{\varepsilon},x}}^{2} }}} \right)$$

(51) $$\theta \!\mid\!{\mib{y}} _{{1\: \colon\: T}} ,{\mib \bikappa } ,{\mib \bipsi } _{{\, {\minus}\, \theta }} {\rm \: \sim\: }N\left( {{{\tilde{\sigma }_{\theta }^{2} \mathop{\sum}\nolimits_{t\, {\equals} \, 1}^T {\left( {\kappa _{t} \, {\minus}\, \kappa _{{t\, {\minus}\, 1}} } \right){\plus}\tilde{\mu }_{\theta } \sigma _{\omega }^{2} } } \over {\tilde{\sigma }_{\theta }^{2} T{\plus}\sigma _{\omega }^{2} }},{{\tilde{\sigma }_{\theta }^{2} \sigma _{\omega }^{2} } \over {\tilde{\sigma }_{\theta }^{2} T{\plus}\sigma _{\omega }^{2} }}} \right)$$

(52) $$\sigma _{{{\varepsilon},x}}^{2} \!\mid\!{\mib{y}} _{{1\: \colon\: T}} ,{\mib \bikappa } ,{\mib \bipsi } _{{\, {\minus}\, \sigma _{{{\varepsilon},x}}^{2} }} {\rm \: \sim\: IG}\left( {\tilde{a}_{{\varepsilon}} {\plus}{{pT} \over 2},\:\tilde{b}_{{\varepsilon}} {\plus}{1 \over 2}\mathop{\sum}\limits_{t\, {\equals} \, 1}^T \left( {y_{{xt}} \, {\minus}\, (\alpha _{x} {\plus}\beta _{x} \kappa _{t} )} \right)^{2} } \right)$$

(53) $$\sigma _{\omega }^{2} \!\mid\!{\mib{y}} _{{1\: \colon\: T}} ,{\mib \bikappa } ,{\mib \bipsi } _{{\, {\minus}\, \sigma _{\omega }^{2} }} {\rm \: \sim\: IG}\left( {\tilde{a}_{\omega } {\plus}{T \over 2},\:\tilde{b}_{\omega } {\plus}{1 \over 2}\mathop{\sum}\limits_{t\, {\equals} \, 1}^T \left( {\kappa _{t} \, {\minus}\, (\kappa _{{t\, {\minus}\, 1}} {\plus}\theta )} \right)^{2} } \right)$$

5.2. Non-linear/non-Gaussian state-space inference

In the case of non-linear/non-Gaussian state-space model dynamics such as the stochastic volatility models of LC3-H2 and the LCSV models, the sampler we develop is based on novel developments of the particle Metropolis–Hastings samplers of Andrieu et al. (Reference Andrieu, Doucet and Holenstein2010) adapted to the stochastic mortality models. In particular, we consider a combination of Rao-Blackwellised Kalman filter and particle filter for the latent state process full posterior conditionals, combined with a combination of Metropolis-within-Gibbs and Gibbs sampling steps for the static model parameters, both embedded within a PMCMC framework. We will illustrate the idea of this methodology for the LCSV model where a stochastic volatility dynamics is included in the latent process for the period effect.

5.2.1. Estimation for the LCSV mortality model

The static parameter vector is denoted as $${\mib \bipsi } \, {\equals} \, \left( {\alpha _{{x_{2} \: \colon\: x_{p} }} ,\beta _{{x_{2} \: \colon\: x_{p} }} ,\theta ,\sigma _{{\varepsilon}}^{2} ,\sigma _{\gamma }^{2} ,\lambda _{1} ,\lambda _{2} ,\gamma _{0} } \right)$$ . Note that we treat γ ₀ as a static parameter and our task is to obtain samples from the joint posterior distribution:

(54) $$\pi \left( {\kappa _{{0\: \colon\: T}} ,\gamma _{{1\: \colon\: T}} ,{\mib \bipsi } \!\mid\!{\mib{y}} _{{1\: \colon\: T}} } \right)$$

In this setting one can try a number of different approaches, the first would be to sample jointly from the full posterior distribution (54) via PMCMC methods to be described below. A second approach would be to combine PMCMC methods within a block-Gibbs-based sampler such as the following sampling scheme, where we apply Gibbs sampling to sample from the full conditional densities:

(55a) $$\pi \left( {\kappa _{{0\: \colon\: T}} \!\mid\!{\mib \bipsi } ,\gamma _{{1\: \colon\: T}} ,{\mib{y}} _{{1\: \colon\: T}} } \right)$$

(55b) $$\pi \left( {{\mib \bipsi } \!\mid\!\kappa _{{0\: \colon\: T}} ,\gamma _{{1\: \colon\: T}} ,{\mib{y}} _{{1\: \colon\: T}} } \right)$$

(55c) $$\pi \left( {\gamma _{{1\: \colon\: T}} \!\mid\!{\mib \bipsi } ,\kappa _{{0\: \colon\: T}} ,{\mib{y}} _{{1\: \colon\: T}} } \right)$$

Note that sampling from (55a) can be achieved by the FFBS procedure described in Algorithm 3, as one can apply Kalman filtering since γ _1:T is assumed to be given. The only difference compared to section 5.1.1 is that the term $$\sigma _{\omega }^{2} $$ is replaced by exp{γ _t } in Kalman filtering.

In the following, we provide details on how to sample from either the full posterior (54) or from full conditionals such as the density in (55c), via PMCMC method. Sampling from the posteriors of the static parameters (55b) is detailed in section 5.2.4.

5.2.2. PMCMC for mortality models

In this section, we explain the generic form of the PMCMC methodology that can be applied in a range of approaches for state-space stochastic mortality models. In general a PMCMC sampling method is a class of MCMC method where SMC algorithm is used as a proposal distribution within a MCMC algorithm. Though this seems trivial, it is actually based on a key observation that by using such a filter within the MCMC, the dimension of the acceptance probability in the Metropolis–Hastings acceptance–rejection stage is significantly reduced and can therefore facilitate much better mixing performance of the resulting Markov chain, reducing variance in estimation, see discussion in detail in Andrieu et al. (Reference Andrieu, Doucet and Holenstein2010).

To bring out the essence of PMCMC, we first discuss a generic approach to sample from a target distribution:

(56) $$\pi \left( {{\mib \biphi } _{{1\: \colon\: T}} ,{\mib \bipsi } \!\mid\!{\bf z}_{{1\: \colon\: T}} } \right)$$

where ϕ _1:T and $$\bipsi $$ are the latent state and static parameters of a general state-space model. Note, the state processes in this context are generally non-linear and potentially non-Gaussian.

From the perspective of obtaining the most efficiently mixing Markov chain to sample from this posterior, the ideal proposal distribution for constructing the Markov chain for $$\left( {{\mib \biphi } '_{{1\: \colon\: T}} ,{\mib \bipsi '} } \right)$$ is easily seen to be given by

(57) $$q\left( {{\mib \bipsi '} \!\mid\!{\mib \bipsi } } \right)p_{{{\mib \bipsi '} }} \left( {{\mib \biphi \prime \!}_{{1\: \colon\: T}}\!\mid\!{\mib{z}} _{{1\: \colon\: T}} } \right)$$

where q( $$\bipsi $$ ′| $$\bipsi $$ ) is a proposal for the parameters and the proposal for the latent state, $$p_{{{\mib \bipsi '} }} \left( {{\mib \biphi \prime } _{{1\: \colon\: T}}{\mib z}_{{1\: \colon\: T}} } \right)$$ , is from the state equation (given $$\bipsi $$ ′). Here ( ϕ _1:T , $$\bipsi $$ ) is the current state at MCMC iteration j−1 and $$\left( {{\mib \biphi \prime}\!_{{1\: \colon\: T}} ,{\mib \bipsi '} } \right)$$ is the proposed next move at MCMC iteration j.

In this ideal case, the acceptance probability of this ideal proposal is given by

(58) $$\alpha \left( {\left( {{\mib \biphi }{{\mib \prime} }\!_{{1\: \colon\: T}} ,{\mib \bipsi '} } \right),\left( {{\mib \biphi } _{{1\: \colon\: T}} ,{\mib \bipsi } } \right)} \right)\, {\equals} \, 1\wedge{{p\left( {{\mib \biphi }{{\mib \prime} }\! _{{1\: \colon\: T}} ,{\mib \bipsi '} \!\mid\!{\mib{z}} _{{1\: \colon\: T}} } \right)q\left( {{\mib \bipsi } \!\mid\!{\mib \bipsi '} } \right)p_{{\mib \bipsi } } \left( {{\mib \biphi } _{{1\: \colon\: T}} \!\mid\!{\mib{z}} _{{1\: \colon\: T}} } \right)} \over {p\left( {{\mib \biphi } _{{1\: \colon\: T}} ,{\mib \bipsi } \!\mid\!{\mib{z}} _{{1\: \colon\: T}} } \right)q\left( {{\mib \bipsi '} \!\mid\!{\mib \bipsi } } \right)p_{{{\mib \bipsi '} }} \left( {{\mib \biphi } _{{1\: \colon\: T}}^{{\mib \prime} } \!\mid\!{\mib{z}} _{{1\: \colon\: T}} } \right)}}$$

(59) $$\, {\equals} \, 1\wedge{{p_{{{\mib \bipsi '} }} \left( {{\mib \biphi } _{{1\: \colon\: T}} \!\mid\!{\bf z}_{{1\: \colon\: T}} } \right)p\left( {{\mib \bipsi '} \!\mid\!{\mib{z}} _{{1\: \colon\: T}} } \right)q\left( {{\mib \bipsi } \!\mid\!{\mib \bipsi '} } \right)p_{{\mib \bipsi } } \left( {{\mib \biphi } _{{1\: \colon\: T}} \!\mid\!{\mib{z}} _{{1\: \colon\: T}} } \right)} \over {p_{{\mib \bipsi } } \left( {{\mib \biphi } _{{1\: \colon\: T}} \!\mid\!{\mib{z}} _{{1\: \colon\: T}} } \right)p\left( {{\mib \bipsi } \!\mid\!{\mib{z}} _{{1\: \colon\: T}} } \right)q\left( {{\mib \bipsi '} \!\mid\!{\mib \bipsi } } \right)p_{{{\mib \bipsi '} }} \left( {{\mib \biphi }{{\mib \prime} }\! _{{1\: \colon\: T}} \!\mid\!{\mib{z}} _{{1\: \colon\: T}} } \right)}}$$

(60) $$\, {\equals} \, 1\wedge{{p\left( {{\mib \bipsi '} \!\mid\!{\mib{z}} _{{1\: \colon\: T}} } \right)q\left( {{\mib \bipsi } \!\mid\!{\mib \bipsi '} } \right)} \over {p\left( {{\mib \bipsi } \!\mid\!{\mib{z}} _{{1\: \colon\: T}} } \right)q\left( {{\mib \bipsi '} \!\mid\!{\mib \bipsi } } \right)}}$$

(61) $$\, {\equals} \, 1\wedge{{p_{{{\mib \bipsi '} }} ({\mib{z}} _{{1\: \colon\: T}} )p({\mib \bipsi '} )q\left( {{\mib \bipsi } \!\mid\!{\mib \bipsi '} } \right)} \over {p_{{\mib \bipsi } } ({\mib{z}} _{{1\: \colon\: T}} )p({\mib \bipsi } )q({\mib \bipsi '} \!\mid\!{\mib \bipsi } )}}$$

where r ₁Λr ₂:=min(r ₁, r ₂). A desirable property of the ideal proposal is that the acceptance probability depends only on the marginal likelihood, together with the prior and proposal for the static parameters. This is optimal in the sense that the dimension of the numerator and denominator is reduced significantly to the static model parameter dimensions, and not including explicitly the path-space latent process dimensions, a reduction of d×T dimensions for a d-dimensional state vector ϕ _t . However, clearly one can never achieve this goal as it requires perfect knowledge of $$p_{{{\mib \bipsi '} }} \left( {{\mib \biphi }{{\mib \prime} }\!_{{1\: \colon\: T}} } \!\mid\!{\mib{z}} _{{1\: \colon\: T}} } \right)$$ as well as the ability to sample this distribution, both of which are unachievable except in the special case of the Linear Gaussian case explained in section 5.1.1.

To circumvent this problem, the particle marginal Metropolis–Hastings sampler (PMMH; Andrieu et al., Reference Andrieu, Doucet and Holenstein2010) applies SMC method to obtain an approximate of the state transition density (which is also the state proposal):

(62) $$\hat{p}_{{{\mib \bipsi '} }} \left( {{\mib \biphi } _{{1\: \colon\: T}} \!\mid\!{\mib{z}} _{{1\: \colon\: T}} } \right)\, {\equals} \, \mathop{\sum}\limits_{i\, {\equals} \, 1}^N w_{T}^{{(i)}} \delta _{{{\mib \biphi } _{{1\: \colon\: T}}^{{(i)}} }} \left( {{\mib \biphi } _{{1\: \colon\: T}} } \right)$$

where $$w_{T}^{{(i)}} $$ is the importance weight, $$\delta _{x} (X)$$ a Dirac mass function centred at X and a proposed next move of the latent state is drawn from this discrete approximate distribution. Moreover, a by-product of a SMC algorithm is the marginal likelihood, $$\hat{p}_{{\mib \bipsi } } ({\mib{z}} _{{1\: \colon\: T}} )$$ , which has the following important property:

Lemma 5.1 A SMC proposal admits as a by-product an unbiased estimator of the marginal likelihood p $$_\bipsi $$ ( z _1:T ) given by

(63) $$\hat{p}_{{\mib \bipsi } } \left( {{\mib{z}} _{{1\: \colon\: T}} } \right)\,\colon\! {\equals} \, \prod\limits_{t\, {\equals} \, 2}^T \hat{p}_{{\mib \bipsi } } \left( {{\mib{z}} _{t} \!\mid\!{\mib{z}} _{{1\: \colon\: t\, {\minus}\, 1}} } \right)$$

where a SMC approximation with N-particles produces, for all t

(64) $$\hat{p}_{{\mib \bipsi } } \left( {{\mib{z}} _{t} \!\mid\!{\mib{z}} _{{1\: \colon\: t\, {\minus}\, 1}} } \right)\, {\equals} \, {1 \over N}\mathop{\sum}\limits_{i\, {\equals} \, 1}^N w_{t}^{{(i)}} $$

which is an unbiased particle estimate of p $$_\bipsi $$ ( z _t | z _1:t−1). This non-trivial unbiasedness was first presented in Del Moral (Reference Del Moral2004) and has since been utilised to great advantage as explained in Chopin et al. (Reference Chopin, Jacob and Papaspiliopoulos2013). In addition, the variance of this estimator typically only grows linearly with T.

The unbiased approximate marginal likelihoods are then used in the acceptance probability (61):

(65) $$\alpha \left( {\left( {{\mib \biphi } _{{1\: \colon\: T}}^{{\mib \prime} } ,{\mib \bipsi '} } \right),\left( {{\mib \biphi } _{{1\: \colon\: T}} ,{\mib \bipsi } } \right)} \right)\, {\equals} \, 1\wedge{{\hat{p}_{{{\mib \bipsi '} }} \left( {{\mib{z}} _{{1\: \colon\: T}} } \right)p\left( {{\mib \bipsi '} } \right)q\left( {{\mib \bipsi } \!\mid\!{\mib \bipsi '} } \right)} \over {\hat{p}_{{\mib \bipsi } } \left( {{\mib{z}} _{{1\: \colon\: T}} } \right)p\left( {\mib \bipsi } \right)q\left( {{\mib \bipsi '} \!\mid\!{\mib \bipsi } } \right)}}$$

Due to the unbiasedness of the estimated marginal likelihood, Andrieu et al. (Reference Andrieu, Doucet and Holenstein2010) show that, even though only SMC approximates are used (with finite number of particles N), the invariant distribution of PMMH is the target distribution π( ϕ _1:T , $$\bipsi $$ | z _1:T ).

To apply PMCMC for an efficient estimation of the LCSV model, we first notice that we can obtain explicitly the posteriors of static parameters via conjugate priors. As a result we are only required to sample from the density π(γ _1:T | $$\bipsi $$ , κ _0:n , y _1:T ), instead of the joint density π(γ _1:T , $$\bipsi $$ |κ _0:n , y _1:T ). It turns out that there is a class of PMCMC algorithm, called particle-independent Metropolis–Hastings sampler (PIMH), which provide a mechanism to sample exactly from π(γ_1:T | $$\bipsi $$ , κ _0:n , y _1:T ).

Our approach to sampling from the joint posterior distribution, π(k _0:T , γ _1:T, $$\bipsi $$ | y _1:T , of the LCSV model is summarised in Algorithm 4.

5.2.3. PIMH: sampling from π(γ _1:T | $$\bipsi $$ , κ _0:n , y _1:T )

We first note that π(γ _1:T | $$\bipsi $$ , κ _0:n , y _1:T =π $$_\bipsi $$ (γ _1:T |k _0:n ) given the structure of the LCSV model. Using an independent proposal density, q $$_\bipsi $$ (γ_1:T |κ_0:n ), in the Metropolis–Hastings algorithm, the acceptance probability is given by

(66) $$\alpha \left( {\gamma '_{{1\: \colon\: T}} ,\gamma _{{1\: \colon\: T}} } \right)\, {\equals} \, 1\wedge{{\pi _{{\bf \psi }} (\gamma{{\mib \prime} }\! _{{1\: \colon\: T}} \!\mid\!\kappa _{{0\: \colon\: n}} )q_{{\bf \psi }} \left( {\gamma _{{1\: \colon\: T}} \!\mid\!\kappa _{{0\: \colon\: n}} } \right)} \over {\pi _{{\bf \psi }} (\gamma _{{1\: \colon\: T}} \!\mid\!\kappa _{{0\: \colon\: n}} )q_{{\bf \psi }} \left( {\gamma _{{1\: \colon\: T}}^{{\mib \prime} } \!\mid\!\kappa _{{0\: \colon\: n}} } \right)}}$$

Ideally, one may take q $$_\bipsi $$ (γ _1:T |κ _0:n )=π $$_\bipsi $$ (γ _1:T |κ _0:n ). However, in most cases such an ideal choice is impossible to sample from and to evaluate. The PIMH sampler proposes instead to use the SMC approximation $$\hat{\pi }_{\psi } (\gamma _{{1\: \colon\: T}} \!\mid\!\kappa _{{0\: \colon\: n}} )$$ as the proposal density and calculate the acceptance probability as

(67) $$\alpha \left( {\gamma{{\mib \prime} }\! _{{1\: \colon\: T}} ,\gamma _{{1\: \colon\: T}} } \right)\, {\equals} \, 1\wedge{{\hat{\pi }_{{\mib \bipsi } } \left( {\kappa _{{0\: \colon\: n}} } \right){\prime} } \over {\hat{\pi }_{{\mib \bipsi } } \left( {\kappa _{{0\: \colon\: n}} } \right)[j\, {\minus}\, 1]}}$$

where $$\hat{\pi }_{{\mib \bipsi } } \left( {\kappa _{{0\: \colon\: n}} } \right){\prime} $$ and $$\hat{\pi }_{{\mib \bipsi } } \left( {\kappa _{{0\: \colon\: n}} } \right)[j\, {\minus}\, 1]$$ are unbiased marginal likelihoods estimated by SMC (see Lemma 5.1) in the current MCMC iteration j and the previous iteration j−1, respectively. It can be shown that the invariant distribution of the PIMH sampler is the target distribution π $$_\bipsi $$ (γ _1:T |κ _0:n ) (Andrieu et al., Reference Andrieu, Doucet and Holenstein2010).

It remains to specify an SMC approximation $$\hat{\pi }_{{\mib \bipsi } } \left( {\gamma _{{1\: \colon\: T}} \!\mid\!\kappa _{{0\: \colon\: n}} } \right)$$ (Appendix B). We use the so-called bootstrap filter, that is, the proposal distribution in the SMC algorithm to draw γ _t is given by the state equation (23c):

(68) $$g_{t} \left( {\gamma _{t} \!\mid\!\gamma _{{1\: \colon\: t\, {\minus}\, 1}} ,\kappa _{{0\: \colon\: t}} } \right)\,\colon\! {\equals} \, \pi _{{\mib \bipsi } } \left( {\gamma _{t} \!\mid\!\gamma _{{t\, {\minus}\, 1}} } \right)$$

Consequently, the importance weight is evaluated as

(69) $$\tilde{w}_{t} \propto\tilde{w}_{{t\, {\minus}\, 1}} \pi _{{\mib \bipsi } } \left( {\kappa _{t} \!\mid\!\gamma _{t} ,\kappa _{{t\, {\minus}\, 1}} } \right)$$

where π(κ _t |γ _t , κ _t−1) is the incremental importance weight. Our approach for sampling from π $$_\bipsi $$ (γ _1:T |κ _0:T ) is summarised in Algorithm 5 (together with Algorithm 6).

5.2.4. Sampling from π( $$\bipsi $$ |κ _0:T , γ _1:T , y _1:T )

We assume the prior for $$\big( {\alpha _{{x_{2} \: \colon\: x_{p} }} ,\,\beta _{{x_{2} \: \colon\: x_{p} }} ,\,\theta ,\,\sigma _{{\varepsilon}}^{2} ,\,\sigma _{\gamma }^{2} ,\,\lambda _{1} ,\,\lambda _{2} ,\,\gamma _{0} } \big)$$ is given by

(75a) $$\alpha _{x} \: \sim\: N\left( {\tilde{\bimu }_{\alpha } ,\tilde{\sigma }_{\alpha }^{2} } \right),\quad \beta _{x} \: \sim\: N\left( {\tilde{\bimu }_{\beta } ,\tilde{\sigma }_{\beta }^{2} } \right),\quad \theta \: \sim\: N\left( {\tilde{\bimu }_{\theta } ,\tilde{\sigma }_{\theta }^{2} } \right),\quad \sigma _{{\varepsilon}}^{2} \: \sim\: {\rm IG}\left( {\tilde{a}_{{\varepsilon}} ,\tilde{b}_{{\varepsilon}} } \right)$$

(75b) $$\sigma _{\gamma }^{2} \: \sim\: {\rm IG}\left( {\tilde{a}_{\gamma } ,\,\tilde{b}_{\gamma } } \right),\:\,\lambda _{1} \: \sim\: N_{{[\, {\minus}\, 1,1]}} \left( {\tilde{\bimu }_{{\lambda _{1} \!}} ,\,\tilde{\sigma }_{{\lambda _{1} }}^{2} } \right),\,\lambda _{2} \: \sim\: N\left( {\tilde{\bimu }_{{\lambda _{2} }\!} ,\,\tilde{\sigma }_{{\lambda _{2} }}^{2} } \right)\:\gamma _{0} \: \sim\: N\left( {\tilde{\bimu }_{{\gamma _{0\!} }} ,\,\tilde{\sigma }_{{\gamma _{0} }}^{2} } \right)$$

where x=x ₂, … , x _p and N _{[−1, 1]} denotes a truncated Gaussian with support [−1, 1]. It is assumed that the priors for all parameters are independent.

Samples from the density π( $$\bipsi $$ |κ _0:T , γ _1:T , Y _1:T ) are obtained by sampling from the following posteriors:

(76) $$\alpha _{x} \!\mid\!{\mib{y}} _{{1\: \colon\: T}} ,\,{\mib{\bikappa }} ,\,{\mib{\bigamma }} ,\,{\mib{\bipsi }} _{{\, {\minus}\, \alpha _{x} }} \: \sim\: N\left( {{{\tilde{\mu }_{\alpha } \sigma _{{\varepsilon}}^{2} {\plus}\tilde{\sigma }_{\alpha }^{2} \mathop{\sum}\nolimits_t {\left( {y_{{xt}} \, {\minus}\, \beta _{x} \kappa _{t} } \right)} } \over {T\tilde{\sigma }_{\alpha }^{2} {\plus}\sigma _{{\varepsilon}}^{2} }},{{\tilde{\sigma }_{\alpha }^{2} \sigma _{{\varepsilon}}^{2} } \over {T\tilde{\sigma }_{\alpha }^{2} {\plus}\sigma _{{\varepsilon}}^{2} }}} \right)$$

(77) $$\beta _{x} \!\mid\!{\mib{y}} _{{1\: \colon\: T}} ,\,{\mib{\kappa }} ,\,{\mib \gamma } ,\,{\mib \bipsi } _{{\, {\minus}\, \beta _{x} }} \: \sim\: N\left( {{{\tilde{\sigma }_{\beta }^{2} \mathop{\sum}\nolimits_t {\left( {y_{{xt}} \, {\minus}\, \alpha _{x} } \right)\kappa _{t} {\plus}\tilde{\mu }_{\beta } \sigma _{{\varepsilon}}^{2} } } \over {\tilde{\sigma }_{\beta }^{2} \mathop{\sum}\nolimits_t {\kappa _{t}^{2} {\plus}\sigma _{{\varepsilon}}^{2} } }},{{\tilde{\sigma }_{\beta }^{2} \sigma _{{\varepsilon}}^{2} } \over {\tilde{\sigma }_{\beta }^{2} \mathop{\sum}\nolimits_t {\kappa _{t}^{2} {\plus}\sigma _{{\varepsilon}}^{2} } }}} \right)$$

(78) $$\theta \!\mid\!{\mib{y}} _{{1\: \colon\: T}} ,\,{\mib \kappa } ,\,{\mib \gamma } ,\,{\mib \bipsi } _{{\, {\minus}\, \theta }} \: \sim\: N\left( {{{\tilde{\mu }_{\theta } \,/\,\tilde{\sigma }_{\theta }^{2} {\plus}\mathop{\sum}\nolimits_t {\left( {\kappa _{t} \, {\minus}\, \kappa _{{t\, {\minus}\, 1}} } \right)\,/\,e^{{\gamma _{t} }} } } \over {1\,/\,\tilde{\sigma }_{\theta }^{2} {\plus}\mathop{\sum}\nolimits_t {1\,/\,e^{{\gamma _{t} }} } }},{1 \over {1\,/\,\tilde{\sigma }_{\theta }^{2} {\plus}\mathop{\sum}\nolimits_t {1\,/\,e^{{\gamma _{t} }} } }}} \right)$$

(79) $$\sigma _{{\varepsilon}}^{2} \!\mid\!{\mib{y}} _{{1\: \colon\: T}} ,\,{\mib{\kappa }} ,\,{\mib \gamma } ,\,{\mib \bipsi } _{{\, {\minus}\, \sigma _{{\varepsilon}}^{2} }} \: \sim\: {\rm IG}\left( {\tilde{a}_{{\varepsilon}} {\plus}{{pT} \over 2},\:\,\tilde{b}_{{\varepsilon}} {\plus}{1 \over 2}\mathop{\sum}\limits_{t\, {\equals} \, 1}^T {\mathop{\sum}\limits_{x\, {\equals} \, 1}^p {\left( {y_{{xt}} \, {\minus}\, \left( {\alpha _{x} {\plus}\beta _{x} \kappa _{t} } \right)} \right)^{2} } } } \right)$$

(80) $$\sigma _{\gamma }^{2} \!\mid\!{\mib{y}} _{{1\: \colon\: T}} ,\,{\mib{\kappa }} ,\,{\mib \gamma } ,\,{\mib \bipsi } _{{\, {\minus}\, \sigma _{\gamma }^{2} }} \: \sim\: {\rm IG}\left( {\tilde{a}_{\gamma } {\plus}{T \over 2},\,\:\tilde{b}_{\gamma } {\plus}{1 \over 2}\mathop{\sum}\limits_{t\, {\equals} \, 1}^T \left( {\gamma _{t} \, {\minus}\, \lambda \gamma _{{t\, {\minus}\, 1}} } \right)^{2} } \right)$$

(81) $$\lambda _{1} \!\mid\!{\mib{y}} _{{1\: \colon\: T}} ,\,{\mib{\kappa }} ,{\mib \gamma } ,{\mib \bipsi } _{{\, {\minus}\, \lambda _{1} }} \: \sim\: N_{{[\, {\minus}\, 1,1]}} \left( {{{\sigma _{\gamma }^{2} \tilde{\mu }_{{\lambda _{1} }} {\plus}\tilde{\sigma }_{{\lambda _{1} }}^{2} \mathop{\sum}\nolimits_t {\gamma _{{t\, {\minus}\, 1}} \gamma _{t} } } \over {\sigma _{\gamma }^{2} {\plus}\tilde{\sigma }_{{\lambda _{1} }}^{2} \mathop{\sum}\nolimits_t {\gamma _{{t\, {\minus}\, 1}}^{2} } }},\,{{\tilde{\sigma }_{{\lambda _{1} }}^{2} \sigma _{\gamma }^{2} } \over {\sigma _{\gamma }^{2} {\plus}\tilde{\sigma }_{{\lambda _{1} }}^{2} \mathop{\sum}\nolimits_t {\gamma _{{t\, {\minus}\, 1}}^{2} } }}} \right)$$

(82) $$\lambda _{2} \!\mid\!{\mib{y}} _{{1\: \colon\: T}} ,\,{\mib \kappa } ,\,{\mib \gamma } ,\,{\mib \bipsi } _{{\, {\minus}\, \lambda _{2} }} \: \sim\: N\left( {{{\sigma _{\gamma }^{2} \tilde{\mu }_{{\lambda _{2} }} {\plus}\tilde{\sigma }_{{\lambda _{2} }}^{2} \mathop{\sum}\nolimits_t {\left( {\gamma _{t} \, {\minus}\, \lambda _{1} \gamma _{{t\, {\minus}\, 1}} } \right)} } \over {\sigma _{\gamma }^{2} {\plus}T\tilde{\sigma }_{{\lambda _{2} }}^{2} }},\,{{\tilde{\sigma }_{{\lambda _{2} }}^{2} \sigma _{\gamma }^{2} } \over {\sigma _{\gamma }^{2} {\plus}T\tilde{\sigma }_{{\lambda _{2} }}^{2} }}} \right)$$

(83) $$\gamma _{0} \!\mid\!{\mib{y}} _{{1\: \colon\: T}} ,{\mib \kappa } ,{\mib \gamma } ,{\mib \bipsi } _{{\, {\minus}\, \gamma _{0} }} \: \sim\: N\left( {{{\sigma _{\gamma }^{2} \tilde{\mu }_{{\gamma _{0} }} {\plus}\tilde{\sigma }_{{\gamma _{0} }}^{2} \lambda \gamma _{1} } \over {\sigma _{\gamma }^{2} {\plus}\tilde{\sigma }_{{\gamma _{0} }}^{2} \lambda ^{2} }},\,{{\tilde{\sigma }_{{\gamma _{0} }}^{2} \sigma _{\gamma }^{2} } \over {\sigma _{\gamma }^{2} {\plus}\tilde{\sigma }_{{\gamma _{0} }}^{2} \lambda ^{2} }}} \right)$$

where the posterior distributions are obtained similarly as in section 5.1.2.

6.1. Data description

The data set consists of Danish male population death rates for 21 age groups (0, 1–4, 5–9, … , 95–99) from year 1835 to 2010 where we fix year 2010 as the end year. Figure 1 displays some of the time series of the log death rates for the Danish male population. It is clear that the multi-dimensional time series exhibit different volatility for different age groups, which justify the introduction of heteroscedasticity into the observation equation as discussed in section 3.1. We also observe that, mainly before 1950, there are periods that the volatility of death rates for some age groups are markedly different. Such a change of volatility in the temporal dimension suggests that stochastic volatility may be present in the underlying time period effect.

Figure 1 Time series of log death rates for Danish male population from year 1835 to 2010.

6.2. Estimation results

In our empirical study we focus on Bayesian inference and forecasting. We assume vague priors so that all inferences are mainly based on the data and the impact of the prior is not material. Taking the LCSV model as an example, we assume κ ₀ $$\sim$$ N(0, 10), α _x $$\sim$$ N(0, 10), β _x $$\sim$$ N(0, 10), θ $$\sim$$ N(0, 10), $$\sigma _{{\varepsilon}}^{2} \: \sim\: {\rm IG}(2.001,0.001)$$ , $$\sigma _{\gamma }^{2} \: \sim\: {\rm IG}(2.001,0.001)$$ , λ ₁ $$\sim$$ N(0, 10), λ ₂ $$\sim$$ N(0, 10) and γ ₀ $$\sim$$ N(0, 10), where x∈{x ₂, … , x _p }. The number of iterations of the Markov chain is 15,000 with 5,000 burn-in. We fix $$\alpha _{{x_{1} }} \, {\equals} \, {1 \over T}\mathop{\sum}\nolimits_{t\, {\equals} \, 1}^T y_{{x_{1} ,t}} $$ and $$\beta _{{x_{1} }} \, {\equals} \, 0.2$$ as an identification constraint; see Remark 3.1 for a discussion for our choice of values here.

Estimated values of the static parameters (except α and β ) for the Danish mortality data (1835–2010) are shown in Table 3. The rest of the estimated parameters and states are displayed in Figure 2. Here, we only show the plots for the LCSV-H model since the corresponding figures obtained from the LC, LC-H and LCSV model are visually similar to the case of the LCSV-H model.

Figure 2 Estimation of (upper panels) α , β and $$\sigma _{{x_{1} \,\colon\,x_{{21}} ,{\epsilon}}}^{2} $$ ; (lower panels) time effect κ _1834:2010, log-volatility γ _1835:2010 and first difference $$\Delta \bar{\kappa }_{t} $$ , for Danish male mortality data (1835–2010) using the LCSV-H model. LCSV-H, Lee–Carter stochastic volatility model with heteroscedasticity; CI, confidence interval.

Table 3 Estimated values of the static parameters (except α and β ) for the Danish male mortality data (1835–2010).

The range in (,) represents 95% credible interval.

LC, Lee–Carter; LC-H, LC model with heteroscedasticity; LCSV, LC stochastic volatility; LCSV-H, LCSV model with heteroscedasticity.

It is evident from Figure 2 that there are periods, namely 1850–1870, 1910–1920, 1930–1950, that the time effect κ exhibits higher volatility compared with other periods. We also observe that κ accelerates markedly downward after 1990 and is relatively smooth in the recent period 1950–2010. The filtering of the log-volatility process γ _1835:2010 (Figure 2) quantifies the volatility level $$\left( {e^{{\gamma _{t} }} } \right)$$ of the time effect and gives further evidence on the stochastic volatility nature of mortality. To see more clearly the phenomenon of changing volatility, we plot the first difference $$\Delta \bar{\kappa }_{t} \, {\equals} \, \bar{\kappa }_{t} \, {\minus}\, \bar{\kappa }_{{t\, {\minus}\, 1}} $$ in Figure 2 for the LCSV-H model, where $$\bar{\kappa }_{t} $$ denotes the posterior mean of κ _t , t=1836, … , 2010. It shows evidently the change of volatility level in the latent process κ _t . The patterns of the estimated log-volatility γ _1835:2010 and the first difference $$\Delta \bar{\kappa }_{t} $$ clearly suggest that it is not appropriate to assume constant volatility $$\left( {\sigma _{\omega }^{2} } \right)$$ for the time effect.

The state-space modelling approach is able to uncover the age-specific heteroscedasticity structure hidden in the Danish mortality time series. Figure 2 reveals that variability is particularly high for the very young and very old age group. Implications of the heteroscedastic structure on forecasting will be discussed in section 6.4.

To investigate the forecasting properties of the stochastic volatility model, we also estimate the models based on calibration periods 1835–1990 and 1950–1990. Figures 3 and 4 show the estimated parameters and states for the LCSV-H model in those periods.

Figure 3 Estimation of (upper panels) α , β and $$\sigma _{{x_{1} \,\colon\,x_{{21}} ,{\epsilon}}}^{2} $$ ; (lower panels) time effect κ _1834:1900, log-volatility γ _1835:1990 and first difference $$\Delta \bar{\kappa }_{t} $$ , for Danish male mortality data (1835–1990) using the LCSV-H model. LCSV-H, Lee–Carter stochastic volatility model with heteroscedasticity; CI, confidence interval.

Figure 4 Estimation of (upper panels) α , β and $$\sigma _{{x_{1} \,\colon\,x_{{21}} ,{\epsilon}}}^{2} $$ ; (lower panels) time effect κ _1949:1990, log-volatility γ _1950:1990 and first difference $$\Delta \bar{\kappa }_{t} $$ , for Danish male mortality data (1950–1990) using the LCSV-H model. LCSV-H, Lee–Carter stochastic volatility model with heteroscedasticity; CI, confidence interval.

6.3. Model assessment

To compare the fit of the models to the data, we apply deviance information criterion (DIC) as a Bayesian measures of model complexity and fit (Spiegelhalter et al., Reference Spiegelhalter, Best, Carlin and van der Linde2002). It is common to assess and compare models with latent variables using conditional DIC (Berg et al., Reference Berg, Meyer and Yu2004; Celeux et al., 2006). Specifically, we use the so-called conditional log-likelihood which is calculated as

(84) $${\rm ln} f\left( {{\mib{y}} _{{1\: \colon\: T}} \!\mid\!{\mib \bipsi } ,\kappa _{{1\: \colon\: T}} } \right)\, {\equals} \, \mathop{\sum}\limits_{x\, {\equals} \, x_{1} }^{x_{p} } \mathop{\sum}\limits_{t\, {\equals} \, 1}^T \left( {\, {\minus}\, {1 \over 2}{\rm ln}\,2\pi \, {\minus}\, {\rm ln}\,\sigma _{{{\varepsilon},x}} \, {\minus}\, {1 \over 2}\left( {{{y_{{x,t}} \, {\minus}\, \left( {\alpha _{x} {\plus}\beta _{x} \kappa _{t} } \right)} \over {\sigma _{{{\varepsilon},x}} }}} \right)^{2} } \right)$$

Note that the likelihood is conditional on parameters that include both static parameters and the latent process κ. Using the conditional log-likelihood function, the deviance is defined as

(85) $$D({\bf \Psi })\, {\equals} \, \, {\minus}\, 2{\rm ln}\,f\left( {{\mib{y}} _{{1\: \colon\: T}} \!\mid\!{\bf \Psi }){\plus}2{\rm ln}\,h({\mib{y}} _{{1\: \colon\: T}} } \right)$$

where $${\bf \Psi } $$ =( $$\bipsi $$ , κ _1:T ) and we assume h( y _1:T )=1 since in the models we consider it plays the role of a constant which is the same for competing models. The effective dimension, p _D , is evaluated as

(86) $$p_{D} \, {\equals} \, \bar{D}({\bf \Psi })\, {\minus}\, D({\bf \bar{\Psi }})$$

where $$\bar{D}({\bf \Psi })$$ and $${\bf \bar{\Psi }}$$ denote, respectively, the mean of D( $${\bf \Psi } $$ ) and the mean of the posterior distribution of $${\bf \Psi } $$ . The conditional DIC is then given by

(87) $${\rm DIC}\,\colon\! {\equals} \, \bar{D}({\bf \Psi }){\plus}p_{D} \, {\equals} \, 2\bar{D}({\bf \Psi })\, {\minus}\, D({\bf \bar{\Psi }})$$

which can be evaluated straightforwardly using MCMC samples.

The DIC values for the models with different calibration periods are shown in Table 4 Footnote ⁶ . The LCSV-H model clearly outperforms other models for all calibration periods considered. The LC-H is the second best; in fact it even outperforms the LCSV model for long calibration periods 1835–2010 and 1835–1990. These results suggest that one is strongly encouraged to include heteroscedasticity structure when modelling Danish mortality data. Note that for long calibration periods, the incorporation of stochastic volatility can markedly improve model fit as can be seen by comparing the LCSV-H model and the LC-H model, as well as the comparison between LCSV model and the LC model.

Table 4 Deviance information criterion of models with different calibration periods.

LC, Lee–Carter; LC-H, LC model with heteroscedasticity; LCSV, LC stochastic volatility; LCSV-H, LCSV model with heteroscedasticity.

6.4. Forecasting

In this section, we investigate the forecasting properties of the mortality models summarised in Table 2 where heteroscedasticity as well as stochastic volatility structures are incorporated. Our analysis is based on the forecasting distributions of (log) death rates and life expectancy. The Bayesian state-space framework allows us to obtain the forecasting distributions using MCMC samples which is shown below.

6.4.1. Death rates

For the LC (LC-H) model, the k-step ahead forecasting distribution of y _T+k , given y _1:T , is given by

(88) $$\pi \left( {{\mib{y}} _{{T{\plus}k}} \!\mid\!{\mib{y}} _{{1\: \colon\: T}} } \right)\, {\equals} \, {\int} {\pi \left( {{\mib{y}} _{{T{\plus}k}} \!\mid\!\kappa _{{T{\plus}k}} ,{\mib \bipsi } } \right)\pi \left( {\kappa _{{T{\plus}k}} \!\mid\!\kappa _{{T{\plus}k\, {\minus}\, 1}} ,{\mib \bipsi } } \right)\,\ldots\,\pi \left( {\kappa _{T} ,{\mib \bipsi } \!\mid\!{\mib{y}} _{{1\: \colon\: T}} } \right)\:d{\mib \bipsi } d\kappa _{{T\: \colon\: T{\plus}k}} } $$

where $$\bipsi $$ is the parameter vector for the LC (LC-H) model; (88) suggests that we can sample recursively to obtain the forecasting distribution, for k≥1, as follows:

(89a) $$\kappa _{{T{\plus}k}}^{{(\ell )}} \: \sim\: N\left( {\kappa _{{T{\plus}k\, {\minus}\, 1}}^{{(\ell )}} {\plus}\theta ^{{(\ell )}} ,\left( {\sigma _{\omega }^{2} } \right)^{{(\ell )}} } \right)$$

(89b) $${\mib{y}} _{{T{\plus}k}}^{{(\ell )}} \: \sim\: N\left( {{\mib{\alpha }} ^{{(\ell )}} {\plus}{\mib{\beta }} ^{{(\ell )}} \kappa _{{T{\plus}k}}^{{(\ell )}} ,\,\Sigma ^{{(\ell )}} } \right)$$

where $$\ell \, {\equals} \, 1,\,\ldots\:,L$$ and L is the number of MCMC iterations after burn-in. Here Σ is a diagonal matrix with $$\sigma _{{{\varepsilon},x}}^{2} $$ on the diagonal for the LC-H model and $$\sigma _{{\varepsilon}}^{2} $$ for the LC model. This procedure generates an estimate of the forecasting distribution.

Similarly, the forecasting distribution of y _T+k , given y _T , for the LCSV (LCSV-H) model is given by

(90) $$\eqalignno{ \pi \left( {{\mib{y}} _{{T{\plus}k}} \!\mid\!{\mib{y}} _{{1\: \colon\: T}} } \right)\, {\equals} \, &#x0026; {\int} {\pi \left( {{\mib{y}} _{{T{\plus}k}} \!\mid\!\kappa _{{T{\plus}k}} ,\,{\mib \bipsi } } \right)\pi \left( {\kappa _{{T{\plus}k}} \!\mid\!\kappa _{{T{\plus}k\, {\minus}\, 1}} ,\,\gamma _{{T{\plus}k}} ,\,{\mib \bipsi } } \right)\,\ldots\,} \cr &#x0026;\pi \left( {\gamma _{{T{\plus}1}} \!\mid\!\gamma _{T} ,\,{\mib \bipsi } } \right)\pi \left( {\kappa _{T} ,\,\gamma _{T} ,\,{\mib \bipsi } \!\mid\!{\mib{y}} _{{1\: \colon\: T}} } \right)d{\mib \bipsi } d\kappa _{{T\: \colon\: T{\plus}k}} d\gamma _{{T\: \colon\: T{\plus}k}} $$

For k≥1, the forecasting distribution can be obtained by sampling recursively

(91a) $$\gamma _{{T{\plus}k}}^{{(\ell )}} \: \sim\: N\left( {\lambda _{1}^{{(\ell )}} \gamma _{{T{\plus}k\, {\minus}\, 1}}^{{(\ell )}} {\plus}\lambda _{2}^{{(\ell )}} ,\,\big( {\sigma _{\gamma }^{2} } \big)^{{(\ell )}} } \right)$$

(91b) $$\kappa _{{T{\plus}k}}^{{(\ell )}} \: \sim\: N\left( {\kappa _{{T{\plus}k\, {\minus}\, 1}}^{{(\ell )}} {\plus}\theta ^{{(\ell )}} ,exp\left\{ {\gamma _{{T{\plus}k}}^{{(\ell )}} } \right\}} \right)$$

(91c) $${\mib{y}} _{{T{\plus}k}}^{{(\ell )}} \: \sim\: N\left( {{\mib \bialpha } ^{{(\ell )}} {\plus}{\mib \bibeta } ^{{(\ell )}} \kappa _{{T{\plus}k}}^{{(\ell )}} ,\,\Sigma ^{{(\ell )}} } \right)$$

where $$\ell \, {\equals} \, 1,\,\ldots\,,L$$ , and Σ is a diagonal matrix with $$\sigma _{{{\varepsilon},x}}^{2} $$ on the diagonal for the LCSV-H model and $$\sigma _{{\varepsilon}}^{2} $$ for the LCSV model.

Figure 5 shows the forecasted log death rates based on the LC-H, LCSV and LCSV-H model, using the LC model as a benchmark. We show age groups 5–9, 35–39, 65–69 and 95–99 as representatives of young, adult, old and very old age. The models are estimated using data for the period 1835–2010 and forecast for 30 years.

Figure 5 (Colour online) 30-year forecasted log death rates (2011–2041) for Danish male population under (left column) LC-H model, (middle column) LCSV model and (right column) LCSV-H model in comparison with LC model. Calibration period: 1835–2010. LC, Lee–Carter; LC-H, LC model with heteroscedasticity; LCSV, LC stochastic volatility; LCSV-H, LCSV model with heteroscedasticity; CI, confidence interval.

The heteroscedasticity structure, from the LC-H model, gives rise to materially larger forecasting intervals for the young and very old age group, while the forecasting interval for the age group 35–39 is narrower than predicted by the LC model. The LCSV model, on the other hand, produces a wider forecasting interval compared to the LC model except for the very old age group. The observed wider forecasting interval is due to the fact that the volatility level is increasing in the last estimation periods and is larger than $$\sigma _{\omega }^{2} $$ estimated in the LC model. Moreover, as the estimated β _x is close to zero at older ages (Figure 2), the impact of the forecasted κ on the prediction of death rates diminished significantly as older ages are considered. The LCSV-H model exhibits similar features of the LC-H and the LCSV model. It is interesting to note that the forecasted means obtained from the different models are very similar and their differences mainly lie in the forecasting interval.

To illustrate further the forecasting property of the LCSV model, we estimate the models for the period 1835–1990 and plot 20-year out-of-sample forecasted log death rates in Figure 6. It turns out the forecasting intervals predicted by the LCSV model tends to be narrower than the LC model, as the estimated $$\sigma _{\omega }^{2} $$ in the LC model is larger than the volatility level at the last estimation period for the LCSV model in this case. Note that the forecasted distributions produced by the LC-H model are biased compared to the benchmark LC model since the fitted rates at the last estimation period, that is year 1990, are different for the LC and LC-H model. This feature is known as jump-off error (Lee & Miller, Reference Lee and Miller2001). One may remove this jump-off bias by forcing the forecasted death rates to start at the actual rates instead of the fitted rates (Bell, Reference Bell1997; Shang et al., Reference Shang, Booth and Hyndman2011). In this paper we do not perform this procedure, however.

Figure 6 (Colour online) 20-year out-of-sample forecasted log death rates for Danish male population under (left column) LC-H model, (middle column) LCSV model and (right column) LCSV-H model in comparison with LC model. Calibration period: 1835–1990. LC, Lee–Carter; LC-H, LC model with heteroscedasticity; LCSV, LC stochastic volatility; LCSV-H, LCSV model with heteroscedasticity; CI, confidence interval.

Figure 7 shows the forecasting distributions of log death rates where we assume a shorter calibration period from 1950 to 1990. For all the models, the estimated β _x for all age groups, except for age groups 0, 1–4 and 5–9, are very close to zero. It is in fact expected since there is no clear downward trend in the observed mortality data besides the first few age groups, during the period 1950–1990. Therefore, there is only small difference between the forecasting distributions produced by the LC model and the LCSV model, except for young age groups. Note that there is a clear change of downward trend for some of the middle age groups for the Danish male mortality data as shown in Figure 7. It results in the out-of-sample data falling out of the lower bound of the 95% credible intervals and its consequence for the forecasting of life expectancy will be discussed in section 6.4.2 and generally in section 6.4.3.

Figure 7 (Colour online) 20-year out-of-sample forecasted log death rates for Danish male population under (left column) LC-H model, (middle column) LCSV model and (right column) LCSV-H model in comparison with LC model. Calibration period: 1950–1990. LC, Lee–Carter; LC-H, LC model with heteroscedasticity; LCSV, LC stochastic volatility; LCSV-H, LCSV model with heteroscedasticity; CI, confidence interval.

By comparing the forecast performance using in-sample data from 1835 to 1990 and from 1950 to 1990 displayed in Figures 6 and 7, we expose the influence that leaving out important historical events, that may affect the mortality rates markedly in a population, can have on the ability to accurately model trend and volatility structures in population dynamics. In particular, we observe that one must be cautious as forecast performance can degrade markedly when important historical events are excluded from the sample as the forecast using data from 1835 to 1990 has clearly outperformed the forecast using only shorter calibration data from 1950 to 1990.

6.4.2. Life expectancy

Using the samples of the forecasted log death rates $$y_{{x,t}}^{{(\ell )}} \, {\equals} \, {\rm ln}\,\hat{m}_{{x,t}}^{{(\ell )}} $$ , where $$\ell \, {\equals} \, 1,\,\ldots\,,L$$ and L is the number of MCMC samples, we can obtain the so-called period life expectancy at different ages by constructing an abridged life table, since we use age group data, as follows (Koissi et al., Reference Koissi, Shapiro and Hognas2006; Yusuf et al., Reference Yusuf, Martins and Swanson2014). We consider age group x∈{0, 1–4, 5–9, … , 95–99} and $$\tilde{x}$$ is defined as the initial age of age group $$x$$ , that is $$\tilde{x}\in\{ 0,1,5,\,\ldots\,,90,95\} $$ . Define $$n_{{\tilde{x}}} $$ as the length of the interval of age group x (corresponds to $$\tilde{x}$$ ) and hence we have n ₀=1, n ₁=4, n ₅=5, … , n ₉₅=5. We then calculate the (crude) death probabilityFootnote ⁷ that a person aged $$\tilde{x}$$ in year t will die in the next $$n_{{\tilde{x}}} $$ years as

(92) $$_{{n_{{\tilde{x}}} }} \hat{q}_{{\tilde{x},t}}^{{(\ell )}} \, {\equals} \, {{n_{{\tilde{x}}} \:\hat{m}_{{x,t}}^{{(\ell )}} } \over {1{\plus}n_{{\tilde{x}}} \left( {1\, {\minus}\, a\left( {\tilde{x},n_{{\tilde{x}}} } \right)} \right)\hat{m}_{{x,t}}^{{(\ell )}} }}$$

where $$a\left( {\tilde{x},n_{{\tilde{x}}} } \right)$$ is the average fraction of the $$n_{{\tilde{x}}} $$ years lived by the people who is initially aged $$\tilde{x}$$ in that interval. Using the assumption that deaths are distributed uniformly in the interval, we set $$a(\tilde{x},n_{{\tilde{x}}} )\, {\equals} \, 0.5$$ for every $$\tilde{x}$$ Footnote ⁸ . The hypothetical number of people alive at age $$\tilde{x}{\plus}n_{{\tilde{x}}} $$ , $$l_{{\tilde{x}{\plus}n_{{\tilde{x}}} ,t}}^{{(\ell )}} $$ , is determined by $$l_{{\tilde{x}{\plus}n_{{\tilde{x}}} ,t}}^{{(\ell )}} \, {\equals} \, l_{{\tilde{x},t}}^{{(\ell )}} \big( {1\, {\minus}\, _{{n_{{\tilde{x}}} }} q_{{\tilde{x},t}}^{{(\ell )}} } \big)$$ , where $$l_{{0,t}}^{{(\ell )}} $$ is assumed to be 100,000. We can then calculate the number of deaths $$_{{n_{{\tilde{x}}} }} d_{{\tilde{x},t}}^{{(\ell )}} \, {\equals} \, l_{{\tilde{x},t}}^{{(\ell )}} \, {\minus}\, l_{{\tilde{x}{\plus}n_{{\tilde{x}}} ,t}}^{{(\ell )}} $$ and the person-years lived $$_{{n_{{\tilde{x}}} }} L_{{\tilde{x},t}}^{{(\ell )}} \, {\equals} \, n_{{\tilde{x}}} \big( {l_{{\tilde{x}{\plus}n_{{\tilde{x}}} ,t}}^{{(\ell )}} {\plus}a(\tilde{x},n_{{\tilde{x}}} ){\times}_{{n_{{\tilde{x}}} }} d_{{\tilde{x},t}}^{{(\ell )}} } \big)$$ . The total future lifetime of the $$l_{{\tilde{x},t}}^{{(\ell )}} $$ persons who attain age $$\tilde{x}$$ is $$T_{{\tilde{x},t}}^{{(\ell )}} \, {\equals} \, \mathop{\sum}\nolimits_{i\geq \tilde{x}} _{{n_{{\tilde{x}}} }} L_{{i,t}}^{{(\ell )}} $$ , where i∈{0, 1, 5, … , 90, 95}. Finally, a sample of the period life expectancy at age $$\tilde{x}$$ is obtained as

(93) $$e_{{\tilde{x},t}}^{{(\ell )}} \, {\equals} \, {{T_{{\tilde{x},t}}^{{(\ell )}} } \mathord{\left/ {\vphantom {{T_{{\tilde{x},t}}^{{(\ell )}} } {l_{{\tilde{x},t}}^{{(\ell )}} }}} \right. \kern-\nulldelimiterspace} {l_{{\tilde{x},t}}^{{(\ell )}} }}$$

and the distributions are obtained in different forecasting year t=T+k, where k≥1.

Remark 6.1 (Period and cohort life expectancy): Period life expectancy assumes there is no trend for future death rates (it is evaluated based on the age-specific death rates in a fixed year t) while cohort life expectancy assumes death rates following the lifetime of a cohort and hence it takes mortality trend into account. For example, to evaluate period life expectancy at age 65 in year t, one needs $$\{ _{5} q_{{65,t}} ,_{5} q_{{70,t}} ,\,\ldots\,,_{5} q_{{95,t}} \} $$ while for cohort life expectancy, $$\{ _{5} q_{{65,t}} ,_{5} q_{{70,t{\plus}5}} ,\,\ldots\,,_{5} q_{{95,t{\plus}30}} \} $$ are used instead. However, the cohort life expectancy for people born in recent years cannot be evaluated using data alone since some of the death rates data are yet to be observed. As a result we focus on period life expectancy so that our forecasts can be compared with the observed data.

Figure 8 shows the forecasted life expectancy at birth, age 65 and age 85 for all the models estimated using data from 1835 to 2010. Interestingly, the forecasted life expectancy at birth is very similar for the LC model and LC-H model for the calibration period 1835–2010. It reflects the fact that forecasting intervals of death rates produced by the LC-H model are wider for some age groups and narrower for other, compared to the LC model, as can be seen from Figure 5. It turns out that these effects almost cancel each other out for this case as death rates are aggregated for all age groups to form the life expectancy at birth distribution. This explanation does not apply to life expectancy at ages 65 and 85, however, since only forecasted death rates for age groups larger than 65 and 85, respectively, are used to obtain the corresponding life expectancy distribution. As the forecasting intervals of death rates generated by the LC-H model tend to be narrower for old age groups compared with the LC model, the interval for the forecasted life expectancy at ages 65 and 85 distribution produced by the LC-H model is observably narrower than the LC model. Note that we only show four different age groups in Figure 5; projected death rates for other age groups are also relevant but are not displayed due to limited space. In contrast to the LC-H model, we can observe from Figure 5 that LCSV model produces generally wider forecasting intervals for death rates at different age groups compared to the LC model. Consequently, the forecasting intervals for life expectancy at various ages generated by the LCSV model is wider than the LC model, as shown in Figure 8. Similarly to the case of death rates forecasting, the LCSV-H model has both the features of the LC-H and LCSV model in terms of life expectancy prediction.

Figure 8 (Colour online) 30-year forecasted life expectancy (2011–2041) at birth, ages 65 and 85 for Danish male population under (left column) LC-H model, (middle column) LCSV model and (right column) LCSV-H model in comparison with LC model. Calibration period: 1835–2010. LC, Lee–Carter; LC-H, LC model with heteroscedasticity; LCSV, LC stochastic volatility; LCSV-H, LCSV model with heteroscedasticity; CI, confidence interval.

Results for forecasted life expectancy using data from 1835 to 1990 are shown in Figure 9. It is clear that the forecasting intervals for life expectancy at different ages produced by the LCSV model are narrower compared to the LC model, which can be traced back to the fact that the LCSV model predicts narrower forecasting intervals for death rates at various age groups than the prediction by the LC model, see Figure 6. Apparently the jump-off bias are quite different for the LC-H (LCSV-H) model and the LC model. Note also that the “cancel out” effect of aggregated death rates (Figure 6) for the life expectancy at birth for the LC-H model is less prominent for the calibration period 1835–1990 than the period 1835–2010.

Figure 9 (Colour online) 20-year out-of-sample forecasted life expectancy (1991–2010) at birth, ages 65 and 85 for Danish male population under (left column) LC-H model, (middle column) LCSV model and (right column) LCSV-H model in comparison with LC model. Calibration period: 1835–1990. LC, Lee–Carter; LC-H, LC model with heteroscedasticity; LCSV, LC stochastic volatility; LCSV-H, LCSV model with heteroscedasticity; CI, confidence interval.

As we use the fitted death rates instead of the observed death rates in the jump-off year (i.e. year 2010), there is a jump-off bias in the forecasted death rates. The forecasted life expectancy at age 65 is particularly sensitive to this jump-off bias. It comes from a sudden decline of death rates for age groups larger than 65 beginning in year 1990, hence a significant increase of life expectancy at age 65 is observed. The jump-off bias is significantly smaller when the calibration period 1835–1990 is considered, see Figure 9. As expected, the forecasted distributions of life expectancy at birth and age 65 are similar for all the models estimated using mortality data from year 1950–1990 (Figure 10). Note that the 95% credible intervals capture poorly the out-of-sample data in this case except for the life expectancy at age 85. It is a consequence of the sudden change of significant downward trend for the death rates observed in the middle age groups of the Danish mortality data starting from around 1990, see Figure 7, as well as the jump-off bias. We discuss about the linear trend assumption and jump-off bias in the next section.

Figure 10 (Colour online) 20-year out-of-sample forecasted life expectancy (1991–2010) at birth, ages 65 and 85 for Danish male population under (left column) LC-H model, (middle column) LCSV model and (right column) LCSV-H model in comparison with LC model. Calibration period: 1950–1990. LC, Lee–Carter; LC-H, LC model with heteroscedasticity; LCSV, LC stochastic volatility; LCSV-H, LCSV model with heteroscedasticity; CI, confidence interval.