Identification and Forecasting in Mortality Models

Abstract

Mortality models often have inbuilt identification issues challenging the statistician. The statistician can choose to work with well-defined freely varying parameters, derived as maximal invariants in this paper, or with ad hoc identified parameters which at first glance seem more intuitive, but which can introduce a number of unnecessary challenges. In this paper we describe the methodological advantages from using the maximal invariant parameterisation and we go through the extra methodological challenges a statistician has to deal with when insisting on working with ad hoc identifications. These challenges are broadly similar in frequentist and in Bayesian setups. We also go through a number of examples from the literature where ad hoc identifications have been preferred in the statistical analyses.

1. Introduction

Mortality models are commonly used in a wide range of fields such as actuarial sciences, epidemiology, and sociology. They are often used in important decisions such as how to deal with unisex legislation in the pension industry; see Ornelas et al. [1] and Jarner and Kryger [2]. However, such models do often have inbuilt identification issues stemming from overparametrisation. While identification issues are omnipresent in statistical modelling, this paper focuses on mortality modelling, where estimated parameters are treated as time series and extrapolated to give forecasts of future mortality. The underlying theme of this paper is to provide strategies of avoiding arbitrariness resulting from the identification process. We suggest two ways forward. First, we can reparametrise the model in terms of a freely varying parameter, which therefore has to be of lower dimension than the original parameter. Secondly, we can work with an identified version of the original parameter as long as we keep track of the consequences of the identification choice. That way we ensure that two researchers making different identification choices get the same statistical inferences and forecasts.

A simple example is the age-period model for an age-period array of mortality rates. It is well-known that the levels of the age- and period-effects cannot be determined from the likelihood representing the overparametrisation of the model. When the estimated age- and period-effects are treated as time series and subjected to plotting and extrapolation, then our approach ensures that the statistical analysis is the same for two researchers identifying the above model in two different ways. Whereas this issue is relatively simple for the age-period model, identification becomes more tricky for complicated models such as the age-period-cohort model and the model of Lee and Carter [3], let alone two-sample situations.

Mortality models are built as a combination of age, period, and cohort-effects, but the likelihood only varies with a surjective function of these time effects. The time effects can be divided into two parts. One part that moves the likelihood function and another part which does not induce variation in the likelihood function. We will argue that all inferences and forecasts should be concerned primarily with the part of the parameter that moves the likelihood function. This does not preclude the researcher from working with the time effects, but it gives some limitations on what can be done. This is important because the motivation and the intuition of mortality models typically originate in the time effects. For instance, in the context of an age-period-cohort model linear trends cannot be identified so time series plots of the time effects need to be invariant to linear trends and extrapolations of time effects must preserve the arbitrary linear trend in the time effects. This applies regardless of whether the identification issue is dealt with in a frequentist manner or by Bayesian methods.

To formalise the discussion slightly return to the age-period example. Denote the predictor for the age-period data array by 𝜇. The age-period model then determines how the predictor 𝜇 varies with a vector 𝜃 summarising age and period effects. That vector is split into two components 𝜉 and 𝜆 so that the predictor only depends on 𝜃 through 𝜉 but not on 𝜆 which cannot be identified by statistical analysis. In the age-period example 𝜉 could reflect the contrasts and the overall level of the predictor 𝜇, whereas 𝜆 reflects the level of the age effect. The more principled solution is then to work exclusively with 𝜉 and simply consider 𝜃 as a motivation rather than the objective of the analysis. Another solution is to ad hoc identify 𝜆 based on a notion of mathematical convenience or based on a particular purpose given the substantive context.

Once an ad hoc identification of 𝜆 is chosen the identification problem appears to go away, because the likelihood analysis can now go through. The reason is that the variation of 𝜃 is now reduced to the variation of 𝜉 precisely because 𝜆 is fixed. Suppose two researchers choose the same likelihood and the same parametrisation of 𝜉 but different ad hoc identifications 𝜆^† and 𝜆^‡. Which of their conclusions will be the same and which will be different? As the likelihood only depends on 𝜉 the fits of the two researchers will be identical. But differences might arise if the statistical inference or forecasting or any other statistical analysis involves 𝜆 in some way.

Indeed, with many extrapolation methods forecasts will be invariant to the choice of 𝜆. But, there will also be extrapolation methods where this is not the case. Examples arise in the age-period-cohort model, where linear trends have to be handled with care.

We will start by analysing linearly parametrised models at a rather general level. We do this with two aspects in mind. First, we need to step back to a point in the analysis before ad hoc identification is made. Secondly, we also want to avoid the discussion of how to choose 𝜉 and 𝜆, which tend to be specific to the mortality model in question. Working at the general level we can focus on the mappings between different parametrisations and the invariance properties coming from these mappings. It is then seen that the parameter 𝜉 arises as a maximal invariant. The general setting also allows the formulation of a series of results discussing different types of ad hoc identification, first in a frequentist fashion and then in a Bayesian fashion.

Subsequently, we will consider the age-period-cohort model in detail, both for one- and two-sample situations. Using the general results it becomes easier to see that a number of popular methods inadvertently include features that are not invariant to ad hoc identification. These include the “intrinsic estimator” advocated by Yang et al. [4], the “mixed model approach” by Yang and Land [5], the Bayesian approach by Berzuini and Clayton [6], and the two-sample analysis by Riebler and Held [7]. Finally, we consider the nonlinearly parametrised model of Lee and Carter [3]. The nonlinearity gives a further complication since the mapping from the time effects to the mortality predictor is nondifferentiable. As it turns out the mortality predictor varies in a smooth space, so the nondifferentiability is avoided by working directly with the mortality predictor instead of the original time effects. Instead, a Lee-Carter application should consider whether a certain matrix has rank of unity or zero. Apart from that the analysis is similar to that of linearly parametrised models. Likewise a theory is given for two-sample situations.

Throughout the paper our concern rests exclusively with the identification problem and the consequences of ad hoc identification for estimation, plots, inference, and forecasting. In practice, important additional concerns are how to choose appropriate models and forecasting methods. We would like to refer to Girosi and King [8], Pitacco et al. [9] for general discussions of these issues, and also to Kuang et al. [10] and Coelho and Nunes [11] for discussions of forecast methods in the light of structural breaks. Instead, the aim of the paper is to present an overall framework that can help streamlining the identification discussion that has appeared in so many papers in so many fields over so many years.

Section 2 of this paper considers standard linear statistical models, which lend themselves to a relative straightforward analysis based on linear algebra. Any ad hoc identification splits the time effect into two components. The first component is an arbitrary component, which is not needed for the identification of the likelihood. The other component is necessary and sufficient to identify the model and hence sufficient for statistical analysis. In Section 3 it is outlined how to analyze the statistical model when the latter component is ad hoc identified. It is argued that this can cause difficulties for estimation, interpretation, and forecast. In Section 4 it is shown that Bayesian analysis shares the same challenges as the frequentist approach. In Sections 5 and 6 we study the two particular examples: the omnipresent age-period-cohort and Lee-Carter mortality models. All proofs are collected in the Appendix.

2. Statistical Models with Linear Parametrisations

In this section we present the identification problem in a linear framework. The problem is solved by analysing the mapping from the original time effect to the predictor which, in turn, leads to standard statistical analysis. In Section 6 we show how these ideas transfer to a nonlinear context. This contrasts with Section 3 in which we illustrate the analytical challenges and inconveniences arising from ad hoc identification.

In Section 2.1 we present the overparametrized linear model for the mortality predictor. The identification problem is defined in Section 2.2 via the likelihood. In an overparametrized linear model two different parameters might produce the same likelihood. In Section 2.3 we analyze the mapping from the overparametrised parameter to the predictor. This mapping enables us to split the overparametrised parameter into two. One arbitrary parameter and one parameter identify the model without being overparametrised. This latter parameter is shown to be a maximal invariant parameter. In Section 2.4 it is demonstrated how any statistical analysis can be based on this maximal invariant parameter alone. In particular we comment that visual data representations, hypothesis testing, and forecasting are simple and well defined. This in turn leads to standard statistical analysis.

The analysis of the linearly parametrised involves projections on linear or affine spaces and on their orthogonal complements. It is therefore convenient to introduce the following notation. A matrix 𝑚 has full column rank if 𝑚^′𝑚 is invertible. In this case the orthogonal complement 𝑚_⟂ is a matrix so 𝑚^′_⟂𝑚=0 and (𝑚,𝑚_⟂) is invertible. Thus, when 𝑚 itself is invertible then 𝑚_⟂ is the empty matrix. It is not difficult to calculate 𝑚_⟂ in practise, an explicit construction of 𝑚_⟂ follows from a singular value decomposition of 𝑚𝑚^′, choosing 𝑚_⟂ as the eigenvectors associated with the zero eigenvalues. Moreover, let 𝑚=𝑚(𝑚^′𝑚)⁻¹ so that 𝑚^′𝑚 is the identity matrix, while 𝑚_⟂=𝑚_⟂(𝑚^′_⟂𝑚_⟂)⁻¹.

2.1. The Model

Think of the time effect 𝜃 as our preferred intuitive, but unidentified parameter, and think of the predictor 𝜇 as some function of 𝜃 specifying the model at hand. In a Poisson type model, where the mean specifies the distribution, 𝜇 could be the log of that mean. Such Poisson models are omnipresent in mortality models. We will often think of 𝜃 as containing some time effects. Often forecasting is carried out simply by isolating and extrapolating such a time effect.

Consider a data vector 𝑌 of dimension 𝑛. This could, for instance, be the vector consisting of the stacked mortality rates for a rectangular age-period array of dimension 𝐼×𝐽 in which case 𝑛=𝐼𝐽. The statistical model for 𝑌 could be a generalized linear model. This involves an appropriately chosen distribution and a link function, which links the expected mortality rate to an 𝑛-dimensional predictor, which is denoted by 𝜇. Taken together this defines a likelihood function 𝖫(𝜇;𝑌).

The model for the predictor 𝜇 is constructed in terms of, for instance, age, period, and cohort time effects. These time effects are summarized in a vector 𝜃, which is of dimension 𝑞<𝑛. Therefore 𝜇 is a surjective function of 𝜃. For the moment the specification of the predictor is assumed linear so that 𝜇=𝐷𝜃for𝜃∈Θ=ℝ^𝑞, (1) for some design matrix 𝐷∈ℝ^𝑛×𝑞. We refer to this specification as the mortality model, while the space Θ is the time effect space. The time effect space is chosen as an unrestricted real space in accordance with the starting point of most mortality analyses.

The parameter space for the likelihood function and therefore for the statistical model is given by the range of variation for the predictor 𝜇; that is, 𝑀=𝜇∈ℝ^𝑛∶𝜇=𝐷𝜃for𝜃∈Θ=ℝ^𝑞. (2) The likelihood function is assumed uniquely identified on this space in the sense that for all pairs of predictors so 𝜇^†≠𝜇^‡ then the likelihood of 𝜇^†, 𝜇^‡ differ; that is, 𝖫𝜇^†;𝑌≠𝖫𝜇^‡;𝑌, (3) for 𝑌 in a set with positive probability.

2.2. The Identification Problem

The identification problem of mortality models arises when the mapping from the time effect space Θ to the parameter space 𝑀 is surjective but not injective. With a linear parametrisation this arises when the design matrix 𝐷 has reduced column rank 𝑝<𝑞 so 𝐷^′𝐷 is singular. In this situation there exists time effects 𝜃^†≠𝜃^‡ with the same likelihood: 𝖫𝐷𝜃^†;𝑌=𝖫𝐷𝜃^‡;𝑌, (4) for all data 𝑌. Then the time effect space Θ is not useful as parameter space for the statistical model.

2.3. Analysing the Mapping 𝜃↦𝜇

When analysing the mapping from our intuitively preferred parametrisation 𝜃 into the linear predictor 𝜇, we will be able to rewrite 𝜃 as a sum of two components: one is a function of the predictor and the other is the arbitrary part varying with 𝜃, but not with the predictor. We provide two methods for analysis.

The first method is to find a basis 𝑋∈ℝ^𝑛×𝑝 with full column rank𝑝 for the design 𝐷. The design matrix of the mortality model can then be expressed as 𝐷=𝑋𝐴^′ for some matrix 𝐴∈ℝ^𝑞×𝑝 with full column rank𝑝. Introduce a new 𝑝-dimensional parameter: 𝜉=𝐴^′𝜃. (5) The parameter space 𝑀 can then be written more parsimoniously as 𝑀=𝜇∈ℝ^𝑛∶𝜇=𝑋𝜉for𝜉∈ℝ^𝑝. (6) The mapping from 𝜉 to 𝜇 is bijective, so the statistical model can just as well be parametrised in terms of 𝜉∈Ξ=ℝ^𝑝.

Alternatively, the identification problem can be expressed through an invariance argument. This argument relates to the parameterization but resembles the classical invariance argument for reduction of data; see Cox and Hinkley [12, page 157]. With a linear parametrisation the argument involves the orthogonal complement to the matrix 𝐴. That is a matrix 𝐴_⟂∈ℝ^{𝑞×(𝑞−𝑝)} which has the properties that 𝐴^′_⟂𝐴=0 and that (𝐴,𝐴_⟂) is invertible. The mortality model (1) is defined by the mapping 𝜃⟼𝜇=𝐷𝜃=𝑋𝐴^′𝜃, (7) from Θ=ℝ^𝑞 to 𝑀. This mapping is surjective in that two different values of 𝜃 may result in the same 𝜇 and therefore the same likelihood. These equivalence classes in the time effect space can be described by the group of transformations 𝑔∶𝜃⟼𝜃+𝐴_⟂𝜁, (8) acting on Θ for arbitrary 𝜁∈ℝ^𝑞−𝑝. Indeed, it holds that 𝜃 and 𝑔(𝜃) will result in the same 𝜇. The mapping (7) is therefore invariant to the group 𝑔. We will argue that the parameter 𝜉=𝐴^′𝜃 is a maximal invariant to the group 𝑔 acting on Θ, which provides a link with (6). It has to be argued that for any 𝜃^†, 𝜃^‡ so that 𝜉^†=𝐴^′𝜃^† equals 𝜉^‡=𝐴^′𝜃^‡ then 𝜃^‡=𝑔(𝜃^†), see Cox and Hinkley [12, page 159]. For this argument use the orthogonal projection identity to write 𝜃=𝐴𝐴^′𝐴⁻¹𝜉+𝐴_⟂𝐴^′_⟂𝐴_⟂⁻¹𝜑; (9) for unique 𝜉=𝐴^′𝜃 and 𝜑=𝐴^′_⟂𝜃. Thus, if 𝐴^′𝜃^‡=𝐴^′𝜃^† then 𝜃^‡=𝑔(𝜃^†) with 𝜁=𝜑^‡−𝜑^†=𝐴^′_⟂(𝜃^‡−𝜃^†).

In applications it can be difficult to find a basis 𝑋 for the design 𝐷. It can be easier to find a group 𝑔 and hence 𝐴_⟂ and then use this information to construct 𝐴 and a candidate basis 𝑋=𝐷𝐴, noting that 𝐷=𝑋𝐴^′. This argument leaves it to be proven that 𝑋 is a basis, or equivalently, that the suggested group 𝑔 actually describes the equivalence classes of the mapping from 𝜃 to 𝜇.

It is useful to note that in the choices of 𝑋, 𝐴 only the spaces spanned by them are unique since 𝑋𝐴^′=𝑋𝑚𝑚⁻¹𝐴^′ for any invertible 𝑚∈ℝ^𝑝×𝑝. Likewise, the maximal invariant 𝜉 is only unique up to bijective transformations. This lack of uniqueness has no impact on the analysis of the likelihood albeit it influences interpretations.

2.4. Statistical Analysis Using the Maximal Invariant Parameter

The statistical model parametrised with the maximal invariant parameter 𝜉 can be analysed by standard statistical techniques. This contrasts to a range of problems that arise when working with an ad hoc identified time effect 𝜃. In the following the relatively simple standard statistical analysis of the model parametrised by 𝜉 is discussed with respect to likelihood theory, interpretation, plots, hypothesis testing, forecasting, and Bayesian analysis. In Sections 3 and 4 we give an overview of the much more complicated theory underpinning models parametrised by the ad hoc identified time effect 𝜃. Age-period-cohort examples follow in Section 5.

2.4.1. Exponential Family Theory

Suppose the likelihood is drawn from a generalized linear model based on an exponential family. Then the model is actually a regular exponential family where the maximal invariant parameter 𝜉 is the canonical parameter since it is freely varying in a real space; see Barndorff-Nielsen [13, page 116]. This opens up for a wealth of convenient statistical properties such as a likelihood equation with a simple expression and explicit conditions for a unique solution. In contrast, ad hoc identified parameters are based on an injective mapping of the canonical parameter 𝜉 into 𝜃; see Sections 3.1 and 3.2. It is then more difficult to fully exploit the exponential family theory.

2.4.2. Interpretation and Plots

The maximal invariant parameter 𝜉 varies freely in ℝ^𝑝. It can therefore be interpreted as the parameter of any standard statistical model. Since 𝜉 is freely varying the coordinates of 𝜉 can be interpreted independently. When 𝜃 is a collection of time effects then 𝜉 can be organised as a collection of time series. Since the coordinates of 𝜉 are freely varying the time series plots of the components of 𝜉 have the usual interpretation of time series. In contrast, ad hoc identified estimators are constrained to a 𝑝-dimensional subspace Θ_𝜆 of Θ=ℝ^𝑞, which is often affine but can be more complicated. A consequence is that plots are complicated to evaluate; see Section 3.4.1.

2.4.3. Hypothesis Testing

Hypotheses are easily formulated and analysed when using the maximal invariant parametrisation. An affine hypothesis that restricts 𝜉 to vary in a 𝑝_𝐻-dimensional affine subspace can be formulated as 𝐻^′𝜉=𝜂 for known matrices 𝐻∈ℝ^{𝑝×(𝑝−𝑝𝐻)}, 𝜂∈ℝ^{𝑝−𝑝𝐻}. This implies a restriction on the predictor 𝜇=𝑋𝜉 of (6). Form the orthogonal complement 𝐻_⟂ and recall the orthogonal projection identity 𝐼_𝑛=𝐻𝐻^′+𝐻_⟂𝐻^′_⟂ so that 𝜇=𝑋𝐻𝐻^′𝜉+𝑋𝐻_⟂𝐻^′_⟂𝜉. Introduce a 𝑝_𝐻-dimensional parameter 𝜑=𝐻^′_⟂𝜉, a design matrix 𝑋_𝐻=𝑋𝐻_⟂, and an offset 𝑍_𝐻=𝑋𝐻𝜂. The restricted parameter space is 𝑀_𝐻=𝜇∈ℝ^𝑛∶𝜇=𝑋_𝐻𝜑+𝑍_𝐻for𝜑∈ℝ^𝑝𝐻. (10) In an exponential family context both the unrestricted model and the restricted model form regular exponential families. A variety of nice properties then follow for the estimators and the test statistics from the exponential family theory. Examples are given in Sections 5.3 and 5.5.3. In contrast, the hypothesis derived from restrictions on ad hoc identified parameters and the resulting degrees of freedom are complicated to analyse; see Section 3.4.2.

