4.1.1 AP model

Under the AP model, we have

(16) \begin{equation} \log \mu_{x, y} = \alpha_x + \kappa_y,\; x = 1,\ldots, n_a,\; y = 1,\ldots, n_y\end{equation}

First, we write the model in the standard form (1) and compute its rank. Let $\boldsymbol{\theta} = (\boldsymbol{\alpha}^\prime,\boldsymbol{\kappa}^\prime)^\prime$ be the vector of regression coefficients. The model matrix is

(17) \begin{equation} \textbf{\textit{X}} = [\textbf{\textit{X}}_a\,{:}\,\textbf{\textit{X}}_y] = [\textbf{1}_{n_y} \otimes \textbf{\textit{I}}_{n_a}\,{:}\,\textbf{\textit{I}}_{n_y} \otimes \textbf{1}_{n_a}],\; n_a\, n_y \times (n_a+n_y)\end{equation}

and $\log \boldsymbol{\mu} = \textbf{\textit{X}} \boldsymbol{\theta}$ . Second, we find a basis for the null space of $\textbf{\textit{X}}$ . Since $\textbf{\textit{X}}$ is $n_a\,n_y \times (n_a+n_y)$ and $r(\textbf{\textit{X}}) = n_a + n_y - 1$ , the dimension of $\mathcal{N}(\textbf{\textit{X}})$ is one. In this simple case, we can just write down a basis vector for $\mathcal{N}(\textbf{\textit{X}})$ . We consider

(18) \begin{equation} \textbf{\textit{n}} = \left(\begin{array}{c} \phantom{-}\textbf{1}_{n_a} \\[6pt] - \textbf{1}_{n_y}\end{array}\right)\end{equation}

and it is easy to check that $\textbf{\textit{X}} \textbf{\textit{n}} = \textbf{0}$ and so $\textbf{\textit{n}}$ spans $\mathcal{N}(\textbf{\textit{X}})$ .

We adopt the following approach to computation. We use two constraint systems: first, a standard constraint system, i.e., one that is found in the literature: for the AP model we take $\sum \kappa_y = 0$ ; second, a random constraint system: we set $\sum_1^{n_a+n_y} u_i\theta_i= 0$ where $u_i,\, i = 1,\ldots, n_a+n_y$ , are realisations of independent uniform variables on [0, 1], i.e., $\mathcal{U}(0,1)$ . Let $\hat {\boldsymbol{\theta}}_s$ and $\hat {\boldsymbol{\theta}}_r$ be the estimates of $\boldsymbol{\theta}$ under the two systems. Then, by our fundamental result (2), $\hat {\boldsymbol{\theta}}_s - \hat {\boldsymbol{\theta}_r} \in \mathcal{N}(\textbf{\textit{X}})$ and so

(19) \begin{equation} \hat {\boldsymbol{\theta}}_s - \hat {\boldsymbol{\theta}}_r = A\left(\begin{array}{c} \phantom{-}\textbf{1}_{n_a} \\[6pt] - \textbf{1}_{n_y}\end{array}\right)\end{equation}

for some scalar A. Equating coefficients, we find that

(20) \begin{equation}\hat {\boldsymbol{\alpha}}_s - \hat {\boldsymbol{\alpha}}_r = A \textbf{1}_{n_a}, \quad \hat {\boldsymbol{\kappa}}_s - \hat {\boldsymbol{\kappa}}_r = - A \textbf{1}_{n_y}\end{equation}

where $\hat {\boldsymbol{\alpha}}_s$ , $\hat {\boldsymbol{\kappa}}_s$ , $\hat {\boldsymbol{\alpha}}_r$ and $\hat {\boldsymbol{\kappa}}_r$ are the components of $\hat{\boldsymbol{\theta}}_s$ and $\hat {\boldsymbol{\theta}_r}$ , respectively. The value of A depends on the particular values of $u_i$ ; in our computation, we found $A = 20.7$ although any value of A is possible since it corresponds to adding A to $\alpha_x$ and subtracting A from $\kappa_y$ in (16). We note that (20) confirms what has been widely observed, namely that the $\boldsymbol{\alpha}$ and $\boldsymbol{\kappa}$ are only estimable up to an additive constant. A careful discussion of what can and cannot be estimated in the AP model can be found in Clayton & Schifflers (Reference Clayton and Schifflers1987a).

Forecasting in the AP model is done by forecasting the $\boldsymbol{\kappa}$ values, keeping the $\boldsymbol{\alpha}$ values fixed at their estimated values and then using equation (16) to forecast the values of $\log \mu$ at each age. Informally, we argue that since $\hat {\boldsymbol{\kappa}}_s$ and $\hat {\boldsymbol{\kappa}}_r$ are parallel with $\hat {\boldsymbol{\kappa}}_s - \hat {\boldsymbol{\kappa}}_r = - A \textbf{1}_{n_y}$ , we require that their forecast values, $\hat {\boldsymbol{\kappa}}_{s, f}$ and $\hat {\boldsymbol{\kappa}}_{r, s}$ say, are also parallel with $\hat {\boldsymbol{\kappa}}_{s, f} - \hat {\boldsymbol{\kappa}}_{r, f} = - A \textbf{1}_{n_f}$ , where $n_f$ is the length of the forecast. Now, the forecast values under both constraint systems will be equal since the change in the $\boldsymbol{\kappa}$ forecast values is exactly compensated for by the change in the $\boldsymbol{\alpha}$ values. This will be achieved with an ARIMA $(p, \delta, q)$ model, $\delta = 1$ or 2, since the fitted means of the forecasts for $\hat {\boldsymbol{\kappa}}_s$ and $\hat {\boldsymbol{\kappa}}_r$ will differ by $-A$ .

This informal argument extends to the case when $\hat {\boldsymbol{\kappa}}_s$ and $\hat {\boldsymbol{\kappa}}_r$ differ by a linear function, as is the case in the age-period-cohort model for example. Here, we want the forecasts to obey the same linear relationship. A formal argument is given in section 4.3.2 when we discuss the age-period-cohort-improvements model.

We make the important remark that not all forecasting models used in the forecasting of mortality lead to forecasts which are invariant with respect to the choice of constraints. For example, if we forecast $\boldsymbol{\kappa}$ with an AR(1) model without a mean, then the forecasts will not be invariant. In this paper, we take the position that if two time series differ by a function, then the forecasts should differ by the same function. For example, in the AP model, $\hat {\boldsymbol{\kappa}}_s$ and $\hat {\boldsymbol{\kappa}}_r$ differ by a constant function, i.e., are parallel; our forecasts should also be parallel. In the APC model, estimates of the period terms and of the cohort terms differ by linear functions. Again, we want our forecasts of the period and cohort terms to obey the same functional relationship. In this way, we will find that not only are fitted values of mortality invariant with respect to the choice of constraints so also are their forecast values.

4.1.2 Smooth AP model

We turn now to the effect of constraints when some model terms are smoothed. In the AP model, the forecast values of $\log \boldsymbol{\mu}$ are more regular by age if the age parameters $\boldsymbol{\alpha}$ are smoothed (see Delwarde et al., Reference Delwarde, Denuit and Eilers2007; Currie, Reference Currie2013). Let $\boldsymbol{\alpha} =\textbf{\textit{B}}_a \textbf{\textit{a}}$ , where $\textbf{\textit{B}}_a$ is a B-spline regression matrix along age; here, $\textbf{\textit{B}}_a$ is $n_a \times c_a$ , where $c_a$ is the number of B-splines in the basis. The model matrix (17) becomes

(21) \begin{equation} \textbf{\textit{X}} = [\textbf{\textit{X}}_a\,{:}\,\textbf{\textit{X}}_y] = [\textbf{1}_{n_y} \otimes \textbf{\textit{B}}_a\,{:}\,\textbf{\textit{I}}_{n_y} \otimes \textbf{1}_{n_a}],\; n_a\,n_y \times (c_a+n_y)\end{equation}

with rank $c_a+n_y-1$ ; the regression coefficients $\boldsymbol{\theta} =(\textbf{\textit{a}}^\prime, \boldsymbol{\kappa}^\prime)^\prime$ . We distinguish between the regression coefficients denoted $\boldsymbol{\theta}$ and the parameters of interest which we denote $\boldsymbol{\theta}^*$ ; here $\boldsymbol{\theta}^* = (\boldsymbol{\alpha}^\prime, \boldsymbol{\kappa}^\prime)^\prime$ . Let $\textbf{\textit{X}}^*$ be the model matrix (17) in the unsmoothed case. We note that

(22) \begin{equation}\textbf{\textit{X}} \boldsymbol{\theta} = \textbf{\textit{X}}^* \boldsymbol{\theta}^*\end{equation}

since $(\textbf{1}_{n_y} \otimes \textbf{\textit{B}}_a) \textbf{\textit{a}} = (\textbf{1}_{n_y} \otimes \textbf{\textit{I}}_{n_a}) \boldsymbol{\alpha}$ . We will see that (22) is a very convenient identity. If there is no smoothing then $\boldsymbol{\theta} = \boldsymbol{\theta}^*$ and $\textbf{\textit{X}} = \textbf{\textit{X}}^*$ .

The penalty matrix is

\begin{equation*} \textbf{\textit{P}} = \mbox{blockdiag}\{\lambda_a \textbf{\textit{D}}^\prime \textbf{\textit{D}}, \textbf{0}\}\end{equation*}

where $\textbf{\textit{D}}$ is the difference matrix of order d acting on the smoothed age coefficients $\textbf{\textit{a}}$ , $\textbf{0}$ is the $n_y \times n_y$ matrix of 0s acting on the unsmoothed year coefficients $\boldsymbol{\kappa}$ and $\lambda_a$ is the smoothing parameter. We compute $\mathcal{N}(\textbf{\textit{X}}^\prime \tilde {\textbf{\textit{W}}} \textbf{\textit{X}} +\textbf{\textit{P}})$ for the model in equation (21). Let

\begin{equation*} \textbf{\textit{n}} = \left(\begin{array}{c} \phantom{-}\textbf{1}_{c_a} \\[6pt] - \textbf{1}_{n_y}\end{array}\right)\end{equation*}

Using the fact that $\textbf{\textit{B}}_a \textbf{1}_{c_a} = \textbf{1}_{n_a}$ , we see that $\textbf{\textit{X}} \textbf{\textit{n}} = \textbf{0}$ , i.e., $\textbf{\textit{n}} \in \mathcal{N}(\textbf{\textit{X}})$ . Furthermore, we have $\textbf{\textit{D}} \textbf{1}_{c_a} = \textbf{0}$ for any order $d > 0$ of the penalty; hence, $\textbf{\textit{P}} \textbf{\textit{n}} = \textbf{0}$ , i.e., $\textbf{\textit{n}} \in \mathcal{N}(\textbf{\textit{P}})$ . Hence, by (13), $\mathcal{N}(\textbf{\textit{X}}^\prime \tilde {\textbf{\textit{W}}} \textbf{\textit{X}}+ \textbf{\textit{P}}) = \mathcal{N}(\textbf{\textit{X}}^\prime \tilde {\textbf{\textit{W}}} \textbf{\textit{X}}) \cap \mathcal{N}(\textbf{\textit{P}}) = \textbf{\textit{n}}$ . Thus, if $\hat {\boldsymbol{\theta}_s}$ and $\hat {\boldsymbol{\theta}_r}$ are any two estimates of $\boldsymbol{\theta}$ , we have

