2.1 Structure of age/period mortality models

An AP mortality model in discrete time is one which assumes that mortality rates can be modelled as a series of terms involving functions of age, x, and period, t Footnote ¹. In the notation of Hunt & Blake (Reference Hunt and Blake2020b), this can be written as

(1) \begin{align}\eta _{x,t} &= \alpha_x + \sum_{i=1}^N \beta^{(i)}_x \kappa^{(i)}_t\end{align}

where $\eta_{x,t}$ is a link function transforming the raw data, $\alpha_x$ is a static function of ageFootnote ², $\kappa^{(i)}_t$ are N period functions governing the evolution of mortality with time and $\beta^{(i)}_x$ are age functions modulating the impact of this change over the age range. This structure does not include any allowance for the lifelong effects of different birth years (called cohort effects) on mortality.

The structure in equation (1) as it is currently written does not require any of the functions to be known in advance of fitting the model to data. As such, it has what we refer to as a “non-parametric” structure. We consider this as the most general form of an AP mortality model and discuss its identifiability issues in section 4. We will also consider the “parametric” case where $\beta^{(i)}_x$ is a parametric function of age, $\beta^{(i)}_x = f^{(i)}(x;\,\theta^{(i)})$ , in section 6, and models which mix parametric and non-parametric age functions in section 7. Whether $\beta^{(i)}_x$ is parametric or non-parametric will affect the interpretation of the model, as discussed in Hunt & Blake (Reference Hunt and Blake2020b), and also lead to subtly different identification issues.

The form given in equation (1) is widely used and lends itself naturally to interpreting the parameters as measuring either an age or a period feature of mortality rates. Alternatively, when analysing this structure, we may find it useful to consider the static age function, $\alpha_x$ , and the age functions, $\beta^{(i)}_x$ , as being column vectors in $\mathbb{R}^X$ instead of functions of age, and the period functions, $\kappa^{(i)}_t$ , as row vectors in $\mathbb{R}^T$ , rather than a time series. Considering the parameters in this context, it is natural to define inner products, $<.,.>$ on $\mathbb{R}^X$ and $\mathbb{R}^T$ , respectively, and use these to compare the different functions. For instance, we could define the “scale” of an age function by taking

\begin{align*}\| \beta^{(i)}_x \| &= \: < \beta^{(i)}_x, \beta^{(i)}_x >\end{align*}

or the “angle”, $\theta$ , between age functions as

\begin{align*}\cos \theta = \frac{<\beta^{(i)}_x,\beta^{(j)}_x>} {\sqrt{\|\beta^{(i)}_x \|\|\beta^{(j)}_x \|}}\end{align*}

We can think of the inner products being the standard Euclidean inner products, i.e., that $<\beta^{(i)}_x, \beta^{(j)}_x > \: = \: \sum_x \beta^{(i)}_x \beta^{(j)}_x$ and $<\kappa^{(i)}_t, \kappa^{(j)}_t > \: = \: \sum_t \kappa^{(i)}_t \kappa^{(j)}_t$ .

If the period parameters, $\kappa^{(i)}_t$ , are interpreted as random variables, then we also see that this standard Euclidean inner product can be interpreted in terms of their sample mean and sample variance

\begin{align*}\bar{\kappa}^{(i)} = \frac{1}{T}\sum_t \kappa^{(i)}_t &= \frac{1}{T}<\kappa^{(i)}_t, 1 > \\\sigma^2_{\kappa^{(i)}} = \frac{1}{T} \sum_t \left( \kappa^{(i)}_t - \bar{\kappa}^{(i)}\right)^2 &= \frac{1}{T}<\kappa^{(i)}_t-\bar{\kappa}^{(i)}, \kappa^{(i)}_t-\bar{\kappa}^{(i)} > \\&= \frac{1}{T} \| \kappa^{(i)}_t-\bar{\kappa}^{(i)} \|\end{align*}

Similarly, we see that the sample correlation between two period functions is given by the angle between them

\begin{align*}\text{Corr}(\kappa^{(i)}_,\kappa^{(j)}_t) &= \frac{\sum_t \left(\kappa^{(i)}_t - \bar{\kappa}^{(i)}\right) \left(\kappa^{(j)}_t - \bar{\kappa}^{(j)}\right)}{\sqrt{\sigma^2_{\kappa^{(i)}} \sigma^2_{\kappa^{(j)}} }} \\&=\frac{<\kappa^{(i)}_t-\bar{\kappa}^{(i)}, \kappa^{(j)}_t-\bar{\kappa}^{(j)}>} {\sqrt{\|\kappa^{(i)}_t-\bar{\kappa}^{(i)} \|\|\kappa^{(j)}_t-\bar{\kappa}^{(j)} \|}} \\&= \cos \theta_{\kappa-\bar{\kappa}}\end{align*}

Consequently, the standard Euclidean inner product has a number of helpful interpretations and is widely usedFootnote ³. However, we could equally reasonably choose other inner products on $\mathbb{R}^X$ and $\mathbb{R}^T$ if these are more convenientFootnote ⁴.

When projecting the period functions using multivariate time series processes, it is helpful to define vectors

\begin{align*}\boldsymbol{\kappa}_t &= \begin{pmatrix} \kappa^{(1)}_t, & \ldots & \kappa^{(N)}_t \end{pmatrix}^\top \\\boldsymbol{\beta}_x &= \begin{pmatrix} \beta^{(1)}_x, & \ldots & \beta^{(N)}_x \end{pmatrix}^\top\end{align*}

The model, therefore, has the vector structure

(2) \begin{align}\eta_{x,t} &= \alpha_x + \boldsymbol{\beta}^\top_x \boldsymbol{\kappa}_t \end{align}

In order to project the model, the vector $\boldsymbol{\kappa}_t$ can be modelled using VARIMA processes. This is considered further in section 9.

We can also construct matrices for the age and period functions as $\beta = \{\beta^{(1)}_x \; \beta^{(2)}_x \ldots \beta^{(N)}_x \}$ and $\kappa = \{ \kappa^{(1)}_t; \kappa^{(2)}_t; \ldots; \kappa^{(N)}_t \}$ and therefore re-write equation (1) in matrix form

(3) \begin{align}H = \alpha \texttt{1}^\top + \beta \kappa\end{align}

where

• H is the $(X \times T)$ matrix of transformed data (i.e. $H = \{ \eta_{x,t} \}$ ),

• $\alpha$ is a $(X \times 1)$ matrix of the static age function,

• $\texttt{1}$ is a $(T \times 1)$ matrix of ones, and

• $\beta$ and $\kappa$ are the $(X \times N)$ matrix and $(N \times T)$ matrix of age and period functions constructed above, respectively.

When expressed in this form, AP models can be analysed through the prism of matrix algebra and linear mathematics. Specifically, we can see that an AP mortality model is a mapping, $\Theta$ , from the space of parameters to the model space, $\mathcal{M}$ , of fitted mortality rates.

(4) \begin{align}\Theta(\alpha_x, \beta^{(i)}_x, \kappa^{(i)}_t) \,:\, \mathbb{R}^X \times \mathbb{R}^{NX} \times \mathbb{R}^{NT} \to \mathcal{M} \subset \mathbb{R}^{X \times T}\end{align}

Analysing AP mortality models as linear transformations can be very useful and is pursued in sections 2.2 and 4 and in Appendix B. However, whilst such an abstraction can be useful for some purposes, it is important to remember that the parameters in the model have specific interpretations, for instance, that the period functions are ordered chronologically, and so the problem of identifiability should not be seen purely as an exercise in linear mathematics.

2.2 Identifiability in age/period models

An AP mortality model cannot, in general, be estimated as it stands. This is because any parameter estimates would not be unique, since equation (3) is not, in general, fully identifiable.

A model is fully identified when all the parameters in it can be uniquely determined by reference to the available data. In contrast, most mortality models are not fully identified – there exist different sets of parameters which will give the same fitted mortality rates and consequently the same goodness of fit. Although this phenomenon is not unique to mortality models, it is very widespread in mortality modelling and has significant implications when we come to project these models.

The models are not fully identifiable because the space of the parameters for the model, $\mathbb{R}^X \times \mathbb{R}^{NX} \times \mathbb{R}^{NT}$ has a higher dimension than that of the model space, $\mathcal{M}$ , as we show later. Therefore, the mapping $\Theta$ in equation (4) cannot be injectiveFootnote ⁵, since we cannot find a one-to-one mapping from a higher dimension space to a lower one. In practice, this means that we can find transformations of the parameters

(5) \begin{align}\{ \alpha_x, \beta^{(i)}_x, \kappa^{(i)}_t \} \to \{ \hat{\alpha}_x, \hat{\beta}^{(i)}_x, \hat{\kappa}^{(i)}_t \}\end{align}