2.4.4. Forecasting

Most often the objective of a mortality study is to forecast the future mortality. In the linear context, 𝜇=𝑋𝜉, this is done by extending the design 𝑋 and by extrapolating 𝜉.

It is usually easy to extend the design 𝑋 into the forecast horizon. This involves the construction of a triangular block matrix with an appropriate number of extra rows corresponding to the data over the forecast horizon as well as extra columns representing the extra parameters that would be needed: 𝑋^ℎ=𝑋𝑋0_ℎ1𝑋^ℎ₂. (11) Extrapolating 𝜉 into a vector 𝜉̃ then gives the forecast 𝜇=𝑋^ℎ₁,𝑋^ℎ₂𝜉𝜉̃. (12) The extrapolation of the parameter 𝜉 can be done as follows. The estimated parameter, or part of it, can be thought of as a time series. Any forecast techniques from the time series literature applied directly to 𝜉 can be used, subject to the usual contextual considerations.

Ad hoc identified time effects can be extrapolated in a similar way; see Section 3.4.3. This may, however, result in avoidable arbitrary effects in the forecast. Necessary and sufficient conditions for this eventuality are given for age-period-cohort models in Section 5.4.3. The practical examples are mainly Bayesian in nature and are discussed next.

2.4.5. Bayesian Analysis

The introduction of the canonical parameter shows that the likelihood, in Bayesian notation, is of the form 𝑝(𝑦∣𝜃)=𝑝(𝑦∣𝜉) where 𝜉 is freely varying. A purist Bayesian analysis can simply introduce a prior on the canonical parameter, 𝑝(𝜉). This is updated in a straight forward way, resulting in the posterior 𝑝(𝜉∣𝑦)=𝑝(𝑦∣𝜉)𝑝(𝜉)/𝑝(𝑦).

In contrast, introducing a prior on ad hoc identified parameters gives various difficulties. Only parts of the prior are updated by the likelihood, so that it becomes unclear which information arises from the data and which information arises from the ad hoc identification. Moreover, avoidable arbitrariness is introduced in the forecast; see Section 4. Introduction of hyperparameters exacerbates the issue. Examples are given in Sections 5.4.4, 5.5.2, and 6.1.6.

3. Working with the Time Effects

In Section 2 we considered the situations where estimation, hypothesis testing a hypothesis, or forecasting is carried out using the canonical parameter. However, there might be situations, where the original time effect parametrisation is preferred, perhaps because it is felt that this parametrisation is particularly helpful in guiding the intuition. This requires ad hoc identification of the time effect. In this section we will guide the considerations a statistician that has to go through when insisting on an analysis based on some nonunique parametrisations. As in Section 2 we focus on linearly parametrised models. Specific examples follow in Sections 5 and 6.

In Section 3.1 ad hoc identification is defined. As an example we consider a least squares estimation problem with collinear regressors in Section 3.2. For the age-period-cohort model reviewed in Section 5 it is common to ad hoc identify in two steps: first identifying levels then the linear trends. We consider such two-step ad hoc identification in Section 3.3. The consequence of ad hoc identification is considered in Section 3.4. Indeed, when forecasting the time effect, we do not want the forecast to depend on the identification scheme. The same applies to graphical visualisation of our data, where the eye may extract patterns that depend on the identification scheme. Likewise, confusion may arise when formulating a hypothesis directly on the time effect parameters.

3.1. Ad Hoc Identification

In this section the time effect parametrisation is considered. An identification scheme has to be introduced when working with the time effects. This may rest on mathematical convenience or it may be chosen for a particular purpose given the substantive context. We therefore call it ad hoc identification. Here we consider a simple identification scheme but turn to a more common two-step identification scheme in Section 3.3.

Once the canonical parameter 𝜉 has been estimated there is often a wish to return to the original time effect 𝜃. The two are linked through the surjective mapping 𝜃⟼𝜉=𝐴^′𝜃, (13) from Θ=ℝ^𝑞 to Ξ=ℝ^𝑝. Indeed, since 𝜉 is constructed as a function of 𝜃 the notation for 𝜉 is often chosen to reflect 𝜃. The canonical parameter 𝜉 does, however, only give partial information about 𝜃. The remaining part, say 𝜆, of 𝜃 will have to be chosen by the researcher and combined with 𝜉.

A linear ad hoc identification of 𝜃 comes about by the researcher choosing a constraint 𝐿^′𝜃=𝜆 (14) for some known 𝜆∈ℝ^𝑞−𝑝 and some matrix 𝐿∈ℝ^{𝑞×(𝑞−𝑝)} chosen so the square matrix (𝐴,𝐿) is invertible. The time effect space Θ is now reduced to an affine subspace Θ_𝜆=𝜃_𝜆∈Θ∶𝐿^′𝜃_𝜆=𝜆. (15) Given 𝜃 we can find 𝜉, 𝜆 through (13) and (14) as (𝜉^′,𝜆^′)^′=(𝐴,𝐿)^′𝜃. At the same time, given values of 𝜉, 𝜆 and the invertibility of (𝐴,𝐿), the ad hoc identified time effect is found through 𝜃_𝜆=𝐴^′𝐿^′−1𝜉𝜆=𝐿_⟂𝐴^′𝐿_⟂⁻¹𝜉+𝐴_⟂𝐿^′𝐴_⟂⁻¹𝜆. (16) In this notation a subindex 𝜆 is introduced to avoid confusion with the time effect 𝜃 in the original mortality model. Indeed, there are now four different parameters in play, namely, the original time effect 𝜃∈Θ, the predictor 𝜇∈𝑀, the maximal invariant parameter 𝜉∈Ξ and the ad hoc identified time effect 𝜃_𝜆∈Θ_𝜆, each of which has a different interpretation. The mapping from 𝜃 to each of 𝜇, 𝜉, and 𝜃_𝜆 is surjective, while there are bijective mappings between the latter three. The interpretations of the time effect 𝜃 and the canonical parameter 𝜉 will inevitably be different. For a start they have different dimensions. Endowing the spaces with Euclidean norms shows that distances in the two spaces Θ and Ξ will be judged differently. The time effect 𝜃 and the ad hoc identified time effect 𝜃_𝜆 will similarly have different interpretations. Although they have the same dimensions the Euclidean norms on Θ and Θ_𝜆 will be rather different. Confusion may arise in the interpretation of a mortality analysis if there is no clear distinction between 𝜃 and 𝜃_𝜆. In addition an unnecessary arbitrariness may arise when making inference on 𝜃_𝜆 or extrapolating 𝜃∈Θ_𝜆. We will return to these issues in Section 3.4.

It is perhaps interesting to note that despite the linear parametrisation the ad hoc identification need not be done in a linear fashion as in (14). Indeed it is common for Poisson models with a log link to ad hoc identify 𝜃 through the original multiplicative scale. That means that the ad hoc identification is done nonlinearly through 𝜉=𝐴^′𝜃_𝜆,𝜆=𝑓𝜃_𝜆. (17)

The fit of the model is unaffected by the ad hoc identification. Indeed the fit is measured in terms of the estimate of the predictor 𝜇=𝐷𝜃_𝜆 where 𝐷=𝑋𝐴^′. Since the identification is made so 𝜉=𝐴^′𝜃_𝜆; the estimated predictor reduces to 𝜇=𝐷𝜃_𝜆=𝑋𝐴^′𝜃_𝜆=𝑋𝜉̂, (18) regardless of the choice of ad hoc identification.

3.2. A Least Squares Example

As an illustration of estimation in the presence of ad hoc identification consider a normal likelihood. Different, but equivalent, expressions can be found depending on the parametrisation. The likelihood of the predictor 𝜇 is 𝖫𝜇,𝜎²;𝑌=2𝜋𝜎^2−𝑛/2−1exp2𝜎²𝑌−𝜇^′𝑌−𝜇for𝜇∈𝑀,𝜎²>0. (19) Rewriting it in terms of the canonical parameter it is 𝖫𝜉,𝜎²;𝑌=2𝜋𝜎^2−𝑛/2−1exp2𝜎²(𝑌−𝑋𝜉)^′(𝑌−𝑋𝜉)for𝜉∈Ξ=ℝ^𝑝,𝜎²>0, (20) while introducing the time effect parameter gives 𝖫𝜃,𝜎²=;𝑌2𝜋𝜎^2−𝑛/2−1exp2𝜎²𝑌−𝑋𝐴^′𝜃^′𝑌−𝑋𝐴^′𝜃for𝜃∈Θ=ℝ^𝑞,𝜎²>0. (21) The likelihood (20) of the canonical parameter 𝜉 can be analysed by the least squares method since the design 𝑋 has full column rank. The maximum likelihood estimator for 𝜉 and the predictor for the data are 𝜉̂=𝑋^′𝑋⁻¹𝑋^′𝑌,𝑌=𝑋𝜉̂=𝑋𝑋^′𝑋⁻¹𝑋^′𝑌. (22) Along with the residual variance this is all the information that is given by the likelihood.

The likelihood (21) of the time effect 𝜃 only depends on 𝜃 through 𝜉=𝐴^′𝜃. The lack of identification means that the maximum likelihood estimator for 𝜃 has an arbitrary element, so that it is a set valued estimator. Based on (16) this can be expressed by Θ_𝜆=𝐿_⟂𝐴^′𝐿_⟂⁻¹𝜉̂+𝐴_⟂𝐿^′𝐴_⟂⁻¹𝜆where𝜃_𝜆∈Θ_𝜆⊂Θ, (23) for any 𝐿 so (𝐴,𝐿) is invertible and for any 𝜆∈ℝ^𝑞−𝑝. The fit, however, remains the same and (18) becomes 𝜇=𝐷𝜃_𝜆=𝑋𝐴^′𝐿_⟂𝐴^′𝐿_⟂⁻¹𝜉̂+𝐴_⟂𝐿^′𝐴_⟂⁻¹𝜆=𝑋𝜉̂=𝑌. (24) In order to compute actual estimates then 𝐿, 𝜆 have to be chosen, which amounts to ad hoc identification. For instance, with the ad hoc identifying restrictions 𝐿=𝐴_⟂ and 𝜆=0 then 𝜃_𝜆 can be thought of as the least squares estimator of 𝑌 on 𝐷 using the Moore-Penrose generalised inverse for the singular matrix 𝐷^′𝐷; see Searle [14, page 212]. See Section 5.4.1 for an example.

3.3. Step-Wise Identification

It is common to ad hoc identify parameter in a step-wise fashion. In the first step the time effect parameter is only partially constrained. The full identification then follows in a second step. An example is given in Section 5.4.1 for an age-period-cohort model in which the levels of the time effects are constrained in the first step leaving the ad hoc identification of the linear trends to the second step.

The first step constraints are affine of the type 𝐶^′𝜃_𝐶=𝜓, (25) for known matrices 𝐶∈ℝ^{𝑞×(𝑞−𝑞𝐶)}, 𝜓∈ℝ^{𝑞−𝑞𝐶}. The constrained time effect space is then Θ_𝐶=𝜃_𝐶∈Θ∶𝐶^′𝜃_𝐶=𝜓. (26) Thereby the 𝑞-dimensional time effect space Θ is reduced to a 𝑞_𝐶-dimensional variation. The properties of this partially ad hoc identified parameter space depends on the rank of the matrix (𝐴,𝐶). If the number of constraints, 𝑞−𝑞_𝐶, is at most equal to the number of unidentified components 𝑞−𝑝, it is possible that (𝐴,𝐶) has full column rank. In that case the constraint implies a partial ad hoc identification without constraining the parameter space 𝑀 of the statistical model. This is shown in Theorem 1; see also Section 5.4.1 for an example, while the proof is given in the Appendix. When (𝐴,𝐶) has reduced rank the parameter space 𝑀 is also constrained; see Section 3.4.2 for a discussion.

Theorem 1. Suppose (𝐴,𝐶) has full column rank. Then the matrix 𝑚=𝐴^′_⟂𝐶∈ℝ^{(𝑞−𝑝)×(𝑞−𝑞𝐶)} has full column rank and the constraint (25) does not constrain the canonical parameter 𝜉 and the predictor 𝜇. Hence, the predictor space remains of the form (2). The equivalence classes in Θ_𝐶 under the mapping 𝜃↦𝜇=𝑋𝐴^′𝜃 are given by the group 𝑔_𝐶∶𝜃⟼𝜃+𝐴_⟂𝑚_⟂𝜁, (27) for arbitrary 𝜁∈ℝ^{𝑞𝐶−𝑝} where 𝑚_⟂∈ℝ^{(𝑞−𝑝)×(𝑞𝐶−𝑝)} is the orthogonal complement of 𝑚. The maximal invariant remains 𝜉=𝐴^′𝜃.

The partial ad hoc identification by (25) implies that any time series analysis of the time effects has to happen relative to the constrained space Θ_𝐶 rather than the space Θ. This is awkward as discussed in Section 3.4 below. It is also considerably more complicated than working with the freely varying canonical parameter 𝜉; see Section 2.4.2.

3.4. Consequences of Ad Hoc Identification

In the following we will look closer at the consequences of working with the ad hoc identified time effect parameter 𝜃 in the context of a linear mortality model of the form 𝜇=𝐷𝜃. We consider the consequences for plotting, hypothesis testing, and forecasting.

3.4.1. Plots of Time Effects

In the mortality model (1) the time effect 𝜃 is the concatenation of age, period, and cohort effects. It seems natural to think of these individual time effects as time series and to plot them against time. As the time effect 𝜃 varies in the unrestricted space Θ=ℝ^𝑞 this maps the 𝑞-vector into unrestricted time series.

Estimates of the time effects are constructed by combining an estimate of 𝜉 with an ad hoc chosen value for 𝜆=𝐿^′𝜃, see (14). The resulting estimate 𝜃_𝜆 is therefore constrained to the space Θ_𝜆⊂Θ. The interpretation of the estimate 𝜃_𝜆 is therefore different from the interpretation of the original time effect 𝜃. Distances on the spaces Θ and Θ_𝜆 are judged differently and the variability of 𝜃_𝜆 is deduced exclusively from 𝜉̂ through (16). The time series components of 𝜃_𝜆 are now restricted through 𝜆=𝐿^′𝜃_𝜆. Plots of the 𝜃_𝜆-time series are therefore interpreted differently from the imagined plots of the original 𝜃-time series and from the plots of the maximal invariant parameter 𝜉 discussed in Section 2.4.2. Indeed, if one were to analyse the estimated 𝜃_𝜆-time series statistically the linear constraint should be taken into account. This is a bit complicated as illustrated below, but it is the consequence of working with the ad hoc identified parameter 𝜃_𝜆 rather than the canonical parameter 𝜉.

Attempts to give intrinsic meaning to 𝜆 will be specific to the index set for the data set at hand. For instance, the requirement that the age effect should be zero on average does not carry over when looking at a subsample or when forecasting. It is not obvious that such an ad hoc identification is any more or less arbitrary than saying that, for instance, the first or the last age effect should have a particular value.

Adding confidence bands to a plot of 𝜃_𝜆 is in itself not difficult. If 𝜉̂ is asymptotically normal with mean 𝜉 and variance Σ, then 𝜃_𝜆 is asymptotically normal with mean 𝜃_𝜆 and variance 𝐿_⟂(𝐴^′𝐿_⟂)⁻¹Σ(𝐿^′_⟂𝐴)⁻¹𝐿^′_⟂. This is a normal distribution on the space Θ_𝜆. The interpretation of these standard errors will therefore be similar to that of 𝜃_𝜆 itself.

Finally, it may be of interest to analyse the estimated 𝜃_𝜆-time series statistically. Denote this time series by 𝑥_𝜆. Its sample space is now Θ_𝜆. A statistical model on Θ_𝜆 can be built as follows. The starting point could be a time series model for unrestricted variables 𝑥 on the sample space Θ. This gives a joint density for 𝑥∈Θ, which can be reduced by marginalisation to a density for 𝑥_𝜆∈Θ_𝜆. Whether one is working with the unrestricted model for 𝑥∈Θ or the restricted model for 𝑥_𝜆∈Θ_𝜆 inferences that are invariant to 𝑔 must be based on those statistics of 𝑥 or 𝑥_𝜆 that are invariant to 𝑔. Thus, inferences must be based on the maximal invariant under 𝑔. For a general overview of invariant reduction see Cox and Hinkley [12, page 175f], whereas Nielsen [15] gives the argument in some detail for an autoregression with a linear trend.

3.4.2. Hypothesis Testing

Having formulated the model in terms of time effects it may be of interest to test the hypothesis that one of these time effects is absent. No identification issues arise when the hypothesis is formulated as a restriction on the canonical parameter 𝜉 as discussed in Section 2.4.3. But one has to be careful when formulating hypotheses in terms of the original time effect. See Sections 5.4.5, 5.5.3, and 5.5.4 for examples.

Affine hypotheses on the time effect are of the form 𝑅^′𝜃_𝑅=𝜌, (28) for known matrices 𝑅∈ℝ^{𝑞×(𝑞−𝑞𝑅)}, 𝜌∈ℝ^{𝑞−𝑞𝑅}. The constrained time effect space is then Θ_𝑅=𝜃_𝑅∈ℝ^𝑞∶𝑅^′𝜃_𝑅=𝜌. (29) To see how the restriction (28) restricts the predictor space 𝑀⊂ℝ^𝑛 recall that the predictor 𝜇 only depends on 𝜃 through 𝜉=𝐴^′𝜃. Thus, the analysis of the restriction (28) depends on the interplay between the matrices 𝐴, 𝑅. Theorem A.3 in Appendix A.3 gives a general result to that effect. It shows that the hypothesis (28) restricts the predictor space 𝑀 to a 𝑝_𝑅-dimensional affine subspace of ℝ^𝑛 in so far as it restricts the canonical parameter 𝜉. In particular, the degrees of freedom of the hypothesis, 𝑝−𝑝_𝑅, may in general be different from the dimension reduction of the time effect parameter, 𝑞−𝑞_𝑅. When this is the case the restriction (28) has an element of ad hoc identifying the time effect.

3.4.3. Forecasts

Forecasts can be made by extrapolating the ad hoc identified time effects 𝜃_𝜆. Two researchers choosing different ad hoc identification schemes, but otherwise making the same analysis, may make different forecasts. This can be avoided if the extrapolation method is chosen with some care.