\begin{equation*} \hat {\boldsymbol{\theta}_s} - \hat {\boldsymbol{\theta}_r} = A \left(\begin{array}{c} \phantom{-}\textbf{1}_{c_a} \\[6pt] - \textbf{1}_{n_y}\end{array}\right)\end{equation*}

for some scalar A. Pre-multiplying both sides by $\mbox{blockdiag}\{\textbf{\textit{B}}_a\,{:}\,\textbf{\textit{I}}_{n_y}\}$ , we find

\begin{equation*} \left(\begin{array}{c} \hat {\boldsymbol{\alpha}}_s \\[6pt] \hat {\boldsymbol{\kappa}}_s\end{array}\right)- \left(\begin{array}{c} \hat {\boldsymbol{\alpha}}_r \\[6pt] \hat {\boldsymbol{\kappa}}_r\end{array}\right) = A \left(\begin{array}{c} \phantom{-}\textbf{1}_{n_a} \\[6pt] - \textbf{1}_{n_y}\end{array}\right)\end{equation*}

exactly as in the unsmoothed case. We conclude that smoothing has no effect on our invariance results for the AP model.

We have presented a mathematical discussion of the AP model. An important practical point is that from a computational point of view, we do not need the detailed understanding provided by this discussion. We can find the rank and effective rank and check the invariance of the fitted and forecast values easily in R (or other computational tool). We will make further comments on the computational aspects of this work as we go.

4.2.1 APC model

Under the APC model, we have

(23) \begin{equation} \log \mu_{x, y} = \alpha_x + \kappa_y + \gamma_{c(x, y)},\; x = 1,\ldots, n_a,\; y = 1,\ldots, n_y \end{equation}

where c(x, y) is the cohort index for age x in year y; with our notation, we have $c(x, y) = n_a - x + y$ . First, we write model (23) in the standard form (1) and compute its rank. Let $\boldsymbol{\theta} = (\boldsymbol{\alpha}^\prime,\boldsymbol{\kappa}^\prime, \boldsymbol{\gamma}^\prime)^\prime$ . The model matrix is

(24) \begin{equation} \textbf{\textit{X}} = [\textbf{\textit{X}}_a\, :\, \textbf{\textit{X}}_y \, :\, \textbf{\textit{X}}_c],\; n_a\,n_y \times (2n_a+2n_y-1)\end{equation}

where $\textbf{\textit{X}}_a$ and $\textbf{\textit{X}}_y$ are defined in (17); the row of $\textbf{\textit{X}}_c$ corresponding to age x and year y contains a one in column $n_a-x+y$ and zeros elsewhere. The rank of $\textbf{\textit{X}}$ is $2n_a+2n_y-4$ and so the dimension of $\mathcal{N}(\textbf{\textit{X}})$ is three. The relationship between the estimates of $\boldsymbol{\theta}$ under different constraint systems is determined by $\mathcal{N}(\textbf{\textit{X}})$ .

We now find a basis $\{\textbf{\textit{n}}_1, \textbf{\textit{n}}_2, \textbf{\textit{n}}_3\}$ for $\mathcal{N}(\textbf{\textit{X}})$ . We note that this basis is not unique, but we can find a basis that clarifies the relationship between the different estimates of $\boldsymbol{\theta}$ . First, we argue that each of $\textbf{\textit{X}}_a$ , $\textbf{\textit{X}}_y$ and $\textbf{\textit{X}}_c$ contains a one in each row and zeros elsewhere. Hence, we can take $\textbf{\textit{n}}_1$ and $\textbf{\textit{n}}_2$ equal to

(25) \begin{equation} \textbf{\textit{n}}_1 = \left(\begin{array}{c} \phantom{-} \textbf{1}_{n_a} \\[6pt] - \textbf{1}_{n_y} \\[6pt] \phantom{-} \textbf{0}_{n_c}\end{array}\right), \quad \textbf{\textit{n}}_2 = \left(\begin{array}{c} \phantom{-} \textbf{1}_{n_a} \\[6pt] \phantom{-} \textbf{0}_{n_y} \\[6pt] - \textbf{1}_{n_c}\end{array}\right)\end{equation}

It remains to find a suitable $\textbf{\textit{n}}_3$ . Some computation is helpful. We fit the APC model under a set of standard constraints, which we take as

(26) \begin{equation} \sum_1^{n_y} \kappa_y = \sum_1^{n_c} \gamma_c = \sum_1^{n_c}c \gamma_c = 0\end{equation}

where c is the cohort index (Cairns et al., Reference Cairns, Blake, Dowd, Coughlan, Epstein, Ong and Balevich2009). An alternative is to weight by the number of times cohort c appears in the data, say, $w_c$ . Thus, we could use

\begin{equation*} \sum_1^{n_y} \kappa_y = \sum_1^{n_c} w_c \gamma_c = \sum_1^{n_c}w_c c \gamma_c = 0\end{equation*}

as in Richards et al. (to appear). Our main point is that, although we will get (slightly) different parameter estimates with these constraint systems, our forecasts of mortality will be identical. We also fit the model with a set of random constraints

\begin{align*} \sum u_{i, j} \theta_j = 0,\;i = 1,2,3,\; j = 1,\ldots, n_a+n_y+n_c\end{align*}

where the $u_{i, j}$ are independent realisations from the $\mathcal{U}(0,1)$ distribution. Let $\hat {\boldsymbol{\theta}}_s$ and $\hat {\boldsymbol{\theta}}_r$ be the maximum likelihood estimates of $\boldsymbol{\theta}$ under the two constraint systems. Let $\hat {\boldsymbol{\alpha}}_s$ , $\hat {\boldsymbol{\kappa}}_s$ and $\hat {\boldsymbol{\gamma}}_s$ be the components of $\hat {\boldsymbol{\theta}}_s$ with a corresponding notation for the components of $\hat{\boldsymbol{\theta}}_r$ . The null space of $\textbf{\textit{X}}$ characterises $\hat {\boldsymbol{\theta}}_s - \hat {\boldsymbol{\theta}}_r$ , so we define

(27) \begin{equation} \Delta \hat {\boldsymbol{\alpha}} = \hat {\boldsymbol{\alpha}}_s - \hat {\boldsymbol{\alpha}}_r, \; \Delta \hat {\boldsymbol{\kappa}} = \hat {\boldsymbol{\kappa}}_s - \hat {\boldsymbol{\kappa}}_r, \; \Delta \hat {\boldsymbol{\gamma}} = \hat {\boldsymbol{\gamma}}_s - \hat {\boldsymbol{\gamma}}_r\end{equation}

Figure 2 is a plot of $\Delta \hat {\boldsymbol{\alpha}}$ , $\Delta \hat {\boldsymbol{\kappa}}$ and $\Delta \hat {\boldsymbol{\gamma}}$ and suggests that these values are linear with slopes C, $-C$ and C, respectively, for some C. Thus,

\begin{eqnarray*} \hat {\boldsymbol{\alpha}}_s - \hat {\boldsymbol{\alpha}}_r &=& a \textbf{1}_{n_a} + b \textbf{\textit{x}}_a \\ \hat {\boldsymbol{\kappa}}_s - \hat {\boldsymbol{\kappa}}_r &=& c \textbf{1}_{n_y} - b \textbf{\textit{x}}_y \\ \hat {\boldsymbol{\gamma}}_s - \hat {\boldsymbol{\gamma}}_r &=& d \textbf{1}_{n_c} + b \textbf{\textit{x}}_c\end{eqnarray*}

for some a, b, c and d. Since $\textbf{\textit{n}}_3$ is defined only up to scale, we may take $b = 1$ and conjecture that $\textbf{\textit{n}}_3$ has the form

(28) \begin{equation} \textbf{\textit{n}}_3 = \left(\begin{array}{c} a \textbf{1}_{n_a} + \textbf{\textit{x}}_a \\[6pt] c \textbf{1}_{n_y} - \textbf{\textit{x}}_y \\[6pt] d \textbf{1}_{n_c} + \textbf{\textit{x}}_c\end{array}\right)\end{equation}

Let $\textbf{\textit{x}}^\prime_1$ be the first row of $\textbf{\textit{X}}$ ; we require $\textbf{\textit{x}}^\prime_1 \textbf{\textit{n}}_3 = 0$ . We have

(29) \begin{equation} \textbf{\textit{x}}^\prime_1 \textbf{\textit{n}}_3 = (a+1) + (c-1) + (d+n_a) = a + c + d + n_a = 0\end{equation}

Here, and subsequently, brackets used in this way indicate the terms for age, year and cohort, respectively, and help to clarify the argument. Any solution of (29) will do. A convenient choice is $a = c = 0$ and $d = -n_a$ ; our candidate $\textbf{\textit{n}}_3$ is

(30) \begin{equation} \textbf{\textit{n}}_3 = \left(\begin{array}{c} \phantom{-} \textbf{\textit{x}}_a \\[6pt] - \textbf{\textit{x}}_y \\[6pt] \textbf{\textit{x}}_c - n_a \textbf{1}_{n_c}\end{array}\right)\end{equation}

We check that $\textbf{\textit{X}} \textbf{\textit{n}}_3 = \textbf{0}$ and conclude that $\{\textbf{\textit{n}}_1, \textbf{\textit{n}}_2, \textbf{\textit{n}}_3\}$ is a basis for $\mathcal{N}(\textbf{\textit{X}})$ , where $\textbf{\textit{n}}_1$ and $\textbf{\textit{n}}_2$ are defined in (25), and $\textbf{\textit{n}}_3$ in (30). Our estimates $\hat {\boldsymbol{\theta}}_s$ and $\hat {\boldsymbol{\theta}}_r$ then satisfy

(31) \begin{equation} \hat {\boldsymbol{\theta}}_s - \hat {\boldsymbol{\theta}}_r =A\left(\begin{array}{c} \phantom{-} \textbf{1}_{n_a} \\[6pt] - \textbf{1}_{n_y} \\[6pt] \phantom{-} \textbf{0}_{n_c}\end{array}\right) +B \left(\begin{array}{c} \phantom{-} \textbf{1}_{n_a} \\[6pt] \phantom{-} \textbf{0}_{n_y} \\[6pt] - \textbf{1}_{n_c}\end{array}\right) +C\left(\begin{array}{c} \phantom{-}\textbf{\textit{x}}_a \\[6pt] - \textbf{\textit{x}}_y \\[6pt] \textbf{\textit{x}}_c -n_a \textbf{1}_{n_c}\end{array}\right)\end{equation}

for some scalars A, B and C; in our example, we found $A =1.42$ , $B = 4.10$ and $C = 0.024$ . We note that (31) confirms what is widely known, namely that the age, period and cohort parameters are only estimable up to linear functions; see for example, Cairns et al. (Reference Cairns, Blake, Dowd, Coughlan, Epstein, Ong and Balevich2009).