such that

(6) \begin{align}\Theta(\alpha_x, \beta^{(i)}_x, \kappa^{(i)}_t) &= \Theta(\hat{\alpha}_x, \hat{\beta}^{(i)}_x, \hat{\kappa}^{(i)}_t) \end{align}

We call the transformations of the parameters which satisfy equation (6) “invariant”, because the fitted mortality rates do not change when they are applied to the parameters. The additional degrees of freedom in these invariant transformations correspond to the additional dimensions of the parameter space compared with the model space.

Because $\{ \alpha_x, \beta^{(i)}_x, \kappa^{(i)}_t \}$ and $\{ \hat{\alpha}_x, \hat{\beta}^{(i)}_x, \hat{\kappa}^{(i)}_t \}$ give identical fitted mortality rates and therefore fit observed data equally well, there is no statistical reason to choose between them. In practice, in order to specify a unique set of parameters, constraints independent of the data are imposed – so called “identifiability constraints”. This has the effect of reducing the number of degrees of freedom from the number of parameters. Mathematically, imposing constraints restricts the original parameter space, $\mathbb{R}^X \times \mathbb{R}^{NX} \times \mathbb{R}^{NT}$ , to a subspace, $\mathcal{P}$ , which has fewer dimensions. The aim is to select a subspace, $\mathcal{P}$ , which has the same dimension as the model space, $\mathcal{M}$ , which allows for a one-to-one mapping between the reduced parameter space and the model space. Reducing the dimension of the parameter space can also be achieved by reparameterising the model in a “maximally invariant” form, as discussed in Appendix B.

It is important to know the number of dimensions of the model space, not only to ensure that our model is uniquely estimated but also because this value is used to penalise the likelihood or deviance functions in measures of the goodness of fit, such as the Bayes Information Criterion. A failure to correctly determine the number of free parameters in a model may therefore distort tests of the goodness of fit, such as those performed in Cairns et al. (Reference Cairns, Blake, Dowd, Coughlan, Epstein, Ong and Balevich2009) and Haberman & Renshaw (Reference Haberman and Renshaw2011), and potentially leads to an incorrect assessment about which model gives a superior fit to data. One specific example of this is discussed in Appendix A.

5.1 Location

Because the location of the period functions is not identifiable, Renshaw & Haberman (Reference Renshaw and Haberman2003b) set their level by imposing $\sum_t \kappa^{(i)}_t = 0$ for $i = 1,2$ . As with the LC model, this gives the static age function the demographic significance of representing “average” mortality rates across the period range of the data. This does not cause any additional issues for the LC2 model, so long as it is imposed via an identifiability constraint on $\kappa_t$ and not by imposing the form of $\alpha_x$ (as discussed in Appendix A).

5.2 Scale

To set the scale of the AP terms, Renshaw & Haberman (Reference Renshaw and Haberman2003b) imposed the constraint $\sum_x \beta^{(i)}_x\,{=}\,1$ for $i = 1,2$ , again, in order to be consistent with the convention established by Lee & Carter (Reference Lee and Carter1992). However, the justification for this normalisation scheme makes most sense under the assumption that $\beta^{(i)}_x \geq 0$ for all x – indeed, this is imposed on the LC model in Haberman & Renshaw (Reference Haberman and Renshaw2009) at the expense of goodness of fit to the data. If $\beta^{(i)}_x \geq 0$ , then $\sum_x \beta^{(i)}_x = 1$ constrains the age function to be in the range [0, 1]. The values of $\beta^{(i)}_x$ , therefore, can be felt to represent a proportion of the factor $\kappa^{(i)}_t$ impacting mortality at age x. In general, however, it may be the case that $\beta^{(i)}_x < 0$ at some ages, especially in models with multiple AP terms. If so, the interpretation of the age functions as measuring the proportion of the change is no longer applicable.

Figure 1 shows the age functions from the LC2 model fitted to data for men in the UKFootnote ¹² with the constraint $\sum_x \beta^{(i)}_x = 1$ for $i=1,2$ . We see that if $\beta^{(i)}_x \leq 0$ for some x, as is the case for the second age function, then the identifiability constraint on the age function no longer limits it to a particular range of values. Indeed, $\beta^{(i)}_{x_1}$ can take arbitrarily high positive values, as long as there exists a correspondingly low $\beta^{(i)}_{x_2}$ (i.e. large magnitude and negative) to compensate. This is in contrast to $\beta^{(1)}_x$ , which is greater than zero for all ages, and hence is comparatively close to zero across the whole age rangeFootnote ¹³. This undermines the rationale for selecting a common normalisation scheme for the age functions, which was to aid comparisons of the relative importance of the different AP terms.

Figure 1

LC2 age functions with $\sum_x \beta^{(i)}_x = 1.$

The identifiability constraint $\sum_x \beta^{(i)}_x = 1$ can also, theoretically, lead to numerical problems when fitting the model to data. In practice, the constraint is imposed by taking the set of parameters generated by the fitting algorithm (which do not have any identifiability constraints imposed) and using the transformation in equation (8) with $b = \sum_x \beta^{(i)}_x$ , i.e., $\hat{\beta}^{(i)}_x = \frac{1}{\sum_\xi \beta^{(i)}_\xi }\beta^{(i)}_x$ . This gives an equivalent set of parameters (with the same fit to the data), but where $\sum_x \hat{\beta}^{(i)}_x = 1$ by construction. If, however, $\sum_x \beta^{(i)}_x = 0$ for whatever reason, this procedure will fail as applying the transformation involves dividing by zero, even if the age function fitted originally by the algorithm is reasonable. While this is unlikely, it is far more common that we find $\sum_x \beta^{(i)}_x \approx 0$ , which will then lead to the revised parameters (with the constraint imposed) being infeasibly large, and which may, in turn, generate problems with the fitting algorithm.

Both of these problems with the normalisation scheme are caused because simple summation over x is not a true norm. A true norm, $\|v\|$ , for a vector space, $\mathcal{V}$ , of a vector, v, is defined by the properties

1. 1. $\|v \| \geq 0 \; \forall v \in \mathcal{V}$ ;

2. 2. $\|v \| = 0 \iff v = 0$ ;

3. 3. $\| a v \| = |a| \| v \| \; \forall a \in \mathbb{R}$ ; and

4. 4. $\| v_1 + v_2 \| \leq \|v_1\| + \| v_2 \|$ .

These properties mean that we can use a true norm to define distances and scales within the vector space and therefore make them useful when specifying a normalisation scheme. However, since $\sum_x \beta^{(i)}_x$ is not a true norm in $\mathbb{R}^X$ , we are not able to use this normalisation scheme to compare scales for the age functions, and cannot assume that $\sum_x \beta^{(i)}_x >0$ in our fitting algorithms when we come to impose the identifiability constraints.

Normalisation schemes using true norms on $\mathbb{R}^X$ , such as $\sum_x |\beta^{(i)}_x| = 1$ or $\sum_x (\beta^{(i)}_x)^2 = 1$ , will not suffer from these issues. When it comes to normalising the fitted age function, a procedure using a true norm for the normalisation scheme will never involve division by zero if the transformation in equation (8) is used with any non-trivial age functions. Therefore, in most circumstances, normalisation schemes based on true norms will be preferableFootnote ¹⁴.

However, we note that normalisation schemes based on true norms are not perfectly identified, since the transformation

(14) \begin{align}\{\hat{\beta}^{(i)}_x, \hat{\kappa}^{(i)}_t \} = \{-\beta^{(i)}_x, -\kappa^{(i)}_t\} \end{align}

is an invariant transformation of the parameters where the new parameters still satisfy the identifiability constraints. In principle, we could solve this by choosing alternative sets of normalisation constraints, for instance

\begin{align*}sign\left(\sum_x \beta^{(i)}_x \right) \sum_x \left( \beta^{(i)}_x \right)^2 = 1\end{align*}

which are still based on using true norms but are not invariant to changing the sign of the age function. However, the specific transformation causing this problem has few practical consequences when fitting the model, since the transformation is not continuous. When fitting the LC or LC2 models using maximum likelihood techniques, for instance, we make small adjustments to the parameters at each iteration and so it is not possible to move smoothly from one set of acceptable parameters to another when fitting the model. In addition, the transformation in equation (14) can be applied to any set of parameters after fitting the model and, hence, can be used to select the sign of the age function based on the judgement of the user when reviewing the fitted parameters.