Following the linear approach outlined in Section 2.4.4 the predictor 𝜇=𝐷𝜃=𝑋𝐴^′𝜃 is forecasted by extending the design 𝐷 into 𝐷^ℎ=𝐷𝐷0_ℎ1𝐷^ℎ₂. (30) Extrapolating the ad hoc identified 𝜃_𝜆 into a vector (𝜃^′_𝜆,𝜃^′_𝜆)^′ then gives the forecast 𝜇=𝐷^ℎ₁,𝐷^ℎ₂𝜃_𝜆𝜃_𝜆=𝐷^ℎ₁𝜃_𝜆+𝐷^ℎ₂𝜃_𝜆. (31) Often both components 𝐷^ℎ₁𝜃_𝜆 and 𝐷^ℎ₂𝜃_𝜆 depend on the ad hoc identification. Nonetheless, these dependencies of ad hoc identification may cancel each other so that the overall forecast 𝜇 is invariant to the ad hoc identification. Such invariance would seem desirable in most applications unless there is strong substantial reason for the ad hoc identification scheme. Necessary and sufficient conditions for invariance are presented for the age-period-cohort model in Section 5.4.3 and for a nonlinear model in Section 6.1.5.

In contrast, these considerations are redundant when working with the canonical parameter, 𝜉; see Section 2.4.4.

4. Bayesian Models and Random Effects Models

Mortality analysis is often carried out using either Bayesian methods or random effects methods. The mortality model is then altered through the introduction of a prior distribution on the parameters. One might think that the identification problems become less of an issue or even disappear. This is not the case since the Bayesian method and the random effects method is based on the mortality likelihood which only depends on the time effect 𝜃 through the maximal invariant parameter 𝜉. Thus, the identification challenges remain. The issue is that a prior on the unidentified part, say 𝜆, of the time effect amounts to an ad hoc identification. Indeed, the conditional prior of 𝜆 given 𝜉 is not updated by the mortality likelihood. A main difference is that a maximum likelihood analysis of the original mortality likelihood usually prompts the researcher when there is an identification issue, whereas both Bayesian methods and random effects methods allow computations to go through despite an identification issue.

In Section 4.1 it is seen that introduction of a conditional prior on 𝜆 given 𝜉 is the Bayesian analogue of ad hoc identification. This leads to the same type of forecasting challenges as in the frequentist settings as is seen in Section 4.2. In Section 4.3 we show how the Bayesian identification issues transfer to random effects models.

4.1. Bayesian Estimation

For Bayesian and random effects models we formulate a likelihood and a prior. Thus, consider a likelihood 𝑝(𝑦∣𝜃)=𝐿(𝜃;𝑦). Replacing 𝜃 by 𝜉, 𝜆 the identification problem implies that 𝑝𝑦∣𝜉,𝜆=𝑝𝑦∣𝜉foralloutcomes𝑦. (32) The prior on 𝜃 is factorised as 𝑝(𝜃)=𝑝(𝜉,𝜆)=𝑝(𝜉)𝑝(𝜆∣𝜉). In the case of Bayesian estimation the following result emerges.

Theorem 2. Suppose the likelihood satisfies (32). Then (i)the predictive distribution does not depend on the conditional prior for 𝜆: 𝑝𝑦=𝑝𝑦∣𝜉𝑝(𝜉)𝑑𝜉; (33)(ii)the posterior satisfies 𝑝=𝑝𝑝(𝜉)𝜉∣𝑦𝑦∣𝜉𝑝𝑦,𝑝𝜆∣𝜉,𝑦=𝑝(𝜆∣𝜉); (34)(iii)the posterior means satisfy 𝐸=𝐸(),𝐸=𝐸()𝑝𝜉∣𝑦𝜉𝑝𝜉∣𝑦𝑑𝜉,𝜆∣𝜉,𝑦=𝐸𝜆∣𝜉𝜆∣𝑦𝜆∣𝜉𝜉∣𝑦𝑑𝜉. (35)

Theorem 2 shows that it suffices to give a prior to 𝜉 and ignore 𝜆 as advocated in Section 2.4.5. Indeed the conditional prior for 𝜆 given 𝜉 is not updated. Theorem 2 appears to be well-known; see Poirier [16, Proposition 2] or Smith [17, Section B].

Due to Theorem 2 the Bayesian analyst faces the complications outlined in Section 3.4. Indeed, suppose that two Bayesian researchers choose the same likelihood 𝑝(𝑥∣𝜉,𝜆)=𝑝(𝑥∣𝜉) and the same prior 𝑝(𝜉) for 𝜉, but different conditional priors for 𝜆 given 𝜉. Their marginal distributions for the data are identical, but any inferences regarding interpretation or forecasting will differ in so far as they involve the unidentified parameter 𝜆. A Bayesian researcher should therefore be cautious with inference related to 𝜆. There will of course be situations where the prior knowledge of 𝜆 given 𝜉 is found to be of substantive relevance. In such situations it seems more fruitful to change the likelihood to include that information.

4.2. Forecasting

Bayesian forecasts involve integrating an extrapolative distribution. This can be done in two ways, either working exclusively with the identified, maximal invariant parameter 𝜉 as in Section 2.4.4, or working with the time effect 𝜃=(𝜉,𝜆) as in Section 3.4.3.

4.2.1. Forecasting Using the Maximal Invariant Parameter

Consider first the case where only the maximal invariant parameter 𝜉 is used. In that case the forecast is computed by sampling from the posterior 𝑝(𝜉∣𝑦) and then extrapolating 𝜇 using the sampled value 𝜉 using some extrapolative methods, say 𝑝(𝜇∣𝜉,𝑦). In combination this gives the forecast 𝑝𝜇∣𝑦=𝑝𝜇∣𝜉,𝑦𝑝𝜉∣𝑦𝑑𝜉. (36)

4.2.2. Forecasting Using the Ad Hoc Identified Time Effect

Consider now forecasts involving the full time effect 𝜃=(𝜉,𝜆). Theorem 2(ii) shows that the posterior satisfies 𝑝(𝜃∣𝑦)=𝑝(𝜉∣𝑦)𝑝(𝜆∣𝜉). The distribution forecast with extrapolation 𝑝(𝜇∣𝜉,𝜆,𝑦) is then 𝑝𝜇∣𝑦=𝑝𝜇∣𝜉,𝜆,𝑦𝑝𝜉∣𝑦𝑝(𝜆∣𝜉)𝑑𝜆𝑑𝜉. (37) The concern is now as follows. Suppose a second researcher chooses the same extrapolative method, likelihood, and prior for 𝜉, but different conditional priors 𝑝^†(𝜆∣𝜉). In general, this will result in a different distribution forecast: 𝑝^†𝜇∣𝑦=𝑝𝜇∣𝜉,𝜆,𝑦𝑝𝜉∣𝑦𝑝^†(𝜆∣𝜉)𝑑𝜆𝑑𝜉. (38) The question is then under which conditions will 𝑝(𝜇∣𝑦)=𝑝^†(𝜇∣𝑦) so that the distribution forecasts are invariant to the choice of conditional prior for 𝜆 given 𝜉? A sufficient condition is that the extrapolation method does not depend on 𝜆 so 𝑝𝜇∣𝜉,𝜆,𝑦=𝑝𝜇∣𝜉,𝑦. (39) Condition (39) could alternatively be expressed as requiring that the forecast 𝑝(𝜇∣𝜃,𝑦)=𝑝(𝜇∣𝜉,𝜆,𝑦) is invariant to the group 𝑔 acting on the time effect space Θ so that 𝑝(𝜇∣𝜉,𝜆,𝑦)=𝑝{𝜇∣𝜉,𝑔(𝜆),𝑦}.

Theorem 3. Suppose that the likelihood satisfies (32) and the priors are probabilities. If the extrapolative distribution does not depend on 𝜆 so (39) holds; then the forecast distribution 𝑝(𝜇∣𝑦) computed in (37) is invariant to the choice of conditional prior for 𝜆 given 𝜉. The forecast then reduces to (36).

To summarise, the identification issues surrounding Bayesian analysis are similar to those outlined in the previous sections. Examples of the problems that can arise are discussed in Sections 5.4.4, 5.5.2, and 6.1.6. There are two solutions to the identification problem. The first is only to formulate a prior on 𝜉; see Section 2.4.5. Incidentally, this is what Bernardo and Smith [18, page 218] do in their discussion of the two-way analysis of variance, albeit without linking it to the considerations of Smith [17]. The prior 𝑝(𝜉) can of course be constructed by formulating a prior on 𝜃 and then reduce it to a prior on 𝜉 by marginalisation so 𝑝(𝜉)=∫𝑝(𝜉,𝜆)𝑑𝜆. The other solution is to work with a prior on 𝜃 but avoid those parts of the posterior that depend on 𝜆.

4.3. Random Effects Models

It is common to combine mortality models with a random effects approach, which effectively forms a new model. An example is given in Section 5.4.6. We consider the consequence of the lack of identification.

The random effect models are typically constructed as follows. Suppose the density of the data 𝑦 given the time effects 𝜃=(𝜉,𝜆) is of the form 𝑝(𝑦∣𝜉,𝜆)=𝑝(𝑦∣𝜉) as before; see (32). A prior 𝑝(𝜃∣𝜓) is chosen that now depends on a parameter 𝜓. The prior can be decomposed as 𝑝(𝜃∣𝜓)=𝑝(𝜉∣𝜓)𝑝(𝜆∣𝜉,𝜓). Theorem 2 implies that the density of the data 𝑦 given 𝜓 is 𝑝𝑦∣𝜓=𝑝𝑦∣𝜉𝑝𝜉∣𝜓𝑑𝜉. (40) This in turn is used to form the random effects likelihood of 𝜓 as 𝐿_RE𝜓∣𝑦=𝑝𝑦∣𝜓. (41) This, effectively, defines a new model. The random effects likelihood only depends on the prior 𝑝(𝜃∣𝜓) through 𝑝(𝜉∣𝜓). Two researchers choosing the same prior 𝑝(𝜉∣𝜓) but different conditional priors 𝑝(𝜆∣𝜉,𝜓) will then get the same random effects likelihood and the same maximum likelihood estimator 𝜓.

In mortality modelling it is common to go one step further and estimate the time effects 𝜃 through the mean of the posterior 𝑝(𝜃∣𝑦,𝜓) evaluated at 𝜓=𝜓. Then the identification problem may show up. Theorem 2 shows that 𝑝=𝑝𝑝𝜓𝜉∣𝜓,𝑦𝑦∣𝜉𝜉∣𝑝𝜓,𝑝𝜓,𝑦∣𝜆∣𝜉,𝜓,𝑦=𝑝𝜆∣𝜉, (42) so that the prior for 𝜉 is updated, while the conditional posterior for 𝜆 given 𝜉 is not updated by the data. Thus, in general the estimate for 𝜃 is based, in part, on a prior which is not updated by the data.

5. Age-Period-Cohort Models

We will now apply the theoretical considerations to analyse the age-period-cohort model. The methodological literature on this model is large and the consequences of the above theory are wide ranging.

In Section 5.1 we present the age-period-cohort model along with the maximal invariant parameter. This maximal invariant parameter is also called the canonical parameter because the age-period-cohort model is usually implemented as an exponential family; see Section 2.4.1. When formulating the model we choose a notation matching the age-period-cohort literature rather than the reserving literature. At the same time the exposition takes it starting point in Kuang et al. [19], but the notation deviates.

The implementation of the canonical parameter depends on the type of data array. In Section 5.2 design matrices are given for age-cohort, age-period, and period-cohort data arrays. While they illustrate interesting differences in the structure for these data arrays, they also provide the basis for an immediate implementation via any generalised linear model software. The age-cohort model is expressed as a hypothesis of the age-period-cohort model in Section 5.3. Time effects and forecasting are considered in Section 5.4, while the two-sample age-period-cohort model is discussed in Section 5.5.

5.1. The Model and the Canonical Parameter

Here the age-period model is set up and a quite general identification result is reported.

Consider data 𝑌_𝑖𝑗 indexed by (𝑖,𝑗)∈ℐ where 𝑖 is the age and 𝑗 is the period. The index set may be a rectangle given by 𝑖=1,…,𝐼 and 𝑗=1,…,𝐽 so that the cohort 𝑘=𝐼−𝑖+𝑗 runs from 1 to 𝐾=𝐼+𝐽−1. More generally, the index set could be a generalized trapezoid where two corners are cut off the rectangle so that the cohort 𝑘 runs from 1+ℎ₁ to 𝐼+𝐽−1−ℎ₂ for some ℎ₁,ℎ₂≥0. The class of generalized trapezoids includes the three types of Lexis diagrams discussed by Keiding [20]. We will return to those special cases below.

The statistical model is defined by the assumption that the variables 𝑌_𝑖𝑗 are independent with an exponential family distribution with predictor 𝜇_𝑖𝑗 given by 𝜇_𝑖𝑗=𝛼_𝑖+𝛽_𝑗+𝛾_𝑘+𝛿for𝑖,𝑗∈ℐ. (43) The time effect 𝜃=(𝛼₁…,𝛼_𝐼,𝛽₁,…,𝛽_𝑗,𝛾_ℎ1+1,…, 𝛾_{𝐼+𝐽−1−ℎ2},𝛿)^′ now varies in some time effect space Θ∈ℝ^𝑞 where 𝑞=𝐼+𝐽+𝐾+1−ℎ₁−ℎ₂.

The model (43) is of the form (1) discussed in Section 2. Specifically, the predictors 𝜇_𝑖𝑗 can be stacked in a vector 𝜇 of dimension 𝑛=dimℐ and written as 𝜇=𝐷𝜃. Thus, the parameter space for the model is of the form 𝑀=(𝜇∈ℝ^𝑛∶𝜇=𝐷𝜃 for 𝜃∈Θ) as outlined in (2). The mapping 𝜃↦𝜇 from Θ to 𝑀 is surjective and the equivalence classes in the time effect space can be described by a group of transformations that are discussed in (8). This group can be represented as 𝛼𝑔∶^𝑖𝛽_𝑗𝛾_𝑘𝛿⟼𝛼_𝑖()𝑑𝛽+𝑎+𝑖−1_𝑗𝑑𝛾+𝑏−𝑗−1_𝑘()𝑑()𝑑+𝑐+𝑘−1𝛿−𝑎−𝑏−𝑐−𝐼−1for𝜃∈Θ, (44) for any 𝑎, 𝑏, 𝑐, and 𝑑. This is of the form (8) with 𝜁=(𝑎,𝑏,𝑐,𝑑)^′ although the definition of the matrix 𝐴 depends on the structure of the index set ℐ.

A first clue for the canonical parametrisation is given by Fienberg and Mason [21] and Clayton and Schiffler [22] who pointed out that, on the multiplicative scale, ratios of relative risks are invariant. On the additive scale this amounts to looking at second differences, such as Δ²𝛼_𝑖=𝛼_𝑖−2𝛼_𝑖−1+𝛼_𝑖−2. A graphical illustration of the double differences is given in Figure 1 (graphics were done using R 3.0.2, see R Development Core Team [23]), which is taken from Miranda et al. [24]. Panel (a) illustrates the interpretations of the formula for Δ²𝛼_𝑖 as follows. Consider the 1970 and 1971 cohorts. In 2010 these have ages 40 and 39, while in 2011 these have ages 41 and 40. Thus, Δ²𝛼₄₁ represents the increase in mortality from ages 40 to 41 in 2011 relative to the increase from ages 39 to age 40 in 2010. An equivalent interpretation is that which represents the increase in mortality from ages 40 to 41 for the 1970 cohort relative to the increase from ages 39 to 40 for the 1971 cohort. In a similar way panels (b) and (c) illustrate the formulas for Δ²𝛽₂₀₁₂ and Δ²𝛾₁₉₇₂.

(a)

(b)

(c)

Figure 1 

Illustration of interpretation of Δ2𝛼41, Δ2𝛽2012, and Δ2𝛾1972.

Kuang et al. [19] introduces a parameter formed by these second differences as well as three entries of the predictor; that is, 𝜉=𝜇_𝑖1𝑗1,𝜇_𝑖2𝑗2,𝜇_𝑖3𝑗3,Δ²𝛼₃,…,Δ²𝛼_𝐼,Δ²𝛽₃,…,Δ²𝛽_𝐽,Δ²𝛾_ℎ1+3,…,Δ²𝛾_𝐾−ℎ2. (45) The parameter 𝜉 varies in the space Ξ=ℝ^𝑝 where 𝑝=𝑞−4. If the three points 𝜇_𝑖1𝑗1, 𝜇_𝑖2𝑗2, and 𝜇_𝑖3𝑗3 are chosen not to be linearly related then they define the levels and the linear trends in the predictor. The formal condition is that a certain determinant defined from the indices is nonzero; that is, 𝑖₂𝑗₃−𝑖₃𝑗₂+𝑖₃𝑗₁−𝑖₁𝑗₃+𝑖₁𝑘₂−𝑖₂𝑘₁≠0. (46)

Theorem 4 (see [19], [25, Corollary 2]). Let 𝜇 satisfy (43). If the condition (46) is satisfied then the parameter 𝜉 of (45) satisfies the following:(i)𝜉 is a function of 𝜃 which is invariant to the group 𝑔 in (44);(ii)𝜇 is a function of 𝜉;(iii)the parametrisation of 𝜇 by 𝜉 is exactly identified in that 𝜉^†≠𝜉^‡⇒𝜇(𝜉^†)≠𝜇(𝜉^‡).

Theorem 4 therefore shows that 𝜉 varies freely in Ξ=ℝ^𝑝. Moreover, 𝜉 is a maximal invariant of the mapping 𝑚 from 𝜃 to 𝜇 under the transformations 𝑔. It should be noted that the choice of maximal invariant is not unique. Indeed, any bijective mapping of 𝜉 can serve as maximal invariant. The choice of 𝜉 is convenient since it becomes the canonical parameter in generalized linear models of the exponential family type.

In itself this theorem does not tell how to express the predictor 𝜇 in terms of the canonical parameter 𝜉. The link depends on the structure of the index set ℐ. The above mentioned paper gives implicit expressions for generalized trapezoid index sets. In the following we report explicit expressions for the 3 principal Lexis diagrams.

5.2. Design Matrices for Lexis Diagrams

The link between the canonical parameter 𝜉 and the predictor 𝜇 is analysed for the 3 principal Lexis diagrams. We start with age-cohort data arrays, which were the focus of attention in Kuang et al. [19]. Such arrays are easiest to analyse because all three time scales increase from the point where 𝑖=𝑗=𝑘=1. As a consequence the results are relatively easier for these arrays.

5.2.1. Age-Cohort Data Arrays

Age-cohort data arrays are rectangular in the age and cohort indices and given by ℐ_ac={(𝑖,𝑘)∶𝑖=1,…,𝐼,𝑘=1,…,𝐾}. (47) Consequently, the period index 𝑗=𝑖+𝑘−1 varies over 𝑗=1,…,𝐽=𝐼+𝐾−1. Keiding [20] refers to this Lexis diagram as the first principal set of death.