Figure 2 Values of $\Delta \hat {\boldsymbol{\alpha}}$ , $\Delta \hat {\boldsymbol{\kappa}}$ and $\Delta \hat {\boldsymbol{\gamma}}$ defined in (27) in the APC model (23).

We use the relationships (31) to investigate forecasting with the APC model. We suppose we forecast with an autoregressive integrated moving average model ARIMA $(p, \delta, q)$ . We know from (31) that both $\hat {\boldsymbol{\kappa}}_s$ and $\hat {\boldsymbol{\kappa}}_r$ , and $\hat {\boldsymbol{\gamma}}_s$ and $\hat {\boldsymbol{\gamma}}_r$ differ by linear functions. We consider two cases: (a) forecasting the cohort effects with a simple random walk and (b) models with $\delta= 1$ or 2.

Case (a): Forecasting $\hat{\boldsymbol{\gamma}}$ with a simple random walk.

The constraints (26) imply that fitted values of $\boldsymbol{\gamma}$ have mean and slope of zero, and hence the simple random walk is a plausible model for forecasting cohort effects. We suppose that both $\hat{\boldsymbol{\gamma}}_s$ and $\hat{\boldsymbol{\gamma}}_r$ are forecast in this way. Then a straightforward computation shows that the resulting forecasts of mortality are not equal, i.e., are not invariant with respect to the choice of constraints. The reason for this is that the forecasts of $\hat{\boldsymbol{\gamma}}_s$ and $\hat{\boldsymbol{\gamma}}_r$ are now parallel and the linear relationship between $\hat{\boldsymbol{\gamma}}_s$ and $\hat{\boldsymbol{\gamma}}_r$ in (31) has been broken.

We can use a simple device to illustrate further forecasting cohort effects with a simple random walk. If we fit an ARIMA(0,1,0) model, i.e., a random walk with drift, to $\hat{\boldsymbol{\gamma}}_s$ , then the estimate of the drift parameter, $\mu_s$ say, is

(32) \begin{equation} \hat \mu_s = {\hat {\gamma_s}(n_c) - \hat {\gamma_s}(1) \over n_c -1}\end{equation}

We replace the standard constraints in (26) with

(33) \begin{equation} \sum_1^{n_y} \kappa_y = \sum_1^{n_c} \gamma_c = \gamma_{n_c} - \gamma_1 = 0\end{equation}

We suppose the resulting estimates are denoted $\hat{\boldsymbol{\alpha}}_z$ $\hat{\boldsymbol{\kappa}}_z$ and $\hat{\boldsymbol{\gamma}}_z$ . Now when we fit the ARIMA(0,1,0) model to $\hat{\boldsymbol{\gamma}}_z$ , the estimate of the drift parameter is constrained to be zero by (33), i.e., we have fitted a simple random walk. We must forecast $\hat{\boldsymbol{\gamma}}_s$ with the same ARIMA(0,1,0) model; here, the estimate of the drift parameter is not zero, but the forecasts of mortality under both constraint systems are equal.

Case (b): Forecasting with an ARIMA( $p, \delta, q$ ) model, $\delta = 1,2$ .

We know from (31) that

(34) \begin{equation} \hat {\boldsymbol{\kappa}}_s - \hat {\boldsymbol{\kappa}}_r = -A\textbf{1}_{n_y} - C\textbf{\textit{x}}_y\end{equation}

Thus $\hat {\boldsymbol{\kappa}}_s$ and $\hat {\boldsymbol{\kappa}}_r$ differ by a linear function. We should expect that any reasonable forecasts of $\hat {\boldsymbol{\kappa}}_s$ and $\hat {\boldsymbol{\kappa}}_r$ obey the same functional relation. We show that this is indeed the case when an ARIMA( $p, \delta, q$ ) model with a mean, $\delta = 1,2$ , is used to forecast. A formal proof of this is given in Appendix C. In summary, we will only consider forecasting models that obey this functional relationship since, as we shall see, this will preserve the invariance of forecasts with respect to the choice of constraints.

Let $\hat {\boldsymbol{\kappa}}_{s, f}$ and $\hat {\boldsymbol{\kappa}}_{r, f}$ denote the forecast values of $\hat {\boldsymbol{\kappa}}_s$ and $\hat {\boldsymbol{\kappa}}_r$ , respectively. It follows from (34) and (C.5) in Appendix C with $a = -A$ and $b = -C$ that

(35) \begin{equation} \hat {\boldsymbol{\kappa}}_{s, f} -\hat {\boldsymbol{\kappa}}_{r, f} = -(A + n_yC) \textbf{1}_{n_f} - C \textbf{\textit{x}}_f\end{equation}

where $\textbf{\textit{x}}_f = (1,\ldots, n_f)^\prime$ . The analogous argument with $\hat {\boldsymbol{\gamma}}_s - \hat {\boldsymbol{\gamma}}_r$ shows that

(36) \begin{equation} \hat {\boldsymbol{\gamma}}_{s, f} -\hat {\boldsymbol{\gamma}}_{r, f} = -(B + (n_y-1)C) \textbf{1}_{n_f} + C \textbf{\textit{x}}_f\end{equation}

where $\hat {\boldsymbol{\gamma}}_{s, f}$ and $\hat {\boldsymbol{\gamma}}_{r, f}$ are the forecast values of $\hat {\boldsymbol{\gamma}}_s$ and $\hat {\boldsymbol{\gamma}}_r$ , respectively.

We can now establish the invariance of the forecast values of $\log\mu$ . Define

\begin{equation*} \hat{ \boldsymbol{\theta}}_{s, f} = (\hat{\boldsymbol{\alpha}}_s^\prime, \hat{\boldsymbol{\kappa}}_s^\prime, \hat{\boldsymbol{\kappa}}_{s, f}^\prime, \hat{\boldsymbol{\gamma}}_s^\prime, \hat{\boldsymbol{\gamma}}_{s, f}^\prime)^\prime\end{equation*}

with a similar definition for $\hat{ \boldsymbol{\theta}}_{r, f}$ . Let $\textbf{\textit{X}}_f$ be the model matrix for the APC model for ages $\textbf{\textit{x}}_a$ and years $(\textbf{\textit{x}}_y^\prime, \textbf{\textit{x}}_{y, f}^\prime)^\prime$ , where $\textbf{\textit{x}}_{y, f} =(n_y + 1, \ldots, n_y + n_f)^\prime$ . We can show that $\textbf{\textit{X}}_f(\hat{ \boldsymbol{\theta}}_{s, f} -\hat{ \boldsymbol{\theta}}_{r, f}) = \textbf{0}$ . We omit the proof and instead give a detailed proof in the next section for a smooth version of the APCI model; the present case follows in the same fashion. We can therefore conclude that the fitted and forecast values of mortality in the APC model are invariant with respect to the choice of constraints when an ARIMA $(p, \delta, q)$ model, $\delta =1$ or 2, with a mean is used to forecast.

A “computer proof” of the above result is simple and for many models of mortality is all that is reasonably available or indeed required. We compute fitted and forecast values under any two constraint systems and check the invariance of the fitted and forecast values. Once the basic code is written, different constraint systems and different forecasting regimes are easily applied. We have used this approach in parallel with our theoretical discussion.

4.2.2 Smooth APC model

Just as in the AP model, forecasting of mortality in the APC model is improved if the age parameters $\boldsymbol{\alpha}$ are smoothed. We set $\boldsymbol{\alpha} = \textbf{\textit{B}}_a \textbf{\textit{a}}$ , where $\textbf{\textit{B}}_a$ is $n_a \times c_a$ , and replace $\textbf{\textit{X}}_a$ in (24) with $\textbf{1}_{n_y} \otimes \textbf{\textit{B}}_a$ , as in (21). The model matrix becomes

(37) \begin{equation} \textbf{\textit{X}} = [\textbf{1}_{n_y} \otimes \textbf{\textit{B}}_a\,{:}\,\textbf{\textit{I}}_{n_y} \otimes \textbf{1}_{n_a}\,{:}\,\textbf{\textit{X}}_c],\; n_an_y \times (c_a+n_a+2n_y-1)\end{equation}

and the penalty matrix is

(38) \begin{equation} \textbf{\textit{P}} = \mbox{blockdiag}\{\lambda_a \textbf{\textit{D}}^\prime \textbf{\textit{D}}, \textbf{0}\}\end{equation}

where $\textbf{\textit{D}}$ is the difference matrix of order d acting on the smoothed age coefficients $\textbf{\textit{a}}$ , $\textbf{0}$ is the $(n_y+n_c) \times(n_y+n_c)$ matrix of 0s acting on the unsmoothed year and cohort coefficients $\boldsymbol{\kappa}$ and $\boldsymbol{\gamma}$ , and $\lambda_a$ is the smoothing parameter. With the P-spline system, it is usual to smooth with a second-order penalty. However, if a first-order penalty is used, we find the effective rank of the model is increased by one; in this case, in order to retain invariance the number of constraints should be reduced to two. We continue with a detailed discussion of the case when a second-order penalty is used.

The regression parameter is $\boldsymbol{\theta} = (\textbf{\textit{a}}^\prime, \boldsymbol{\kappa}^\prime, \boldsymbol{\gamma}^\prime)^\prime$ . It is easy to see that

(39) \begin{equation} \textbf{\textit{n}}_1 = \left(\begin{array}{c} \phantom{-}\textbf{1}_{c_a} \\[6pt] -\textbf{1}_{n_y} \\[6pt] \phantom{-}\textbf{0}_{n_c} \end{array}\right), \quad \textbf{\textit{n}}_2 = \left(\begin{array}{c} \phantom{-}\textbf{1}_{c_a} \\[6pt] \phantom{-}\textbf{0}_{n_y} \\[6pt] -\textbf{1}_{n_c} \end{array}\right)\end{equation}

are in $\mathcal{N}(\textbf{\textit{X}})$ . We use an extension of the argument used in the unsmoothed case to find a suitable third basis vector, $\textbf{\textit{n}}_3$ . We fit with the standard constraints (26) and also with some random constraints but with $\lambda_a = 0$ (since we are interested in determining a basis for $\mathcal{N}(\textbf{\textit{X}})$ ). Corresponding to (27), we define

(40) \begin{equation} \Delta \hat {\textbf{\textit{a}}} = \hat {\textbf{\textit{a}}}_s - \hat {\textbf{\textit{a}}}_r, \; \Delta \hat {\boldsymbol{\kappa}} = \hat {\boldsymbol{\kappa}}_s - \hat {\boldsymbol{\kappa}}_r, \; \Delta \hat {\boldsymbol{\gamma}} = \hat {\boldsymbol{\gamma}}_s - \hat {\boldsymbol{\gamma}}_r\end{equation}

In our example, we used a knot spacing of $\Delta_a = 5$ . From Figure 3, we note that $\Delta \hat {\boldsymbol{\kappa}}$ and $\Delta \hat {\boldsymbol{\gamma}}$ are linear with slopes 0.1754 and $-0.1754$ , respectively, while $\Delta \hat {\boldsymbol{\alpha}}$ is linear with slope $-0.8768 \approx -\Delta_a 0.1754$ . Hence, as in (28), we conjecture that $\textbf{\textit{n}}_3$ has the form