To illustrate this, consider the age functions shown in Figure 2 which fit the LC2 model to the same data as in Figure 1 with the normalisation scheme $\sum_x |\beta^{(i)}_x| = 1$ . This normalisation scheme gives a model with exactly the same fit to the data, but the estimated parameters for the age and period functions are now of the same order of magnitudeFootnote ¹⁵, which may make this model easier to project. We also avoid the possibility of any computational problems when imposing the identifiability constraint, since the divisor, $\sum_x |\beta^{(i)}_x|$ , will not be zero for any non-trivial age function.

Figure 2

LC2 age functions with $\sum_x |\beta^{(i)}_x| = 1.$

5.3 Rotation

We established in section 4 that $N(N+1)$ constraints were necessary to restrict the parameters in a general AP mortality model with non-parametric age functions, due to the number of free parameters in the transformations in equations (11) and (12). In the context of the LC2 model, this means that we would require six identifiability constraints. However, only four identifiability constraints (two on the level of the two period functions, two on the normalisation of the two age functions) were described in Renshaw & Haberman (Reference Renshaw and Haberman2003b). We, therefore, have an additional two invariant transformations of the parameters which give the same fit to data and which satisfy the constraints already explicitly imposed by Renshaw & Haberman (Reference Renshaw and Haberman2003b). These can be written as

(15) \begin{align}\begin{pmatrix} \hat{\beta}^{(1)}_x \\[4pt] \hat{\beta}^{(2)}_x \end{pmatrix} &= \left(\begin{array}{c@{\quad}c} \theta & 1-\theta \\[4pt] 0 & 1 \end{array}\right) \begin{pmatrix} \beta^{(1)}_x \\[4pt] \beta^{(2)}_x \end{pmatrix} \notag \\[4pt]\begin{pmatrix} \hat{\kappa}^{(1)}_t \\[4pt] \hat{\kappa}^{(2)}_t \end{pmatrix} &= \frac{1}{\theta} \left(\begin{array}{c@{\quad}c} 1 & \theta-1 \\[4pt] 0 & \theta \end{array}\right)\begin{pmatrix} \kappa^{(1)}_t \\[4pt] \kappa^{(2)}_t \end{pmatrix} \end{align}

and

(16) \begin{align}\begin{pmatrix} \hat{\beta}^{(1)}_x \\[4pt] \hat{\beta}^{(2)}_x \end{pmatrix} &= \left(\begin{array}{c@{\quad}c} 1 & 0 \\[4pt] 1-\phi & \phi \end{array}\right) \begin{pmatrix} \beta^{(1)}_x \\[4pt] \beta^{(2)}_x \end{pmatrix} \notag \\[4pt]\begin{pmatrix} \hat{\kappa}^{(1)}_t \\[4pt] \hat{\kappa}^{(2)}_t \end{pmatrix} &= \frac{1}{\phi} \left(\begin{array}{c@{\quad}c} \phi & 0 \\[4pt] \phi-1 & 1 \end{array}\right)\begin{pmatrix} \kappa^{(1)}_t \\[4pt] \kappa^{(2)}_t \end{pmatrix}\\[-15pt]\nonumber \end{align}

These transformations can be thought of as “rotations” of the AP functions, because they change the angle between age and period functions, but the normalisation scheme $\sum_x \hat{\beta}^{(i)}_x = 1$ still holds. They also clearly illustrate that we have an additional two degrees of freedom, given by the free parameters $\theta$ and $\phi$ , which do not change the fitted mortality rates but which should be used to impose two more identifiability constraints on the model.

This does not necessarily mean that the model in Renshaw & Haberman (Reference Renshaw and Haberman2003b) was poorly identified, however. Although the authors did not explicitly acknowledge the existence of these additional identifiability constraints, their use of singular value decomposition to fit the model imposed them implicitly. By taking singular values (or equivalently, principal components), age and period functions are selected so that $\sum_t \kappa^{(i)}_t \kappa^{(j)}_t = 0$ and $\sum_x \beta^{(i)}_x \beta^{(j)}_x = 0$ for $i \neq j$ . We call such age and period functions “orthogonal” to each other as the angle between them defined earlier using the standard inner product will be $\frac{\pi}{2}$ . This implicit imposition of additional identifiability constraints leads to a fully identified model.

If alternative fitting methods are used, such as maximum likelihood (e.g. in Brouhns et al. (Reference Brouhns, Denuit and Vermunt2002a)) or minimal deviance (e.g. in Renshaw & Haberman (Reference Renshaw and Haberman2003a)), then these constraints must be imposed explicitly in order to obtain a fully identified model. To impose these orthogonality constraints for a general set of LC2 parameters, we would therefore need to solve $\sum_t \hat{\kappa}^{(i)}_t \hat{\kappa}^{(j)}_t = 0$ and $\sum_x \hat{\beta}^{(i)}_x \hat{\beta}^{(j)}_x = 0$ with the transformed parameters defined by equations (15) and (16) in order to find $\theta$ and $\phi$ .

We also note the special case where $A= \left(\begin{array}{c@{\quad}c} 0 & 1 \\[1pt] 1 & 0 \end{array}\right)$ (i.e. $\theta = \phi = 1$ when equations (15) and (16) are composed), which relates to the transformation

(17) \begin{align}\{ \hat{\beta}^{(1)}_x, \hat{\kappa}^{(1)}_t, \hat{\beta}^{(2)}_x, \hat{\kappa}^{(2)}_t \} = \{ \beta^{(2)}_x, \kappa^{(2)}_t, \beta^{(1)}_x, \kappa^{(1)}_t \} \end{align}

This is an invariant transformation of the parameters where the new parameters still satisfy the identifiability constraints. However, it amounts to simply re-labelling the AP terms and arises because the identifiability constraints are the same for all AP terms. Similar to the case in equation (14), this situation could, in principle, be solved by using different identifiability constraints for the different AP terms, for instance

\begin{align*}\sum_x | \beta^{(1)}_x | = 1 \\\sum_x \left( \beta^{(2)}_x \right)^2 = 1\end{align*}

which breaks the symmetry between the different AP terms and, thus, prevents them being relabelled. However, as with equation (14), the transformation in equation (17) has few practical consequences, since it is not continuous and so it is not possible to move smoothly from one set of acceptable parameters to another when fitting the model. Furthermore, using different identifiability constraints for the different AP terms conflicts with a desire for their scale to be comparable with each other and, hence, we do not believe that this issue is important in practice.

If maximum likelihood methods are used to estimate the parameters in a model, it is useful that these estimators are independent of each other. This helps to give more efficient fitting algorithms for estimation and is also useful when allowing for parameter uncertainty using the technique of Brouhns et al. (Reference Brouhns, Denuit and Vermunt2002b) discussed in section 8. Assuming the canonical link function is used as discussed in Hunt & Blake (Reference Hunt and Blake2020b), the independence of the estimators can be assessed by consideration of the information matrix for the different parameters

\begin{align*}I(\beta^{(i)}_x, \beta^{(j)}_x) &= \mathbb{E}\left[\frac{\partial^2 \mathcal{L}}{\partial \beta^{(i)}_x \partial \beta^{(j)}_x} \right] \\[3pt]&= - \sum_t \mathbb{V}ar(D_{x,t}) \kappa^{(i)}_t \kappa^{(j)}_t \\[3pt]I(\kappa^{(i)}_t, \kappa^{(j)}_t) &= \mathbb{E}\left[\frac{\partial^2 \mathcal{L}}{\partial \kappa^{(i)}_t \partial \kappa^{(j)}_t} \right] \\[3pt]&= - \sum_x \mathbb{V}ar(D_{x,t}) \beta^{(i)}_x \beta^{(j)}_x\end{align*}

Therefore, we see that orthogonal age and period functions are independent of each other if $\mathbb{V}ar(D_{x,t})$ is constant across ages and years. This assumption is implicitly made when using singular value decomposition or principal components analysis to estimate parameters. However, the assumption is not consistent with the use of the Poisson or binomial distribution for death counts, as discussed in Hunt & Blake (Reference Hunt and Blake2020b). Under these distributions, the variance of death counts depends upon the exposure to risk at different ages, which changes considerably over different ages and years and is more realistic in practice.

In principle, we could impose independent parameter estimates using the transformations in equations (15) and (16) with carefully selected values of $\theta$ and $\phi$ to obtain an equivalent set of parameters. Doing so would simply be choosing an alternative (but equally valid) set of identifiability constraints. However, in practice, this would mean constraints that are both more difficult to impose than the traditional orthogonality constraints using the Euclidean inner product, and which lose the connection between the inner product and the sample moments of $\kappa^{(i)}_t$ . In practice, imposing $\sum_t \kappa^{(i)}_t \kappa^{(j)}_t = 0$ and $\sum_x \beta^{(i)}_x \beta^{(j)}_x = 0$ for $i \neq j$ to obtain orthogonal age and period functions is a convenient and useful set of identifiability constraints.