Age-cohort arrays are in particular used for reserving in general insurance. In that situation, only the triangle 1≤𝑖,𝑗,𝑘≤𝐼 is observed. The issue is to forecast the other triangle in the square 1≤𝑖,𝑘≤𝐼. In the reserving literature these triangles are referred to as the upper and lower triangles, since the cohort axis has reverse order. The two-factor age-cohort model for triangular age-cohort arrays is known as the chain-ladder model; see England and Verrall [26] for an overview. Zehnwirth [27] introduced an age-period-cohort model for such triangular arrays. The identification issue is analysed in detail in Kuang et al. [19, 25]. Subsequently, Kuang et al. [28] analysed the Poisson likelihood, while Kuang et al. [10] give an empirical analysis focusing on forecasting.

The age-period-cohort model for the age-cohort arrays is parametrised by 𝜇_𝑖𝑘=𝛼_𝑖+𝛽_{𝑖+𝑘−1}+𝛾_𝑘+𝛿for𝑖,𝑘∈ℐ_ac. (48) The time effect 𝜃=(𝛼₁,…,𝛼_𝐼,𝛽₁,…,𝛽_𝐽,𝛾₁,…,𝛾_𝐾,𝛿)^′ now varies in Θ=ℝ^2(𝐼+𝐾).

The design matrix linking the canonical parameter 𝜉 in (45) and the predictor 𝜇 is essentially an identity linking the two parameters. A natural choice of the three levels points to the predictors that are 𝜇₁₁, 𝜇₁₂, and 𝜇₂₁. We then get the representation 𝜇_𝑖𝑘=𝜇₁₁+(𝑖−1)𝜇₂₁−𝜇₁₁+(𝑘−1)𝜇₁₂−𝜇₁₁+^𝑖ℓ_ℓ=3ℎ=3Δ²𝛼_ℎ+^𝑗ℓ_ℓ=3ℎ=3Δ²𝛽_ℎ+^𝑘ℓ_ℓ=3ℎ=3Δ²𝛾_ℎ, (49) with the convention that empty sums are zero, and recalling that second differences are defined as Δ²𝛼_𝑖=𝛼_𝑖−2𝛼_𝑖−1+𝛼_𝑖−2 so that ∑^𝑖_ℎ=3Δ²𝛼_ℎ=Δ𝛼_𝑖−Δ𝛼₂ and ∑^𝑖_ℓ=3∑^ℓ_ℎ=3Δ²𝛼_ℎ=𝛼_𝑖−𝛼₁−(𝑖−1)Δ𝛼₂.

The identity (49) is crucial to the understanding of the age-period-cohort model. It shows that the predictor has a single level expressed as 𝜇₁₁, which in turn satisfies 𝜇₁₁=𝛼₁+𝛽₁+𝛾₁+𝛿. The level 𝜇₁₁ is therefore estimable, but the individual levels 𝛼₁, 𝛽₁, 𝛾₁, and 𝛿 are not identifiable from the model. Further, the model has two linear trends, here expressed with slopes 𝜇₂₁−𝜇₁₁ and 𝜇₁₂−𝜇₁₁ in terms of the age and cohort indices. These slopes can be expressed as 𝜇₂₁−𝜇₁₁=Δ𝛼₂+Δ𝛽₂ and 𝜇₁₂−𝜇₁₁=Δ𝛽₂+Δ𝛾₂. They are estimable, but the individual slopes Δ𝛼₂, Δ𝛽₂, and Δ𝛾₂ are not identifiable.

The design matrix now follows from the identity (49) so that the predictor satisfies 𝜇=𝑋𝜉, where 𝜉=𝜇₁₁,𝜇₂₁−𝜇₁₁,𝜇₁₂−𝜇₁₁,Δ²𝛼₃,…,Δ²𝛼_𝐼,Δ²𝛽₃,…,Δ²𝛽_𝐽,Δ²𝛾₃,…,Δ²𝛾_𝐾^′,𝑋_𝑖𝑘=1,(𝑖−1),(𝑘−1),ℎ(𝑖,3),…,ℎ(𝑖,𝐼),ℎ(𝑗,3),…,ℎ(𝑗,𝐽),ℎ(𝑘,3),…,ℎ(𝑘,𝐾)^′, (50)(51) where 𝜉∈ℝ^𝑝, where 𝑝=2(𝐼+𝐾−2) and ℎ(𝑡,𝑠)=max(𝑡−𝑠+1,0).

The identification relies on Theorem 4, which can be specialised to age-cohort arrays as follows.

Theorem 5 (see [19, Theorem 1]). Let 𝜇 satisfy (48). The parameter 𝜉 of (50) satisfies the following:(i)𝜉 is a function of 𝜃 which is invariant to the group 𝑔 in (44);(ii)𝜇 is a function of 𝜉, because of (49);(iii)the parametrisation of 𝜇 by 𝜉 is exactly identified in that 𝜉^†≠𝜉^‡⇒𝜇(𝜉^†)≠𝜇(𝜉^‡).

Theorem 5 in turn implies that the parameter 𝜉 varies freely in Ξ=ℝ^𝑝, while the design matrix 𝑋 given by (51) has full column rank. Originally, the more general Theorem 4 was proved as a corollary to Theorem 5.

5.2.2. Age-Period Arrays

An age-period data array is rectangular in the age and cohort indices and given by ℐ_ap=𝑖,𝑗∶𝑖=1,…,𝐼,𝑗=1,…,𝐽. (52) Consequently, the cohort index 𝑘=𝑗−𝑖+𝐼 varies over 𝑘=1,…,𝐾=𝐼+𝐽−1. Keiding [20] refers to this Lexis diagram as the third principal set of death.

Age-period arrays are commonly used in epidemiology, in mortality analysis, and in sociology. The analysis of identification issue is largely similar to that of age-cohort arrays. However, the representation of the predictor 𝜇 in terms of 𝜉 differs in an intriguing way, because the third time index, the cohort 𝑘, is the difference of the other two indices.

The age-period-cohort model for the age-period arrays is parametrised by 𝜇_𝑖𝑗=𝛼_𝑖+𝛽_𝑗+𝛾_{𝑗−𝑖+𝐼}+𝛿for𝑖,𝑗∈ℐ_ap. (53) The time effect 𝜃=(𝛼₁,…,𝛼_𝐼,𝛽₁,…,𝛽_𝐽,𝛾₁,…,𝛾_𝐾,𝛿)^′ now varies in Θ=ℝ^2(𝐼+𝐽). A representation of the predictor 𝜇 in terms of the canonical parameter 𝜉 is now 𝜇_𝑖𝑗=𝜇_𝐼1+(𝑖−𝐼)𝜇_𝐼1−𝜇_𝐼−1,1+𝑗−1𝜇_𝐼2−𝜇_𝐼1+^𝐼−2_ℓ=𝑖^𝐼−2_ℎ=ℓΔ²𝛼_ℎ+^𝑗ℓ_ℓ=3ℎ=3Δ²𝛽_ℎ+^{𝑗−𝑖+𝐼ℓ}_ℓ=3ℎ=3Δ²𝛾_ℎ+2. (54) The representation (54) differs from that of (49) in a subtle way. The three reference points for the levels of the predictor are chosen in the corner 𝑖=𝐼, 𝑗=1. From this corner period and cohort indices increase, while age decreases. Hence, the age double differences Δ²𝛼_𝑖 are now cumulated backwards. This phenomenon arises because the cohort index is the difference of the principal indices of age and period, whereas for the age-cohort array the period index is the sum of the principal indices of age and cohort. The predictor is now 𝜇=𝑋𝜉 where, with ℎ(𝑡,𝑠)=max(𝑡−𝑠+1,0), 𝜉=𝜇_𝐼1,𝜇_𝐼1−𝜇_𝐼−1,1,𝜇_𝐼2−𝜇_𝐼1,Δ²𝛼₃,…,Δ²𝛼_𝐼,Δ²𝛽₃,…,Δ²𝛽_𝐽,Δ²𝛾₃,…,Δ²𝛾_𝐾^′,𝑋_𝑖𝑗=1,𝑖−𝐼,𝑗−1,ℎ(1,𝑖),…,ℎ(𝐼−2,𝑖),𝑗,3,…,ℎ𝑗,𝐽,ℎ𝑗−𝑖+𝐼,3,…,ℎ(𝑗−𝑖+𝐼,𝐾)^′. (55)(56) The identification relies on Theorem 4. It is specialised to age-period arrays as follows.

Theorem 6 (see [24, Theorem 4.1]). Let 𝜇 satisfy (53). The parameter 𝜉 of (55) satisfies the following:(i)𝜉 is a function of 𝜃 which is invariant to the group 𝑔 in (44);(ii)𝜇 is a function of 𝜉, because of (54);(iii)the parametrisation of 𝜇 by 𝜉 is exactly identified in that 𝜉^†≠𝜉^‡⇒𝜇(𝜉^†)≠𝜇(𝜉^‡).

The group of transformations in (44) can be specialised as 𝛼𝑔∶^𝑖𝛽_𝑗𝛾_{𝑖−𝑖+𝐼}𝛿⟼𝛼_𝑖𝛽+𝑎+𝑑𝑖_𝑗𝛾+𝑏−𝑑𝑗_{𝑗−𝑖+𝐼}+𝑐+𝑑𝑗−𝑖+𝐼𝛿−𝑎−𝑏−𝑐−𝑑𝐼for𝜃∈Θ=ℝ^2(𝐼+𝐽); (57) see, for instance, Carstensen [29]. This is of the form (8) with 𝜁=(𝑎,𝑏,𝑐,𝑑)^′ and 𝐴^′_⟂=^{11⋯1−111⋯1−1}_{11⋯1−112⋯𝐼−1−2⋯−𝐽12⋯𝐾−𝐼}. (58)

5.2.3. Period-Cohort Arrays

An age-period data arrays is rectangular in the age and cohort indices and given by ℐ_pc=𝑗,𝑘∶𝑗=1,…,𝐽,𝑘=1,…,𝐾. (59) Consequently, the age index 𝑖=𝑗−𝑘+𝐾 varies over 𝑖=1,…,𝐼=𝐽+𝐾−1. Keiding [20] refers to this Lexis diagram as the second principal set of death. Age-period arrays are commonly used in prospective cohort studies in epidemiology and in sociology. The analysis is similar to that of age-period arrays when swapping the role of age and cohort.

The age-period-cohort model for the age-cohort arrays is parametrised by 𝜇_𝑗𝑘=𝛼_{𝑗−𝑘+1}+𝛽_𝑗+𝛾_𝑘+𝛿for𝑗,𝑘∈ℐ_ap. (60) The time effect 𝜃=(𝛼₁,…,𝛼_𝐼,𝛽₁,…,𝛽_𝐽,𝛾₁,…,𝛾_𝐾,𝛿)^′ now varies in Θ=ℝ^2(𝐽+𝐾). A representation of the predictor 𝜇 in terms of the canonical parameter 𝜉 is now 𝜇_𝑗𝑘=𝜇_1𝐾+𝑗−1𝜇_2𝐾−𝜇_1𝐾+(𝑘−𝐾)𝜇_1𝐾−𝜇_1,𝐾−1+^{𝑗−𝑘+1ℓ}_ℓ=3ℎ=3Δ²𝛼_ℎ+^𝑗ℓ_ℓ=3ℎ=3Δ²𝛽_ℎ+^𝐾−2_ℓ=𝑖^𝐾−2_ℎ=ℓΔ²𝛾_ℎ+2. (61) Thus, the canonical parameter and the design matrix are given by 𝜉=𝜇_1𝐾,𝜇_2𝐾−𝜇_1𝐾,𝜇_1𝐾−𝜇_1,𝐾−1,Δ²𝛼₃,…,Δ²𝛼_𝐼,Δ²𝛽₃,…,Δ²𝛽_𝐽,Δ²𝛾₃,…,Δ²𝛾_𝐾^′,𝑋_𝑗𝑘=1,𝑗−1,𝑘−𝐾,ℎ𝑗−𝑘+1,3,…,ℎ𝑗−𝑘+1,𝐼,ℎ(𝑗,3),…,ℎ(𝑗,𝐽),ℎ(1,𝑘),…,ℎ(𝐾−2,𝑘)^′. (62)(63) In parallel with Theorem 6 we then have the following identification result.

Theorem 7. Let 𝜇 satisfy (60). The parameter 𝜉 of (62) satisfies the following:(i)𝜉 is a function of 𝜃 which is invariant to the group 𝑔 in (44);(ii)𝜇 is a function of 𝜉, because of (61);(iii)the parametrisation of 𝜇 by 𝜉 is exactly identified in that 𝜉^†≠𝜉^‡⇒𝜇(𝜉^†)≠𝜇(𝜉^‡).

5.3. Expressing the Age-Cohort Model as a Hypothesis

It is often of interest to test the absence of the period effect. An application to analysing asbestos related mortality can be found in Miranda et al. [24].

The hypothesis is that 𝛽₁=⋯=𝛽_𝐽, when expressed in terms of the time effect parameters. The restricted model is given by, with 𝑘=𝑗−𝑖+𝐼, 𝜇^ac_𝑖𝑗=𝛼_𝑖+𝛾_{𝑗−𝑖+𝐼}+𝛿for𝑖,𝑗∈ℐ_ap. (64)

The identification problem simplifies to a question of determining the levels of 𝛼_𝑖 and 𝛾_𝑘. Therefore the (log) relative risk parameters Δ𝛼_𝑖 are identified as pointed out by Clayton and Schifflers [30]. In this model the cohort index is present and keeps the difference of the principal age and period indices. Therefore the representation of the predictor involves backward cumulated age differences as before but with a subtle change of sign, so that (54) reduces to 𝜇^ac_𝑖𝑗=𝜇_𝐼1−^𝐼−1_ℓ=𝑖Δ𝛼_ℓ+1+^𝑘_ℓ=2Δ𝛾_ℓ. (65) As a consequence the canonical parameter and the design reduce to 𝜇^ac_𝑖𝑗=𝑋^ac_𝑖𝑗𝜉^ac_𝑖𝑗, where 𝑋^ac_𝑖𝑗=1,−1_(1≥𝑖),…,−1_{(𝐼−1≥𝑖)},1_(𝑘≥2),…,1_{(𝑘≥𝐾)},𝜉^ac=𝜇_𝐼1,Δ𝛼₂,…,Δ𝛼_𝐼,Δ𝛾₂,…,Δ𝛾_𝐾^′. (66)(67) Miranda et al. [24, Theorem 4.2] establish an identification result similar to Theorem 6.

The age-cohort model can also be formulated as a hypothesis on the maximal invariant 𝜉 in the age-period-cohort model following Section 2.4.3. The period effects Δ²𝛽_𝑗 are set to zero through 𝐻^′𝜉=0, where 𝐻^′=(0,𝐼_𝐽−2,0). Applying this to the expression for 𝜉 in (55) gives 𝜉_𝐻=𝐻^′_⟂𝜉=𝜇_𝐼1,Δ𝛼_𝐼−Δ𝛾₂,Δ𝛾₂,Δ²𝛼₃,…,Δ²𝛼_𝐼,Δ²𝛾₃,…,Δ²𝛾_𝐾, (68) since in the absence of period effects; then 𝜇_𝐼1−𝜇_𝐼−1,1=Δ𝛼_𝐼−Δ𝛾₂ and 𝜇_𝐼2−𝜇_𝐼1=Δ𝛾₂. The double differences cumulate to first differences through ∑^𝐼_𝑖=3Δ²𝛼_𝑖=Δ𝛼_𝐼−Δ𝛼₂, so the above expression 𝜉_𝐻 is seen to be a linear transformation of 𝜉^ac in (67). In other words the age-cohort model arises from the age-period-cohort model by restricting the maximal invariant parameter.

5.4. Working with the Time Effect

There is a large literature seeking to identify the original time effects 𝛼_𝑖, 𝛽_𝑗, and 𝛾_𝑘 of the age-period-cohort model from the predictor. Here we look closer at some of those ad hoc identification proposals.

5.4.1. Ad Hoc Identification of Levels

For the age-period-cohort model it is popular to impose ad hoc identifications in two steps of the type discussed in Section 3.3. Here the first step is concerned with the level of the time effects and the second step is concerned with the linear trend. Examples are given in Sections 5.4.2 and 5.5.4.

A common first step ad hoc identification is to require that ^𝐼_𝑖=1𝛼_𝑖=^𝐽_𝑗=1𝛽_𝑗=^𝐾_𝑘=1𝛾_𝑘=0. (69) This ad hoc identification is specific to the chosen data range. For instance, the constraint ∑^𝐼_𝑖=1𝛼_𝑖=0 is not easily transferable to a different data set drawn from the same population but with a different set of age groups. This aspect would have to be kept in mind if a substantive motivation was to be found for this constraint. Other ad hoc identification schemes such as 𝛼_𝐼=𝛽_𝐽=𝛾_𝐾=0 have similar problems.

The constraint (69) is a special case of affine constraints of the form 𝐶^′𝜃_𝐶=𝜓 discussed in Section 3.3. The involved dimensions are 𝑞=2(𝐼+𝐽) and 𝑝=𝑞−4, while the number of constrains is 𝑞−𝑞_𝐶=3. The matrix 𝐶^′∈ℝ^{(𝑞−𝑞𝐶)×𝑞} is given by the top left {3×(𝑞−1)}-block of 𝐴^′_⟂ in (58) padded with a column of zeros, while 𝜓∈ℝ³ is given by 𝜓=0. Theorem 1 shows that 𝑚=𝐴^′_⟂𝐶∈ℝ^{(𝑞−𝑝)×(𝑞−𝑞𝑐)} has full rank. Indeed, 𝑚 and its orthogonal complement are given by 𝑚=𝐴^′_⟂𝐼̆̆𝑘̆𝐶=𝐼000𝐽000𝐾𝑖−𝐽𝑗𝐾,𝑚_⟂=−𝑖̆𝑗̆−𝑘̆1, (70) where, for instance, 𝑖̆=𝐼⁻¹∑^𝐼_𝑖=1𝑖=(𝐼+1)/2. Thus, the constrained group of equivalence classes (27) is 𝑔_𝐶∶𝛼_𝑖𝛽_𝑗𝛾_𝑘𝛿⟼𝛼_𝑖𝑖̆𝛽+𝑑𝑖−_𝑗𝑗̆𝛾−𝑑𝑗−_𝑘𝑘̆𝛿+𝑑𝑘−for𝜃∈Θ_𝐶. (71)

5.4.2. Ad Hoc Identification of Slopes: The “Intrinsic” Estimator

The “intrinsic” estimator is a popular estimator in the sociology literature; see Yang et al. [4] and see also O’Brien [31, 32] and Fu et al. [33] for a recent discussion of its merits. It has its roots in a suggestion by Kupper et al. [34], with an early critique given by Holford [35].