\begin{equation*} \textbf{\textit{n}}_3 = \left(\begin{array}{c} a \textbf{1}_{c_a} + \Delta_a \textbf{\textit{x}}_{c_a} \\[6pt] c \textbf{1}_{n_y} - \textbf{\textit{x}}_y \\[6pt] d \textbf{1}_{n_c} + \textbf{\textit{x}}_c\end{array}\right)\end{equation*}

where $\textbf{\textit{x}}_{c_a} = (1,2,3,\ldots, c_a)^\prime$ . Let $\textbf{\textit{x}}^\prime_1$ be the first row of $\textbf{\textit{X}}$ ; we set $\textbf{\textit{x}}^\prime_1 \textbf{\textit{n}}_3 = 0$ . We recall that in the computation of $\textbf{\textit{B}}_a$ , we place a knot at the first age. With this assumption, the first row of $\textbf{\textit{B}}_a$ is $(\frac{1}{6},\frac{2}{3},\frac{1}{6}, {\boldsymbol{0}}^{\prime}_{c_{a}-3})$ . We have

\begin{equation*} \textbf{\textit{x}}^\prime_1 \textbf{\textit{n}}_3 = \left(a+\Delta_a\!\left(\mbox{${1\over6}$} + \mbox{${4\over3}$} + \mbox{${3\over6}$}\right)\right) + (c-1) + (d+n_a) = a + c + d + n_a + 2\Delta_a - 1 = 0.\end{equation*}

A convenient solution is $a = c = 0$ and $d = 1 - n_a - 2\Delta_a$ ; our candidate $\textbf{\textit{n}}_3$ is

(41) \begin{equation} \textbf{\textit{n}}_3 = \left(\begin{array}{c} \Delta_a \textbf{\textit{x}}_{c_a} \\[6pt] - \textbf{\textit{x}}_y \\[6pt] \textbf{\textit{x}}_c - \omega_c \textbf{1}_{n_c}\end{array}\right)\end{equation}

where $\omega_c = -d = n_a + 2\Delta_a -1$ . We now check that $\textbf{\textit{X}} \textbf{\textit{n}}_3 = \textbf{0}$ . We conclude that $\{\textbf{\textit{n}}_1, \textbf{\textit{n}}_2, \textbf{\textit{n}}_3\}$ is a basis for $\mathcal{N}(\textbf{\textit{X}})$ , where $\textbf{\textit{n}}_1$ and $\textbf{\textit{n}}_2$ are defined in (39) and $\textbf{\textit{n}}_3$ in (41). Hence,

(42) \begin{equation}\hat {\boldsymbol{\theta}}_s - \hat {\boldsymbol{\theta}}_r = A \left(\begin{array}{c} \phantom{-}\textbf{1}_{c_a} \\[6pt] -\textbf{1}_{n_y} \\[6pt] \phantom{-}\textbf{0}_{n_c} \end{array}\right) + B \left(\begin{array}{c} \phantom{-}\textbf{1}_{c_a} \\[6pt] \phantom{-}\textbf{0}_{n_y} \\[6pt] -\textbf{1}_{n_c} \end{array}\right) + C \left(\begin{array}{c} \Delta_a \textbf{\textit{x}}_{c_a} \\[6pt] -\textbf{\textit{x}}_y \\[6pt] \textbf{\textit{x}}_c - \omega_c \textbf{1}_{n_c} \end{array}\right)\end{equation}

for some A, B and C; in our example, we found $A = 5.298$ , $B =2.644$ and $C = -0.1754$ . Pre-multiplying (42) through by $\mbox{blockdiag}\{\textbf{\textit{B}}_a, \textbf{\textit{I}}_{n_y}, \textbf{\textit{I}}_{n_c}\}$ and using $\textbf{\textit{B}}_a \textbf{1}_{c_a} = \textbf{1}_{n_a}$ and $\textbf{\textit{B}}_a \Delta_a \textbf{\textit{x}}_{c_a}= \textbf{\textit{x}}_{n_a} + \omega_a \textbf{1}_{n_a}$ , $\omega_a = 2\Delta_a - 1$ , we find

(43) \begin{equation}\begin{array}{c} \phantom{0} \\[6pt] \hat{\boldsymbol{\theta}}_s^* \\[6pt] \phantom{0}\end{array} -\begin{array}{c} \phantom{0} \\[6pt] \hat{\boldsymbol{\theta}}_r^* \\[6pt] \phantom{0}\end{array} =A \left(\begin{array}{c} \phantom{-}\textbf{1}_{n_a} \\[6pt] -\textbf{1}_{n_y} \\[6pt] \phantom{-}\textbf{0}_{n_c} \end{array}\right) + B \left(\begin{array}{c} \phantom{-}\textbf{1}_{n_a} \\[6pt] \phantom{-}\textbf{0}_{n_y} \\[6pt] -\textbf{1}_{n_c} \end{array}\right) + C \left(\begin{array}{c} \textbf{\textit{x}}_a + \omega_a \textbf{1}_{n_a} \\[6pt] -\textbf{\textit{x}}_y \\[6pt] \textbf{\textit{x}}_c - \omega_c \textbf{1}_{n_c} \end{array}\right)\end{equation}

where $\boldsymbol{\theta}^* = (\boldsymbol{\alpha}^\prime, \boldsymbol{\kappa}^\prime, \boldsymbol{\gamma}^\prime)^\prime$ denotes the parameters of interest.

Figure 3 Values of $\Delta \hat {\textbf{\textit{a}}}$ , $\Delta \hat {\boldsymbol{\kappa}}$ and $\Delta \hat {\boldsymbol{\gamma}}$ defined in (40) in the smooth APC model (37) but with $\lambda_a = 0$ .

We show that forecasting with an ARIMA( $p, \delta, q)$ , $\delta = 1$ or 2, is invariant with respect to the choice of constraints. We use the same notation as in the unsmoothed case. Comparing (31) and (43), we see that the middle rows are identical. Hence, the forecasts of $\boldsymbol{\kappa}$ under two different constraint systems in the smoothed case obey the same relation (35) as in the unsmoothed case, i.e.,

\begin{equation*} \hat {\boldsymbol{\kappa}}_{s, f} - \hat {\boldsymbol{\kappa}}_{r, f} = -(A + n_y C) \textbf{1}_{n_f} - C \textbf{\textit{x}}_f\end{equation*}

Forecasts of $\boldsymbol{\gamma}$ obey the following

\begin{eqnarray*} \nonumber \hat {\boldsymbol{\gamma}}_{s, f} - \hat {\boldsymbol{\gamma}}_{r, f} &=& \left(\hat {\boldsymbol{\gamma}}_s(n_y) - \hat {\boldsymbol{\gamma}}_r(n_y)\right) \textbf{1}_{n_f} + C\textbf{\textit{x}}_{n_f} \\ &=& (\!-B +(n_c - \omega_c)C) \textbf{1}_{n_f} + C \textbf{\textit{x}}_{n_f} \quad \mbox{by (\linkref{disp43}{43})} \\ &=& (\!-B +(n_y - 2\Delta_a)C) \textbf{1}_{n_f} + C \textbf{\textit{x}}_{n_f}\end{eqnarray*}

It now follows that $\textbf{\textit{X}}_f^* (\hat{\boldsymbol{\theta}}_{s, f}^{\raisebox{-3pt}{\hbox{\fontsize{8}{10}\selectfont$\ast$}}} - \hat {\boldsymbol{\theta}}_{r, f}^{{\raisebox{-3pt}{\hbox{\fontsize{8}{10}\selectfont$\ast$}}}}) = \textbf{0}$ , where $\textbf{\textit{X}}_f^*$ is the model matrix (24) in the unsmoothed case; the proof is a special case of that given for the smooth APCI model. We have established the invariance of the fitted and forecast values with respect to the choice of constraint system when an ARIMA( $p, \delta, q$ ) model, $\delta =1$ or 2, is used.

4.3.1 APCI model

The APCI model is

(44) \begin{equation} \log\,\mu_{x, y} = \alpha_x + \kappa_y + \gamma_{c(x, y)} + \beta_x(\bar y - y),\;x=1,\ldots,n_a,\;y = 1, \ldots, n_y\end{equation}

where $\bar y = (n_y+1)/2$ (recall the year vector is $\textbf{\textit{x}}_y =(1, \ldots, n_y)^\prime$ ). We note that (44) differs in a minor way from the CMI’s formulation; we have reversed the sign of the $\beta_x$ terms, which now correspond to the $\beta_x$ terms in the Lee–Carter model. The thinking behind the APCI model is as follows: linearise the $\beta_x \kappa_y$ term in the Lee–Carter model and add the cohort term $\gamma_{c(x, y)}$ . This gives a model which takes into account any cohort effects and allows the forecast of the year effect to depend on age. The APCI model is a GLM with model matrix

(45) \begin{equation} \textbf{\textit{X}} = [\textbf{\textit{X}}_a\, :\, \textbf{\textit{X}}_y \, :\, \textbf{\textit{X}}_c \, :\, \textbf{\textit{X}}_b],\; n_an_y \times (3n_a + 2n_y - 1)\end{equation}

where $\textbf{\textit{X}}_a$ and $\textbf{\textit{X}}_y$ are defined in (17), $\textbf{\textit{X}}_c$ in (24) and $\textbf{\textit{X}}_b = (\bar y \textbf{1}_{n_y} - \textbf{\textit{x}}_y)\otimes \textbf{\textit{I}}_{n_a}$ . We define $\boldsymbol{\theta} = (\boldsymbol{\alpha}^\prime, \boldsymbol{\kappa}^\prime, \boldsymbol{\gamma}^\prime, \boldsymbol{\beta}^\prime)^\prime$ and the model has the form (1). The model matrix has $3n_a+ 2n_y - 1$ columns and rank $3n_a + 2n_y - 6$ , so five constraints are required to enable $\boldsymbol{\theta}$ to be estimated uniquely. The null space of $\textbf{\textit{X}}$ is spanned by

(46) \begin{equation} \left(\begin{array}{c} \phantom{-}\textbf{1}_{n_a} \\[6pt] -\textbf{1}_{n_y} \\[6pt] \phantom{-}\textbf{0}_{n_c} \\[6pt] \phantom{-}\textbf{0}_{n_a} \end{array}\right),\; \left(\begin{array}{c} \phantom{-}\textbf{1}_{n_a} \\[6pt] \phantom{-}\textbf{0}_{n_y} \\[6pt] -\textbf{1}_{n_c} \\[6pt] \phantom{-}\textbf{0}_{n_a} \end{array} \right),\; \left( \begin{array}{c} \phantom{-}\textbf{\textit{x}}_a \\[6pt] -\textbf{\textit{x}}_y \\[6pt] \textbf{\textit{x}}_c - n_a \textbf{1}_{n_c} \\[6pt] \phantom{-} \textbf{0}_{n_a} \end{array} \right),\; \left(\begin{array}{c} \textbf{0}_{n_a} \\[6pt] \textbf{\textit{x}}_y - \bar y \textbf{1}_{n_y} \\[6pt] \textbf{0}_{n_c} \\[6pt] \textbf{1}_{n_a} \end{array} \right),\; \left(\begin{array}{c} Q \textbf{\textit{x}}_a + \textbf{\textit{x}}_a^2 \\[6pt] 2n_a \textbf{\textit{x}}_y + \textbf{\textit{x}}_y^2 \\[6pt] n_a^2 \textbf{1}_{n_c} - \textbf{\textit{x}}_c^2 \\[6pt] 2 \textbf{\textit{x}}_a \end{array}\right)\end{equation}