Whichever set of constraints is imposed on the angles between different period functions, the most important thing is, however, to impose some form of constraint. A failure to do so may result in the fitting routine failing to converge or, alternatively, the fitting routine may give model parameters which depend upon the initial parameter estimates used in the algorithm. Similarly, the angles between different age functions must also be constrained in order to fully identify the model. This has implications for estimated parameter uncertainty, as discussed in section 8.

We noted in section 2 that the correlation between two different period functions depends on the angle between them. This means that we see that the correlations we find between period functions from our fitted parameters depends only on the identifiability constraints chosen, and so are not meaningful. For instance, the constraint $\sum_t \kappa^{(i)}_t \kappa^{(j)}_t = 0$ imposes independence on the period functions over the historical range of the data when they are considered as time series. Figure 3 shows period functions for the LC2 model fitted to the same data as above, but with two different constraints on the angles between them. In Figure 3(a), the period functions are orthogonal whereas, in Figure 3(b), they have a correlation of 75%Footnote ¹⁶. However, both sets of parameters give identical fits to the historical data. This will have important consequences when we come to project the model in section 9.

Figure 3

Period functions from the LC2 model.

In situations such as Renshaw & Haberman (Reference Renshaw and Haberman2003b), where orthogonality constraints on the AP terms have been imposed implicitly by the fitting mechanism, we believe that it is important to recognise and state them clearly. Not only will this clarify which features of graphs of the age and period terms are meaningful but it also ensures that we assess the dimension of $\mathcal{P}$ (i.e. the number of degrees of freedom in the model) correctly. This is important when assessing the goodness of fit for the model.

As an example of this, in Haberman & Renshaw (Reference Haberman and Renshaw2011), the LC2 model is compared against other mortality models using various measures including the Akaike Information Criterion, Bayes Information Criterion and Hannan-Quinn Criterion. All of these measures use the number of degrees of freedom (i.e. $dim(\mathcal{P}$ )) of the model to penalise the log-likelihood. By failing to explicitly state the orthogonality constraints placed on the AP terms in the LC2 model and, therefore, failing to include them in the count of restrictions placed upon the model parameters, the study overestimates the number of degrees of freedom in the model. This excessively penalises the LC2 model relative to its comparators.

Using the invariant transformations to impose orthogonality on the age and period functions generalises naturally to more complicated models with $N > 2$ . Identifiability in-sample in a model with non-parametric AP terms is therefore not problematic if fitting methods based on singular value decomposition or principal component analysis are used (except for setting the locations of the $\kappa^{(i)}_t$ and the scale for the $\beta^{(i)}_x$ by imposing an appropriate normalisation scheme).

6.1 Location

We noted in section 4 that the transformation in equation (12) does not change the form of the age functions. Accordingly, it can still be applied to change the levels of the period parameters in exactly the same manner as described in section 4, whilst leaving the fitted mortality rates and the functional forms of the age functions unchanged. The period functions in models with parametric age functions therefore still have unidentified locations, and so we still need to impose levels on the period parameters in exactly the same manner as we did in section 5. Most users of such models impose $\sum_t \kappa^{(i)}_t = 0$ , consistent with the choice made for models with non-parametric age functions and with a similar interpretation. However, for models which have a specific form of the static age function imposed a priori, this is not necessary, as discussed in Appendix A.

6.2 Scale

We see that the transformation in equation (11) takes linear combinations of the old age and period functions in order to create new AP terms. Therefore, these transformations will change the form of the age functions in a model with parametric age functions. Since the form of the age functions defines the model being used, the transformations in equation (11) cannot be used in models with parametric age functions.

In section 5, we saw that these transformations were useful in models with non-parametric age functions when it came to imposing a normalisation scheme on the age functions and orthogonalising them with respect to each other. This was beneficial as it enabled comparability and near-independence between different AP terms. It is therefore desirable to also achieve the same properties for models with parametric age functions.

We also see that although using the transformations in equation (11) in models with parametric age functions gives different age functions (and therefore different models), they do not affect the fitted mortality rates: all the models obtained by using these transformations are equivalent in the sense defined above. It therefore makes sense to choose, from the set of models equivalent to the one we are interested in, a model with age functions which have the desirable properties of possessing a standard normalisation scheme and being orthogonal. We discuss how this can be done in this section and section 6.3, respectively.

Most mortality models with parametric age functions have the age functions defined in their simplest and most natural form. However, choosing definitions for their simplicity rather than for desirable statistical properties, such as having a common normalisation scheme, can lead to issues when comparing age and period terms within the same model and between different models. We show this below for the CBDX and exponential models in equations (18) and (19), respectively. However, for each of these models, we also show how this issue can be resolved by using alternative definitions of the age functions to give models which have far more comparable age and period terms.

First, let us consider how a common normalisation scheme for the age functions can be achieved in the CBDX model in equation (18). In the LC2 model, Renshaw & Haberman (Reference Renshaw and Haberman2003b) imposed the normalisation scheme $\sum_x \beta^{(i)}_x = 1$ on the age functions in the model, using the transformations in equation (11). In contrast, the age functions in equation (18) already have defined scales, i.e., $\sum_x\, f^{(1)}(x) = \sum_x 1 = X$ and $\sum_x\, f^{(2)}(x) = \sum_x (x-\bar{x}) = 0$ .

However, these defined scales cause problems when it comes to comparing the AP terms. The most important of these issues is that the scale of $f^{(2)}(x)$ is zero, which is not sensible for a functions which is not identically equal to zero. This is a consequence of using a normalisation scheme which is not based on using a true norm. In section 5, we saw that a more sensible choice of normalisation scheme was to use $\sum_x |\beta^{(i)}_x|$ to define the scales of the age functions. Using this for the CBDX model, we find $\sum_x |\,f^{(1)}(x)| = \sum_x 1 = X$ and $\sum_x |\,f^{(2)}(x) |= \sum_x |(x-\bar{x})| = 0.25X^2$ if X is even or $0.25(X-1)(X+1)$ if X is odd.

However, this fails to resolve the second problem, which is that different scales are defined for each of the AP terms, i.e., the scale of the first age function is proportional to the number of ages, X, whilst the scale of the second in proportional to $X^2$ . This makes comparisons difficult, both between the CBDX and LC2 models and between the first and second AP terms within the CBDX model. The differing scales of the corresponding period functions can also lead to numerical problems when we try to project them using multivariate methods, as discussed in section 9.

To ensure that the age functions have the same scale, we need to define a model equivalent to that in equation (18) where the age functions have this property. Trivially, we see that the model

(20) \begin{align}\eta_{x,t} &= \alpha_x + \frac{1}{X}\kappa^{(1)}_t + \frac{4(x-\bar{x})}{X^2} \kappa^{(2)}_t \end{align}

(assuming X is even) is equivalent to the model in equation (18)Footnote ¹⁸. All that differs between the models in equations (18) and (20) is the precise definition of the age functions, although the age functions in both models have the same functional form (i.e. a constant and a linear function of age, x). In addition, we see that in the model in equation (20), $\sum_x |\,f^{(i)}(x) |= 1$ for both age functions. In particular, this has the advantage of greater comparability between the AP terms.

To illustrate the impact of ensuring that the age functions have a common normalisation scheme, Figure 4 shows the period functions from the two CBDX models in equations (18) and (20), fitted to the same data as used for the LC2 model in section 5, with both the original and the revised normalisation schemes. We see that the magnitude of the different period functions fitted with the original model in equation (18) differs enormouslyFootnote ¹⁹. This can be a problem as most numerical algorithms for analysing time series are optimised to work best on series of comparable orders of magnitude. In contrast, the revised CBDX model in equation (20) gives period functions of comparable magnitudeFootnote ²⁰. The common scale also means that it is easier to compare these period functions with those in Figure 3 from the LC2 model.

Figure 4

Period functions from the CBDX model.

Turning now to the exponential model in equation (19), we find similar issues for the normalisation scheme of the age functions. In the exponential model, $\sum_x |\,f^{(1)}(x)| = X$ as before for the CBDX model, which can be dealt with in exactly the same manner. In addition, $\sum_x |\,f^{(2)}(x;\,\lambda)| = \sum_x e^{-\lambda x} = \frac{e^{-\lambda}(1-e^{-\lambda (X+1)})}{ 1-e^{-\lambda}} \approx \frac{e^{-\lambda}}{1-e^{-\lambda}}$ for the second age function. Not only will this be different from the scale of the first AP term, but the scale is a function of the free parameter $\lambda$ . Since $\lambda$ varies during the fitting process, this will alter the scale of $f^{(2)}(x;\,\lambda)$ . Hence, $\lambda$ will be trying to fulfil two purposes simultaneously: first, describing the shape of the age function and second, determining its scale, i.e., the relative importance of the AP term. This confusion of different purposes can cause numerical instability in most fitting algorithms, which may be one reason why age functions with free parameters have not been commonly used in practice.