The “intrinsic” estimator is defined in two steps. In the first step, the levels are identified by the ad hoc constraint (69). Three of the 𝜃-coordinates are then dropped; that is 𝛼_𝐼, 𝛽_𝐽, and 𝛾_𝐾 are dropped. In a second step the linear trend is ad hoc identified using a Moore-Penrose inverse as in (23).

We can analyse these steps using the developed framework. The first step identifies the levels by the ad hoc constraint (69), which is a constraint of the form 𝐶^′𝜃=0 for the 𝐶 discussed in Section 5.4.1. This 𝜃 is defined on Θ_𝐶 which is a linear subspace with a dimension deficiency of 3. Introduce a selection matrix 𝑆_⟂∈ℝ^{𝑞×(𝑞−3)} that selects all coordinates of 𝜃 except 𝛼_𝐼, 𝛽_𝐽, and 𝛾_𝐾. Thus 𝑆_⟂ arises as a 𝑞-dimensional with 3 columns deleted corresponding to 𝛼_𝐼, 𝛽_𝐽, and 𝛾_𝐾. This is chosen so that (𝐶,𝑆_⟂) is invertible. Then 𝑆^′_⟂𝜃 is freely varying in that 𝑆^′_⟂Θ_𝐶=ℝ^𝑞−3. The skew projection identity 𝐼_𝑞=𝑆(𝐶^′𝑆)⁻¹𝐶^′+𝐶_⟂(𝑆^′_⟂𝐶_⟂)⁻¹𝑆^′_⟂ and the constraint 𝐶^′𝜃=0 then implies that 𝜃=𝐶_𝑆𝜗 where 𝐶_𝑆⟂=𝐶_⟂(𝑆^′_⟂𝐶_⟂)⁻¹ and 𝜗=𝑆^′_⟂𝜃∈ℝ^𝑞−3. Note that while 𝐶_𝑆⟂ depends on 𝑆_⟂ and 𝐶_⟂, it does not depend on the normalisation of 𝐶_⟂, since we can replace 𝐶_⟂ by 𝐶_⟂𝑚 for arbitrary invertible matrices 𝑚∈ℝ^{(𝑞−3)×(𝑞−3)}. This implies that 𝐶_𝑆 is a function of 𝑆_⟂ and 𝐶. The predictor 𝜇 is now parametrised by 𝜇=𝑋𝐴^′𝜃=𝑋𝐴^′_𝐶𝜗 with 𝐴^′_𝐶=𝐴^′𝐶_𝑆. This corresponds to equation 5 of Yang et al. [4] who use the notation 𝑋 and 𝑏 for 𝑋𝐴^′_𝐶 and 𝜗, respectively.

In the second step the linear trend is ad hoc identified through a time effect parameter of the form (23) with 𝐴, 𝜃 replaced by 𝐴_𝐶, 𝜗 so that 𝜃_ad.hoc=𝐶_⟂𝜗_ad.hoc where 𝜗_ad.hoc=𝐿_⟂(𝐴^′_𝐶𝐿_⟂)⁻¹𝜉+(𝐴_𝐶)_⟂{𝐿^′(𝐴_𝐶)_⟂}⁻¹𝜆 for some scalar 𝜆 and some matrix 𝐿_⟂∈ℝ^{(𝑞−3)×𝑝}.

The “intrinsic” estimator is ad hoc identified through the choices 𝜆=0 and 𝐿_⟂=𝐴_𝐶, while 𝐶 is chosen by (69). It therefore estimates an “intrinsic” parameter: 𝜃_intrinsic=𝐶_𝑆⟂𝐶^′_𝑆⟂𝐴𝐴^′𝐶_𝑆⟂𝐶^′_𝑆⟂𝐴⁻¹𝜉, (72) which depends on the choices of 𝑆^′_⟂, 𝐶, and 𝐴_⟂. However, since we can replace 𝐶_⟂ by 𝐶_⟂𝑚 for arbitrary invertible matrices 𝑚∈ℝ^{(𝑞−3)×(𝑞−3)} without changing 𝜃_intrinsic the expression 𝜃_intrinsic does not depend on the normalisation of 𝐶_⟂. The “intrinsic” parameter satisfies the following result.

Theorem 8. The “intrinsic” parameter is an injective mapping of the canonical parameter 𝜉∈ℝ^𝑝 into a 𝑝=𝑞−4 dimensional linear subspace Θ_intrinsic of Θ=ℝ^𝑞. The “intrinsic” time effect space is a 𝑝-dimensional linear subspace of ℝ^𝑞 of the form Θ_intrinsic=𝜃∈ℝ^𝑞∶𝜃=𝐶_𝑆⟂𝐶^′_𝑆⟂𝐴𝐴^′𝐶_𝑆⟂𝐶^′_𝑆⟂𝐴⁻¹𝜉𝑓𝑜𝑟𝜉∈ℝ^𝑝,=𝜃∈ℝ^𝑞∶𝐶^′𝜃=0,𝑤^′𝐶^′_⟂𝑆_⟂𝑆^′_⟂𝐶_⟂𝐶^′_⟂𝜃=0, (73) where 𝑤∈ℝ^𝑞−3 is uniquely defined up to a scale by 𝑤^′𝐶^′_⟂𝐴=0.

Theorem 8 implies that the “intrinsic” parameter should be interpreted as an object varying in the linear subspace Θ_intrinsic rather than in the unrestricted time effect space Θ=ℝ^𝑞. As outlined in Section 3.4 this has consequences for the interpretation of plots of the time effects, hypothesis testing, and forecasts. A consequence of this argument is that different choices of 𝐶, 𝑆_⟂, 𝐿, and 𝜆 would lead to other ad hoc identified parameters varying in other affine subspaces of Θ. In other words, the “intrinsic” estimator carries the cost of working with the somewhat complicated linear subspace Θ_intrinsic. This effort may be worthwhile if the particular choice of 𝐶, 𝑆_⟂, 𝐿, and 𝜆 can be made on substantive grounds.

5.4.3. Forecasting

Forecasting of future mortality rates involves an extrapolation of the time parameters. In Section 2.4.4 it was argued that ad hoc identification may introduce an undesirable arbitrariness in the forecast. When working exclusively with the canonical parameter 𝜉 this arbitrariness is avoided. It is, however, also possible to work with ad hoc identified time effects under specific circumstances that we characterise here for age-period arrays. This builds on the theory developed in Kuang et al. [25] for age-cohort data arrays.

In the context of an age-period data array ℐ_ap it is often of interest to forecast ℎ periods ahead. Suppose it is of interest to forecast the mortality at age 𝑖 in period 𝐽+ℎ, so that the cohort is 𝑘=𝐼+𝐽+ℎ−𝑖. This requires an extrapolation of the period effect. If the cohort index is sufficiently large, that is, 𝑘>𝐾, then the cohort effect needs to be extrapolated too. Thus, there are two forecast index arrays of interest: 𝒥_ap=,𝒦𝑖,𝑗∶𝑖=1,…,𝐼;𝑗=𝐽+1,…,𝐽+ℎ;𝑘≤𝐾_ap=𝑖,𝑗∶𝑖=1,…,𝐼;𝑗=𝐽+1,…,𝐽+ℎ;𝑘>𝐾. (74) Figure 2 illustrates these forecast index arrays.

Figure 2 

ℐap is the data array. 𝒥ap,1 is the forecast array where only period parameters need to be extrapolated. 𝒥2 is the forecast array where both period and cohort parameters need to be extrapolated. Cohorts are indicated by dashed lines.

Identification plays a role when extrapolating the estimates obtained on the data array ℐ_ap. The identification issues can be ignored if the investigator simply extrapolates Δ²𝛽_𝑗 and Δ²𝛾_𝑘. In the context of ad hoc identified time effects arbitrary linear trends are introduced in the model. The forecast of the predictor 𝜇_{𝑖,𝐽+ℎ} is invariant to these if and only if the chosen extrapolation method for 𝛽_𝑗, 𝛾_𝑘 preserves these linear trends so that they can cancel with the arbitrary linear trend in 𝛼_𝑖. The next result gives a precise formulation of this statement. It applies both to point forecasts and distribution forecasts.

Theorem 9. Consider the predictor 𝜇_𝑖𝑗 for 𝑖,𝑗∈ℐ_ap as given in (53). Suppose the time effects 𝛼_𝑖, 𝛽_𝑗, and 𝛾_𝑘 are ad hoc identified. Consider the class of ℎ periods-ahead forecasts over 𝒥_ap constructed as 𝜇_{𝑖,𝐽+ℎ}=𝛼_𝑖+𝛽_𝐽+ℎ+𝛾_{𝐼+𝐽+ℎ−𝑖}+𝛿, where 𝛽_𝐽+ℎ+𝛾_{𝐼+𝐽+ℎ−𝑖} is a function of the ad hoc identified estimate 𝜃. Let 𝑔 be the group (57). Invariance of the forecast 𝜇_{𝑖,𝐽+ℎ} with respect to the ad hoc identification is equivalent to either of the following:(i)the extrapolation method for period and cohort effects is linear trend-preserving: 𝛽_𝐽+ℎ𝜃+𝛾_{𝐼+𝐽+ℎ−𝑖}𝜃=𝛽_𝐽+ℎ𝑔𝜃(+𝛾−𝑏+𝑑𝐽+ℎ)_{𝐼+𝐽+ℎ−𝑖}𝑔𝜃−𝑐−𝑑(𝐼+𝐽+ℎ−𝑖)∀𝑏,𝑐,𝑑∈ℝ; (75)(ii)functions 𝑓_𝛽, 𝑓_𝛾 exist so that with 𝜉̂_𝛽,𝛾=(Δ²𝛽₃,…,Δ²𝛽_𝐽,Δ²𝛾₃,…,Δ²𝛾_𝐾)^′; then 𝛽_𝐽+ℎ𝜃+𝛾_{𝐼+𝐽+ℎ−𝑖}𝜃=𝛽_𝐽+ℎΔ𝛽_𝐽+𝑓_𝛽𝜉̂+𝛾_𝐾+(ℎ−𝑖+1)Δ𝛾_𝐾+𝑓_𝛾𝜉̂. (76)

To illustrate the use of Theorem 9 consider the extrapolation methods 𝛽_𝐽+ℎ=𝛽_𝐽 and Δ𝛽_𝐽+ℎ=Δ𝛽_𝐽. The first forecast is a random walk forecast and it is seen to violate (ii). The second forecast is a cumulated random walk and satisfies (ii). The reason is that 𝛽_𝐽+ℎ=𝛽_𝐽+∑^ℎ_ℓ=1Δ𝛽_𝐽+ℓ. Since Δ𝛽_𝐽+ℓ=Δ𝛽_𝐽, then 𝛽_𝐽+ℎ=𝛽_𝐽+ℎΔ𝛽_𝐽. Further examples of forecasts that are linear trend-preserving as well as some which are not are given Kuang et al. [25, Table 1].

Kuang, Nielsen, and Nielsen [10] apply this to reserving data organised in an age-cohort array ℐ_ac and discuss the issue of robustification of forecast with respect to structural breaks at the forecast origin. Miranda et al. [24] give an application to asbestos related mortality using an age-period array ℐ_ap.

5.4.4. Bayesian Ad Hoc Identification Using a Dynamic Prior

A Bayesian ad hoc identification using a dynamic prior does not solve the identification problem as discussed in Section 4 and the same care has to be exercised to avoid the problems outlined in Section 3.4. Berzuini and Clayton [6] suggest such an ad hoc identification approach. On page 831 they write “Identificability problems may be solved by imposing an arbitrary linear constraint on the log-linear trend components of age, period and cohort effects. Happily, such an arbitrary constraint has no effect on the predictions of the model.” The previous analysis suggests that this is far from innocent.

The Berzuini-Clayton suggestion is to ad hoc identify the model (53) through 𝛼_𝑖=𝛼₁+𝛼₂𝑖+^𝑖ℓ_ℓ=3ℎ=3Δ²𝛼_ℎ,𝛽_𝑗=𝛽₁+𝛽₂𝑗+^𝑗ℓ_ℓ=3ℎ=3Δ²𝛽_ℎ,𝛾_𝑘=𝛾₁+𝛾₂𝑘+^𝑘ℓ_ℓ=3ℎ=3Δ²𝛾_ℎ,𝛿=0. (77) A dynamic prior is chosen so that the double differences Δ²𝛼_𝑖, Δ²𝛽_𝑗, and Δ²𝛾_𝑘 are independent zero mean normal with variances 𝜙=(𝜎²_𝛼,𝜎²_𝛽,𝜎²_𝛾) that have 𝜒²-type prior. The purpose of this is in part to facilitate extrapolations Δ²𝛼_𝑖, Δ²𝛽_𝑗, and Δ²𝛾_𝑘 for 𝑖>𝐼, 𝑗>𝐽, and 𝑘>𝐾, which is done through further draws from normal distributions. The level/trend effects 𝜃_level=(𝛼₁,𝛼₂,𝛽₁,𝛽₂,𝛾₁,𝛾₂)^′ have independent uniform priors on some large intervals.

We will analyse the Berzuini-Clayton model as applied to an age-period data array ℐ_ap. Decompose the canonical parameter 𝜉 from (54) into two parts: the slope and level parameters, say 𝜉_𝜇=(𝜇_𝐼1,𝜇_𝐼1−𝜇_𝐼−1,1,𝜇_𝐼2−𝜇_𝐼1)^′, and the collection of double differences, say 𝜉_Δ. The assumed prior for 𝜉_Δ is a simple collection of independent normal distributions with variances 𝜙. The assumed prior for 𝜉_𝜇 is a linear combination of not only the independent uniform variables 𝜃_level, but also on 𝜉_Δ, since the age double differences Δ²𝛼_𝑖 are cumulated backwards in (54), but forwards in (77). Thus, the prior for 𝜉=(𝜉^′_𝜇,𝜉^′_Δ)^′ depends on the 𝜃_level construction.

We get a hyper-parameter 𝜆_hyper=(𝜆,𝜙), where 𝜆 is some three-dimensional ad hoc identified level/trend effect dependent on 𝜃_level,𝜉_Δ. We will argue that the ad hoc identified level/trend effect 𝜆 will wash out in the Berzuini-Clayton model. However, the level/trend parameter 𝜉_𝜇 is a function of the 𝜃_level construction that is tailored to the ad hoc identification. That construction remains in the analysis.

In the presentation of the posterior Berzuini and Clayton are careful only to consider the double differences 𝜉_Δ and stay clear of the ad hoc identified level/trend effect 𝜃_level. Theorem 2 yields the posterior 𝑝(𝜉∣𝑦)=𝑝(𝑦∣𝜉)𝑝(𝜉)/𝑝(𝑦). Thus, the marginal posterior for the double differences is 𝑝(𝜉_Δ∣𝑦)=∫𝑝(𝑦∣𝜉_Δ,𝜉_𝜇)𝑝(𝜉_Δ,𝜉_𝜇)𝑑𝜉_𝜇/𝑝(𝑦). This links 𝜉_Δ to 𝜉_𝜇 and in turn to the 𝜃_level construction.

The extrapolative method is based on double differences so it only depends on 𝜆_hyper through 𝜙 due to Theorem 9 and the subsequent discussion. Thus, the extrapolative method is of the form 𝑝(𝜇∣𝜉,𝜆_hyper,𝑦)=𝑝(𝜇∣𝜉,𝜙,𝑦). By construction it does not reduce to 𝑝(𝜇∣𝜉,𝑦) so that condition (39) for Theorem 3 is not satisfied. The distribution forecast is of the form 𝑝𝜇∣𝑦=𝑝𝜇∣𝜉,𝜙,𝑦𝑝𝜉∣𝑦𝑝𝜙∣𝜉𝑑𝜙𝑑𝜉, (78) which, apart from depending on the 𝜃_level construction, also depends on the conditional prior 𝑝(𝜙∣𝜉), which is not updated by the likelihood.

In summary, it appears that the Berzuini-Clayton analysis depends on the 𝜃_level construction as well as the conditional prior 𝑝(𝜙∣𝜉). The dependence on the 𝜃_level construction could of course be addressed by introducing priors directly on 𝜉_𝜇, which in turn would be updated by the likelihood. Since the conditional prior 𝑝(𝜙∣𝜉) cannot be updated by the likelihood that its sole justification rests on the substantial context.

5.4.5. A Functional Form Hypothesis

It is instructive to consider functional form restrictions on the time effects. Such hypotheses can be analysed using the results outlined in Section 3.4.2. As an example restrict the age effect to be quadratic in a similar way to Yang and Land [5] so that 𝛼_𝑖=𝜎₀+𝜎₁𝑖+𝜎₂𝑖²for𝑖=1,…,𝐼. (79)

This restriction on the time effect can be analysed by writing it on the form 𝑅^′𝜃=𝜌, see (28), and then applying Theorem A.3. Alternatively, in this particular case, we can show that the restriction actually only affects the ad hoc identified time effect through the canonical parameter, so a simpler analysis can be made.

A quadratic polynomial has constant second order derivative. Therefore the restriction (79) implies Δ²𝛼_𝑖=2𝜎₂for𝑖=3,…,𝐼. (80) This expression has one free parameter. Thus, it is useful to consider the third order difference: Δ³𝛼_𝑖=Δ²𝛼_𝑖−Δ²𝛼_𝑖−1=0for𝑖=4,…,𝐼. (81) This gives 𝐼−3 linear restrictions on the canonical parameter. The age time effect 𝛼_𝑖 then has three remaining parameters, say 𝛼₁, 𝛼₂, and 𝛼₃. These are freely varying since the parameters 𝜎₀, 𝜎₁, and 𝜎₂ are freely varying.

If the constraint is imposed directly on the canonical parameter, the restricted model is a regular exponential family with the advantages outlined in Section 2.4. However, if the analysis is done with the time effect the levels and trend will have to be ad hoc identified while bearing in mind the issues discussed above.

5.4.6. The “Hierarchical Age-Period Cohort Regression Model”

In some cases a random effects approach can be used to get an overview of the many parameters of the age-period model. When applied to the time effects this implies an ad hoc identification. An example is the “hierarchical age-period cohort regression model” by Yang and Land [5]. In that paper the age effect is given a quadratic structure, but that does not have to be the case. The model is then given by 𝛼_𝑖=𝜎₀+𝜎₁𝑖+𝜎₂𝑖²,𝛽=_𝑗^𝖣𝖭0,𝜎²_𝛽,𝛾=_𝑘^𝖣𝖭0,𝜎²_𝛾,𝛿=0. (82) Since random effects are only introduced for some of the time effects, the analysis of Section 4.3 has to modified in a similar way to the analysis in Section 5.4.4.