where $Q = -2(n_a + \bar y)$ . The first three basis vectors in (46) are the same as the basis vectors for the APC model with $\textbf{0}_{n_a}$ appended. The fourth basis vector can be deduced by observing that

\begin{equation*} \kappa_y + (\bar y - y) \beta_x = \kappa_y - (\bar y - y) \omega + (\bar y - y)(\beta_x + \omega) = \kappa_y^* + (\bar y - y) \beta_x^*\end{equation*}

for any scalar $\omega$ . The fifth basis vector can be found with the same method as used to find the third basis vector for the APC model. We denote the basis vectors in (46) by $\textbf{\textit{n}}_i,\,i=1,\ldots,5$ .

We comment briefly on the CMI’s approach to identifiability. The CMI (2016a, 7.3) gives the following transformation of the parameters which leaves the values of $\log \boldsymbol{\mu}$ unchanged. We give this transformation in our notation and with the CMI’s $\beta_x$ replaced by ours.

\begin{multline} \nonumber \left(\begin{array}{c} \Delta \alpha_x \\[6pt] \Delta \kappa_y \\[6pt] \Delta \gamma_{y-x} \\[6pt] \Delta \beta_x \end{array}\right)= \theta_1 \left(\begin{array}{c} \phantom{-} 1 \\[6pt] \phantom{-} 0 \\[6pt] -1 \\[6pt] \phantom{-} 0 \end{array}\right)+ \theta_2 \left(\begin{array}{c} - (x - \bar x) \\[6pt] y - \bar y \\[6pt] (x - \bar x) - (y - \bar y) \\[6pt] 0 \end{array}\right)+ \\[6pt] \theta_3 \left(\begin{array}{c} (x - \bar x)^2 \\[6pt] (y - \bar y)^2 \\[6pt] - \left((x - \bar x) - (y - \bar y)\right)^2 \\[6pt] 2(x - \bar x) \end{array}\right) +\theta_4 \left(\begin{array}{c} \phantom{-} 1 \\[6pt] -1 \\[6pt] \phantom{-} 0 \\[6pt] \phantom{-} 0 \end{array}\right)+ \theta_5 \left(\begin{array}{c} 0 \\[6pt] -(y - \bar y) \\[6pt] 0 \\[6pt] (y - \bar y) \end{array}\right)\end{multline}

Writing the transformation in this way, we see that it is equivalent to a particular null basis for $\textbf{\textit{X}}$ , namely,

(47) \begin{align} \left(\begin{array}{c} \phantom{-}\textbf{1}_{n_a} \\[6pt] \phantom{-}\textbf{0}_{n_y} \\[6pt] -\textbf{1}_{n_c} \\[6pt] \phantom{-}\textbf{0}_{n_a} \end{array}\right),\; \left( \begin{array}{c} -(\textbf{\textit{x}}_a - \bar a \textbf{1}_{n_a}) \\[6pt] \textbf{\textit{x}}_y - \bar y \textbf{1}_{n_y} \\[6pt] -(\textbf{\textit{x}}_c - \bar c \textbf{1}_{n_c}) \\[6pt] \phantom{-} \textbf{0}_{n_a} \end{array} \right),\; \left(\begin{array}{c} (\textbf{\textit{x}}_a - \bar a \textbf{1}_{n_a})^2 \\[6pt] (\textbf{\textit{x}}_y - \bar y \textbf{1}_{n_y})^2 \\[6pt] -(\textbf{\textit{x}}_c - \bar c \textbf{1}_{n_c})^2 \\[6pt] 2(\textbf{\textit{x}}_a - \bar a \textbf{1}_{n_a}) \end{array} \right),\; \left(\begin{array}{c} \phantom{-}\textbf{1}_{n_a} \\[6pt] -\textbf{1}_{n_y} \\[6pt] \phantom{-} \textbf{0}_{n_c} \\[6pt] \phantom{-}\textbf{0}_{n_a} \end{array} \right),\; \left(\begin{array}{c} \textbf{0}_{n_a} \\[6pt] \textbf{\textit{x}}_y - \bar y \textbf{1}_{n_y} \\[6pt] \textbf{0}_{n_c} \\[6pt] \textbf{1}_{n_a} \end{array}\right)\end{align}

where $\bar a$ is the mean age, $\bar y$ is the mean year and $\bar c= n_a - \bar a + \bar y$ . Three of the basis vectors are the same as ours; the other two have the same form. In particular, the third basis vector in (47) has quadratic and linear functions in the equivalent positions to (46).

We fit the model under the standard constraints (CMI, 2016a)

(48) \begin{equation} \sum \kappa_y = \sum \gamma_c = \sum c \gamma_c = \sum c^2 \gamma_c = \sum y \kappa_y = 0\end{equation}

where c is the cohort index, $c = 1, \ldots, n_c$ ; again, we could use weighted constraints as in Richards et al. (to appear). We also use the random constraints

(49) \begin{equation} \sum u_{i, j} \theta_j = 0, \, i = 1,\ldots, 5,\, j = 1,\ldots, 3n_a+2n_y -1\end{equation}

where the $u_{i, j}$ are independent $\mathcal{U}(0,1)$ . With our usual notation for the estimates of $\boldsymbol{\theta}$ under the standard and random constraints, we have

\begin{equation*} \hat {\boldsymbol{\theta}}_s - \hat {\boldsymbol{\theta}}_r = A \textbf{\textit{n}}_1 + B \textbf{\textit{n}}_2 + C \textbf{\textit{n}}_3 + D \textbf{\textit{n}}_4 + E \textbf{\textit{n}}_5\end{equation*}

for some A, B, C, D and E. Equating coefficients in (46), we immediately find that $\hat {\boldsymbol{\alpha}}_s -\hat {\boldsymbol{\alpha}}_r$ , $\hat {\boldsymbol{\kappa}}_s - \hat {\boldsymbol{\kappa}}_r$ and $\hat {\boldsymbol{\gamma}}_s - \hat {\boldsymbol{\gamma}}_r$ are quadratic functions of $\textbf{\textit{x}}_a$ , $\textbf{\textit{x}}_y$ and $\textbf{\textit{x}}_c$ , respectively, with E, E and $-E$ as the coefficients of the quadratic terms; furthermore, $\hat {\boldsymbol{\beta}}_s - \hat {\boldsymbol{\beta}}_r$ is linear with slope 2E. These relationships have implications for invariance when it comes to forecasting.

Figure 4 confirms the relationships among the coefficients. Since $\Delta \hat {\boldsymbol{\alpha}} = \hat {\boldsymbol{\alpha}}_s -\hat {\boldsymbol{\alpha}}_r$ is quadratic with the coefficient of the quadratic term equal to E, first differences of $\Delta \hat {\boldsymbol{\alpha}}$ , $D(\Delta \hat {\boldsymbol{\alpha}})$ , will be linear with slope 2E as shown in Figure 4; the same remark applies to the first differences $D(\Delta \hat {\boldsymbol{\kappa}})$ and $D(\Delta \hat {\boldsymbol{\gamma}})$ . The slope of $\Delta \hat {\boldsymbol{\beta}}$ is 2E. In our example, we found $A = -624.3$ , $B = -70.02$ , $C=18.24$ , $D = 20.32$ and $E = 0.002385$ .

Figure 4 Plot of first differences, denoted D, of $\Delta \hat {\boldsymbol{\alpha}} = \hat {\boldsymbol{\alpha}}_s - \hat {\boldsymbol{\alpha}}_r$ , $\Delta \hat {\boldsymbol{\kappa}} = \hat {\boldsymbol{\kappa}}_s - \hat {\boldsymbol{\kappa}}_r$ and $\Delta \hat {\boldsymbol{\gamma}} = \hat {\boldsymbol{\gamma}}_s - \hat {\boldsymbol{\gamma}}_r$ ; plot of $\Delta \hat {\boldsymbol{\beta}} = \hat {\boldsymbol{\beta}}_s- \hat {\boldsymbol{\beta}}_r$ in the APCI model (44).

Figure 5 Left: plot of $\hat {\boldsymbol{\alpha}}_s$ , $\hat {\boldsymbol{\beta}}_s$ ; right: plot of $\hat {\boldsymbol{\alpha}}_r$ and $\hat {\boldsymbol{\beta}}_r$ in the APCI model (44).

Figure 5 is a plot of the estimated $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$ under our two constraint systems; the estimates are wildly different both in shape and range. Under the standard constraint, the estimate of $\boldsymbol{\alpha}$ has a familiar look, while the estimate of $\boldsymbol{\beta}$ is very close to the estimate of $\boldsymbol{\beta}$ in the Lee–Carter model; see Figure 8 for the corresponding plot. However, Figure 5 serves to emphasise our main point: it is how the coefficients join forces that matter, not the individual coefficients themselves.

The CMI smooth both $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$ in their treatment of the APCI model, and we will focus on this case. First, we discuss briefly the unsmoothed case. We forecast with our two constraint systems, namely standard and random. The quadratic nature of the fifth basis vector in the null space of the APCI model (see (46)) implies that a first-order ARIMA model will not give invariant forecasts. In this paper, we concentrate on forecasting models which lead to invariant forecasts so we do not pursue this further.

4.3.2 Smooth APCI model

The case for smoothing both $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$ is strong: forecasts will be much less prone to irregular behaviour by age (see Delwarde et al., Reference Delwarde, Denuit and Eilers2007; Currie, Reference Currie2013). We will not smooth either $\boldsymbol{\kappa}$ or $\boldsymbol{\gamma}$ , since we are interested in stochastic forecasts of these parameters. This is in contrast to the CMI which in addition to smoothing $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$ also smooth $\boldsymbol{\kappa}$ or $\boldsymbol{\gamma}$ ; the CMI’s approach is designed to facilitate deterministic targeting of future mortality. We do not comment on this debate here; see Booth & Tickle (Reference Booth and Tickle2008), Richards et al. (to appear).

Let $\textbf{\textit{B}}_a,\,n_a \times c_a$ , be a regression matrix evaluated over a basis of cubic B-splines for age. We use the same regression matrix to smooth both $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$ and set $\boldsymbol{\alpha} = \textbf{\textit{B}}_a \textbf{\textit{a}}$ and $\boldsymbol{\beta} = \textbf{\textit{B}}_a \textbf{\textit{b}}$ . The model matrix is now

(50) \begin{equation} \textbf{\textit{X}} = [\textbf{1}_{n_y} \otimes \textbf{\textit{B}}_a\, : \,\textbf{\textit{I}}_{n_y} \otimes \textbf{1}_{n_a}\, : \,\textbf{\textit{X}}_c\, : \,(\bar y \textbf{1}_{n_y} - \textbf{\textit{x}}_y) \otimes \textbf{\textit{B}}_a]\end{equation}