For the CBDX model, we obtained a common normalisation scheme for the age functions by choosing slightly different definitions for the age functions, i.e., we defined alternative age functions which were equal to the original ones, but rescaled by $\sum_x |\,f^{(i)}(x)|$ . For the exponential model, we do the same thing, to obtain

(21) \begin{align}\eta_{x,t} &= \alpha_x + \frac{1}{X}\kappa^{(1)}_t + \frac{ 1-e^{-\lambda}}{e^{-\lambda}(1-e^{-\lambda (X+1)})} e^{-\lambda x} \kappa^{(2)}_t \end{align}

The only difference in this case is that the second age function is rescaled by a function of the free parameter, $\lambda$ , rather than a constant in the case of the CBDX model. Again, we see that the age functions have the same functional forms (a constant and an exponential function of age) as before, but with the normalisation scheme $\sum_x |\,f^{(2)}(x;\,\lambda)| = 1 \; \forall \lambda$ as $\lambda$ is varied when fitting the model. This contrasts with the model in equation (19) and ensures that both age functions have the same normalisation scheme and so are more comparable.

We call age functions such as the revised $f^{(2)}(x;\,\lambda)$ in equation (21) “self-normalising”, as they have the property that our desired normalisation scheme is imposed automatically for all values of the free parameters in the age function (i.e. $\sum_x |\,f^{(i)}(x;\, \theta^{(i)})| = 1 \; \forall \theta^{(i)}$ ). Self-normalisation is an important and useful property. Most importantly, the common normalisation scheme allows for comparability between different age functions (potentially with very different functional forms) in a model, independent of their shape. Furthermore, by allowing the value of the free parameter to describe the shape of the age function, without impacting the scale of the AP term, we find that self-normalising age functions are considerably more robust (in the sense of being likely to converge) and stable to small changes in the data. For this reason, the age functions used in the “toolkit” in the Appendix of Hunt & Blake (Reference Hunt and Blake2014) are all self-normalising with respect to the normalisation scheme $|\,f^{(i)}(x;\,\theta^{(i)})|\,{=}\,1$ Footnote ²¹. However, the trade-off is that the numerical routines are significantly more complicated to implement and may need to be written specially for the specific circumstances, rather than adapted from “off-the-shelf” statistical packagesFootnote ²².

In summary, we see that, when the age functions in a mortality model are defined parametrically, a common normalisation scheme for all of them can be achieved by defining the age functions carefully. For more sophisticated age functions involving free parameters estimated from the data, this means defining age functions which are self-normalising, so that the normalisation scheme holds for all values of these parameters as they are varied during the fitting procedure.

6.3 Rotation

In section 6.2, we saw that for models with parametric age functions, we could ensure that the age functions had the same normalisation scheme by carefully defining them to have this property when we specified the model. The same is also true if we want our age functions to be orthogonal to each other.

Again, similar to section 6.2, we start from the fact that most mortality models have their age functions defined in the simplest form, such as in equations (18) and (19). These simple forms are not, necessarily, orthogonal. However, we can define equivalent models where the age functions are orthogonal. Unlike the case of ensuring a common normalisation scheme, however, we will see that orthogonality between age functions is not always a desirable property and may conflict with other desirable properties, such as the terms in the model having distinct demographic significance. Therefore, the choice of whether to define orthogonal age functions or not will depend upon the model in question and the aims of the model user.

For example, consider the CBDX model of equation (18) before normalisation. The model already has orthogonal age functions, since $\sum_x\, f^{(1)}(x) f^{(2)}(x) = \sum_x (x-\bar{x}) = 0$ . However, we could also consider an equivalent model, with simpler definitions of the age functions of the form

(22) \begin{align}\eta_{x,t} &= \alpha_x + \kappa^{(1)}_t + x \kappa^{(2)}_t \end{align}

This model is more similar to the form of the original CBD model proposed in Cairns et al. (Reference Cairns, Blake and Dowd2006a). However, we observe that the age functions are not orthogonal, i.e., $\sum_x\, f^{(1)}(x) f^{(2)}(x) = \sum_x x= \frac{1}{2}X(X+1)$ . It is easy to see that models in equations (18) and (22) are equivalent, in that they give the same fitted mortality rates and are linked through a transformation of the form in equation (11). The form of the age functions in equation (18) was introduced in Cairns et al. (Reference Cairns, Blake, Dowd, Coughlan, Epstein, Ong and Balevich2009) and, in practice, has proved far more popular than the simpler age functions in equation (22), in part because it is more robust to fit to data due to the parameter estimates for the period functions being nearly independent of each other. Consequently, we see that defining orthogonal age functions can be desirable, even if it comes at the expense of a slightly more complicated definition of the age functions.

The age functions in the CBDX model are of constant and linear form, i.e., polynomials of order zero and one, respectively. Defining orthogonal age functions, as in equation (18), has not changed this form, merely selected the first two members of the orthogonal family of polynomials, i.e., the Legendre polynomialsFootnote ²³. The orthogonal age functions in equation (18) have the same demographic significance as the simpler age functions in equation (22), but the additional desirable property of orthogonality. Generalising this, we see that choosing orthogonal age functions does not change their form and hence does not affect their demographic significance when the age functions come from the same functional family (e.g. polynomials).

However, this is not the case when the age functions come from different functional families. We see this by considering the exponential model once more. To define orthogonal age functions for this model, we could select a model equivalent to that in equation (19) with orthogonal age functions, namely,

\begin{align*}f^{(2)}(x;\,\lambda) &= e^{-\lambda x} - \frac{e^{-\lambda}(1-e^{-\lambda (X+1)})}{ 1-e^{-\lambda}}\end{align*}

We see that the age functions in this model are orthogonal as $\sum_x\, f^{(1)}(x) f^{(2)}(x;\,\lambda) = 0 \; \forall \lambda$ . This revised model is equivalent to that in equation (19), as it gives the same fitted mortality rates and the two models are linked by a transformation of the form of that in equation (11).

However, it is likely that we originally selected an exponential function for its demographic significance (e.g. a mortality effect which decreases rapidly with age, such as that associated with the relatively high rate of infant mortality). The redefined $f^{(2)}(x;\,\lambda)$ will not possess this demographic significance, as it will start positive and then tend rapidly to a negative constant. This lack of demographic significance is unlikely to be desirable. Therefore, orthogonal age functions can conflict with a desire for each AP term to have distinct demographic significance for models with parametric age functions coming from different functional families.

In summary, we find that orthogonality between age functions makes most sense when the age functions come from the same family, such as polynomials, and therefore can be orthogonalised easily. For models with very different functional forms for the age functions, orthogonalisation is unlikely to be desirable as it will conflict with a desire to give each AP term distinct demographic significance.

8.1 Parameter uncertainty

Having obtained a set of parameters by fitting a model to data with some set of arbitrary identifiability constraints, it is common to investigate the degree of uncertainty associated with these estimated parameters. A number of techniques have been developed to do this, for instance

• using the asymptotic normality of parameters estimated by maximum likelihood methods, as in Brouhns et al. (Reference Brouhns, Denuit and Vermunt2002b);

• using a “semi-parametric” bootstrap based on Poisson (or binomial) death counts, as in Brouhns et al. (Reference Brouhns, Denuit and Van Keilegom2005);

• using a residual bootstrapping method, such as that developed in Koissi et al. (Reference Koissi, Shapiro and Hognas2006) or the more complicated techniques discussed in D’Amato et al. (Reference D’Amato, Lorenzo, Haberman, Russolillo and Sibillo2011) and Debón et al. (Reference Debón, Montes, Mateu, Porcu and Bevilacqua2008), (Reference Debón, Martinez-Ruiz, Montes, Martínez-Ruiz and Montes2010), and

• using Bayesian Markov chain Monte Carlo (MCMC) methods, as in Czado et al. (Reference Czado, Delwarde and Denuit2005).

All of these techniques were developed for the LC model, as the simplest and most widely used mortality model. In the following section, we follow this convention and implicitly assume that we are dealing with the LC model. However, in principle, they could all be used with any other AP mortality model.