From (80) it is seen that the model restricts Δ²𝛼_𝑖=2𝜎₂. Thus, divide the canonical parameter 𝜉 into three elements: the slope and level parameters, say 𝜉_𝜇=(𝜇_𝐼1,𝜇_𝐼1−𝜇_𝐼−1,1,𝜇_𝐼2−𝜇_𝐼1)^′, the age-double differences 𝜉_𝛼=(Δ²𝛼₃,…,Δ²𝛼_𝐼), and the remaining double differences 𝜉_𝛽,𝛾. Here 𝜉_𝛼 is restricted by the hypothesis 2𝜎₂ and 𝜉_𝛽,𝛾 is linear function of the normal random effects, while 𝜉_𝜇 is a three-dimensional linear function of 𝜎₂ and of the six-dimensional object 𝜈=(𝜎₀,𝜎₁,𝛽₁,𝛽₂,𝛾₁,𝛾₂)^′. This leaves a three-dimensional ad hoc identified level/slope parameter 𝜆 which is also a function of 𝜈 but not entering the likelihood. Let 𝜓=(𝜎₀,𝜎₁,𝜎₂,𝜎²_𝛽,𝜎²_𝛾).

The random effects likelihood are constructed in three steps. First, we have the usual age-period-cohort likelihood 𝑝(𝑦∣𝜉). Secondly, the random effects distribution for 𝜉_𝜇, 𝜉_𝛽,𝛾, and 𝜆 is multivariate normal, while 𝜉_𝛼 is deterministic function of 𝜓. Thus, decompose the prior as 𝑝(𝜉_𝜇,𝜉_𝛽,𝛾,𝜆∣𝜓)=𝑝(𝜉_𝜇,𝜉_𝛽,𝛾∣𝜓)𝑝(𝜆∣𝜉_𝜇,𝜉_𝛽,𝛾,𝜓). Thirdly, following Section 4.3 the random effects likelihood will not depend on 𝑝(𝜆∣𝜉_𝜇,𝜉_𝛽,𝛾,𝜓) and it is given by 𝑝=𝑝𝑦∣𝜓𝑦∣𝜉_𝛼,𝜉_𝜇,𝜉_𝛽,𝛾𝑝𝜉_𝛼,𝜉_𝜇,𝜉_𝛽,𝛾∣𝜓𝑑𝜉_𝜇,𝜉_𝛽,𝛾. (83) The prior 𝑝(𝜆∣𝜉_𝜇,𝜉_𝛽,𝛾,𝜓) is not updated by the data. Plots and inferences based on the posterior 𝑝(𝜃∣𝜓,𝑦) will then suffer from the ad hoc identification issues outlined in Section 3.4.

5.5. A Two-Sample Age-Period-Cohort Model

When confronted with two samples for women and for men it may be of interest to apply the age-period-cohort model (43) to each of the samples and impose that some of the time effects are the same across samples. The models for samples 𝑟=1,2 are 𝜇_𝑖𝑗𝑟=𝛼_𝑖𝑟+𝛽_𝑗𝑟+𝛾_𝑘𝑟+𝛿_𝑟for𝑖,𝑗∈ℐ,𝑟=1,2. (84) The time effect 𝜃=(…,𝛼_𝑖𝑟,𝛽_𝑗𝑟,𝛾_𝑘𝑟,𝛿_𝑟,…)^′ now varies in Θ=ℝ^𝑞 where 𝑞=4(𝐼+𝐽).

5.5.1. Analysis of the Unrestricted Two-Sample Model

The unrestricted two-sample model is simply analysed as two copies of the one sample model of Section 5.1. The time effects of each copy are only defined up to linear trends. The group of transformations characterizing the identification problem combines two copies of the one sample group (44). The maximal invariant parameter is 𝜉=(𝜉^′₁,𝜉^′₂)^′∈ℝ^𝑝 where 𝑝=4(𝐼+𝐽−2) and each of 𝜉_𝑟 are of the form (45). The benefits of Section 2 hold when working with that parameter.

5.5.2. Bayesian Ad Hoc Identification Using a Dynamic Model

An application of the unrestricted two-sample model can be found in Cairns et al. [36]. The two samples are the population of England and Wales and the subpopulation of assured lives, so the substantive question is whether there is a selection effect for the assured lives. A Bayesian model with dynamic prior is used. It shares some features with the Berzuini and Clayton [6] model discussed in Section 5.4.4 although the details of the ad hoc identification of the levels and slopes are slightly different. When it comes to forecasting the extrapolative method appears to depend on the ad hoc identified parameter as well as the hyperparameters. This complicates the analysis of the forecast relatively the discussion in Section 5.4.4.

5.5.3. The Hypothesis of Common Period Parameters

The two-sample model allows the possibility for adding cross-sample restrictions on the parameters. As an example we consider the hypothesis of common period parameters.

Working with the canonical parameter the hypothesis is Δ²𝛽_𝑖1=Δ²𝛽_𝑖2for𝑗=3,…,𝐽. (85) This is a simple linear restriction as that discussed in Section 2.4.3. It is readily seen that the degrees of freedom of the hypothesis are 𝑝−𝑝_𝐻=𝐽−2 so the dimension of the restricted model is 𝑝_𝐻=4𝐼+3𝐽−6. The canonical parameter under the hypothesis is then 𝜉_𝐻=…,𝜇_11𝑟,𝜇_21𝑟−𝜇_11𝑟,𝜇_12𝑟−𝜇_11𝑟,Δ²𝛼_𝑖𝑟,Δ²𝛽_𝑗,Δ²𝛾_𝑘𝑟,…^′. (86)

The same result arises when writing the hypothesis in terms of time effects so that 𝛽_𝑗1=𝛽_𝑗2for𝑗=1,…,𝐽. (87) Such hypotheses on the time effect were discussed in Section 3.4.2. It can be analysed using the general result in Theorem A.3. However, we will take the simpler route of arguing that this only restricts the canonical parameter given a hypothesis of the type (85). The argument relies on noting that analysing the restriction for the predictors 𝜇_𝑖𝑗1 and 𝜇_𝑖𝑗2 is equivalent to analysing the restriction for the predictors 𝜇_𝑖𝑗1 and 𝜇_𝑖𝑗2−𝜇_𝑖𝑗1, where the cross-sample differenced predictor is of the form 𝜇_𝑖𝑗2−𝜇_𝑖𝑗1=𝛼_𝑖2−𝛼_𝑖1+𝛽_𝑗2−𝛽_𝑗1+𝛾_𝑘2−𝛾_𝑘1+𝛿₂−𝛿₁. (88) Now, the restricted model for the cross-sample differenced predictor 𝜇_𝑖𝑗2−𝜇_𝑖𝑗1 is an age-cohort model: 𝜇_𝑖𝑗2−𝜇_𝑖𝑗1=𝛼_𝑖2−𝛼_𝑖1+𝛾_𝑘2−𝛾_𝑘1+𝛿₂−𝛿₁. (89) Following the analysis of Section 5.3 the (87) therefore implies the 𝐽−2 linear restrictions given by (85). At the same time the predictor for the first sample 𝜇_𝑖𝑗1 is left unrestricted by (87). In summary, the restrictions (85) and (87) are equivalent.

The restriction has an interesting implication for the interpretation of the involved double differences. For the unrestricted model it was found that only plain double differences, such as Δ²𝛼_𝑗𝑟, are identified. Under the restriction the cross-sample differenced predictor is of age-cohort form (89) so also the cross-double differences Δ(𝛼_𝑖2−𝛼_𝑖1) and Δ(𝛾_𝑘2−𝛾_𝑘1) are identified.

5.5.4. Step-Wise Ad Hoc Identification under the Hypothesis

The analysis of Riebler and Held [7] finds that the difference 𝛼_𝑖2−𝛼_𝑖1 is identified under the hypothesis (85). This is not consistent with the above analysis showing that the cross-sample differenced predictor is an age-cohort model under the hypothesis, for which levels such as 𝛼_𝑖2−𝛼_𝑖1 are identified.

The apparent difference comes about because Riebler and Held follow a step-wise identification approach along the lines of Sections 3.3 and 5.4.1. In a first step the time effects 𝛼_𝑖𝑟, 𝛽_𝑗𝑟, and 𝛾_𝑘𝑟 are constrained to have zero-sums as in (69). In a second step the slopes are ad hoc identified using a Bayesian approach similar to that of Berzuini and Clayton [6]; see Sections 4 and 5.4.4 for a discussion of the consequences.

The identification in the first step implies that 𝛼_𝑖2−𝛼_𝑖1 has a zero sum. Under the hypothesis (85) this is exactly what is needed to ad hoc identify the levels in the age-cohort model (89). In other words a different level identification in the first step leads to a different level for the difference 𝛼_𝑖2−𝛼_𝑖1.

6. Models with Nonlinear Parametrisations

Some additional issues arise when looking at models with nonlinear parametrisations. A prominent example is the mortality model proposed by Lee and Carter [3] and which is the current benchmark in mortality studies done by government agencies and pension funds. For this model the time effect space Θ has a nondifferentiability which can actually be avoided by working directly with the parameter space 𝑀.

We analyze the Lee-Carter model in Section 6.1. In Section 6.2 we turn to a two-sample problem where some additional difficulties can arise when forecasting.

6.1. The Lee-Carter Model

The mortality model proposed by Lee and Carter [3] has predictor of the form 𝜇_𝑖𝑗=𝛼_𝑖+𝛽_𝑖𝜅_𝑗for𝑖,𝑗∈ℐ_ap. (90) The time effects 𝜃=(𝛼₁,…,𝛼_𝐼,𝛽₁,…,𝛽_𝐼,𝜅₁,…,𝜅_𝐽) vary in Θ=ℝ^2𝐼+𝐽.

Lee and Carter pointed towards two identification issues of the model. If 𝛼, 𝛽, and 𝜅 are one solution to (90), then 𝛼−𝛽𝑐,𝛽,𝜅+𝑐 is also a solution for any scalar 𝑐, just as 𝛼, 𝛽/𝑑, and 𝜅𝑑 are a solution for any 𝑑≠0. Consequently, they proposed the ad hoc identification: _𝑖𝛽_𝑖=1,_𝑗𝜅_𝑗=0. (91) This is, however, not the full story about the identification issues. To get at this we follow the outline from the linear parametrised models and start by finding the parameter space for the predictor 𝜇.

6.1.1. The Parameter Space

We start by finding the predictor space 𝑀. Write the model in matrix form. Let 𝜇 denote the 𝐼×𝐽-matrix of 𝜇_𝑖𝑗. Then 𝜇=𝛼𝜄^′+𝛽𝜅^′, (92) where 𝛼, 𝛽, and 𝜅 are vectors concatenating 𝛼_𝑖, 𝛽_𝑖, and 𝜅_𝑗 and where 𝜄=(1,…,1)^′∈ℝ^𝐽. Postmultiply by the projection identity 𝐼_𝐽=𝜄𝜄^′+𝜄_⟂𝜄^′_⟂ to get 𝜇=𝛼𝜄^′+𝛽𝜅^′𝜄𝜄^′+𝜄_⟂𝜄^′_⟂=𝛼+𝛽𝜅^′𝜄𝜄^′+𝛽𝜅^′𝜄_⟂𝜄^′_⟂, (93) where the orthogonal complement 𝜄_⟂ can be chosen so that 𝜄^′_⟂𝜅=(Δ𝜅₂,…,Δ𝜅_𝐽)^′ but could also be chosen otherwise. Equation (93) shows that the model is composed of two matrices with rank one. Thus, the parameter space is given by 𝑀=𝜇∈ℝ^𝐼×𝐽∶𝜇=𝛾𝜄^′+𝛿𝜄^′_⟂for𝛾,𝛿∈ℝ^𝐼×ℝ^{𝐼×(𝐽−1)}sorank(𝛿)≤1. (94) Note that 𝑀 does not depend on the normalisation of 𝜄_⟂ since 𝛿 is freely varying. The space 𝑀 is a manifold since the space of matrices 𝛿 with an upper bound to the rank is a manifold as opposed to the space where 𝛿 has rank of unity. This space can be parametrised parsimoniously by 𝜉=𝛾,𝛿where𝛾=𝛼+𝛽𝜅^′𝜄,𝛿=𝛽𝜅^′𝜄_⟂, (95) varying in the manifold Ξ=𝛾,𝛿∈ℝ^𝐼×ℝ^{𝐼×(𝐽−1)}∶rank(𝛿)≤1. (96) The 𝜉 is the candidate for the maximal invariant describing the equivalence classes of the mapping from the time effect 𝜃 to the predictor 𝜇.

The next step is to analyse the time effect space Θ. It is convenient to decompose 𝑀 into two disjoint sets depending on the rank of 𝛿. These sets are 𝑀₁=𝜇∈ℝ^𝐼×𝐽∶𝜇=𝛾𝜄^′+𝛿𝜄^′_⟂for𝛾,𝛿∈ℝ^𝐼×ℝ^{𝐼×(𝐽−1)}(𝛿),𝑀sorank=1₀=𝜇∈ℝ^𝐼×𝐽∶𝜇=𝛾𝜄^′for𝛾∈ℝ^𝐼. (97) Correspondingly, the time effect space Θ can be decomposed into two disjoint sets: Θ₁=𝜃∈Θ∶∃𝑖,𝑗so𝛽_𝑖𝜅_𝑗≠𝛽_𝑖𝜅_𝐽,Θ₀=𝜃∈Θ∶∀𝑖,𝑗so𝛽_𝑖𝜅_𝑗=𝛽_𝑖𝜅_𝐽. (98) Note that 𝛿=0 if and only if 𝜃∈Θ₀. Consider first the time effect space Θ₁, which is implicitly what Lee and Carter had in mind. The mapping 𝜃↦𝜇 on Θ₁ to 𝑀 is invariant to the group of transformations: 𝑔₁∶𝛼_𝑖𝛽_𝑖𝜅_𝑗⟼𝛼_𝑖+𝛽_𝑖𝑐𝛽_𝑖𝑑𝜅_𝑗𝑑,−𝑐 (99) acting on Θ₁ for all 𝑐∈ℝ and all 𝑑≠0. The parameter 𝜉=(𝛾^′,𝛿^′)^′ is invariant under 𝑔₁ acting on Θ₁. Now, consider the space Θ₀ with deficient rank. Then 𝛼_𝑖,𝛽_𝑖, and 𝜅_𝑗 map into 𝛼_𝑖+𝜑_𝑖 where 𝜑_𝑖=𝛽_𝑖𝜅_𝑗 is constant in 𝑗, so that 𝛿=𝛽𝜅^′𝜄_⟂=0. This mapping is invariant to the group of transformations: 𝑔₀∶𝛼_𝑖𝛽_𝑖𝜅_𝑗⟼𝛼_𝑖+𝑎_𝑖𝛽_𝑖𝜅_𝑗−𝑎_𝑖, (100) acting on Θ₀ for all (𝑎₁,…,𝑎_𝐼)^′∈ℝ^𝐼.

Theorem 10. Let 𝜇∈𝑀. The parameter 𝜉∈Ξ of (95) satisfies the following:(i)𝜉 is a function of 𝜃∈Θ which is invariant to the groups 𝑔₀, 𝑔₁ in (99) and (100);(ii)𝜇 is a function of 𝜉;(iii)the parametrisation of 𝜇 by 𝜉 is exactly identified in the sense that 𝜉^†≠𝜉^‡⇒𝜇(𝜉^†)≠𝜇(𝜉^‡).

Theorem 10 shows that 𝜉 varies freely on the space Ξ and it gives a unique parametrisation of 𝜇. As a function of 𝜃 it is invariant to 𝑔₀,𝑔₁; hence it is a maximal invariant.

It is interesting to compare the properties of the spaces 𝑀, Ξ, and Θ. The spaces 𝑀 and Ξ are spaces of matrices with deficient rank. These are smooth spaces, but they are not vector spaces since the sum of matrices with rank one may have rank larger than one. In contrast Θ is a vector space. The mapping from Θ to 𝑀 will inevitably be nondifferentiable. This nondifferentiability is avoided by working directly with 𝑀. Likewise, in a Bayesian setting it would seem more difficult to introduce a meaningful prior of Θ with its nondifferentiability than on 𝑀.

6.1.2. Maximum Likelihood Estimation

The maximum likelihood estimator for 𝜉 can be derived analytically in the normal case.

Consider a situation where the data array is of age-period form so 𝑌_𝑖𝑗 for (𝑖,𝑗)∈ℐ_ap. Suppose 𝑌_𝑖𝑗 are independent normal with mean 𝜇_𝑖𝑗 and variance 𝜎². Organise the data in a matrix 𝑌. Then the log likelihood is of the form ℓ𝜇,𝜎²;𝑌=−𝐼𝐽2log2𝜋𝜎²−12𝜎²tr𝑌−𝜇𝑌−𝜇^′. (101) The maximum likelihood estimator is of the following form. Subsequently, this is related to the estimator suggested by Lee and Carter.

Theorem 11. For a normal age-period array parametrised by (94) the maximum likelihood estimators are 𝛾=𝑌𝜄𝜄^′𝜄⁻¹,𝛿=svd₁𝑌𝜄_⟂𝜄^′_⟂𝜄_⟂⁻¹𝜄^′_⟂𝜄_⟂𝜄^′_⟂𝜄_⟂⁻¹, (102) where svd₁(⋅) is the singular value decomposition truncated to one factor.

Thus, 𝛾 is estimated by the row-averages of the data matrix, while 𝛿 is estimated by the singular value decomposition of the row-wise demeaned data matrix.

6.1.3. Estimation of Ad Hoc Identified Time Effects

The ad hoc identification (91) gives a time effect 𝜃_𝜆 varying in a 2𝐼+𝐽−2 dimensional affine subspace of Θ=ℝ^2𝐼+𝐽. The ad hoc identified 𝜃_𝜆 can now be expressed in terms of the maximal invariant parameter 𝜉 using (95). In the case where 𝛿≠0 then it has singular value decomposition 𝛿=𝛿_𝐿𝛿_𝑆𝛿^′_𝑅 for two vectors 𝛿_𝐿∈ℝ^𝐼 and 𝛿_𝑅∈ℝ^𝐽−1 so 𝛿^′_𝐿𝛿_𝐿=1 and 𝛿_𝑅𝛿^′_𝑅=1, while 𝛿_𝑆>0 is a positive scale. The ad hoc identification of Lee and Carter then gives 𝛼_𝜆=𝛾,𝛽_𝜆=𝛿_𝐿𝜄^′𝛿_𝐿⁻¹,𝜅_𝜆=𝜄_⟂𝛿_𝑅𝜄^′𝛿_𝐿𝛿_𝑆. (103) Inserting the maximum likelihood estimators from Theorem 11 yields the estimators proposed by Lee and Carter. However, the disentangling of the singular values and singular vectors of 𝛿 is done by the ad hoc identification 𝛽^′𝜄=1 and 𝜅^′𝜄=0. These estimators are therefore specific to the considered data array and data set in parallel with the discussion in Sections 3.2 and 5.4.1.

6.1.4. Consequences of the Possible Rank Deficiency