where $\textbf{\textit{X}}_c$ is defined below (24). The model matrix is $n_a\,n_y \times (2c_a + n_y + n_c)$ with $r(\textbf{\textit{X}}) = 2c_a + n_y +n_c -5$ , so $\mathcal{N}(\textbf{\textit{X}})$ has dimension five. We fit model (50) but without smoothing, i.e., with the smoothing parameters set to zero. Let $\boldsymbol{\theta} = (\textbf{\textit{a}}^\prime,\boldsymbol{\kappa}^\prime, \boldsymbol{\gamma}^\prime, \textbf{\textit{b}}^\prime)^\prime$ be the vector of regression coefficients and $\boldsymbol{\theta}^* = (\boldsymbol{\alpha}^\prime, \boldsymbol{\kappa}^\prime, \boldsymbol{\gamma}^\prime, \boldsymbol{\beta}^\prime)^\prime$ denote the parameters of interest. We denote the five vectors

(51) \begin{align} \left(\begin{array}{c} \phantom{-}\textbf{1}_{c_a} \\[6pt] -\textbf{1}_{n_y} \\[6pt] \phantom{-}\textbf{0}_{n_c} \\[6pt] \phantom{-}\textbf{0}_{c_a} \end{array}\right),\; \left(\begin{array}{c} \phantom{-}\textbf{1}_{c_a} \\[6pt] \phantom{-}\textbf{0}_{n_y} \\[6pt] -\textbf{1}_{n_c} \\[6pt] \phantom{-}\textbf{0}_{c_a} \end{array} \right),\; \left( \begin{array}{c} \Delta_a \textbf{\textit{x}}_{c_a} \\[6pt] -\textbf{\textit{x}}_y \\[6pt] \textbf{\textit{x}}_c - \omega_c \textbf{1}_{n_c} \\[6pt] \phantom{-} \textbf{0}_{c_a} \end{array} \right),\; \left(\begin{array}{c} \textbf{0}_{c_a} \\[6pt] \textbf{\textit{x}}_y - \bar y \textbf{1}_{n_y} \\[6pt] \textbf{0}_{n_c} \\[6pt] \textbf{1}_{c_a} \end{array} \right),\; \left(\begin{array}{c} Q_a(\textbf{\textit{x}}_{c_a}) \\[6pt] Q_\kappa(\textbf{\textit{x}}_{n_y}) \\[6pt] Q_\gamma(\textbf{\textit{x}}_{n_c}) \\[6pt] L_b(\textbf{\textit{x}}_{c_a}) \end{array} \right)\end{align}

by $\textbf{\textit{n}}_1, \ldots,\textbf{\textit{n}}_5$ , respectively; here $Q_a(\cdot)$ , $Q_\kappa(\cdot)$ and $Q_\gamma(\cdot)$ are quadratic functions and $L_b(\cdot)$ is a linear function. We can check that $\textbf{\textit{n}}_1,\ldots,\textbf{\textit{n}}_4$ are all in $\mathcal{N}(\textbf{\textit{X}})$ . Figure 6 is a plot of differences in estimates under the standard and random constraints (48) and (49) and suggests that the fifth basis vector in $\mathcal{N}(\textbf{\textit{X}})$ has the form $\textbf{\textit{n}}_5$ . Additionally, the form of the fifth basis in (46) in the unsmoothed case supports this idea. We will not need an exact expression for $\textbf{\textit{n}}_5$ .

Figure 6 Plot of first differences, denoted D, of $\Delta \hat {\textbf{\textit{a}}} = \hat {\textbf{\textit{a}}}_s - \hat {\textbf{\textit{a}}}_r$ , $\Delta \hat {\boldsymbol{\kappa}} = \hat {\boldsymbol{\kappa}}_s - \hat {\boldsymbol{\kappa}}_r$ and $\Delta \hat {\boldsymbol{\gamma}} = \hat {\boldsymbol{\gamma}}_s - \hat {\boldsymbol{\gamma}}_r$ ; plot of $\Delta \hat {\textbf{\textit{b}}} = \hat {\textbf{\textit{b}}}_s - \hat {\textbf{\textit{b}}}_r$ in APCI model with the smooth model matrix (50) but with $\lambda_a = \lambda_b = 0$ .

Let $\textbf{\textit{D}}_a$ and $\textbf{\textit{D}}_b$ be difference matrices of orders $d_a$ and $d_b$ applied to the coefficients $\textbf{\textit{a}}$ and $\textbf{\textit{b}}$ , respectively. The penalty matrix is

\begin{equation*} \textbf{\textit{P}} = \mbox{blockdiag}\{\lambda_a \textbf{\textit{D}}_a^\prime \textbf{\textit{D}}_a,\, \textbf{0},\, \lambda_b \textbf{\textit{D}}_b^\prime \textbf{\textit{D}}_b\}\end{equation*}

where $\textbf{0}$ is a square matrix of 0s of size $n_y+n_c$ . The null space of the penalised model is given by $\mathcal{N}(\textbf{\textit{X}}) \cap \mathcal{N}(\textbf{\textit{P}})$ ; see (13). Now $\textbf{\textit{P}} \textbf{\textit{n}}_i = \textbf{0}$ for $i = 1,\ldots, 4$ for $d_a \ge 2$ and $d_b \ge 2$ ; for $\textbf{\textit{P}} \textbf{\textit{n}}_5 = \textbf{0}$ , we need $d_a \ge 3$ and $d_b \ge 2$ . It is usual to smooth with difference matrices of order two; in this case $\textbf{\textit{n}}_5\notin \mathcal{N}(\textbf{\textit{P}})$ . We continue with the case $d_a = d_b = 2$ , in which case the null space of the model is given by the first four vectors in (51). Direct computation of $r(\textbf{\textit{X}}^\prime \tilde {\textbf{\textit{W}}} \textbf{\textit{X}} + \textbf{\textit{P}})$ confirms that the effective rank of the smooth APCI model with $d_a = d_b = 2$ is $2c_a + n_y +n_c - 4$ . We conclude that $\{\textbf{\textit{n}}_1, \ldots,\textbf{\textit{n}}_4\}$ is a basis for the null space of the smooth model.

We make a brief comment on the effect of the choice of order of the penalty. We take $d_a =\break d_b = 2$ . Under the P-spline system of smoothing, the rank of the model is in effect $r(\textbf{\textit{X}}^\prime \textbf{\textit{X}} +\textbf{\textit{P}})$ , while the rank of the model without smoothing is $r(\textbf{\textit{X}}^\prime \textbf{\textit{X}})$ ; denote these ranks by $r_s$ and r, respectively. Usually, $r = r_s$ ; indeed we are not aware that the phenomenon $r_s >r$ has been previously observed. Certainly, the implicit constraint in the smooth model can lead to very different estimates of the parameters. We argue that since $\boldsymbol{\alpha}$ , $\boldsymbol{\kappa}$ , $\boldsymbol{\gamma}$ and $\boldsymbol{\beta}$ are not identifiable, we should not be concerned about these differences. The fitted force of mortalities (which are identifiable) under the two models will be very close.

We fit the model with four standard constraints (in (48) we drop the constraint $\sum c^2 \gamma_c= 0$ ) and four random constraints. Let $\hat {\boldsymbol{\theta}_s}$ and $\hat {\boldsymbol{\theta}_r}$ be the estimates of $\boldsymbol{\theta}$ as usual. Then, from (51), we have

(52) \begin{equation}\hat {\boldsymbol{\theta}}_s - \hat {\boldsymbol{\theta}}_r = A \left(\begin{array}{c} \phantom{-}\textbf{1}_{c_a} \\[6pt] -\textbf{1}_{n_y} \\[6pt] \phantom{-}\textbf{0}_{n_c} \\[6pt] \phantom{-} \textbf{0}_{c_a} \end{array}\right) + B \left(\begin{array}{c} \phantom{-}\textbf{1}_{c_a} \\[6pt] \phantom{-}\textbf{0}_{n_y} \\[6pt] -\textbf{1}_{n_c} \\[6pt] \phantom{-} \textbf{0}_{c_a} \end{array}\right) + C \left(\begin{array}{c} \Delta_a \textbf{\textit{x}}_{c_a} \\[6pt] -\textbf{\textit{x}}_y \\[6pt] \textbf{\textit{x}}_c - \omega_c \textbf{1}_{n_c} \\[6pt] \textbf{0}_{c_a} \end{array}\right) + D \left(\begin{array}{c} \textbf{0}_{c_a} \\[6pt] \textbf{\textit{x}}_y - \bar y \textbf{1}_{n_y} \\[6pt] \textbf{0}_{n_c} \\[6pt] \textbf{1}_{c_a} \end{array}\right)\end{equation}

for some A, B, C and D; here, as before, $\omega_c = n_a + 2\Delta_a - 1$ . In our example, we found $A = 5.742$ , $B = -0.701$ , $C= -0.0343$ and $D = -0.523$ . Pre-multiplying (52) through by $\mbox{blockdiag}\{\textbf{\textit{B}}_a, \textbf{\textit{I}}_{n_y}, \textbf{\textit{I}}_{n_c}, \textbf{\textit{B}}_a\}$ and using $\textbf{\textit{B}}_a \textbf{1}_{c_a} = \textbf{1}_{n_a}$ and $\textbf{\textit{B}}_a\Delta_a \textbf{\textit{x}}_{c_a} = \textbf{\textit{x}}_{n_a} + \omega_a \textbf{1}_{n_a}$ , we find

(53) \begin{align}\hat {\boldsymbol{\theta}}_s^* - \hat {\boldsymbol{\theta}}_r^* =A \left(\begin{array}{c} \phantom{-}\textbf{1}_{n_a} \\[6pt] -\textbf{1}_{n_y} \\[6pt] \phantom{-}\textbf{0}_{n_c} \\[6pt] \phantom{-} \textbf{0}_{n_a} \end{array}\right) + B \left(\begin{array}{c} \phantom{-}\textbf{1}_{n_a} \\[6pt] \phantom{-}\textbf{0}_{n_y} \\[6pt] -\textbf{1}_{n_c} \\[6pt] \phantom{-}\textbf{0}_{n_a} \end{array}\right) + C \left(\begin{array}{c} \textbf{\textit{x}}_a + \omega_a \textbf{1}_{n_a} \\[6pt] -\textbf{\textit{x}}_y \\[6pt] \textbf{\textit{x}}_c - \omega_c \textbf{1}_{n_c} \\[6pt] \textbf{0}_{n_a} \end{array}\right) + D \left(\begin{array}{c} \textbf{0}_{n_a} \\[6pt] \textbf{\textit{x}}_y - \bar y \textbf{1}_{n_y} \\[6pt] \textbf{0}_{n_c} \\[6pt] \textbf{1}_{n_a} \end{array}\right)\end{align}

where as before $\omega_a = 2 \Delta_a - 1$ . It follows immediately that

(54) \begin{eqnarray} \hat {\boldsymbol{\alpha}}_s - \hat {\boldsymbol{\alpha}}_r &=& C \textbf{\textit{x}}_a + (A + B + w_a C) \textbf{1}_{n_a} \end{eqnarray}

(55) \begin{eqnarray}\hat {\boldsymbol{\kappa}}_s - \hat {\boldsymbol{\kappa}}_r &=& (D- C) \textbf{\textit{x}}_y - (A + D \bar y) \textbf{1}_{n_y} \end{eqnarray}