The first three of these methods have been tested and compared in Renshaw & Haberman (Reference Renshaw and Haberman2008) and all four were compared in Li (Reference Li2014). It is important that any conclusions drawn from them do not depend upon the arbitrary identifiability constraints imposed in the model. Since the fitted mortality rates do not change under the invariant transformations of the model, their variability due to parameter uncertainty should not depend on the identifiability constraints imposed either. Appropriate methods for determining parameter uncertainty should ensure this. Two users of a mortality model, using the same data and method for investigating parameter uncertainty, but using different (but equally valid) identifiability constraints should find the same degree of variability of mortality rates under parameter uncertainty.

It is therefore desirable to start from the difference between the observed and fitted mortality rates, since this will be independent of the identifiability constraints chosen from them model and ensure that our results are consistent with observations. For instance, in Brouhns et al. (Reference Brouhns, Denuit and Van Keilegom2005), Poisson-distributed random death counts were generated at each age and yearFootnote ²⁴. The distribution of the bootstrapped death counts is therefore unaffected by which identifiability constraints are imposed. Likewise, the fitting residuals used in Koissi et al. (Reference Koissi, Shapiro and Hognas2006) depend only on the actual and fitted death counts and thus not on the identifiability constraints used in fitting the model. Therefore, estimates of the impact of parameter uncertainty on observable quantities, such as fitted mortality rates or life expectancies, will be independent of the arbitrary identifiability constraints.

However, estimates for the variability of the model parameters will still only be valid conditional on the chosen set of identifiability constraints. For instance, imposing the constraint $\kappa^{(i)}_1 = 0$ in a model will mean that $\kappa^{(i)}_1$ will trivially not show any variability using the Brouhns et al. (Reference Brouhns, Denuit and Van Keilegom2005) or Koissi et al. (Reference Koissi, Shapiro and Hognas2006) methods, but this will not be the case for other choices of constraints. Therefore, the observed parameter uncertainty should be seen only in the context of the identifiability constraints applied.

It is also important to ensure that the model is fully identified when using these bootstrapping approaches. If the model is not fully identified, we may observe spurious variation in the parameters which does not lead to real variability in the fitted mortality rates. This is of most practical relevance with the orthogonality constraints for models such as the LC2 model in equation (13), as these are often overlooked if maximum likelihood or minimum deviance techniques are used to fit the model.

The alternative approach to starting from the difference between observed and expected mortality rates is to consider the distribution of the model parameters directly. However, methods which generate new samples of parameters directly, such as the asymptotic method of Brouhns et al. (Reference Brouhns, Denuit and Vermunt2002b) or the Bayesian techniques of Czado et al. (Reference Czado, Delwarde and Denuit2005), must be used with considerably more care.

First, consider the asymptotic method of Brouhns et al. (Reference Brouhns, Denuit and Vermunt2002b). This assumes that the variation of the maximum likelihood parameters is given by the information matrix (i.e. the second derivative of the log-likelihood, $\mathcal{L}$ , with respect to the model parameters evaluated at the selected parameter estimates). The first thing to note here is that, in order to identify the model, the likelihood being maximised is the constrained likelihood. Starting from the forms of the likelihood function in Hunt & Blake (Reference Hunt and Blake2020b), this means that we use Lagrangian multipliers to impose the constraints. For example, to impose the Lee & Carter (Reference Lee and Carter1992) model constraints involves adjusting the likelihood function by

\begin{align*}\mathcal{L}(d_{x,t};\,\{\alpha,\beta,\kappa\}) &\to \mathcal{L}(d_{x,t};\,\{\alpha,\beta,\kappa\}) - \lambda_1 \sum_t \kappa_t - \lambda_2 \left(1-\sum_x \beta_x\right)\end{align*}

Therefore, the information matrix is explicitly dependent upon the identifiability constraints imposed. For instance, we can see this by considering the second derivative of the likelihood with respect to the age function $\beta_x$

\begin{align*}\frac{\partial^2 \mathcal{L}}{\partial (\beta_x)^2} &= -\sum_t \mathbb{V}ar(D_{x,t}) (\kappa_t)^2\end{align*}

if we use the canonical link function, as discussed in Hunt & Blake (Reference Hunt and Blake2020b). If we apply the transformation in equation (9), $\beta_x$ is unchanged. However, we have

\begin{align*}\frac{\partial^2 \mathcal{L}}{\partial (\hat{\beta}_x)^2} &= -\sum_t \mathbb{V}ar(D_{x,t}) (\hat{\kappa}_t)^2 \\&= -\sum_t \mathbb{V}ar(D_{x,t}) (\kappa_t + b)^2 \\&= \frac{\partial^2 \mathcal{L}}{\partial (\beta_x)^2} +2b\sum_t \mathbb{V}ar(D_{x,t})\kappa_t - b^2 \sum_t \mathbb{V}ar(D_{x,t})\end{align*}

In this case, the form of the information matrix with respect to $\beta_x$ has changed under a transformation which did not change $\beta_x$ itself. This needs to be taken into consideration carefully and may explain the variation in the uncertainty in the fitted mortality rates observed in Renshaw & Haberman (Reference Renshaw and Haberman2008) when the identifiability constraints are altered.

Next, we consider Bayesian techniques, such as MCMC. As discussed in Nielsen & Nielsen (Reference Nielsen and Nielsen2014), these can often appear to solve identifiability issues but in fact confuse and disguise them. The use of Bayesian methods often involves consideration of the posterior distribution, $\pi$ , of the parameters given by

\begin{align*}\ln(\pi(\{\alpha,\beta,\kappa\})) &= \mathcal{L}(d_{x,t};\,\{\alpha,\beta,\kappa\}) + \ln(\phi(\{\alpha,\beta,\kappa\})) + \text{constant}\end{align*}

where $\phi$ is the prior distribution for the parameters. The log-likelihood function, $\mathcal{L}(d_{x,t};\,\{\alpha,\beta,\kappa\})$ , is unchanged by the invariant transformations of the model parameters and so does not depend upon the chosen identifiability constraints. However, in general, the prior distribution $\phi$ will change under these transformations, unless it is very carefully chosen. This, in turn, means that the posterior distribution will also vary under the invariant transformations of the model, and so will depend implicitly on any identifiability constraints imposed.

A poorly chosen set of priors implicitly imposes a set of identifiability constraints upon the model. For example, a prior distribution that assumes $\kappa^{(i)}_t$ follows an AR(1) process around zero implicitly imposes a level on the period parameters. These implicit constraints may conflict with the explicit constraints subsequently imposed (such as a subsequent choice of the level of $\kappa^{(i)}_t$ ). Even when there are no conflicts, this implicit selection of identifiability constraints is opaque, and it is not clear which features of the posterior distribution are meaningful and which are mere artefacts of the identifiability scheme implicit in the prior.

We therefore recommend that the prior distribution of the model parameters, $\phi$ , is selected so that it is unchanged by the invariant transformations of the model. This enables a single set of identifiability constraints to be imposed upon the model without internal conflicts, with these constraints being clear and transparent to all other model users, and with the posterior distribution being independent of the arbitrary choice of identifiability constraints (just as the likelihood is).

8.2 Hypothesis testing

Identifiability issues also have important consequences if hypothesis testing on the parameters is performed. In general, hypotheses cannot be tested on the parameter values directly, since they depend upon the identifiability constraints. For instance, testing the hypothesis $\kappa_T = 0$ in the LC model is meaningless, since we can impose $\kappa_T = 0$ (or any other value) by our choice of identifiability constraint. We might be tempted to find combinations of the parameters which are invariant to the transformations of the parameters and test hypotheses based on these. For instance, we may wish to test the hypothesis that mortality is declining faster at age $x_1$ than at age $x_2$ using the LC model. To do this, we might note that the expected value of $B \equiv \frac{\beta_{x_1}}{\beta_{x_2}}$ is invariant under the transformation in equation (11) and so does not depend on the identifiability constraints, making it a suitable candidate for hypothesis testing. However, we would have to take care when using a statistic such as this, since it will be undefined in the case $\beta_{x_2} = 0$ , which could not be known before the model is fitted to data. In general, therefore, any tests of hypotheses should be performed on observable quantities such as the fitted mortality rates rather than the model parameters.

Direct hypothesis testing of the parameters in an AP model is not often performed in the literature, and therefore this discussion may appear to be of theoretical interest only. However, it is common to use a variety of statistical tests when determining the time series properties of the period functions. For instance, in Lee & Carter (Reference Lee and Carter1992) and Cairns et al. (Reference Cairns, Blake, Dowd, Coughlan, Epstein and Khalaf-Allah2011), Box–Jenkins methods were used to determine the preferred time series process for the period functions of different models. Based on the conclusions above, in many cases, the results of these statistical tests will depend on arbitrary choices made when identifying the model. The properties typically tested, such as stationarity, lagged dependence and cross correlation, will affect our projected mortality rates and so are matters of great practical importance. We should therefore treat with extreme caution the results of any such analysis. This subject is dealt with further in section 9.