The parameter space 𝑀 was split into spaces 𝑀₁ and 𝑀₀ depending on the rank of 𝛿. The space 𝑀₀ is a Lebesgue null set relative to 𝑀. Broadly speaking, there are two consequences of the possible rank deficiency. The first consequence is an estimation problem arising in the vicinity of 𝑀₀. The second consequence is that the usual normal asymptotic distribution theory does not apply in the vicinity of 𝑀₁. Whether this becomes a problem in practice depends on the data. One solution is to ensure that the time effect really is present when using the Lee-Carter model.

Investigate whether the time effects are present amounts to estimating the rank of 𝛿. For a given data set two Lee-Carter models can be estimated. The first model with predictor space 𝑀 is the unrestricted model in which rank(𝛿)≤1. The second model has predictor space 𝑀₀ so 𝛿=0. Twice the difference of the likelihood values gives a likelihood ratio test statistic which is asymptotically 𝜒². If the smaller model, 𝑀₀, is accepted this is used in subsequent analysis. However, if the smaller model 𝑀₀ is rejected then it is likely that the predictor is not located in the vicinity of 𝑀₀ and it is then safe to work with the predictor space 𝑀₁.

The consistency of this step-wise procedure is discussed in a cointegration context by Johansen [37, Section 12]. Even when this procedure points towards working with the parameter space 𝑀₁ the rank deficiency may still affect inference under 𝑀₁. Analysis of simple canonical correlation models suggests that inference under 𝑀₁ will be nearly similar if the distance to 𝑀₀ is sufficiently large. A problem is that the distribution for the test statistic will have poor finite sample properties when the parameter value is close to 𝑀₀. A simple way to get around this problem is to test for 𝑀₀ using a test with lower level than the conventional level. A more complicated way to address this is to employ a finite sample correction when seeking to test for 𝑀₀. See Nielsen [38, 39] for further discussion of this issue in the context of simple canonical correlation models.

The rank deficiency issue is typically not encountered in a standard Lee-Carter analysis. The reason is that the analysis is typically applied to data where there is a marked improvement in mortality rates over time. A Lee-Carter analysis could however run into trouble if it were applied to data without a strong calendar effect. The issue becomes more pertinent when extending Lee-Carter model with a cohort component such as 𝜇_𝑖𝑗=𝛽_𝑖⁽¹⁾+𝛽_𝑖⁽²⁾𝜅_𝑗⁽²⁾+𝛽_𝑖⁽³⁾𝛾⁽³⁾_{𝑗−𝑖+𝐼}; (104) see Renshaw and Haberman [40]. If the cohort effect is modest the latter matrix is nearly rank deficient and the likelihood will be nearly flat in certain directions. This is presumably the reason for the estimation problem noted by Cairns et al. [41].

6.1.5. Forecasting

The purpose of Lee-Carter model is usually to forecast future mortality. This issue is considered for the model with parameter space 𝑀₁. The standard approach is to extrapolate 𝜅, ad hoc identified through, for instance, 𝜅^′𝜄=0. The ℎ-step ahead extrapolation of 𝜅_𝐽+ℎ based on some forecast methods is denoted by 𝜅_𝐽+ℎ(𝜅). Combined with the estimates 𝛼_𝑖, 𝛽_𝑖 this gives the overall forecast 𝜇_{𝑖,𝐽+ℎ}𝜃=𝛼_𝑖+𝛽_𝑖𝜅_𝐽+ℎ(𝜅). (105) The identification question is then for which extrapolation methods this equals 𝜇_{𝑖,𝐽+ℎ}𝑔₁𝜃=𝛼_𝑖+𝛽_𝑖𝑐+𝛽_𝑖𝑑𝜅_𝐽+ℎ{(𝜅−𝑐)𝑑}. (106) The condition for avoiding adverse impact of the ad hoc identification is as follows.

Theorem 12. Let 𝜇∈𝑀₂. The forecast 𝜇_{𝑖,𝐽+ℎ} in (105) is invariant to ad hoc identification if and only if the extrapolation method for the period effect is location-scale preserving: 𝜅_𝐽+ℎ{(𝜅−𝑐𝜄)𝑑}=𝜅_𝐽+ℎ(𝜅)−𝑐𝑑∀𝑐∈ℝ,∀𝑑≠0. (107)

The default forecast method in the literature is a random walk with a drift, which was the preferred forecast of Lee and Carter [3]. This is given by 𝜅_𝐽+ℎ=𝜅_𝐽+ℎ−1+𝜈_𝑐+𝜀_ℎ, (108) with estimates 𝜈_𝑐=(𝐽−1)⁻¹∑^𝐽_𝑗=2(𝜅_𝑗−𝜅_𝑗−1) and normal errors 𝜀_ℎ with mean zero and estimated variance 𝜎²(𝜅)=(𝐽−2)⁻¹∑^𝐽_𝑗=2(𝜅_𝑗−𝜅_𝑗−1−𝜈_𝑐)². This extrapolation method is location-scale preserving as required in Theorem 12. It is even linear trend preserving. Other valid forecasts are a random walk without intercept as given by the equation 𝜅_𝐽+ℎ=𝜅_𝐽+ℎ−1+𝜀_ℎ, or an autoregression given by 𝜅_𝐽+ℎ=𝜌𝜅_𝐽+ℎ−1+𝜈_𝑐+𝜀_ℎ.

An alternative approach to forecasting would consider the predictor of the model for a particular age ground, say 𝑖. That predictor is 𝜇_𝑖=𝑒^′_𝑖(𝛾𝜄^′+𝛿𝜄^′_⟂), where 𝑒_𝑖 is the 𝑖th unit vector. From this we can generate forecasts 𝜇_{𝑖,𝐽+ℎ} using any time series method. The resulting forecast will in general depend on 𝜅 as well as 𝛼_𝑖, 𝛽_𝑖 and it is therefore more general than the forecasts discussed in Theorem 12, which only depends on 𝜅. The forecast for another age group, say 𝑖^†, should be the same up to a linear transformation dictated by the Lee-Carter structure. Thus, the ℎ-step ahead forecasts for the entire array are 𝜇_𝐽+ℎ=𝛿𝜄^′_⟂𝑒_̆𝚤𝑒^′_𝑖𝛿𝜄^′_⟂𝑒_̆𝚤𝜇_{𝑖,𝐽+ℎ}−𝑒^′_𝑖𝛾𝜄^′+𝛾𝜄^′, (109) for an index 𝚤̆ is chosen so that 𝑒^′_𝑖𝛿𝜄^′_⟂𝑒_̆𝚤≠0.

6.1.6. Bayesian Ad Hoc Identification Using a Dynamic Model

A Bayesian model with dynamic specification of the prior has been suggested by Pedroza [42]. Dynamic priors are presented for the time effects 𝜃=(𝜉,𝜆) involving a hyper parameter 𝜙. The ad hoc identification (91) is imposed so that analysis is made for an ad hoc identified time effect 𝜃_𝜆.

Pedroza presents posteriors for 𝜃_𝜆. When evaluating this posterior one should bear in mind that the conditional prior 𝑝(𝜆∣𝜉) is not updated by the data; see Theorem 2. The presented extrapolative method does not depend on 𝜆. Even so, the forecast will depend on conditional prior 𝑝(𝜙∣𝜉) which is not updated by the data; see Theorem 3.

6.2. The Two-Sample Lee-Carter Model

We now turn to applications of the Lee-Carter model in two-sample problems. Suppose two samples are for women and men. One approach would be to fit separate Lee-Carter models to the two datasets. These Lee-Carter models are of the form 𝜇_𝑖𝑗𝑟=𝛼_𝑖𝑟+𝛽_𝑖𝑟𝜅_𝑗𝑟for𝑖,𝑗∈ℐ_ap,𝑟=1,2. (110) The objective is now to extrapolate the period effects 𝜅_𝑗𝑟. Extrapolating the two models separately using separate random walks is often seen to be volatile, so methods that seek to combine information from both estimated series 𝜅_𝑗𝑟 are sought after. The next result describes the invariance problem in forecasting.

Theorem 13. Let 𝜇_𝑟∈𝑀_2𝑟 for 𝑟=1,2. The forecast 𝜇_{𝑖,𝐽+ℎ,1} for sample 𝑟=1 is invariant to ad hoc identification if the extrapolation method 𝜅_𝐽+ℎ,1 preserves location/scale for sample 1, but is invariant to location and scale for sample 2. That is for all 𝑐₁,𝑐₂∈ℝ and all 𝑑₁,𝑑₂≠0; then 𝜅_𝐽+ℎ,1𝑑₁𝜅₁−𝑐₁,𝑑₂𝜅₂−𝑐₂=𝑑₁𝜅_𝐽+ℎ,1𝜅₁,𝜅₂−𝑐₁. (111)

For one sample the standard forecasting technique appears to be the random walk with a drift as in (108). For the two-sample problem a suggestion could be that women and men should share a common random walk with a drift but deviate from this by a stationary process. In econometrics this idea is referred to as cointegration as proposed by Engle and Granger [43]; see also Johansen [37] for a likelihood based vector autoregressive approach. It is tempting to require that the calendar effects should cointegrate with coefficients of unity, so 𝜅_𝑗1−𝜅_𝑗2 should be stationary. However, that apparently intuitive choice violates Theorem 13 because the locations and scales of 𝜅_𝑗𝑟 are different and arbitrary.

There are two fixes to this problem. The first solution is to work directly with the mortality predictors 𝜇_𝑖𝑗𝑟 for an arbitrary age group 𝑖 as outlined for the one-sample case in connection with (109). Since no identification is involved it is permitted to impose that 𝜇_𝑖𝑗1 and 𝜇_𝑖𝑗2 cointegrate with coefficients of unity. The forecast for age group 𝑖 is then carried over to other age groups. The second solution is to work with the estimated series 𝜅_𝑗𝑟 but estimate the cointegrating coefficients from the data. In other words, the cointegrating relation 𝜅_𝑗1−𝜑𝜅_𝑗2−𝜓 should be zero mean, stationary, with coefficients 𝜑, 𝜓 estimated from the data. This can, for instance, be done by Johansen’s approach for a bivariate vector autoregression; see Hendry and Nielsen [44, Section 17].

7. Conclusion

Ad hoc identification is intimately linked to interpretation, inference, numerical analysis, and forecasting. The ad hoc identification will often introduce an arbitrary element in the statistical analysis, whether it is based on frequentist or Bayesian methods. This arbitrary element is entirely avoidable and is in our view best avoided unless there is a clear substantial motivation for ad hoc identification. For decades there has been a debate over how it is best to ad hoc identify mortality models. Our proposal is to bypass this discussion by analysing the surjective mapping between the unidentified time effect parameter and the predictor of the model and then deduce a maximal invariant parametrisation. In our experience there are typically two substantial benefits. First, it simplifies estimation and other statistical computations which is what we have focused on here. Secondly and perhaps more importantly, it helps to focus the substantial question that gives rise to the analysis in the first place.

The issue of dealing with two time scales also occurs in other statistical models, such as the Cox regression model; see Cabrera et al. [45] for a recent application. In future research it would be interesting to consider whether the analysis presented here has any bearing on that problem.

Appendix

A. Proofs

A.1. Some Linear Algebra Results

Lemma A.1. Consider 𝐴∈ℝ^𝑞×𝑝 and 𝐶∈ℝ^{𝑞×(𝑞−𝑞𝐶)} so 𝑝,𝑞_𝐶≤𝑞. Suppose they have full column rank. Then the following statements are equivalent for some 𝑝_𝐶≤𝑝∶(i)(𝐴,𝐶)∈ℝ^{𝑞×(𝑝+𝑞−𝑞𝐶)} has rank 𝑝_𝐶+𝑞−𝑞_𝐶;(ii)𝐴^′_⟂𝐶∈ℝ^{(𝑞−𝑝)×(𝑞−𝑞𝐶)} has rank 𝑝_𝐶+𝑞−𝑞_𝐶−𝑝;(iii)𝐶^′_⟂𝐴∈ℝ^{𝑞𝐶×𝑝} has rank 𝑝_𝐶.

Proof of Lemma A.1. (i)⇔(ii) Premultiply the matrix (𝐴,𝐶) with the invertible matrix (𝐴,𝐴_⟂)^′ to get the identity 𝐴^′𝐴^′_⟂()=𝐼𝐴,𝐶^𝑝𝐴^′𝐶0𝐴^′_⟂𝐶=𝐼_𝑝00𝐴^′_⟂𝐶𝑀𝐼where𝑀=^𝑝𝐴^′𝐶0𝐼_{𝑝+𝑞−𝑞𝐶}. (A.1) Since the first matrix (𝐴,𝐴_⟂)^′ and the last matrix 𝑀 have full rank then rank(𝐴,𝐶)=rank𝐼_𝑝+rank𝐴^′_⟂𝐶=𝑝+rank𝐴^′_⟂𝐶. (A.2) (i)⇔(iii) Swap the roles of 𝐴, 𝐶 so rank(𝐴,𝐶)=𝑞−𝑞_𝐶+rank(𝐶^′_⟂𝐴).

Lemma A.2. Consider 𝐴∈ℝ^𝑞×𝑝 and 𝐶∈ℝ^{𝑞×(𝑞−𝑞𝐶)} so 𝑝≤𝑞_𝐶≤𝑞. Suppose (𝐴,𝐶)∈ℝ^{𝑞×(𝑝+𝑞−𝑞𝐶)} has full column rank 𝑝+𝑞−𝑞_𝐶. Then 𝑚=𝐴^′_⟂𝐶∈ℝ^{(𝑞−𝑝)×(𝑞−𝑞𝐶)} has full column rank and (𝐴,𝐶) has orthogonal complement given by (𝐴,𝐶)_⟂=𝐴_⟂𝑚_⟂ where 𝑚_⟂∈ℝ^{(𝑞−𝑝)×(𝑞𝐶−𝑝)}.

Proof of Lemma A.2. Since (𝐴,𝐶) has full column rank Lemma A.1 (i, ii) implies that 𝑚 has full column rank. Since 𝑚^′_⟂𝐴^′_⟂(𝐴,𝐶)=(𝑚^′_⟂𝐴^′_⟂𝐴,𝑚^′_⟂𝑚)=0 we argue that (𝐴,𝐶,𝐴_⟂𝑚_⟂) is invertible. Premultiply by the invertible matrix (𝐴,𝐴_⟂)^′ to get 𝐴^′𝐴^′_⟂𝐴,𝐶,𝐴_⟂𝑚_⟂=𝐼_𝑝𝐴^′𝐶00𝑚𝐴^′_⟂𝐴_⟂𝑚_⟂. (A.3) The matrix (𝑚,𝐴^′_⟂𝐴_⟂𝑚_⟂) has full rank. Indeed, its inverse is 𝑚,𝐴^′_⟂𝐴_⟂𝑚_⟂⁻¹=𝐴^′_⟂𝐴_⟂⁻¹𝑚𝑚^′𝐴^′_⟂𝐴_⟂⁻¹𝑚⁻¹,𝑚_⟂𝑚^′_⟂𝐴^′_⟂𝐴_⟂𝑚_⟂^−1′. (A.4) Thus, the block triangular matrix (A.3) and, hence, (𝐴,𝐶,𝐴_⟂𝑚_⟂), have full rank.

A.2. Proofs of Main Theorems

Proof of Theorem 1. Since (𝐴,𝐶) has full column rank then Lemma A.2 shows it has orthogonal complement 𝐴_⟂𝑚_⟂ so that (𝐴,𝐶,𝐴_⟂𝑚_⟂) is invertible. The orthogonal projection identity shows 𝜃=𝐴𝐴,𝐶^′𝐶^′𝜃+𝐴_⟂𝑚_⟂𝐴_⟂𝑚_⟂^′𝜃=𝜉𝜓𝐴,𝐶+𝐴_⟂𝑚_⟂𝜁, (A.5) by the constraint 𝐶^′𝜃=𝜓 and the definitions 𝐴^′𝜃=𝜉 and (𝐴_⟂𝑚_⟂)^′𝜃=𝜁. This defines the constrained time effect space Θ_𝐶. Consider now the mapping 𝜃↦𝜇=𝑋𝐴^′𝜃. Premultiply the above expression for 𝜃 by 𝐴^′=(𝐼_𝑝,0)(𝐴,𝐶)^′ to get 𝐴^′𝜃=𝜉 so 𝜇=𝑋𝜉. Thus, since 𝜓 is fixed, the equivalence classes in Θ_𝐶 are given by 𝑔_𝐶∶𝜃↦𝐴_⟂𝑚_⟂𝜁, with 𝜉=𝐴^′𝜃 as a maximal invariant.

Proof of Theorem 2. (i) With the likelihood (32) so 𝑝(𝑦∣𝜉,𝜆)=𝑝(𝑦∣𝜉) then 𝑝𝑦=𝑝𝑝()𝑝()=𝑝𝑝(𝜉)𝑦∣𝜉,𝜆𝜉,𝜆𝑑𝜆𝑑𝜉=𝑝𝑦∣𝜉𝜉,𝜆𝑑𝜆𝑑𝜉𝑦∣𝜉𝑑𝜉. (A.6)
(ii) By Bayes formula and the likelihood (32) then 𝑝=𝑝𝑝(𝜉)𝜉∣𝑦𝑦∣𝜉𝑝𝑦,𝑝=𝑝𝑝()𝜆∣𝜉,𝑦𝑦∣𝜉,𝜆𝜆∣𝜉𝑝().𝑦∣𝜉=𝑝𝜆∣𝜉 (A.7)
(iii) The posterior means are 𝖤=𝖤()𝖤=𝜆𝑝=()𝜉∣𝑦𝜉𝑝𝜉∣𝑦𝑑𝜉,𝜆∣𝜉,𝑦=𝜆𝑝𝜆∣𝜉,𝑦𝑑𝜆=𝜆𝑝𝜆∣𝜉𝑑𝜆,𝜆∣𝑦𝜆𝑝𝜆∣𝑦𝑑𝜆=𝜆,𝜉∣𝑦𝑑𝜉𝑑𝜆𝑝𝜉∣𝑦𝜆𝑝𝜆∣𝜉𝑑𝜆𝑑𝜉, (A.8) noting that 𝑝(𝜆∣𝜉,𝑦)=𝑝(𝜆∣𝜉).

Proof of Theorem 3. Consider the expressions in (37) and (38); that is, 𝑝𝜇∣𝑦=𝑝𝜇∣𝜉,𝜆,𝑦𝑝𝜉∣𝑦𝑝(𝜆∣𝜉)𝑑𝜆𝑑𝜉 (A.9) and a similar expression involving 𝑝^†. The question is when they are identical. Assuming 𝑝(𝜇∣𝜉,𝜆,𝑦)=𝑝(𝜇∣𝜉,𝑦) as in (39) the expression reduces to 𝑝=𝑝𝑝𝑝()𝑝𝜇∣𝑦𝜇∣𝜉,𝑦𝜉∣𝑦𝜆∣𝜉𝑑𝜆𝑑𝜉=𝑝𝜇∣𝜉,𝑦𝜉∣𝑦𝑑𝜉, (A.10) since the conditional prior integrates to unity. The same applies for the expression involving 𝑝^†(𝜆∣𝜉).