(56) \begin{eqnarray} \hat {\boldsymbol{\gamma}}_s - \hat {\boldsymbol{\gamma}}_r &=& C \textbf{\textit{x}}_c - (B + w_c C) \textbf{1}_{n_c} \end{eqnarray}

(57) \begin{eqnarray} \hspace*{-50pt}\hat {\boldsymbol{\beta}}_s - \hat {\boldsymbol{\beta}}_r &=& D \textbf{1}_{n_a}\hspace*{10pt}\end{eqnarray}

We now turn to forecasting with the smooth APCI model. We use the same notation as in the APC model. Let $n_f$ be the length of the forecast and $\textbf{\textit{x}}_f = (1, \ldots, n_f)^\prime$ . Let $\hat {\boldsymbol{\kappa}}_{s, f}$ and $\hat {\boldsymbol{\kappa}}_{r, f}$ be the forecast values of $\hat {\boldsymbol{\kappa}}_s$ and $\hat {\boldsymbol{\kappa}}_r$ , and $\hat {\boldsymbol{\gamma}}_{s, f}$ and $\hat {\boldsymbol{\gamma}}_{r, f}$ be the forecast values of $\hat {\boldsymbol{\gamma}}_s$ and $\hat {\boldsymbol{\gamma}}_r$ , respectively. We know from (55) and (56) that $\hat {\boldsymbol{\kappa}}_s - \hat {\boldsymbol{\kappa}}_r$ and $\hat {\boldsymbol{\gamma}}_s - \hat {\boldsymbol{\gamma}}_r$ are both linear so, exactly as in the APC model, we can show that

(58) \begin{eqnarray} \hat {\boldsymbol{\kappa}}_{s, f} - \hat {\boldsymbol{\kappa}}_{r, f} &=& (n_y (D - C) - (A + D \bar y)) \textbf{1}_{n_f} + (D - C) \textbf{\textit{x}}_f \end{eqnarray}

(59) \begin{eqnarray} \hat {\boldsymbol{\gamma}}_{s, f} - \hat {\boldsymbol{\gamma}}_{r, f} &=& ((n_c - w_c)C - B) \textbf{1}_{n_f} + C \textbf{\textit{x}}_f\end{eqnarray}

We now establish the invariance of the forecast values of $\log\boldsymbol{\mu}$ . Define

\begin{equation*} \hat {\boldsymbol{\theta}}_{s, f}^* = (\hat {\boldsymbol{\alpha}}_s^\prime, \hat {\boldsymbol{\kappa}}_s^\prime, \hat {\boldsymbol{\kappa}}_{s, f}^\prime, \hat {\boldsymbol{\gamma}}_s^\prime, \hat {\boldsymbol{\gamma}}_{s, f}^\prime, \hat {\boldsymbol{\beta}}_s^\prime)^\prime\end{equation*}

with a similar definition for $\hat {\boldsymbol{\theta}}_{r, f}^*$ . Let $\textbf{\textit{X}}_f$ be the model matrix for the smooth APCI model for ages $\textbf{\textit{x}}_a$ and years $(\textbf{\textit{x}}_y^\prime, \textbf{\textit{x}}_{y, f}^\prime)^\prime$ , where $\textbf{\textit{x}}_{y, f} = (n_y + 1, \ldots, n_y + n_f)^\prime$ ; let $\textbf{\textit{X}}_f^*$ be the corresponding model matrix for the unsmoothed model. We wish to show that $\textbf{\textit{X}}_f (\hat {\boldsymbol{\theta}}_{s, f} - \hat {\boldsymbol{\theta}}_{r, f}) =\textbf{0}$ . We observe that $\textbf{\textit{X}}_f \hat {\boldsymbol{\theta}}_{s, f} = \textbf{\textit{X}}_f^*\hat {\boldsymbol{\theta}}^*_{s, f}$ , since $(\textbf{1}_{n_y} \otimes \textbf{\textit{B}}_a) \hat{\textbf{\textit{a}}} = (\textbf{1}_{n_y} \otimes \textbf{\textit{I}}_{n_a}) \hat {\boldsymbol{\alpha}}$ and $((\bar y \textbf{1}_{n_y} - \textbf{\textit{x}}_y) \otimes \textbf{\textit{B}}_a) \hat {\textbf{\textit{b}}} =((\bar y \textbf{1}_{n_y} - \textbf{\textit{x}}_y) \otimes \textbf{\textit{I}}_{n_a}) \hat{\boldsymbol{\beta}}$ .

Figure 7 Three regions for fitted and forecast values of $\log \boldsymbol{\mu}$ for ages $1,\ldots,55$ , years $1,\ldots,45$ , and a 10-year-ahead forecast.

There are three cases to consider, as shown in Figure 7. In region 1, the fitted values depend on the estimated values of $\boldsymbol{\alpha}$ , $\boldsymbol{\kappa}$ , $\boldsymbol{\gamma}$ and $\boldsymbol{\beta}$ ; in region 2, the forecast values depend on the estimated values of $\boldsymbol{\alpha}$ , $\boldsymbol{\gamma}$ and $\boldsymbol{\beta}$ , and the forecast values of $\boldsymbol{\kappa}$ ; and in region 3, the forecast values depend on the estimated values of $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$ , and the forecast values of $\boldsymbol{\kappa}$ and $\boldsymbol{\gamma}$ . We see from (54), (55), (56), (57), (58) and (59) that

\begin{equation*} \textbf{\textit{X}}^*_f (\hat {\boldsymbol{\theta}}^*_{s, f} - \hat {\boldsymbol{\theta}}^*_{r, f}) = A \textbf{\textit{v}}_1 + B \textbf{\textit{v}}_2 + C \textbf{\textit{v}}_3 + D \textbf{\textit{v}}_4\end{equation*}

for some vectors $\textbf{\textit{v}}_1$ , $\textbf{\textit{v}}_2$ , $\textbf{\textit{v}}_3$ and $\textbf{\textit{v}}_4$ . We show that $\textbf{\textit{v}}_1 = \textbf{\textit{v}}_2 = \textbf{\textit{v}}_3 = \textbf{\textit{v}}_4 = \textbf{0}$ ; this will establish invariance of fitted and forecast values.

It is convenient to use a double-suffix notation to identify the entries in $\textbf{\textit{v}}_1$ , $\textbf{\textit{v}}_2$ , $\textbf{\textit{v}}_3$ and $\textbf{\textit{v}}_4$ . Thus, cell (i, j) is associated with $\textbf{\textit{v}}_1(i, j)$ and corresponds to the $n_a(j-1) + i$ entry in $\textbf{\textit{v}}_1$ ; we have a similar correspondence for $\textbf{\textit{v}}_2$ , $\textbf{\textit{v}}_3$ and $\textbf{\textit{v}}_4$ . We check invariance region by region.

Consider region 1. Suppose cell (i, j) lies in region 1 with cohort index $n_a - i + j$ . We pick out the terms in A, B, C and D from (54), (55), (56) and (57) as appropriate. The bracket notation $(\;)$ indicates which terms in $\textbf{\textit{X}}^*_f (\hat {\boldsymbol{\theta}}^*_{s, f} - \hat {\boldsymbol{\theta}}^*_{r, f})$ are used, i.e., $\boldsymbol{\alpha}$ , $\boldsymbol{\kappa}$ , $\boldsymbol{\gamma}$ and $\boldsymbol{\beta}$ in turn,

\begin{eqnarray*} A &:& (1) + (-1) + (0) + (0) = 0 \\[4pt] B &:& (1) + (0) + (-1) + (0) = 0 \\[4pt] C &:& (i + \omega_a) + (-j) + (n_a - i + j - \omega_c) + (0) = 0 \\[4pt] D &:& (0) + (j - \bar y) + (0) + (\bar y - j) = 0\end{eqnarray*}

and we have verified invariance in region 1. Of course, we already knew this from the general result on fitted values.

In region 2, we consider the j-ahead forecast for $\boldsymbol{\kappa}$ . This is cell $(i, n_y + j)$ with cohort index $n_a - i + n_y + j$ . Using (54), (56), (57) and (58), we find

\begin{eqnarray*} A &:& (1) + (-1) + (0) + (0) = 0 \\[4pt] B &:& (1) + (0) + (-1) + (0) = 0 \\[4pt] C &:& (i + \omega_a) + (- n_ y - j) + (n_a - i + n_y + j - \omega_c) + (0) = 0 \\[4pt] D &:& (0) + (n_y - \bar y + j) + (0) + (\bar y - n_y - j) = 0\end{eqnarray*}

as required. Finally, in region 3, we consider the j-ahead forecast for $\boldsymbol{\gamma}$ . For age i, we have the $(i+j-1)$ -ahead forecast for $\boldsymbol{\kappa}$ , i.e., cell $(i, n_y + i + j -1)$ . We use (54), (57), (58) and (59) and find

\begin{eqnarray*} A &:& (1) + (-1) + (0) + (0) = 0 \\[4pt] B &:& (1) + (0) + (-1) + (0) = 0 \\[4pt] C &:& (i + \omega_a) + (- n_ y - i - j + 1) + (n_c - \omega_c + j) + (0) = 0 \\[4pt] D &:& (0) + (n_y - \bar y + i + j - 1) + (0) + (\bar y - n_y - i - j + 1) = 0\end{eqnarray*}

We conclude that the fitted and forecast values are invariant with respect to the choice of constraints when an ARIMA $(p, \delta, q)$ model with fitted mean, $\delta =1$ or 2, is used to forecast. We note that the proof of invariance in the smooth APC model is a special case of the above. In (52), we omit the final row and column, and we have reduced (52) to (42).

4.4.1 LC model

The LC model is

(60) \begin{equation} \log \mu_{x, y} = \alpha_x + \beta_x \kappa_y,\, x = 1,\ldots, n_a,\, y = 1,\ldots, n_y\end{equation}

We require one location and one scale constraint to make the parameters estimable. We use

(61) \begin{equation} \sum \kappa_y = 0,\; \sum \beta_x = 1\end{equation}

as in Lee & Carter (Reference Lee and Carter1992). The following transformation leaves the fitted values of $\log \boldsymbol{\mu}$ unchanged for any values of A and B.

(62) \begin{equation} \alpha_x + \beta_x \kappa_y \mapsto (\alpha_x + A\beta_x) + (\beta_x/B)(B \kappa_y - AB) = \alpha_x^* + \beta_x^* \kappa_y^*\end{equation}

Let $\hat {\boldsymbol{\alpha}}_s$ , $\hat {\boldsymbol{\beta}}_s$ and $\hat {\boldsymbol{\kappa}}_s$ be the maximum likelihood estimates of $\boldsymbol{\alpha}_s$ , $\boldsymbol{\beta}_s$ and $\boldsymbol{\kappa}_s$ under the constraints (61) with the usual corresponding notation for the estimates $\hat {\boldsymbol{\alpha}}_r$ , $\hat {\boldsymbol{\beta}}_r$ and $\hat {\boldsymbol{\kappa}}_r$ under some random constraint system. It follows from (62) that