In summary, not only do our estimates of the parameters of an AP model depend on the identifiability constraints when fitting the model, so do our estimates of the uncertainty attached to those parameter estimates. We should therefore avoid testing hypotheses on these parameter estimates, as our results will be dependent on the arbitrary identification scheme imposed. In general, methods of estimating parameter uncertainty which use bootstrapping techniques on the fitted mortality rates, which are independent of our choice of identifiability constraints, are likely to be preferred over methods which target the parameters directly. We must still ensure, however, that our models are fully identified when testing parameter uncertainty, as the parameters in a poorly identified model may show spurious differences in ways which do not affect the variability of the fitted mortality rates.

9.1 Models with non-parametric age functions

Consider the case of projecting an AP mortality model with non-parametric age functions, which has been fitted using data over the period [1, T] to give mortality rates at time $\tau > T$ . From equation (2), we could write this as

\begin{align*}\eta_{x,\tau} &= \alpha_x + \boldsymbol{\beta}^\top_x \boldsymbol{\kappa}_\tau\end{align*}

We can also see that the projected mortality rates for the future are unchanged by the use of the invariant transformations of the parameters in equations (12) and (11), just as the fitted mortality rates were for the past, i.e.,

\begin{align*}\eta_{x,\tau} &= \hat{\alpha}_x + \boldsymbol{\hat{\beta}}^\top_x \boldsymbol{\hat{\kappa}}_\tau\end{align*}

where

\begin{align*}\hat{\boldsymbol{\kappa}}_\tau &= A \boldsymbol{\kappa} _\tau + B \\\hat{\boldsymbol{\beta}}^\top_x &= \boldsymbol{\beta}^\top_x A^{-1} \\\hat{\alpha}_x &= \alpha_x - \boldsymbol{\beta}^\top_x A^{-1}B\end{align*}

Unlike the fitted parameters, however, the projected $\boldsymbol{\kappa}_\tau$ will be some random variable, whose distribution is a function of the fitted parameters, i.e., $\boldsymbol{\kappa}_\tau = P_\kappa(\tau;\,\{\boldsymbol{\kappa}\})$ . We said previously that we should use the same method of projection for all sets of parameters as a first step in ensuring that the projected mortality rates do not depend upon the identifiability constraints. However, for different identifiability constraints, these processes will be estimated from different sets of fitted parameters, e.g., if we use $P_\kappa(\tau;\,\{\boldsymbol{\kappa}\})$ to project the untransformed period parameters, we must use $P_\kappa(\tau;\,\{\boldsymbol{\hat{\kappa}}\})$ to project the transformed period parameters. If we combine this with the invariance of the projected mortality rates, we have

(27) \begin{align}\alpha_x + \boldsymbol{\beta}^\top_x P_\kappa(\tau;\,\{\boldsymbol{\kappa}\}) &= \hat{\alpha}_x + \boldsymbol{\hat{\beta}}^\top _x P_\kappa(\tau;\,\{\boldsymbol{\hat{\kappa}}\}) \notag \\[3pt]\alpha_x + \boldsymbol{\beta}^\top _x P_\kappa(\tau;\,\{\boldsymbol{\kappa}\}) &= \alpha_x - \boldsymbol{\beta}^\top_x A^{-1}B + \boldsymbol{\beta}^\top_x A^{-1} P_\kappa(\tau;\,\{A\boldsymbol{\kappa} + B\}) \notag \\[3pt]\boldsymbol{\beta}^\top_x P_\kappa(\tau;\,\{\boldsymbol{\kappa}\}) &= \boldsymbol{\beta}^\top_x A^{-1} \left[P_\kappa (\tau;\,\{A\boldsymbol{\kappa} + B\}) - B \right] \notag \\[3pt]P_\kappa(\tau;\,\{\boldsymbol{\kappa}\}) &= A^{-1} \left[P_\kappa(\tau;\,\{A\boldsymbol{\kappa} + B\}) - B \right] \notag \\[3pt]P_\kappa (\tau;\,\{A\boldsymbol{\kappa} + B\}) &= A P_\kappa (\tau;\,\{\boldsymbol{\kappa}\}) + B \end{align}

for general $\boldsymbol{\beta}_x$ , i.e., that the time series processes we use to project the period functions are location and scale preserving. This is also discussed in Nielsen & Nielsen (Reference Nielsen and Nielsen2014).

One common practice is to use univariate time series processes to project the period functions, on the grounds that they are uncorrelated over the historical sample. For example, in Hyndman & Ullah (Reference Hyndman and Ullah2007), p. 4948, when considering the selection of suitable time series processes for projecting a model with non-parametric age functions, it was statedFootnote ²⁶

This logic was reiterated in Hyndman et al. (Reference Hyndman, Booth and Yasmeen2013) for a related model, as “There is no need to consider vector models because the $\kappa^{(i)}_{t}$ coefficients are all uncorrelated by construction”.

However, we saw in section 5 that the lack of correlation between the different period functions is a product of the choice of identifiability constraints, and that we could find alternative parameters which gave identical fitted mortality rates which had non-zero correlation. Choosing univariate time series processes will therefore not give well-identified projections, but instead will give projected mortality rates which are dependent upon the identifiability constraints chosen.

The first conclusion we can draw is that we should always use multivariate processes to project mortality models with more than one AP term. Using a multivariate framework allows us to consider the period functions together and so encourages a unified approach to modelling them, rather than focusing on each period function separately. It also allows the invariant transformations in equations (11) and (12) to be applied to the time series processes directly to check whether they are well-identified.

The use of multivariate processes means that the order of integration of each of the time series processes should be the same. We should only consider the stationarity of the vector process as a whole, rather than of its individual components. It is common practice to use the highest order of integration for any of the individual period functions (usually first order) as the order of integration for all of them to avoid identification issues.

We can see this by taking a general multivariate time series process for $\boldsymbol{\kappa}_t$ from the class of VARIMA(p,d,q) processes

(28) \begin{align}\Delta^d \boldsymbol{\kappa}_t &= \boldsymbol{\mu} + \sum_{s=1}^p \Phi_s \Delta^d \boldsymbol{\kappa}_{t-s} + \sum_{r=0}^{q} \Psi_r \boldsymbol{\epsilon}_{t-r} \end{align}

and applying the transformations in equation (11) and (12) to give

\begin{align*}\Delta^d \boldsymbol{\hat{\kappa}}_t &= A\boldsymbol{\mu} + \Delta^d B - \sum_{s=1}^p A \Phi_s A^{-1} \Delta^d B + \sum_{s=1}^p A \Phi_s A^{-1} \Delta^d \boldsymbol{\hat{\kappa}}_{t-s} + \sum_{r=0}^{q} A\Psi_r \boldsymbol{\epsilon}_{t-r} \\&= \boldsymbol{\hat{\mu}} + \sum_{s=1}^p \hat{\Phi}_s \Delta^d \boldsymbol{\hat{\kappa}}_{t-s} + \sum_{r=0}^{q} \hat{\Psi}_r \boldsymbol{\hat{\epsilon}}_{t-r}\end{align*}

We therefore see that all general VARIMA(p, d, q) processes are location and scale invariant in the sense of equation (27), and so are well identified.

However, we also see from this that any specific structure we impose a priori on $\boldsymbol{\mu}$ , $\Phi_s$ and $\Psi_r$ will not be invariant under these transformations. Our second conclusion is, therefore, that we should not assume any pre-specified locations, scales or correlations between our period functions by assuming a prior structure for the matrices governing the time series processes that drives them.

In practice, in order to be invariant to transformations of the form in equation (11), we should always allow for the possibility of both cross-lags between the time series and contemporaneous correlations between the innovations, even if these are not evident from inspection of the fitted time series. In situations where our arbitrary identification constraints set some of these time series parameters to zero, this will emerge naturally from their estimation and do not need to be imposed by the model user.

Finally, we observe that all VARIMA time series models are invariant to simple rescalings of the period functions, i.e., using the transformation in equation (11), the matrix A being diagonal. Therefore, all time series processes are invariant under alternative choices of normalisation scheme. However, having a consistent scale for all period functions is desirable as it assists with the numerical estimation of the time series parameters.

In summary, the use of multivariate time series processes means that we should not treat the period functions differently when projecting them, as the invariant transformation in equation (11) means that the AP terms are interchangeable, which, in turn, means that we can rotate them without changing the fit to data or the demographic significance of any of the parameters.

9.2 Projecting the LC2 model