Proof of Theorem 5. Similar to the proof of Kuang et al. [19, Theorem 1], albeit for a rectangular instead of a triangular data array.

Proof of Theorem 8. Recall 𝜃_intrisic=𝐶_𝑆⟂𝐶^′_𝑆⟂𝐴(𝐴^′𝐶_𝑆⟂𝐶^′_𝑆⟂𝐴)⁻¹𝜉 where 𝐶_𝑆⟂=𝐶_⟂(𝑆^′_⟂𝐶_⟂)⁻¹ as defined in (72). Premultiply by 𝐴^′ and 𝐶^′ to see that 𝐴^′𝜃_intrisic=𝜉 and 𝐶^′𝜃_intrisic. This does, however, not describe the full variation of 𝜃_intrisic since (𝐴,𝐶) is not a square matrix. We must extend the matrix (𝐴,𝐶) with columns so that it is square and invertible.
We find a vector 𝑤∈ℝ^𝑞−3 so that (𝐴,𝐶,𝐶_⟂𝑤) is invertible. Recall 𝐴∈ℝ^𝑞×𝑝 where 𝑞=𝑝+4. The matrix 𝐶∈ℝ^𝑞×3 is chosen so that (𝐴,𝐶) has full column rank 𝑞−1=𝑝+3 as discussed in Section 5.4.1. Apply Lemma A.2, swapping the role of 𝐴,𝐶, to see 𝑤_⟂=𝐶^′_⟂𝐴, say, has full column rank. Then (𝐴,𝐶) has orthogonal complement 𝐶_⟂𝑤.
We show that for any invertible matrix 𝑀∈ℝ^𝑞×𝑞 then the 1×2 block matrix {(𝐴,𝐶),𝑀𝐶_⟂𝑤} is invertible. To see this hold, premultiply by {(𝐴,𝐶),𝐶_⟂𝑤}^′ to see that an invertible upper triangular matrix arises.
To analyse the properties of 𝜃_intrisic it suffices to analyse {(𝐴,𝐶),𝑀𝐶_⟂𝑤}^′𝜃_intrisic, since there is a bijective mapping between the two. Choose 𝑀=𝐶_⟂(𝐶^′_⟂𝑆_⟂)(𝑆^′_⟂𝐶_⟂)𝐶^′_⟂+𝐶𝐶^′. Then it holds that 𝐶^′_⟂𝑀^′=(𝐶^′_⟂𝑆_⟂)(𝑆^′_⟂𝐶_⟂)𝐶^′_⟂ so that 𝐶^′_⟂𝑀^′𝐶_𝑆⟂𝐶^′_𝑆⟂=𝐶^′_⟂ and therefore 𝑤^′𝐶^′_⟂𝑀^′𝜃_intrisic=𝑤^′𝐶^′_⟂𝐴𝐴^′𝐶_𝑆⟂𝐶^′_𝑆⟂𝐴⁻¹𝜉=0, (A.11) since 𝑤^′𝐶^′_⟂𝐴=0 by construction. Thus, it holds that 𝑤^′(𝐶^′_⟂𝑆_⟂)(𝑆^′_⟂𝐶_⟂)𝐶^′_⟂𝜃_intrisic=0 as required.

Proof of Theorem 9. This is a generalisation of the proof of Kuang et al. [25, Theorems 1,2]. Let 𝑘_ℎ=𝐽+ℎ+𝐼−𝑖−𝐾=ℎ−𝑖−1.
(i) Recall the group 𝑔 in (57). Then (i) follows by comparing the equations 𝜇_{𝑖,𝐽+ℎ}𝜃=𝛼_𝑖+𝛽_𝐽+ℎ𝜃+𝛾_𝐾+𝑘ℎ𝜃+𝜇𝛿,_{𝑖,𝐽+ℎ}𝑔𝜃=𝛼_𝑖+𝑎+𝑑𝑖+𝛽_𝐽+ℎ𝑔𝜃+𝛾_𝐾+𝑘ℎ𝑔𝜃+𝛿−𝑎−𝑏−𝑐−𝑑𝐼. (A.12) (ii) As in Section 3.1 there is a bijective mapping from 𝜃 to 𝜆, 𝜉, where 𝜉 is invariant to 𝑔, but 𝜆 is not. The choice of 𝜆 is not important. Any extrapolations of 𝛽_𝑗, 𝛾_𝑘 can then be written in the form 𝛽_𝐽+ℎ𝜃=𝛽_𝐽+ℎΔ𝛽_𝐽+𝐹_𝛽𝜉̂,𝛾𝜆,_𝐾+𝑘ℎ𝜃=𝛾_𝐾+𝑘_ℎΔ𝛾_𝐾+𝐹_𝛾𝜆,𝜉̂, (A.13) for some functions 𝐹_𝛽, 𝐹_𝛾. Applying the group 𝑔 it follows 𝛽_𝐽+ℎ𝑔𝜃=𝛽_𝐽+𝑏−𝑑𝐽+ℎΔ𝛽_𝐽−𝑑+𝐹_𝛽𝑔𝜆,𝜉̂,𝛾_𝐾+𝑘ℎ𝑔𝜃=𝛾_𝐾+𝑐+𝑑𝐾+𝑘_ℎΔ𝛾_𝐾+𝑑+𝐹_𝛾𝑔𝜆,𝜉̂. (A.14) Due to (i) it must hold 𝐹(𝜆,𝜉)=𝐹_𝛽(𝜆,𝜉)+𝐹_𝛾(𝜆,𝜉) and must equal 𝐹{𝑔(𝜆),𝜉}=𝐹_𝛽{𝑔(𝜆),𝜉}+𝐹_𝛾{𝑔(𝜆),𝜉}. The function 𝐹 must then be constant in the first argument. This must apply in the forecast region 𝒥_ap where 𝛾 is not extrapolated. Therefore, it must also hold that the 𝐹_𝛽(𝜆,𝜉)=𝐹_𝛽{𝑔(𝜆),𝜉}, and in turn that 𝐹_𝛾(𝜆,𝜉)=𝐹_𝛾{𝑔(𝜆),𝜉}. Conversely, if the functions 𝐹_𝛽, 𝐹_𝛾 are constant in 𝜆 then the forecast is invariant to 𝑔.

Proof of Theorem 10. (i) Equation (95) shows that 𝜉 is a function of 𝜃.
(ii) Equation (93) shows that 𝜇=𝛾𝜄^′+𝛿𝜄^′_⟂ is a function of 𝜉.
(iii) The decomposition 𝜇=𝜇𝜄𝜄^′+𝜇𝜄^′_⟂𝜄_⟂ shows that there is a one-one mapping between 𝜇 and (𝜇𝜄,𝜇𝜄^′_⟂). In turn, (93) shows that (𝜇𝜄,𝜇𝜄^′_⟂)=(𝛾,𝛿). Thus if (𝛾^†,𝛿^†)≠(𝛾^‡,𝛿^‡) then 𝜇^†≠𝜇^‡.

Proof of Theorem 11. Rewrite the trace term using the identity 𝐼_𝐽=𝜄𝜄^′+𝜄_⟂𝜄^′_⟂ to get 𝒯=tr𝑌−𝜇𝑌−𝜇^′=tr𝑌−𝜇𝜄𝜄^′𝑌−𝜇^′+tr𝑌−𝜇𝜄_⟂𝜄^′_⟂𝑌−𝜇^′. (A.15) By (94) then 𝜇=𝛾𝜄^′+𝛿𝜄^′_⟂ so that 𝒯=𝒯₁+𝒯₂ where 𝒯₁=tr𝑌−𝛾𝜄^′𝜄𝜄^′𝑌−𝛾𝜄^′′=𝜄^′𝜄tr𝑌𝜄−𝛾𝑌𝜄−𝛾^′,𝒯₂=tr𝑌−𝛿𝜄^′_⟂𝜄_⟂𝜄^′_⟂𝑌−𝛿𝜄^′_⟂^′=tr𝑌𝜄_⟂𝜄^′_⟂−𝛿𝜄^′_⟂𝑌𝜄_⟂𝜄^′_⟂−𝛿𝜄^′_⟂^′. (A.16) The term 𝒯₁ has a minimum of zero if and only if 𝑌𝜄=𝛾. In 𝒯₂ replace for a moment 𝛿𝜄^′_⟂ by a matrix 𝜙 with rank of at most one. The altered 𝒯₂ is minimised when 𝜙 is the singular value decomposition of 𝑌𝜄_⟂𝜄^′_⟂ truncated to rank one due; see Golub and van Loan [46, Theorem 2.5.2]. That singular value decomposition has the property that it is zero when multiplied by 𝜄. Therefore it is also the minimiser of the original problem.

Proof of Theorem 12. It follows by comparing (105) and (106).

Proof of Theorem 13. Write 𝜇_{𝑖,𝐽+ℎ,1}(𝜃₁,𝜃₂)=𝛼_𝑖1+𝛽_𝑖1𝜅_𝐽+ℎ,1(𝜅₁,𝜅₂). This is equals to 𝜇_{𝑖,𝐽+ℎ,1}{𝑔₂(𝜃₁),𝑔₂(𝜃₂)}=(𝛼_𝑖1+𝛽_𝑖1𝑐₁)+(𝛽_𝑖1/𝑑₁)𝜅_𝐽+ℎ,1{𝑑₁(𝜅₁−𝑐₁),𝑑₂(𝜅₂−𝑐₂)} under the given condition.

A.3. A Further Result on Time Effect Hypotheses

Consider the restriction 𝑅^′𝜃_𝑅=𝜌 of (28). The following result holds.

Theorem A.3. Consider the restriction 𝑅^′𝜃_𝑅=𝜌 (A.17) of (28) where 𝐴∈ℝ^𝑞×𝑝 and 𝑅∈ℝ^{𝑞×(𝑞−𝑞𝑅)} so 𝑝,𝑞_𝑅≤𝑞. Suppose 𝐴 and 𝑅 have full column rank. Then, for some 𝑝_𝑅≤min(𝑝,𝑞_𝑅) it holds 𝑝_𝑅=rank(𝑅^′_⟂𝐴)=rank(𝐴,𝑅)−(𝑞−𝑞_𝑅). Then write 𝑅^′_⟂𝐴=𝑎𝑏^′ for some matrices 𝑎∈ℝ^{𝑞𝑅×𝑝𝑅} and 𝑏∈ℝ^{𝑝×𝑝𝑅} with full column rank. The hypothesis (28) restricts the canonical parameter affinely through 𝑏^′_⟂𝜉=𝑏^′_⟂𝐴^′𝑅𝜌, (A.18) so that the degrees of freedom of the restriction is 𝑝−𝑝_𝑅. Introduce the parameters 𝜑₁=𝑎^′𝑅^′_⟂𝜃_𝑅,𝜑₂=𝑏^′𝐴^′𝜉=𝜑₁+𝑏^′𝐴^′𝑅𝜌. (A.19) Then, the predictor 𝜇 can be written as the 𝑝_𝑅-dimensional affine subspace of the form 𝑀_𝑅=𝜇∈ℝ^𝑛∶𝜇=𝑋𝑏𝜑₁+𝑋𝐴^′𝑅𝜌𝑓𝑜𝑟𝜑₁∈ℝ^𝑝𝑅=𝜇∈ℝ^𝑛∶𝜇=𝑋𝑏𝜑₂+𝑋𝑏_⟂𝑏^′_⟂𝐴^′𝑅𝜌𝑓𝑜𝑟𝜑₂∈ℝ^𝑝𝑅. (A.20)(A.21) The mapping 𝜃_𝑅↦𝜇=𝑋𝐴^′𝜃_𝑅 on 𝜃_𝑅∈Θ_𝑅 is invariant with respect to the group 𝑔_𝑅∶𝜃_𝑅↦𝜃_𝑅+𝑅_⟂𝑎_⟂𝜁₃ with 𝜑₁ and 𝜑₂ as maximal invariants.
A special case arises if the restriction combines a restriction on 𝐴 with ad hoc identification. That is, if 𝑅=(𝐴𝑑,𝐶) where 𝑑∈ℝ^{𝑝×(𝑝−𝑝𝑅)} and 𝐶∈ℝ^{𝑞×{(𝑞−𝑝)−(𝑞𝑅−𝑝𝑅)}} so (𝐴,𝐶) has full column rank. Then 𝑏=𝑑_⟂ so that the restriction (A.18) reduces to 𝑑^′𝜉=(𝐼_{𝑝−𝑝𝑅},0)𝜌.

Proof of Theorem A.3. Apply Lemma A.1 (i, iii) to see that rank(𝐴,𝑅)=𝑝_𝑅+𝑞−𝑞_𝑅 is equivalent to rank(𝑅^′_⟂𝐴)=𝑝_𝑅. We can then write 𝑅^′_⟂𝐴=𝑎𝑏^′ for some matrices 𝑎∈ℝ^{𝑞𝑅×𝑝𝑅} and 𝑏∈ℝ^{𝑝×𝑝𝑅} with full column rank.
Exploit the projection identity 𝐼_𝑞=𝑅_⟂𝑅^′_⟂+𝑅𝑅^′ to get 𝜉=𝐴^′𝜃_𝑅=𝐴^′𝑅_⟂𝑅^′_⟂𝜃_𝑅+𝐴^′𝑅𝑅^′𝜃_𝑅. Insert 𝐴^′𝑅_⟂=𝑏𝑎^′ and 𝑅^′𝜃_𝑅=𝜌 to get 𝜉=𝑏𝑎^′𝑅^′_⟂𝜃_𝑅+𝐴^′𝑅𝜌 and therefore 𝜃⟼𝜇=𝑋𝐴^′𝜃=𝑋𝑏𝑎^′𝑅^′_⟂𝜃_𝑅+𝑋𝐴^′𝑅𝜌. (A.22) Inserting the orthogonal projection 𝐼_𝑝=𝑏𝑏^′+𝑏_⟂𝑏^′_⟂ this can also be written as 𝜃⟼𝜇=𝑋𝐴^′𝜃=𝑋𝑏𝑎^′𝑅^′_⟂𝜃_𝑅+𝑏^′𝐴^′𝑅𝜌+𝑋𝑏_⟂𝑏^′_⟂𝐴^′𝑅𝜌. (A.23) Noting that 𝜑₁=𝑎^′𝑅^′_⟂𝜃_𝑅 is freely varying so is 𝜑₂=𝑎^′𝑅^′_⟂𝜃_𝑅+𝑏^′𝐴^′𝑅𝜌. To rewrite 𝜑₂ premultiply the first term by 𝐼_𝑝𝑅=𝑏^′𝑏 to get 𝜑₂=𝑏^′(𝑏𝑎^′𝑅^′_⟂𝜃_𝑅+𝐴^′𝑅𝜌). Then insert 𝜌=𝑅^′𝜃_𝑅 and 𝑏𝑎^′=𝐴^′𝑅_⟂ to get 𝜑₂=𝑏^′𝐴^′(𝑅𝑅^′_⟂+𝑅𝑅^′)𝜃_𝑅=𝑏^′𝐴^′𝜃_𝑅=𝑏^′𝐴^′𝜉. This gives the space 𝑀_𝑅 in (A.21).
To derive the reduced group 𝑔_𝑅 choose a 𝜃_𝑅∈Θ_𝑅 so that 𝑅^′𝜃_𝑅=𝜌. Any 𝜃^†_𝑅∈Θ can be written as 𝜃^†_𝑅=𝜃_𝑅+𝑅𝜁₁+𝑅_⟂𝑎𝜁₂+𝑅_⟂𝑎_⟂𝜁₃ since (𝑅,𝑅_⟂) and (𝑎,𝑎_⟂) have full rank. Since 𝜃_𝑅∈Θ_𝑅 then 𝜃^†_𝑅∈Θ_𝑅 if and only if 𝜁₁=0. We now consider whether 𝜃_𝑅, 𝜃^†_𝑅 are equivalent with respect to the restricted mapping (A.22). It holds 𝜃^†_𝑅⟼𝜇^†=𝑋𝑏𝑎^′𝑅^′_⟂𝜃_𝑅+𝑅_⟂𝑎𝜁₂+𝑅_⟂𝑎_⟂𝜁₃+𝑋𝐴^′𝑅𝜌. (A.24) Noting that 𝜇=𝑋𝑏𝑎^′𝑅^′_⟂𝜃_𝑅+𝑋𝐴^′𝑅𝜌 and 𝑎^′𝑅^′_⟂𝑅_⟂𝑎=𝐼_𝑝𝑅, while 𝑎^′𝑅^′_⟂𝑅_⟂𝑎_⟂=0, then 𝜇^†=𝜇+𝑋𝑏𝜁₂ which reduces to 𝜇 if and only if 𝜁₂=0. Thus, the mapping 𝜃_𝑅↦𝜇 on 𝜃_𝑅∈Θ_𝑅 is invariant with respect to the group 𝑔_𝑅∶𝜃_𝑅↦𝜃_𝑅+𝑅_⟂𝑎_⟂𝜁₃ with 𝜑₁ and 𝜑₂ as maximal invariants.
Now, suppose 𝑅=(𝐴𝑑,𝐶). Then Lemma A.2 shows that 𝑚=𝐴^′_⟂𝐶 has full column rank, while (𝐴,𝐶) has orthogonal complement 𝐴_⟂𝑚_⟂. Thus, 𝑅 has orthogonal complement 𝑅_⟂=𝑑𝐴,𝐶^⟂0,(𝐴,𝐶)_⟂=𝐼𝐴,𝐶^𝑝0𝑑_⟂,𝐴_⟂𝑚_⟂. (A.25) In particular, it holds 𝐴^′𝑅_⟂=(𝐼_𝑝,0)(𝐴,𝐶)^′𝑅_⟂=(𝑑_⟂,0)=𝑑_⟂(𝐼_𝑝𝑅,0), which is denoted 𝑏𝑎^′ in the general case, so that 𝑏=𝑑_⟂. Consider the restriction (A.18). Here, 𝑏^′_⟂𝜉=𝑑^′𝜉, while 𝑏^′_⟂𝐴^′=𝑑^′𝐴^′=(𝐼_{𝑝−𝑝𝑅},0)(𝐴𝑑,𝐶)^′=(𝐼_{𝑝−𝑝𝑅},0)𝑅^′ so that 𝑏^′_⟂𝐴^′𝑅𝜌=(𝐼_{𝑝−𝑝𝑅},0)𝜌.

Conflict of Interests

The authors do not have any conflict of interests.

Acknowledgment

Comments from Andrew Hunt and María Dolores Martínez Miranda are gratefully acknowledged.