(63) \begin{eqnarray} \hspace*{-16pt}\hat {\boldsymbol{\alpha}}_s &=& \hat {\boldsymbol{\alpha}}_r + A \hat {\boldsymbol{\beta}}_r \end{eqnarray}

(64) \begin{eqnarray} \hspace*{-33pt}\hat {\boldsymbol{\beta}}_s &=& \hat {\boldsymbol{\beta}}_r/B \end{eqnarray}

(65) \begin{eqnarray} \hat {\boldsymbol{\kappa}}_s &=& B \hat {\boldsymbol{\kappa}}_r - AB \textbf{1}_{n_y}\end{eqnarray}

for some scalars A and B.

Brouhns et al. (Reference Brouhns, Denuit and Vermunt2002) used maximum likelihood to estimate the parameters. We also use maximum likelihood but use (15); this is a full Newton–Raphson scheme which allows estimation with both smoothing and constraints. Following Currie (Reference Currie2013), we consider two coupled GLMs:

(66) \begin{align} \mbox{GLM1:}\, \log \boldsymbol{\mu} = \textbf{1}_{n_y} \otimes \tilde {\boldsymbol{\alpha}} + \textbf{\textit{X}}_1 \boldsymbol{\beta},\; \textbf{\textit{X}}_1 = \tilde {\boldsymbol{\kappa}} \otimes \textbf{\textit{I}}_{n_a} \end{align}

(67) \begin{align} \hspace{6.5pt}\mbox{GLM2:}\,\log \boldsymbol{\mu} = \textbf{\textit{X}}_2 \boldsymbol{\theta} , \; \textbf{\textit{X}}_2 = [\textbf{1}_{n_y} \otimes \textbf{\textit{I}}_{n_a} \, :\ \textbf{\textit{I}}_{n_y} \otimes \tilde {\boldsymbol{\beta}}]\end{align}

where $\boldsymbol{\theta} = (\boldsymbol{\alpha}^\prime, \boldsymbol{\kappa}^\prime)^\prime$ ; here $\tilde{\boldsymbol{\alpha}}$ , $\tilde {\boldsymbol{\beta}}$ and $\tilde{\boldsymbol{\kappa}}$ represent current estimates of $\boldsymbol{\alpha}$ , $\boldsymbol{\beta}$ and $\boldsymbol{\kappa}$ , respectively. Both GLM1 and GLM2 are in the class defined in equation (1) but note that in GLM1 the offset is $\log \textbf{\textit{e}} + \textbf{1}_{n_y} \otimes \tilde {\boldsymbol{\alpha}}$ . In GLM1, we estimate $\boldsymbol{\beta}$ for current values of $\boldsymbol{\alpha}$ and $\boldsymbol{\kappa}$ , while in GLM2 we estimate $\boldsymbol{\alpha}$ and $\boldsymbol{\kappa}$ for given values of $\boldsymbol{\beta}$ . We now iterate between GLM1 and GLM2 until convergence.

We fit the model with (a) the standard constraints (61) and (b) the random constraints in GLM1 on $\boldsymbol{\beta}$ , i.e., $\sum_1^{n_a} u_i \beta_i = 1$ , and the random constraints in GLM2 on $\boldsymbol{\theta} = (\boldsymbol{\alpha}^\prime, \boldsymbol{\kappa})^\prime$ , i.e., $\sum_1^{n_a+n_y} u_i\theta_i = 0$ ; the $u_i$ are independent $\mathcal{U}(0,\,1)$ .

The upper panels and the lower left panel of Figure 8 show the estimates of $\boldsymbol{\alpha}$ , $\boldsymbol{\kappa}$ and $\boldsymbol{\beta}$ , respectively, under the standard and random constraints. In our example, we found $A = 5.907$ , $B = 1.983$ and $AB = 11.71$ . The parameter estimates in Figure 8 satisfy (63), (64) and (65) with these values of A and B. The lower right panel shows the single (invariant) estimate of log mortality at age 65.

Figure 8 Parameter estimates in the LC model (60) under standard and random constraints: $\hat {\boldsymbol{\alpha}}_s = \hat {\boldsymbol{\alpha}}_r + A \hat {\boldsymbol{\beta}}_r$ , $\hat {\boldsymbol{\beta}}_s = \hat {\boldsymbol{\beta}}_r/B$ , $\hat {\boldsymbol{\kappa}}_s = B \hat {\boldsymbol{\kappa}}_r - AB \textbf{1}_{n_y}$ ; $A = 5.907$ , $B= 1.983$ and $AB = 11.71$ . Fitted invariant $\log \mu$ for age 65.

Forecasting in the LC model is particularly straightforward. As usual, we forecast with an ARIMA(p, $\delta$ , q) model, $\delta=1$ or 2. Let $\hat {\boldsymbol{\kappa}}_{s, f}$ and $\hat {\boldsymbol{\kappa}}_{r, f}$ denote some forecast values of $\hat {\boldsymbol{\kappa}}_s$ and $\hat {\boldsymbol{\kappa}}_r$ , respectively. We have from (65) $\hat{\boldsymbol{\kappa}}_s = B\hat {\boldsymbol{\kappa}}_r - AB \textbf{1}_{n_y}$ ; it follows immediately from Appendix C that $\hat {\boldsymbol{\kappa}}_{s, f} = B\hat {\boldsymbol{\kappa}}_{r, f} - AB \textbf{1}_{n_f}$ , where $n_f$ is the number of years in the forecast. Thus, the forecasts of $\boldsymbol{\kappa}$ satisfy (65), and hence forecasts of mortality are invariant with respect to the choice of constraint system.

4.4.2 Smooth LC model

One unsatisfactory feature of the Lee–Carter model, particularly for actuaries, is that irregularities in the estimate of $\boldsymbol{\beta}$ lead to irregular forecasts of mortality at individual ages and can even lead to forecasts at adjacent ages crossing over; see Currie (Reference Currie2013) for examples of this behaviour. Delwarde et al. (Reference Delwarde, Denuit and Eilers2007) addressed this problem by smoothing the estimate of $\boldsymbol{\beta}$ with the method of P-splines. They set $\boldsymbol{\beta} = \textbf{\textit{B}}_a \textbf{\textit{b}}$ , where $\textbf{\textit{B}}_a$ is a regression matrix evaluated on a basis of B-splines along age. The model matrix for GLM1 in (66) becomes

\begin{equation*} \mbox{GLM1:} \quad \log \boldsymbol{\mu} = \textbf{1}_{n_y} \otimes \tilde {\boldsymbol{\alpha}} + \textbf{\textit{X}}_1 \textbf{\textit{b}},\; \textbf{\textit{X}}_1 = \tilde {\boldsymbol{\kappa}} \otimes \textbf{\textit{B}}_a\end{equation*}

However, $[\tilde {\boldsymbol{\kappa}} \otimes \textbf{\textit{B}}_a] \textbf{\textit{b}} = [\tilde {\boldsymbol{\kappa}} \otimes \textbf{\textit{I}}_{n_a}] \boldsymbol{\beta}$ , i.e., we are back in the model definition (66). The model matrix for GLM2 is not affected by the smoothing of $\boldsymbol{\beta}$ . We deduce that (63), (64) and (65) hold in the smooth case too. We conclude that invariance of fitted and forecast values holds for the smooth model. Numerical work supports this conclusion.

4.5.1 CBD model with cohort effects

Cairns et al. (Reference Cairns, Blake and Dowd2006) introduced the model $\log \mu_{x, y} =\kappa_y^{(1)} + (x-\bar x)\kappa_y^{(2)}$ which is often referred to as the CBD model. Cairns et al. (Reference Cairns, Blake, Dowd, Coughlan, Epstein, Ong and Balevich2009) modified this model with the addition of cohort effects:

(68) \begin{equation} \log \mu_{x, y} = \kappa_y^{(1)} + \kappa_y^{(2)}(x-\bar x) + \gamma_{c(x, y)},\; x = 1,\ldots, n_a,\, y =1,\ldots,n_y\end{equation}

where $\bar x$ is the mean age. The model matrix is

(69) \begin{equation} \textbf{\textit{X}} = [\textbf{\textit{I}}_{n_y} \otimes \textbf{1}_{n_a}\ :\ \textbf{\textit{I}}_{n_y} \otimes (\textbf{\textit{x}}_a - \bar x \textbf{1}_{n_a})\ :\ \textbf{\textit{X}}_c].\end{equation}

The coefficient vector is $\boldsymbol{\theta} = ({\boldsymbol{\kappa}^{(1)}}^\prime,{\boldsymbol{\kappa}^{(2)}}^\prime, \boldsymbol{\gamma}^\prime)^\prime$ with length $n_a+ 3n_y - 1$ . The rank of $\textbf{\textit{X}}$ is $n_a + 3n_y - 3$ , so two constraints are required to give unique parameter estimates. We use $\sum \gamma_c = \sum c \gamma_c = 0$ as our standard constraints (Cairns et al., Reference Cairns, Blake, Dowd, Coughlan, Epstein, Ong and Balevich2009) and two sets of random constraints on $\boldsymbol{\theta}$ . We compute $\mathcal{N}(\textbf{\textit{X}})$ , the null space of $\textbf{\textit{X}}$ , and conclude that

(70) \begin{eqnarray} \boldsymbol{\kappa}^{(1)}_r &=& \boldsymbol{\kappa}^{(1)}_s + [A + (\bar x - n_x)B] \textbf{1}_{n_y} - B \textbf{\textit{x}}_y \nonumber \\ \boldsymbol{\kappa}^{(2)}_r &=& \boldsymbol{\kappa}^{(2)}_s + B \textbf{1}_{n_y} \end{eqnarray}

\begin{align} \hspace*{-65pt}\boldsymbol{\gamma}_r = \boldsymbol{\gamma}_s - A \textbf{1}_c + B \textbf{\textit{x}}_c \nonumber\end{align}

where we have used our usual notation for the estimates under the standard and random constraints. Once again we see that the parameters are linearly related, so provided we choose a forecasting method which preserves these relationships we can conclude that the forecasts of mortality are invariant with respect to the choice of constraints.

4.5.2 Reduced Plat model

Plat (Reference Plat2009) proposed the model

(71) \begin{equation} \log \mu_{x, y} = \alpha_x + \kappa_y^{(1)} + (x-\bar x)\kappa_y^{(2)} + \gamma_{c(x, y)},\; x = 1,\ldots, n_a,\, y =1,\ldots,n_y\end{equation}

Hunt & Blake (Reference Hunt and Blake2020b) refer to this model as the reduced Plat model. We see that it is formally equivalent to the unsmoothed APCI model (44) so the dimension of its null space, $\mathcal{N}(\textbf{\textit{X}})$ , is five, as is readily verified directly. In particular, $\mathcal{N}(\textbf{\textit{X}})$ has a quadratic term, as in (46). This means that an ARIMA(p,1,q) model will not lead to invariant forecasts. We can smooth the age term with a second-order penalty; this will reduce the dimension of the null space from five to four, just as in APCI model, removing the quadratic term from $\mathcal{N}(\textbf{\textit{X}})$ . Once more, we are back in the linear case. We will comment further on this increase in the effective rank of the model in our concluding remarks.