As a practical example of this, consider projecting the LC2 model in section 5. Tests on the fitted time series processes from Figure 3(a) show that they are uncorrelated, which is a direct result of the identifiability constraint $\sum_t \kappa^{(1)}_t \kappa^{(2)}_t = 0$ . However, we saw that the model period functions given in Figure 3(b) had a correlation of 75%, but gave exactly the same fitted mortality rates. We should, therefore, use multivariate processes for both set of parameters.

Testing these parameters for stationarity, we find that both of the period functions in Figure 3 are non-stationary. We would, therefore, be justified in using a multivariate random walk for both sets of period functions (i.e. those from both Figure 3(a) and from Figure 3(b)).

We can see directly that this time series process is well identified, since if

\begin{align*}\boldsymbol{\kappa}_t &= \boldsymbol{\kappa}_{t-1} + \boldsymbol{\mu} + \boldsymbol{\epsilon}_t\end{align*}

then

\begin{align*}\boldsymbol{\hat{\kappa}}_t &= \boldsymbol{\hat{\kappa}}_{t-1} + A \boldsymbol{\mu} + A\boldsymbol{\epsilon}_t\end{align*}

after applying the transformations in equations (12) and (11). We see that integrated time series are unchanged by changes in the level of the period functions, and so are automatically invariant to the transformation in equation (12).

At this point, it is also worth noting an important side effect of imposing orthogonality on the period functions in the LC2 model. $\kappa^{(1)}_t$ is usually found to be linear to quite a good approximation; so much so that this was called the “universal pattern of mortality decline” in Tuljapurkar et al. (Reference Tuljapurkar, Li and Boe2000). By construction, therefore, $\kappa^{(2)}_t$ cannot be roughly linear if we impose orthogonality, which makes projecting it trickier. We believe that this could be one of the reasons why the LC2 model is not more widely used, despite being a natural extension of the classic LC model. Often, the second term appears quadratic to quite a good approximationFootnote ²⁷. Various authors (such as Renshaw & Haberman (Reference Renshaw and Haberman2003b) and Yang et al. (Reference Yang, Yue and Huang2010)) have suggested using break points or “hinges” in order to continue to use linear projection processes. However, this is a case of selecting a time series process specifically because of a feature of the period functions that is present solely because of the particular identifiability constraints imposed, and therefore the resulting projections will not be well identified.

Using a multivariate random walk with drift for the time series processes in Figures 3(a) and (b) gives the projected $\kappa^{(2)}_t$ period functions in Figure 5(a)Footnote ²⁸. While these projections appear quite different, the projected mortality rates from them at age 65, shown in Figure 5(b), are identical, thereby demonstrating that we have, indeed, chosen a well-identified projection method for the LC2 model.

Figure 5

Projections from the LC2 model.

9.3 Models with parametric age functions

In section 6, it was shown that models with parametric age functions have subtly different identifiability issues when fitting them to data to those with non-parametric age functions. This is due to the transformations in equation (11) not being allowed, since they changed the definition of the age functions and hence gave a different, but equivalent, model. However, we saw that this meant we could select between equivalent models, which had different definitions of the age functions, but gave identical fitted mortality rates. This was done in order to choose models with desirable properties such as a common normalisation scheme and orthogonal age functions. These subtle differences are also present when projecting the model.

First, the transformations in equation (12) are used to impose a level on the period functions through identifiability constraints in models with parametric age functions in exactly the same manner as for models with non-parametric age functions. Consequently, we need to ensure that the time series processes used to project the period functions are identifiable under changes in location in exactly the same way as described for non-parametric AP terms above. This means either using integrated time series processes or allowing for mean reversion to a non-zero level.

However, the transformations in equation (11) are not needed in models with parametric age functions, since applying them would fundamentally change the model. Since we cannot normalise the age functions during the fitting process, we must instead define normalised (or self-normalising) age functions in advance. We cannot impose orthogonality on the age functions, although we could define orthogonal age functions a priori.

In addition, we cannot impose orthogonality on the period functions, as was done for the LC2 model, and therefore the period functions in models with parametric age functions will be correlated in general. This means that it is natural to project the period functions in such models using multivariate time series processes, just as we should in models with non-parametric age functions. However, because the transformations in equation (11) are not applicable in models with parametric age functions, if we use a VARIMA(p, d, q) time series process for the period functions, as in equation (28), we only have to ensure that the time series process is invariant to the transformation in equation (12). To do this, we substitute the transformed parameters, $\hat{\boldsymbol{\kappa}}_t = \boldsymbol{\kappa}_y + B$ , into equation (28) to find

\begin{align*}\Delta^d \boldsymbol{\hat{\kappa}}_t &= \boldsymbol{\mu} + \Delta^d B - \sum_{s=1}^p \Phi_s \Delta^d B + \sum_{s=1}^p \Phi_s \Delta^d \boldsymbol{\hat{\kappa}}_{t-s} + \sum_{r=0}^{q} \Psi_r \boldsymbol{\epsilon}_{t-r} \\&= \boldsymbol{\hat{\mu}} + \sum_{s=1}^p \Phi_s \Delta^d \boldsymbol{\hat{\kappa}}_{t-s} + \sum_{r=0}^{q} \Psi_r \boldsymbol{\epsilon}_{t-r}\end{align*}

Although the drift term, $\boldsymbol{\mu}$ , has changed as a result of this transformation the matrices $\Phi_s$ and $\Psi_r$ have not. Consequently, we see that any structure we impose a priori upon the moving average and autocorrelation of the time series process is also unchanged by changes in the identifiability constraints in models with parametric age functions. This means that, in theory, it is possible to give each term distinct structure, such as different orders of integration or numbers of lags. This may be felt to be desirable if doing so gives projections with greater demographic significance.

For example, consider the exponential model in equation (19). In this, we interpret $\kappa^{(2)}_t$ as representing the component of mortality change specific to very young ages, in excess of the changes in general mortality rates governed by $\kappa^{(1)}_t$ . If we had a strong prior belief that these should mean-revert to a natural level (for instance, because we believed that infants should not receive systematically better or worse medical care than the general population), we might chose to allow our subjective demographic significance for the term to overrule a purely statistical evaluation of the time series process in this case. Because we do not use the transformation in equation (11) to enforce a constraint when fitting the model, we do not have to ensure that our projection process is robust to its application when the model is projected.

We may also feel that such a restriction will give projected mortality rates with greater biological reasonableness. For example, we may have biological reasons for believing that infant mortality rates should always be higher than those for young children at age five, say. However, using a non-stationary time series process for $\kappa^{(2)}_t$ allows there to be scenarios with non-zero probability where this is violated, and therefore we might wish to use a stationary time series process for $\kappa^{(2)}_t$ to avoid any scenarios felt to be biologically unreasonableFootnote ²⁹.

However, such arguments ignore the fact that, for any model with parametric age functions, there are a range of equivalent models which give identical fitted mortality rates and so, ideally, should be projected using the same time series processes to give identical projected mortality ratesFootnote ³⁰. There may also be features, such as changes in trend, which are present in the period functions for one model but absent in an equivalent model, and so are not objective features of the data. Since these equivalent models are linked by the transformation in equation (11), it is still highly desirable to use general VARIMA processes, with no a priori structure placed on them, just as for models with non-parametric age functions.

In practice, it is not often that the demographic significance of a term in an AP mortality model leads to specific requirements about how it should be projected. For instance, while we may seek to rule out any possibility of mortality rates being lower at birth than at age five in the exponential model, this is highly unlikely to occur even if non-stationary time series processes are used for $\kappa^{(2)}_t$ , since it is inconsistent with the historical data. We therefore recommend that general, well-identified, multivariate VARIMA processes are used to project the period functions in models with parametric age functions, unless these are shown experimentally to give biologically implausible projected mortality ratesFootnote ³¹.

9.4 Summary

In summary, we can say that in order to obtain projections which are well-identified from an AP model, we need to work backwards from our desire for time series processes which do not change form under the invariant transformations in equations (11) and (12). This means that we should always use multivariate time series processes, as these support a unified approach to projection and allow us to check identifiability easily.

Identifiability also means, in general, that we should not treat the different period functions differently. In practice, this means assuming as little structure a priori for the time series processes as possible and using the same order of integration for each period function. In models with parametric age functions, however, there may be conflicts between achieving this and the biological reasonableness of the projected mortality rates. Treating the different period functions in the same manner is still highly desirable, however, as it avoids using different processes to project equivalent models, and often emerges naturally out of a statistical analysis of the fitted period functions. These conclusions are summarised in Table 1 below.

Table 1.

Requirements for identifiable projection methods in AP mortality models.