Hong Li
email: hong.li@nankai.edu.cn

Yang Lu

Coherent Forecasting of Mortality Rates: A Sparse Vector-Autoregression Approach

Keywords: cohort effect, spatial co-integration, penalized least square

This paper proposes a spatial-temporal autoregressive model for the mortality surface, where mortality rates of each age depend on the historical values of itself (temporality) and the neighbouring ages (spatiality). The mortality dynamics is formulated as a large, first order vector autoregressive model which encompasses standard factor models such as the Lee and Carter (1992) model. Sparsity and smoothness constraints are then introduced, based on the idea that the nearer the two ages, the more important the dependence between mortalities at these ages. Our model has several novelties. First, it ensures that in the long-run, mortality rates at different ages do not diverge. Second, it provides a natural explanation of the so-called cohort effect without identifiability difficulties. Third, the model is easily extended to the multiple-population case in a coherent way. Finally, the model is associated with a closed form, non-parametric estimation method: the penalized least square, which ensures spatial smoothness of the age-dependent parameters.

Using US and UK mortality data, we find that our model produces reasonable projected mortality profile in the long-run, as well as satisfying short-run out-ofsample forecast performance.

Introduction

Forecasting mortality rates has become an important task for insurance companies, pension funds, as well as policy makers, due to the continuous longevity improvement around the world. In most developing countries, the life expectancy increases at an average rate of 0.25 years per annum, and it is estimated that a three-year increase of the life expectancy spells a 50 % percent increase of the cost for various retirement systems1 .

Nevertheless, despite great recent advances in the literature of mortality modelling, several critical issues remain partially unsolved. First, most models decompose mortality improvements into the product of latent mortality trend factors, which are assumed to be non-stationary, and the corresponding age-specific loadings of these trends. As a consequence, the projected mortality rates at different ages will diverge in the longrun, unless all age-specific sensitivity coefficients are identical. For example, [START_REF] Li | Extending the Lee-Carter Method to Model the Rotation of Age Patterns of Mortality Decline for Long-Term Projections[END_REF] note that the difference in the estimated age-specific loadings in the Lee-Carter model leads to increasingly large proportional differences in the projected death rates at adjacent ages. In particular, the projected mortality rates are implausibly low for infant and younger ages relative to older ages, which are inconsistent with the common belief that the mortality profile should vary smoothly and continuously across age. To ensure non-diverging mortality projection in the long-run, [START_REF] Li | Extending the Lee-Carter Method to Model the Rotation of Age Patterns of Mortality Decline for Long-Term Projections[END_REF] propose to gradually rotate the age-specific loadings, in order for them be the same for the majority of ages at some point in the future. The rotation approach, however, depends on expert judgements rather than be data driven. The issue of diverging projected mortality rates applies also to other mortality models such as the Cairns-Blake-Dowd model [START_REF] Cairns | A Two-Factor Model for Stochastic Mortality with Parameter Uncertainty: Theory and Calibration[END_REF]] and the Age-period-cohort model [START_REF] Renshaw | A Cohort-Based Extension to the Lee-Carter Model for Mortality Reduction Factors[END_REF]]. In recent years, modelling mortality rates of multiple populations has gained greater attentions [see, e.g., [START_REF] Li | Coherent Mortality Forecasts for a Group of Populations: An Extension of the Lee-Carter Method[END_REF]; [START_REF] Dowd | A Gravity Model of Mortality Rates for Two Related Populations[END_REF]; [START_REF] Zhou | Modeling Period Effects in Multi-Population Mortality Models: Applications to Solvency II[END_REF]; [START_REF] Hyndman | Coherent Mortality Forecasting: the Product-Ratio Method with Functional Time Series Models[END_REF]]. Most multi-population models satisfy a population coherence assumption, i.e., the projected mortality of different populations for the same age will not diverge in the long-run. However, in these multi-population models, projected mortality rates at different age may still diverge.

Second, although empirical works usually document a significantly better fit when a cohort effect is included [see e.g. [START_REF] Cairns | A Quantitative Comparison of Stochastic Mortality Models Using Data from England and Wales and the United States[END_REF]] along with the period effect, models with age, period, and cohort (APC) effects usually suffer from difficulties in terms of identifiability and interpretability [see e.g. [START_REF] Heckman | Using Longitudinal Data to Estimate Age, Period and Cohort Effects in Earnings Equations[END_REF]; [START_REF] Kuang | Identification of the Age-period-Cohort Model and the Extended Chain-Ladder Model[END_REF]; [START_REF] Hunt | Robustness and Convergence in the Lee-Carter Model with Cohort Effects[END_REF]]. Roughly speaking, different APC models can have the same in-sample fit, but produce substantially different mortality forecasts. This raises doubts on the reliability of such models for forecasting purpose.

Finally, while the widely used factor-based mortality models have rather intuitive structures, the statistical properties of the associated estimation method are not sufficiently understood. Indeed, they are typically non standard due to the non-stationarity of the mortality processes, which renders standard consistency results of, say, maximum likelihood estimation, invalid. In particular, recently, [START_REF] Leng | Inference Pitfalls in Lee-Carter Model for Forecasting Mortality[END_REF] show that the two-stage estimation method of the Lee-Carter model is inconsistent in some cases [see also [START_REF] Callot | Deterministic and Stochastic Trends in the Lee-Carter Mortality Model[END_REF] for a discussion].

This paper addresses the above-mentioned issues within the framework of spatialtemporal autoregressive (STAR) process. The terminology "spatial" is borrowed from the spatial statistics and econometrics literature, which studies, roughly speaking, how time series of quantities observed at various places are correlated. While our paper does not take into account the physical distance between different countries in our analysis, our model has a similar spirit, in the sense that the spatial dependence in the space of different ages is explored. Thus the narrow interpretation of "spatiality" in our paper is the dependence of mortality between neighbouring ages. More precisely, our model emphasizes the similarities of mortality rates at neighbouring ages in terms of: 1) the long-run behaviour; ii) the short-term shock they receive; iii) the smoothness of the age-dependent parameters. Our starting point is a constrained, first-order Vector Autoregressive (VAR) model for the (log) mortality rates at different ages. The parsimony is achieved via appropriate sparsity constraints, which leads to a spatial dependence pattern, in the sense that the closer the two ages, the stronger the dependence between the mortality rates at these ages. Our model has several advantages. First, the constraints we introduce concern both the autoregressive matrix, which captures the Granger causality between mortality rates at different ages, and the variance-covariance matrix of the residual, which captures the instantaneous causality, that is the dependence between residuals 2 . This improves the standard mortality literature, which usually omits the latter. Second, the model implies co-integrated mortality rates at different ages, with the same co-integration vector (1, -1). In other words, the difference of mortality between any two different ages is stationary. Moreover, unlike common factor models, the long-run co-movement constraint is accompanied by a flexible short-term joint dynamics. Indeed, thanks to the spatial dependence pattern of our model, the long-run error correction mechanism takes effect gradually across the age range, which guarantees the flexibility of the short-term dynamics. Third, our model provides a natural interpretation of the so-called cohort effect, without introducing, in an ad hoc way, a separate cohort-specific factor, as is done in the existing literature. Fourth, the model can be readily extended to the simultaneous modelling of two or more populations, such as males and females in a same country, or a group of countries with similar socio-economic characteristics, when the inter-population dependence is introduced. Our model ensures that the (log-) mortality differences between any two ages in any two populations is stationary, with a age-and population-specific equilibrium mean. Finally, in terms of the statistical inference, in order to ensure smoothness of the age-dependent parameters of the model, we propose a two-stage, penalized least square (PLS) estimator. This approach is essentially non-parametric, and contributes to the growing literature on the use of flexible, functional time series to forecast mortality [see e.g. [START_REF] Hyndman | Coherent Mortality Forecasting: the Product-Ratio Method with Functional Time Series Models[END_REF]; [START_REF] Li | A Semiparametric Panel Approach to Mortality Modeling[END_REF]]. We derive the closed form formula of the PLS estimate, and show how to determine the optimal choice of the penalization term via cross-validation.

Our model is applied to US and UK mortality data, and is compared to the benchmark models of [START_REF] Lee | Modeling and Forecasting US Mortality[END_REF] and [START_REF] Li | Coherent Mortality Forecasts for a Group of Populations: An Extension of the Lee-Carter Method[END_REF]. In the empirical analysis, we find that our model outperforms these benchmarks in terms of forecast error, and is able captures the fact that the mortality improvements during the last decades has been faster than before.

The rest of the paper is organized as follows. Section 2 presents the model. Section 3 provides a simple extension of the model for multiple populations. Section 4 discusses the statistical inference. Section 5 applies the model to US and UK populations. Section 6 concludes. Technical details and large figures are gathered in Appendix.

2 See Appendix 4 for the definitions of the Granger causality and the instantaneous causality.

A spatial-temporal autoregressive (STAR) model

The standard literature

The aim of this subsection is not to provide a systematic review of the literature on mortality models [see e.g. [START_REF] Cairns | Modelling and Management of Mortality Risk: a Review[END_REF]; [START_REF] Barrieu | Understanding, Modelling and Managing Longevity Risk: Key Issues and Main Challenges[END_REF]], but rather to present the mainstream methodology, and discuss some of the downsides that our paper seeks to improve.

The Lee-Carter model. Let us denote by m i,t the crude death rate at age i and date t, and y i,t the log transformed death rate, y i,t = log m i,t . For expository purpose we re-index the observable ages by i = 1, 2, ..., I, and the observable dates by t = 1, 2..., T .

In a seminal paper, [START_REF] Lee | Modeling and Forecasting US Mortality[END_REF] postulate that the (log) mortality rates at different ages are captured by a common factor, multiplied by age-specific sensitivity coefficients with respect to the common trend. More precisely, they assume that:

y i,t = a i + b i k t + i,t , ∀i, t, (2.1)
where a i is the time average of y i,t 's, k t is the common factor, and b i is the age-specific sensitivity coefficient for age i. One inconvenience of this approach is that the mortality forecast at different ages diverge in the long run, unless the sensitivity coefficients b i are equal for all ages. As a result, the projected mortality profile would be increasingly discontinuous. For example, assume b i < b i for some i = i . Then the ratio

y i,t y i ,t
increases unlimitedly when t goes to infinity.

Recently, [START_REF] Li | Extending the Lee-Carter Method to Model the Rotation of Age Patterns of Mortality Decline for Long-Term Projections[END_REF] suggest to gradually rotate the sensitivity coefficients, so that these exposures would be equal for the majority of age at some point in the future.

Nevertheless, this approach is rather ad hoc and does not allow a clear justification.

Moreover, it is not clear how to generalize this approach to multi-population mortality models, and/or to take into account the cohort effect.

Models with cohort effect. The demographic literature has documented the existence of a cohort effect that is unexplained by typical age-period based mortality models, such as the Lee-Carter model. For instance, the so-called UK golden cohort(s), who are born between 1925 and 1934, have enjoyed a faster mortality improvement compared to neighbouring cohorts, and this is widely credited to the better childhood health, as well as the war time rationing in their adolescent3 . This phenomenon has motivated various extensions of the Lee-Carter model. For instance4 , [START_REF] Renshaw | A Cohort-Based Extension to the Lee-Carter Model for Mortality Reduction Factors[END_REF] propose:

y i,t = a i + b i k t + c i γ t-i + i,t , (2.2)
where γ t-i captures the cohort effect. One disadvantage of such extensions, despite their gain in terms of goodness-of-fit, is that they suffer from the identifiability difficulty, due to the simple relation that Calendar year = Birth date (that is the cohort Id) + current age of the cohort [see e.g. [START_REF] Heckman | Using Longitudinal Data to Estimate Age, Period and Cohort Effects in Earnings Equations[END_REF], [START_REF] Kuang | Identification of the Age-period-Cohort Model and the Extended Chain-Ladder Model[END_REF], [START_REF] Hunt | Robustness and Convergence in the Lee-Carter Model with Cohort Effects[END_REF]]. Roughly speaking, there exist more than one set of observationally equivalent parameters. This raises doubts of the suitability of such models for forecasting purpose, as the projected mortality rates depend on the (arbitrarily chosen) extra identification constraints.

Another downside of factor models is that due to the non-stationarity of the mortality rates, the consistency of maximum likelihood estimate on factor based models is a tricky issue. For instance, [START_REF] Leng | Inference Pitfalls in Lee-Carter Model for Forecasting Mortality[END_REF] show that, in some cases, the two-stage estimate of the Lee-Carter model is inconsistent.

The vector-autoregressive model and the Granger causality

In this paper, we propose a (sparse) vector autoregressive (VAR) model for the process y t = (y 1,t , y 2,t , ..., y I,t ). More precisely, let us assume that:

y 1,t+1 = y 1,t + m 1 + 1,t , (2.3) y 2,t+1 = (1 -α 2 )y 2,t + α 2 y 1,t + m 2 + 2,t , (2.4)
for the two initial ages, and (2.5) for i = 2, ..., I -1 and t = 1, ..., T -1. In Equation (2.5), α i+1 and β i+1 are positive parameters that are smaller than 1, and the vector residual process ( i,t ) t is a strong Gaussian white noise. That is, it does not feature temporal dependence and follows a

y i+1,t+1 = (1 -α i+1 -β i+1 )y i+1,t + α i+1 y i,t + β i+1 y i-1,t + m i+1 + i+1,t ,
Gaussian distribution N (0, Σ).

The necessity to specify different dynamic equations for (y 2,t ) and (y 1,t ) is due to the lack of observation of mortality rates for ages i ≤ 0. However, we will show later that in our model, for each i, process (y i,t ) t has a random walk dynamics. This, as a consequence, rationalises the specification (2.3) and (2.4).

The term α i+1 y i,t on the RHS of (2.5) captures the cohort effect, that is the persistent effect of the mortality shocks on the same cohort, born in year t -i: the larger α i+1 , the more persistent the impact of the shocks. The term (1 -α i+1 -β i+1 )y i+1,t can be either interpreted as the "learning" effect of the cohort t-i from the (older) neighbouring cohort t -i -1; or as the period effect, that is, it captures the serial correlation of the mortality for the fixed age i. Finally, the term β i+1 y i-1,t captures the "learning" effect of the cohort t -i from the (younger) neighbouring cohort born in t -i + 1.

Thus the Granger causality between the mortality rates of different cohorts only exists between neighbouring cohorts, and within the same cohort. As we will show later, this assumption generalizes the common factor model [see e.g. [START_REF] Lee | Modeling and Forecasting US Mortality[END_REF]],

while still keeping parsimony of the VAR model. Finally, the assumption that the three coefficients, (1 -α i+1 -β i+1 ), α i+1 , and β i+1 sum up to one ensures that the differences between different component processes of (y t ) are stationary [see Section 2.2.1].

Under this specification, the vector process (y 1,t , ..., y I,t ) t has a (large, but sparse) first-order VAR representation:

             y 1,t+1
y 2,t+1 y 3,t+1 . . .

y I,t+1              =              1 0 0 . . . . . . α 2 1 -α 2 0 . . . . . . β 3 α 3 1 -α 3 -β 3 0 . . . . . . . . . . . . . . . . . . . . . 0 β I α I (1 -α I -β I )              R              y 1,t y 2,t y 3,t . . . y I,t              +              m 1,t m 2,t m 3,t . . . m I,t              +              1,t+1 2,t+1 3,t+1
. . . Therefore, our model provides a natural interpretation of the cohort effect, which arises due to the persistence of shocks received by the same cohort during the early life of this cohort. Moreover, due to the learning effect between different cohorts, the cohort effect is expected to be similar for neighbouring cohorts. To summarize, our framework permits correlation not only between cohort effects of different cohorts, but also between the cohort effects and past mortality shocks. These dependencies are typically ignored in standard mortality models, which specify the period and cohort factors independently, and usually do not specify the evolution of the cohort effect.

I,t+1              , ( 2
In mortality modelling, VAR models (or similar vector error correction models) have been suggested by, for example, [START_REF] Lazar | A Multivariate Time Series Approach to Projected Life Tables[END_REF] and [START_REF] Salhi | Basis Risk Modelling: a Co-integration Based Approach[END_REF].

These authors propose to fit unconstrained, large dimensional vector autoregressive models to the time series of mortality rates. However, given the limited sample size of most mortality data, this approach suffers from the curse of dimensionality [see, e.g. [START_REF] Litterman | Forecasting with Bayesian Vector Autoregressions-Five Years of Experience[END_REF]; [START_REF] Tibshirani | Regression Shrinkage and Selection via the Lasso[END_REF] for discussions]. Our constrained model mitigates this problem by incorporating only a small number of lagged regressors that are easily interpretable.

A VAR model has two major ingredients: the Granger causality matrix R in (2.6), which determines the long-run co-movement, or co-integration of the mortality rates, as well as the variance-covariance matrix of the residual, which captures the short-term dependence of the mortality rates. In the mortality literature, these two dependencies are usually captured by different models. For instance, [START_REF] Biffis | A Bidimensional Approach to Mortality Risk[END_REF]; [START_REF] Debón | Modelling Residuals Dependence in Dynamic Life Tables: A Geostatistical Approach[END_REF]; [START_REF] Mavros | Stochastic Mortality Modelling: Key Drivers and Dependent Residuals[END_REF]; [START_REF] Chuliá | Modeling Longevity Risk With Generalized Dynamic Factor Models and Vine Copulae[END_REF] propose to decompose the mortality rates into a "central tendency" process, and a residual process. The first is assumed to follow a standard common factor mortality model such as the Lee-Carter, or the CBD model, whereas the residual is specified using similar spatial constraints. Our approach has several advantages compared to this conventional approach. First, factor models generically do not guarantee the coherence of different age-specific mortality processes (y i,t ) t . Second, such a decomposition is not necessarily unique. In particular, it is shown in the following examples that two of the most popular common factor models, that are the Lee-Carter model and the CBD model, can be represented as special examples of our constrained VAR model.

Example 1: Consider the case where all α i , for i = 2, 3, ..., I and β i , for i = 3, ..., I are equal to zero. Then we have a pure age-period model for the time series (y i+1 , t) t , that is:

y i,t+1 = y i,t + m i + i,t+1 , ∀i ≥ 1.
In particular, if the residuals ( i,t ) satisfies the reduced rank condition:

i,t = m i t , (2.7)
where sequence ( t ) does not depend on i, then we get

y i,t = y i,1 + m i κ t ,
where

κ t = (t -1 + t τ =1 i,τ
) is a random walk with drift. In other words we have obtained the [START_REF] Lee | Modeling and Forecasting US Mortality[END_REF] model.

If, instead of equation (2.7), residuals can be decomposed into:

i,t = 1,t + i 2,t
(2.8) then we have

y i,t = y i,1 + κ 1,t + iκ 2,t ,
which is CBD-type model [see e.g. [START_REF] Cairns | A Quantitative Comparison of Stochastic Mortality Models Using Data from England and Wales and the United States[END_REF]]. Thus we have shown that both the Lee-Carter model and the CBD model can be expressed as a special case of our VAR model. This motivates our modelling strategy compared to the decomposition approach of [START_REF] Debón | Modelling Residuals Dependence in Dynamic Life Tables: A Geostatistical Approach[END_REF]; [START_REF] Mavros | Stochastic Mortality Modelling: Key Drivers and Dependent Residuals[END_REF].

Example 2: Let us now consider the case where α i = 1 for each age i = 2, 3, ..., I, and all β i are zero, for i = 3, ..., I. Then we get:

y i+1,t+1 = y i,t + m i+1 + i+1,t+1 , (2.9)
or equivalently:

y i,t = y 1,t-i + t-1 τ =1 m i+1-τ + t-1 τ =1 i+1-τ,t+1-τ
Thus y i+1,t+1 only depends on the past residuals via i+1-τ,t+1-τ , which belong to the single birth cohort t -i. Thus we have an age-cohort model. Compared to the other stochastic models with cohort effect, for example, [START_REF] Renshaw | A Cohort-Based Extension to the Lee-Carter Model for Mortality Reduction Factors[END_REF], in Model (2.9), the cohort effect evolves as the time increases, instead of being fixed at the birth of the cohort.

Analysis of co-movement

The introduction of the spatial dependence structure has an effect of "smoothing" the vector of mortality rates at a same, fixed time. Let us now discuss its implications for the long-run co-movement of the mortality rates at different ages. Roughly speaking, we want to ensure that in the future, mortality rates at different ages do not diverge.

This is an important criteria for a reasonable mortality model, given that longevity risk is usually a very long-term risk.

Such a non-divergence property is called co-integration [see e.g. [START_REF] Engle | Co-integration and Error Correction: Representation, Estimation, and Testing[END_REF]] and has received much attention in the mortality literature [see Gaille and Sherris different ages i and j, processes (y i,t ) t and (y j,t ) t are co-integrated if: i) both processes are integrated of order one; ii) there exists a constant β i,j such that process y i,t -β i,j y j,t (2.10) is stationary.

In this subsection we will prove a stronger co-integration result:

Proposition 1. Indeed under the specification of Section 2, different component processes (y i,t ) t and (y j,t ) t are co-integrated, with co-integration vector (1, -1).

Proof. First, equation (2.3) indicates that for the lowest age i = 1, process (y 1,t ) t is a random walk with drift m 1 . Then we have, for age i = 2:

y 2,t+1 -y 1,t+1 = (1 -α 2 )(y 2,t -y 1,t ) + m 2 -m 1 + 2,t+1 -1,t+1 ,
in other words process (y 2,t -y 1,t ) is stationary. Similarly, we have:

y 3,t+1 -y 2,t+1 =(1 -α 3 -β 3 )y 3,t + α 3 y 2,t + β 3 y 1,t -(1 -α 2 )y 2,t -α 2 y 1,t + m 3 -m 2 + 3,t+1 -2,t+1 =(1 -α 3 -β 3 )(y 3,t -y 2,t ) -(β 3 -α 2 )(y 2,t -y 1,t ) + m 3 -m 2 + 3,t+1 -2,t+1 . Therefore, 1-(1-α 3 -β 3 )L (y 3,t -y 2,t
) is stationary, where L is the lag operator. Since

(1 -α 3 -β 3 ) lies between 0 and 1, (y 3,t -y 2,t ) is also stationary for each i. By induction, we can show that each process (y i+1,t -y i,t ) is stationary. Thus all the processes (y i,t ) t follow a random walk with the same drift m 1 .

The co-integration property of our model ensures the projected mortality curves are non-divergent. At the same time, component processes are allowed to have flexible shortterm dynamics, especially between remote ages, since it takes h periods for the impact of a shock i,t to reach another cohort that is h years older (resp. younger).

While it is possible to conduct formal test of the assumption on the co-integration rank on the vector of mortality rates (I -1 in model (2.3) to (2.5)) using Johansen's general methodology [see e.g. [START_REF] Lazar | A Multivariate Time Series Approach to Projected Life Tables[END_REF]], such tests are unlikely to have a high power due to data limitation. In particular, a typical dataset may contain 100 ages and fewer than 70 years of obervations. Therefore, while [START_REF] Lazar | A Multivariate Time Series Approach to Projected Life Tables[END_REF] found a co-integration rank of only I -3 using an unconstrained co-integration test, we believe that the coherent property implied by our constrained model is valuable, since it can generate reasonable projections of future mortality rates. This remark also applies to other models with a similar spirit of imposing coherence, such as the Li-Lee model [START_REF] Li | Coherent Mortality Forecasts for a Group of Populations: An Extension of the Lee-Carter Method[END_REF]] and the gravity model [START_REF] Dowd | A Gravity Model of Mortality Rates for Two Related Populations[END_REF]].

The spatial dependence of the residuals

Let us now discuss the specification of the residual ( t ), which captures the instantaneous causality between different component processes (y i,t ) t . The standard literature assumes independence between its different components. To illustrate the necessity of accounting for this dependence, we first simulate paths of the mortality processes which share the same causality matrix, but have different variance-covariance matrix Σ of the residuals.

We will see that quite different patterns emerge when we alter the specification of this latter matrix. We then introduce a spatial autoregressive model for the joint distribution of the residual process.

An illustration by simulation

Let us now report, in a simulation experiment, the role of the cross-sectional dependence structure of the residual process ( t ) on the dynamics of the observed mortality rates.

For expository purpose, we only plot three component processes y t = (y 1,t , y 2,t , y 3,t ).

Example 2 continued. Let us first re-consider Example 2, which has a pure agecohort representation. We assume that α 2 = α 3 = 1, and β 2 = β 3 = 0, and m 1 = 0.02, m 2 = m 3 = 0.06. Firstly, we assume that 1,t , 2,t , 3,t are i.i.d. normal with N (0, 0.05 2 ).

We plot in Figure 2 the path of y 1,t , y 2,t and y 3,t , for periods t = 1, ..., 100. We can observe that it takes one (resp. two) periods for large shocks on (y 1,t ) to propagate to (y 2,t ) t (resp. (y 3,t ) t ). We then keep the values of α 2 , α 3 , β 2 , β 3 but assume that the residuals are identical. The resulting simulated paths of y t are plotted in Figure 3. We can see that the identical residuals make the co-movement between the three processes much more pronounced than the independent residuals shown in Figure 2. A similar case as Example 1. Let us now illustrate the impact of the residual structure on a model with both a period effect and a cohort effect. Figure 4 plots the simulated paths of y t , when the residuals are independent, but α 2 , α 3 , β 2 , β 3 are close to 0. For illustration, we have chosen the value 0.05 for all of these four parameters.

The reason of not choosing the exact value of zero is that in this limiting case, the co-integration relationship between the three time series no longer holds (except when the reduced rank condition of Example 1 is satisfied) and we will have instead three independent random walks. We can remark that compared to the Figure 2, the co-movement is much less pronounced. This is expected, since the inter-cohort transmission term α i , β i are all close to zero. In Figure 5, we keep the same assumption as in Figure 4, except that we assume the components of the residual are identical. Similar to the age-cohort representation, identical residuals lead to a much more pronounced co-movement between the three components than under the assumption of independent residuals.

The model

Let us now propose a specification of the residual processes. There are two approaches to account for the dependence among residuals. The first one is to use parametric models with a small number of parameters to characterize the dependence structure.

For instance, [START_REF] Biffis | A Bidimensional Approach to Mortality Risk[END_REF] and [START_REF] Debón | Modelling Residuals Dependence in Dynamic Life Tables: A Geostatistical Approach[END_REF] propose a spatial random field model in which the covariance E[ i,t j,t ] is only a function of the distance |i -j|, but not of ages i and j. Such assumptions may, however, be too restrictive, as residuals at high ages tend to be much more erratic. Another approach is to leave the variance-covariance matrix completely unconstrained [see e.g. [START_REF] Mavros | Stochastic Mortality Modelling: Key Drivers and Dependent Residuals[END_REF]]. The downside of the latter is the curse of dimensionality 5 , since the number of coefficients to estimate is equal to

I(I+1)

2

, which is quite large so long as I reaches 30 or above.

As a trade-off between these two approaches, we propose a constrained spatial autoregressive model. More precisely, let us consider the following simultaneous equation 5 A recent paper by [START_REF] Chuliá | Modeling Longevity Risk With Generalized Dynamic Factor Models and Vine Copulae[END_REF] do not specify the variance-covariance matrix, but the joint copula function, using a vine structure. Although this specification leads also to O(I) parameters, the structure of the vine, such as which nodes are connected and which are not, should be chosen in a somehow ad hoc way.

I,t = c I I-1,t + η I,t I-1,t = a I-1 I,t + c I-1 I-2,t + η I-1,t I-2,t = a I-2 I-1,t + c I-2 I-3,t + η I-1,t
(2.11)

• • • • • • 1,t = a 1 2,t + η 1,t ,
where η t = (η 1,t , η 2,t , ..., η i,t ) is a Gaussian white noise whose components are independent and have variance σ 2 i for each i, and parameters a i , c i are nonnegative and satisfy the stationary condition:

a i + c i < 1, ∀i = 1, 2, ..., I.
(2.12)

In other words, we assume that each component i,t can be predicted by its neighbouring terms i+1,t and i-1,t alone. This is a first-order, (non-homogeneous) spatial autoregressive model, and is the analogy of the (spatial-temporal) autoregressive model specified in (2.3) to (2.5). A direct consequence is that given two components i,t and j,t , the larger |i -j|, the looser the association between them7 . This specification is more flexible than the parametric approach of [START_REF] Biffis | A Bidimensional Approach to Mortality Risk[END_REF] and [START_REF] Debón | Modelling Residuals Dependence in Dynamic Life Tables: A Geostatistical Approach[END_REF], in the sense that parameters a i , c i are age-dependant and can capture a variety of patterns. In the empirical application, we also propose a non-parametric smoothing technique, which allows these parameters to feature smooth variation across different ages. At the same time, this specification involves 2 × (I -1) parameters, which is significantly smaller than I(I+1)

2

, that is the number of parameters, when the matrix Σ is left unconstrained.

In matrix terms, we have:

             1 -c I 0 . . . . . . -a I-1 1 -c I-1 . . . . . . . . . . . . . . . . . . . . . . . . 0 -a 2 1 -c 2 . . . . . . 0 -a 1 1              M t = η t , (2.13)
where matrix M is invertible under assumption (2.12) 8 ã ȂĆ Therefore, the previous simultaneous system can be rewritten in the reduced form:

t = M -1 η t ,
and the variance-covariance matrix of t is equal to:

Σ = E[ t t ] = (M -1 ) Diag(σ 2 1 , σ 2 2 , ..., σ 2 I )M -1 .
Note that although only elements on the principal diagonal and the upper and lower sub-diagonals of M are non null, in general all the elements of Σ are non null. However, for suitable values of a i , c i and σ i , we can expect that the larger |i -j|, the smaller the covariance E[ i,t j,t ].

A multi-population extension

The aim of this section is to show how the spatial-temporal autoregressive model introduced in Section 2 can be extended to the simultaneous analysis of several populations.

Such joint models are quickly gaining popularity in the literature due to: i) The need for multi-national insurers to evaluate the financial benefits of a geographically diversified portfolio of policyholders [see e.g. [START_REF] Hunt | Modelling Longevity Bonds: Analysing the Swiss Re Kortis Bond[END_REF]]; ii) The need for policy mak-

8 Indeed matrix M is diagonally dominant. Note that this condition is sufficient but not necessary. Nevertheless, we have chosen to impose it since i) it is the stationary condition in spatial autoregressive models, when the cross-sectional dimension I goes to infinity. Roughly speaking, it means that neighbouring observations i+1,t and i-1,t can only partially predict the value of i,t. ii) In the empirical application, we can check that this condition is satisfied.

ers to conduct demographic forecasting on a, say, region-by-region basis. For expository purpose, in this section, we first discuss the two-population case, and generalize the idea to situations with more than two populations.

Let us denote by y 1,i,t and y 2,i,t the log of the mortality rates of the two populations 1 and 2, and assume the following joint dynamic specification:

y j,i+1,t+1 = ρ j,i+1 (1 -α j,i+1 )y j,i+1,t + α j,i+1 y j,i,t + (1 -ρ j,i+1 )y -j,i+1,t + m j,i+1 + j,i+1,t , (3.1)
and

y j,1,t+1 = ρ j,1 y j,1,t + (1 -ρ j,1 )y -j,1,t + m j,1 + j,1,t , (3.2)
with j = 1, 2, and -j = 1 when j = 2 and vice versa, (m j,i ) i -s are constants, and ( j,i,t ) t are Gaussian white noises.

Thus compared to the single population model (2.5), the learning effect from younger cohort, that is β j,i+1 y j,i-1,t , is replaced by a learning effect from the other population 9 .

Let us discuss some special cases. Example 1: if parameters ρ 1,i+1 , ρ 2,i+1 are equal to 1 for all values of i, then there is no cross causality in Equation (3.1). In this case we get two separate spatial-temporal autoregressive models 10 .

Example 2: if, say, ρ 1,i+1 = 1 = ρ 2,i+1 for all i, then Population 1 is the dominant population and Granger causes the mortality rates of population 2, whereas Population 2 is a smaller population such as a sub-population of Population 1. This corresponds, roughly speaking, to the gravity model [see [START_REF] Cairns | Bayesian Stochastic Mortality Modelling for Two Populations[END_REF]; [START_REF] Dowd | A Gravity Model of Mortality Rates for Two Related Populations[END_REF]].

In general, when ρ 1,i+1 , ρ 2,i+1 are within (0, 1) for all i, these two parameters capture the relative importance of one population's own past with respect to the past of the other population, when determining the future mortality rates.

9 It would clearly be possible to retain also the learning effect from the younger cohort. In this case we would have a model of the form, say:

y1,i+1,t+1 = ρ1,i+1 (1-α1,i+1-β1,i+1)y1,i+1,t+α1,i+1y1,i,t+β1,i+1y1,i-1,t +(1-ρ1,i+1)y2,i+1,t+m1,i+1+ 1,i+1,t,
however, empirically we found that our current version has nearly the same prediction power as its extension.

10 However, the residuals of the two populations can still be correlated.

and

y j,1,t+1 = J k=1 ρ j,k,1 y j,k,t + m j,1 + j,1,t , (3.5)
where coefficients ρ j,k,i are positive and sum up to 1, that is, J k=1 ρ j,k,i+1 = 1 for each population j and age x + 1.

Let us now discuss the co-integration in the multi-population case. The following proposition says that the mortality rates are still co-integrated.

Proposition 3. The difference of any two components of the (I × J)-dimensional process (y j,i,t , j = 1, ..., J, i = 1, 2, ..., I) is stationary.

Proof. See Appendix 1.2. These three Propositions [1, 2 and 3] constitute the core contribution of our paper.

They show that in order to ensure the co-integration of different mortality processes, it is no necessary use a common factor model such as Lee-Carter, or an error correction model [such as the gravity model [START_REF] Cairns | Bayesian Stochastic Mortality Modelling for Two Populations[END_REF]]. Rather, the family of VAR models offer much more flexibility, even with a rather sparse autoregressive matrix.

Statistical inference

In practice, the value of I lies between 40 and 100, thus it is computationally cumbersome to conduct maximum likelihood estimation. In this section we propose a two-stage penalized-least square approach. In the first stage, parameters of the autoregression function, α i , β i , and m i , are estimated. These estimates are then plugged into the initial equation to recover the empirical residuals, which are used in the second stage to estimate the variance-covariance matrix. We begin with the simple, benchmark Ordinary Least Square (OLS) estimation, and then generalize to the Penalized Least Square (PLS) estimation, which provides smooth curves of the parameters along the age dimension.

The two-stage approach is similar to the method of [START_REF] Lee | Modeling and Forecasting US Mortality[END_REF], in the sense that both are motivated by the existence of a closed form estimate 11 . Although a one-stage estimation would be more efficient, it falls well beyond the scope of the present paper since it poses much more computational challenges due to the large number of parameters.

The benchmark OLS of the Granger causality parameters

Let us first consider a single population. We rewrite (2.5) as a stationary autoregression:

y i+1,t+1 -y i+1,t = α i+1 (y i,t -y i+1,t )+β i+1 (y i-1,t -y i+1,t )+m i+1 + i+1,t , ∀i ≥ 2, (4.1)
and

y 1,t+1 -y 1,t = m 1 + 1,t (4.2) y 2,t+1 -y 2,t = α 2 (y 1,t -y 2,t ) + m 2 + 2,t . (4.3)
Both sides of these equations are stationary, since process (y i+1,t -y i,t ) is stationary for each i. For each i, the pair of parameters (m i+1 , ρ i+1 ) can be estimated separately, since it only appears in one equation. Such a strategy is called Seemingly Unrelated Regression (SUR) [see e.g. [START_REF] Zellner | An Efficient Method of Estimating Seemingly Unrelated Regressions and Tests for Aggregation Bias[END_REF]], in the sense that each OLS is conducted without taking into account the contemporaneous dependence of the residuals ( i,t ). Under mild conditions, this benchmark OLS is asymptotically consistent when the sample size T goes to infinity.

PLS estimation of the Granger causality parameters

The previous estimates of age-dependent parameters, α i , β i and m i are obtained without any smoothness constraint. Such raw estimates are not satisfactory. Indeed, given the relative small sample size of most mortality data (generically, T is of the order of 30-100),

there is a significant risk that the OLS over-fits the data, and leads to erratic curves for these age-dependent parameters.

In this subsection we propose a (spatially) penalized least square estimation method, also called Tikhonov (or L 2 ) regularization [see e.g. [START_REF] Tsybakov | Introduction to Nonparametric Estimation[END_REF]], to ensure the smoothness of the curves α i , β i and m i . The idea is to add a penalty term to the sum of residuals, i.e., the error function of the OLS estimation, in order to penalize abrupt changes of these parameters across different ages. The penalty term is quadratic in the parameters of the model, which leads to an explicit solution for the optimization problem. More precisely, we consider the following penalized least square estimator:

min L 1 (α i ) i=2,...,I , (β i ) i=3,...,I , (m i ) i=1,...,I (4.4)

where

L 1 = I-1 i=2 T -1 t=1 i+1,t+1 (α i+1 , β i+1 , m i+1 ) 2 + T -1 t=1 2,t+1 (α 2 , β 2 , m 2 ) 2 + T -1 t=1 1,t+1 (m 1 ) 2 (4.5) + λ α I-1 i=3 (α i+1 -α i ) 2 + λ β I-1 i=3 (β i+1 -β i ) 2 + λ m I-1 i=3 (m i+1 -m i ) 2 , (4.6)
with λ ρ , λ β , and λ m being known and nonnegative parameters. The larger these parameters, the smoother the curves of α i , β i and m i . In the special case where λ ρ = λ β = λ m = 0, the unique solution of this minimization problem is the OLS estimator.

The right hand side of (4.5) is the sum of squared errors for each age and year, which are obtained using equations (2.3) to (2.5). Moreover, (4.6) corresponds to the penalty terms. Due to the boundary effect, the penalty does not concern the two lowest ages i = 1 and i = 2.

This penalization scheme is similar to [START_REF] Delwarde | Smoothing the Lee-Carter and Poisson Log-bilinear Models for Mortality Forecasting: A Penalized Log-likelihood Approach[END_REF], who add a quadratic penalty term to the log-likelihood function, to estimate the (Poisson) Lee-Carter model.

Then they propose a numerical, Newton-Raphson based algorithm to estimate the values of the parameters. In our least-square based optimization problem, the objective error function L is quadratic in each parameter α i , β i and m i , which leads to a closed form solution. This result is summarized in the following proposition.

Proposition 4. Estimates of α i , β i and m i can be obtained in closed form formula by solving a linear system.

Proof. See Appendix.

PLS estimation of the variance-covariance matrix

Once we obtain the estimate of the causality matrix, we can recover the empirical values of the residuals, which we use to estimate the variance-covariance matrix. Similar as in the previous case, we use the following PLS estimator:

min L 2 (a i ) i=1,...,I-1 , (c i ) i=2,...,I (4.7)
where

L 2 = I-1 i=2 T t=1 η i,t (a i , c i ) 2 + T t=1 η 1,t (a 1 ) 2 + T t=1 η I,t (c I ) 2 (4.8) + λ a I-1 i=2 (a i+1 -a i ) 2 + λ c I-1 i=2 (c i+1 -c i ) 2 . (4.9)
The right hand side of (4.8) is the sum of squared errors from the dynamics of the residuals, which are obtained from (2.11). Again, λ a and λ c are known and nonnegative parameters. When these two parameters are equal to zero, the model can be estimated by seemingly unrelated ordinary least square, on an equation by equation basis. In the general case, we can show that [see Appendix 1] the estimates of a i , c i allow for closed form expressions.

The two-population model

Similarly, we rewrite equation (3.1) into:

y j,i+1,t+1 -y j,i+1,t = -α j,i+1 ρ j,i+1 (y j,i+1,t -y j,i,t ) + (1 -ρ j,i+1 )(y -j,i+1,t -y j,i+1,t )
+m j,i+1 + j,i+1,t , (4.10) and y j,1,t+1 -y j,1,t = (1 -ρ j,1 )(y -j,1,t -y j,1,t ) + m j,1 + j,1,t for j = 1, 2. From Equation (4.10), we see that if we know ρ j,i and α j,i+1 ρ j,i+1 for j = 1, 2, we can in turn compute α j,i+1 . Furthermore, if we treat α j,i+1 as a separate estimator in the penalized estimation, the parameters are not quadratic in the objective function (with a form similar to (4.5) to (4.6)), and a closed form representation of the estimator is not available. Therefore, to preserve the quadratic property of the estimators, we denote by αj,i = α j,i+1 ρ j,i+1 . Moreover, the estimation and the smoothing is done with respect to αj,i instead of α j,i+1 . In this case, the objective function can be written as min L((α j,i ) j=1,2;i=2,...,I , (m j,i ) j=1,2;i=1,...,I , (ρ j,i ) j=1,2;i=1,...,I ), (4.11)

where

L = 2 j=1 I-1 i=1 T -1 t=1 y j,i+1,t+1 -y j,i+1,t -αj,i+1 (y j,i,t -y j,i+1,t ) -(1 -ρ j,i+1 )(y -j,i+1,t -y j,i+1,t ) -m j,i+1 2 + 2 j=1 T -1 t=1 y j,1,t+1 -y j,1,t -(1 -ρ j,1 )(y -j,1,t -y j,1,t ) -m j,1 2 + 2 j=1 λ j,α I-1 i=2 (α j,i+1 -αj,i ) 2 + 2 j=1 λ j,m I-1 i=2 (m j,i+1 -m j,i ) 2 + 2 j=1 λ j,ρ I-1 i=2 (ρ j,i+1 -ρ j,i ) 2 , (4.12)
where -j means 1 when j = 2 and 2 when j = 1. Again, the estimate is obtained by solving a linear system, which is detailed in the Appendix.

The parameters from the variance-covariance matrix are estimated in a way similar to the single-population case as well. In particular, let ãj,i,t = a j,i,t θ j,i,t and cj,i,t = c j,i,t θ j,i,t for all j, i, and t, the objective function is given by min L 2 (ã j,i ) j=1,2;i=2,...,I , (c j,i ) j=1,2;i=2,...,I , (θ j,i ) j=1,2;i=2,...,I , (4.13)

where

L 2 = 2 j=1 I-1 i=2 T t=1 j,i,t -ãj,i j,i+1,t -ci j,i-1,t -(1 -θ j,i ) -j,i,t ) 2 + 2 j=1 T t=1 j,1,t -ãj,1 j,2,t -(1 -θ j,2 -j,2,t 2 + 2 j=1 T t=1 j,I,t -cj,I j,I-1,t -(1 -θ j,I,t ) -j,I,t 2 + 2 j=1 λ j,a I-1 i=2 (ã j,i+1 -ãj,i ) 2 + 2 j=1 λ j,c I-1 i=2 (c j,i+1 -cj,i ) 2 + 2 j=1 λ j,θ I-1 i=2 ( θj,i+1 -θj,i ) 2 .
4.5 The optimal choice of the smoothing parameter: a cross-validation approach Clearly, the PLS estimate depends on the choice of the nonnegative smoothing parameters:

λ 1 = (λ α , λ β , λ m ), λ 2 = (λ a , λ c ).
When they are all zero, that is when no smoothing constraint is imposed, the agedependant parameters are under-smoothed. On the other hand, when they tend to infinity, these parameters become over-smoothed. In the following, we propose a crossvalidation method to choose the optimal smoothing parameters.

The basic idea is to estimate the model, when observations related to one birth cohort are omitted. The resulting estimates are then used to forecast the mortality rates for this omitted birth cohort. The forecasting practices is done for all birth cohorts in the dataset, and the optimal values of λ 1 and λ 2 are the ones which generate the minimal overall forecasting error for all cohorts. Let us first consider the optimal choice of λ 1 .

The choice of λ 1 proceeds as follows. Denote by G the set of all birth cohorts included in the dataset. For each g ∈ G, we estimate equation (4.4) by leaving out all residuals i,t corresponding to the same cohort g, that is, whose indices satisfy t -i = g. This method is called leave-one-out, and the literature has also proposed other similar cross-validation methods such as the h-block method [see e.g. [START_REF] Burman | A Cross-Validatory Method for Dependent Data[END_REF]; [START_REF] Racine | Feasible Cross-Validatory Model Selection for General Stationary Processes[END_REF]].

It is also possible to leave out the observation of one single year. However, given the emphasize of cohort effect in our model, the cohort-based cross validation method seems to be more intuitive12 .

We denote the resulting estimate by θ(-g, λ 1 ), where -g indicates the omission of birth cohort g. Then we compute the out-of-sample forecast of the y i,t -s related to the birth cohort g, which are collected in a column vector, ŷ(-g, λ), and the corresponding forecast error:

ˆ (-g, λ) := y(g) -ŷ(-g, λ), (4.14) where y(g) is a column vector containing all observations y i,t with (i, t) ∈ g. The dimension of y(g) and ŷ(-g, λ) depends on g. The optimal λ 1 is obtained by minimizing the following measure of the forecast error:

λ 1,opt = arg min e 1 (λ 1 ) = arg min g∈G ˆ (-g, λ) ˆ (-g, λ). (4.15)
Once the value of λ 1 is fixed, we fit the model (2.3) -(2.5) to the whole dataset, and obtain the empirical values of the residuals i,t -s. For each g ∈ G, we then minimize (4.7) while leaving out all η i,t -s with (i, t) ∈ g. Then we compute the out-of-sample forecast of i,t related to g, which we collect in a column vector ˆ (-g, λ 2 ),13 and calculate the corresponding forecast error:

η(-g, λ 2 ) := (g) -ˆ (-g, λ 2 ).
Similarly, (g) is a column vector containing the empirical values of i,t -s related to birth cohort g. The optimal value of λ 2 is the one which minimizes the following measure of forecast error:

λ 2,opt = arg min e 2 (λ 2 ) = arg min g∈G η(-t, λ 2 ) η(-t, λ 2 ). (4.16)
For the two-population model, the optimal smoothing parameters are obtained in an analogous manner.

Empirical application

The mortality data used in this study is downloaded from the Human Mortality Database.14

We focus on the uni-sex populations of the US and of England & Wales (henceforth UK).

We use the logarithm of the single-age crude death rates for ages 0 to 99, and from 1950 to 2013 for both populations. We first fit the single population model described in Section 2.2 and 2.3 to the two populations separately, and then fit the two-population model in Section 3.

The single population case

Let us first report the results of the single-population model. Figure 11 and 12 display the PLS-smoothed and the non-smoothed (with all λ = 0) estimates for the two populations.

For the PLS estimate, the optimal values of λ for the two populations are reported in Table 1. From the figures, we see that the parameters estimated from PLS are much smoother than the non-smoothed estimates for both populations. Moreover, for the smoothed estimates, the conditions α i + β i < 1 and a i + c i < 1 hold for all ages and both populations.

Country λ α λ β λ m λ a λ c US 0.18 0.79 0.01 0.08 0.20 UK 0.42 0.79 1.10 0.43 1.21

Table 1: The optimal values of λ for the US and the UK population in the single population case.

The two-population case

Now we apply the two-population model discussed in Section 3 to the US and the UK population at the same time. Figure 13 and 14 display the smoothed and the original parameter estimates for the product and ratio measures, respectively, in the two-population case. For the PLS estimate, the optimal values of λ for the two populations are shown in Table 2. The optimal values of λ for the average and difference mortality rates in the two-population case.

Country λ α λ ρ λ m λ a λ c λ θ US 0.

Forecasting performance

Let us now evaluate the out-of-sample forecasting performance of our model compared to the benchmarks, that are the Lee-Carter model [START_REF] Lee | Modeling and Forecasting US Mortality[END_REF]] in the single population case, and the Li-Lee model [START_REF] Li | Coherent Mortality Forecasts for a Group of Populations: An Extension of the Lee-Carter Method[END_REF]] in the two-population case.

We first fit the single-and two-population STAR model and the benchmark models to data up to year 2000 for each population, and forecast the future mortality rates from year 2001 to 2013. The forecasting results are compared in terms of their respective root mean squared forecast error (RMSFE). More precisely, if we denote by y j,i,t the corresponding forecast, where j = 1 for US and j = 2 for UK, then the RMSFE for population i is defined as:

RM SF E(j) = 1 I(U -û) U u=û I i=1 (y j,i,u -y j,i,u ) 2 , ( 5.1)
where û and U are year 2001 and 2013, respectively. The RMSFE ratio measures the average of the forecasting error for each year u and age x for each population.

Figure 6 and 7 show the RMSFE for both populations in the single-population and the two-population case, respectively. These RMSFE's are functions of the forecast horizon. We see that the STAR model provides more accurate forecasts than the respective benchmarks for all forecast horizons. Therefore, the STAR model appears to better capture the recent mortality trends. In order to understand this superiority, let us compare the forecast of the mortality rates we obtained from the two models. Figure 8 displays the out-of-sample forecasts of the the age-average of y i,t -s from the STAR model and the Lee-Carter model for both populations. Let us now examine the long-term forecast of the STAR models. In Figure 9 we project the period life expectancies at birth for the next 100 years using both the singlepopulation and the two-population STAR models. We can see that the historical increase of life expectancy is faster for the UK population. Hence, the projected increase of life expectancy from the single-population model is also faster for the UK population.

Due to the coherence property, the projected increase of life expectancy from the twopopulation model should be similar for the two populations in the long-run. As a result, 15 The possibility of introducing correlated random walk has been recently examined by [START_REF] Leng | Inference Pitfalls in Lee-Carter Model for Forecasting Mortality[END_REF]. They show that introducing such a correlation leads to potentially biased estimate of the Lee-Carter model. This is due to the non-stationary factor model structure of the Lee-Carter model. On the other hand, our model, which is a Vector Autoregression, does not suffer from this inconsistency. In particular, under mild conditions, the estimates are all asymptotically consistent, when T goes to infinity.

compared to the independent projections, the two-population model leads to, in the long-run, higher projected life expectancies for the US population, and lower projected life expectancies for the UK population. Nevertheless, as we can see, the forecasts of the single-population and two-population STAR model are quite similar in the short term, and the discrepancy only appears after 2050. This demonstrates the ability of the two-population STAR model to reconcile the long-term convergence with the short-term individual dynamics of the two populations.

Conclusion

This paper has proposed a spatial-temporal autoregressive (STAR) approach for mortality modelling. The model is based on the idea that neighbouring ages have higher degree of interdependence, both in terms of the Granger causality and the instantaneously causality. It provides more flexibility than the benchmark factor model, while at the same time ensures suitable long-run properties of the mortality processes. A closed form, penalized least square estimator has been proposed to ensure the smoothness of the age-dependent parameters both in the Granger causality matrix and the variancecovariance matrix. Finally, the STAR model can be easily extended to model multiple unit circle |z| ≥ 1. This is satisfied since:

| ρ 1,2 (1 -α 1,2 ) -z ρ 2,2 (1 -α 1,2 ) -z | ≥ | 1 -ρ 1,2 (1 -α 1,2 ) 1 -ρ 2,2 (1 -α 2,2 ) | ≥ (1 -ρ 1,2 )(1 -ρ 2,2 )
for all z such that |z| ≥ 1.

Similarly, by induction we can show that process (y 1,i+1,t+1 -y 1,i,t+1 , y 2,i+1,t+1 -y 2,i,t+1 ) is stationary.

A.2 Proof of Proposition 3

The proof of Proposition 3 is similar as in the bivariate case. Let us first consider the vector of mortalities at initial ages Y 1,t = (y 1,1,t , y 2,1,t , ..., y J,1,t ). Equation (3.5) says that Y 1,t follows a VAR(1) model, with autoregressive matrix

R 1 = (ρ j,k,1 ). (A.2)
Since all its entries are positive16 , and the sum of each row is equal to 1, this is a stochastic matrix. By Perron-Frobenius's theorem [see e.g. [START_REF] Seneta | Non-Negative Matrices and Markov Chains[END_REF]], it has a unique eigenvalue on the unit circle, with all other eigenvalues smaller than 1 in modulus.

Thus vector (y 1,1,t , y 2,1,t , ..., y J,1,t ) is co-integrated and the rank of co-integration is J -1.

Let us write the Jordan-Chevalley-Dunford decomposition of the autoregressive matrix R 1 :

R 1 = P -1 DP,
where P is a real invertible matrix, and D is triangular superior. Since by definition the sum of each row of R 1 = (ρ j,k,1 ) is equal to one, namely R 1 e = e, where e is the unitary vector, we can deduce that 1 is the eigenvalue of R 1 . Therefore, without loss of generality, we can assume that d 1,1 is equal to 1. Then we can rewrite (3.5) into:

P Y t = DP Y t-1 + m 1 + t .
Thus the 2, 3, .., J-th components of vector P Y t are all stationary. They constitute the linear span of co-integration vectors of Y t , which is of dimension J -1.

Let us now use the condition R 1 e = e, that is D(P e) = P e. In other words P e is a non zero (since P is invertible) eigenvector associated with the eigenvalue 1 of D. Since D is triangular superior, and 1 is its simple eigenvalue, we deduce that the 2, 3, .., J-th components of P e are all zero. Therefore, we have found J -1 co-integration vector, which are the 2, 3, .., J-th components of matrix P . They are linearly independent (since P is invertible), and share the common property that all the entries sum up to zero. Since the vector space of vectors whose components sum up to zero is of dimension J -1, we conclude that any vector with exactly one entry with value 1, one entry with value -1, and J -2 entries with value 0, belongs to this vector space. Thus each of the two mortality rates y j,1,t and y j ,1,t are co-integrated.

Let us now consider the mortality rates at all ages. First, we remark that the vector Z t = (y j,i,t ) 1≤j≤J,1≤i≤I forms a large VAR system with autoregressive matrix R 2 that is, by definition, row stochastic R 2 e 2 = e 2 , where e 2 is the I × J dimensional unitary vector.

Let us show that all eigenvalues of this large matrix is smaller than 1 in modulus.

To this end, let us write R 2 in blocks:

R 2 =          A 1 0 • • • • • • * A 2 • • • • • • * * A 3 • • • • • • • • • • • • • • •          . (A.3)
This is a matrix of blocks with I rows and I columns, corresponding to the I different ages, and each entry, such as A 1 , is itself a matrix of dimension J × J. Moreover, as a block matrix, R 2 is triangular superior, in the sense that in equation (A.3), all the J × J dimensional block matrices above the diagonal are zero. This is due to the fact that in our multi-population mortality specification, the mortality at age i depends only on past mortality at age i or lower.

We can also remark that the block matrix A 1 is nothing else but the matrix R 1 we considered in equation (A.2). In particular, A 1 has 1 as a simple eigenvalue, and all other eigenvalues are smaller than 1 in modulus.

In order to prove that all other eigenvalues of R 2 are smaller than 1 in modulus, we remark that since R 2 is (block) triangular superior, the characteristic polynomial det(R 2 -xId) can be factorized into:

det(R 2 -xId) = det(A 1 -xId)det(A 2 -xId) • • • det(A I -xId). (A.4)
Thus the spectrum of R 2 is the union of the spectrum of the I different J × J matrices A 1 , ..., A I . It can be easily shown that A 2 , ..., A I all have nonnegative entries only, and the sum of each row is strictly smaller than 117 . Hence they do not have eigenvalues with modulus smaller than 1. Therefore all the roots of det(A i -xId) are larger than unity in modulus, for i = 2, ..., I. As a consequence, we conclude that 1 is the single, simple eigenvalue with modulus non smaller than 1 of R 2 . Thus the vector (y j,i,t ) 1≤j≤J,1≤i≤I has (I × J -1) linearly independent co-integration relationship. Then by the same argument as for the vector of mortality at initial age (y j,1,t ) 1≤j≤J , we can deduce that there exists (I × J -1) co-integration vector whose sums of components are zero, and then, all I × J dimensional vector of the type (0, ..., 1, 0, ..., -1, 0, ...) are co-integration vectors.

B Appendix 2: Closed form formula of the PLS estimate PLS estimation of the causality matrix Let us solve the minimization problem (4.4). The minimiser is a solution of the following equations:

0 = ∂L 1 ∂α i+1 = -2 T -1 t=1 (y i,t -y i+1,t ) y i+1,t+1 -y i+1,t -α i+1 (y i,t -y i+1,t ) -β i+1 (y i-1,t -y i+1,t ) -m i+1 + 2λ α 2α i+1 -α i -α i+2 , 0 = ∂L 1 ∂β i+1 = -2 T -1 t=1 (y i-1,t -y i+1,t ) y i+1,t+1 -y i+1,t -α i+1 (y i,t -y i+1,t ) -β i+1 (y i-1,t -y i+1,t ) -m i+1 + 2λ β 2β i+1 -β i -β i+2 , 0 = ∂L 1 ∂m i+1 = -2 T -1 t=1 y i+1,t+1 -y i+1,t -α i+1 (y i,t -y i+1,t ) -β i+1 (y i-1,t -y i+1,t ) -m i+1 + 2λ m 2m i+1 -m i -m i+2 ,
for i = 2, ..., I -1, and for the boundary ages:

0 = ∂L 1 ∂α 2 = -2 T -1 t=1 (y 1,t -y 2,t ) y 2,t+1 -y 2,t -α 2 (y 1,t -y 2,t ) -m 2 , 0 = ∂L 1 ∂m 2 = -2 T -1 t=1 y 2,t+1 -y 2,t -α 2 (y 1,t -y 2,t ) -m 2 , 0 = ∂L 1 ∂m 1 = -2 T -1 t=1 y I,t+1 -y I,t -m 1 .
Thus parameters m 1 , m 2 , and α 2 are estimated directly from the linear system composed of the three boundary conditions. This is equivalent to solve the OLS for the regression models (2.3) and (2.4).

The rest of the parameters satisfy the linear system in (α i , β i , m i ) i=2,...,I-1 . This system has 3 × (I -2) unknown parameters and as many equations. Thus the solution of this linear system is unique, whenever the determinant of the associated matrix is non null. But this determinant is polynomial in λ α , λ β and λ m , and when λ m = λ α = λ β = 0, the solution of the system is unique (and is, by definition, the OLS estimator). Therefore, the determinant is non null when λ m = λ α = λ β = 0. As a consequence, the determinant is non null for almost all λ α , λ β , λ m ,18 and the unique minimizer of the optimization problem can be computed in closed form, by solving the linear system.

PLS estimation of the variance-covariance matrix Let us now solve the minimization problem (4.7). The minimizer satisfies the following linear system:

0 = ∂L 2 ∂a i = -2 T t=1 i+1,t i,t -a i i+1,t -c i i-1,t + 2λ a 2a i -a i-1 -a i+1 , 0 = ∂L 2 ∂c i = -2 T t=1 i-1,t i,t -a i i+1,t -c i i-1,t + 2λ c 2c i -c i-1 -c i+1 , (B.1)
for all ages i = 2, ..., I -1, as well as:

0 = ∂L 2 ∂a 1 = -2 T t=1 i+1,t ( 1,t -a 1 2,t ) 0 = ∂L 2 ∂c I = -2 T t=1 i+1,t ( I,t -c I I-1,t )
for the two boundary ages. As in the previous subsection, parameters a 1 and c I are estimated using OLS, whereas the rest of the parameters (a i , c i ) i=2,...,I-1 are obtained by solving a linear system (B.1). By construction, the determinant of the latter is polynomial in λ a and λ c , and is non null when λ a = λ c = 0. Therefore, it is non null for almost surely all values λ a and λ c .

The two-population case Let us now solve the minimization problem (4.11). The minimizer satisfies the following linear system:

0 = ∂L ∂ αj,i+1 = -2
T -1 t=1 (y j,i,t -y j,i+1,t ) y j,i+1,t+1 -y j,i+1,t -αj,i+1 (y j,i,t -y j,i+1,t )

-(1 -ρ j,i+1 )(y -j,i+1,t -y j,i+1,t ) -m j,i+1 + 2λ j,α 2 αj,i+1 -αj,iαj,i+2 ,

0 = ∂L ∂m i+1 = -2
T -1 t=1 y j,i+1,t+1 -y j,i+1,t -αj,i+1 (y j,i,t -y j,i+1,t )

-(1 -ρ j,i+1 )(y -j,i+1,t -y j,i+1,t ) -m j,i+1 + 2λ j,m 2m j,i+1 -m j,i -m j,i+2 , 0 = ∂L ∂ρ j,i+1

= -2

T -1 t=1 (y j,i+1,t -y -j,i+1,t ) y j,i+1,t+1 -y j,i+1,t -αj,i+1 (y j,i,t -y j,i+1,t )

-(1 -ρ j,i+1 )(y -j,i+1,t -y j,i+1,t ) -m j,i+1 + 2λ j,β 2ρ j,i+1 -ρ j,i -ρ j,i+2 , for j = 1, 2 and i = 1, ..., I -1. For age 1, the linear system is given by

0 = ∂L ∂m j,1 = -2
T -1 t=1 y j,1,t+1 -y j,1,t -(1 -ρ j,1 )(y -j,1,t -y j,1,t ) -m j,1 , 0 = ∂L ∂ρ j,1 = -2

T -1 t=1 (y j,1,t -y -j,1,t ) y j,1,t+1 -y j,1,t -(1 -ρ j,1 )(y -j,1,t -y j,1,t ) -m j,1 , for j = 1, 2.

D Appendix 4: Granger and instantaneous causality

The aim of this appendix is to review the notions of causality introduced in the (discrete time) literature. The first concept is the directional causality, or Granger's causality, which has been first introduced by [START_REF] Granger | Investigating Causal Relations by Econometric Models and Cross-Spectral Methods[END_REF].

Definition 1 (Granger's (non)causality). Let F = (F t ) and Z = (Z t ) be two discrete time processes, then the process F does not Granger cause the process Z if and only if:

Z t ⊥ ⊥ F t-1 | Z t-1 , ∀t > 0, (D.1)
where the symbol • ⊥ ⊥ • | • means conditional independence.

Similarly, if the noncausality condition is not satisfied, then we say there is Granger causality from F towards Z.

Thus the process F does not Granger cause the process Z, if, at any date, the current value of Z t is independent of the lagged values of F given the lagged values of Z. Equivalently the best predictor of any nonlinear transformation g(Z t ) of Z t given the whole information set (F t-1 , Z t-1 ) does not depend on the past history of F .

In the context of a VAR model of (F t , Z t ), say:

   F t Z t    =   

(

  2011); Li and Hardy (2011); Yang and Wang (2013); Li et al. (2013); Hyndman et al. (2013); Zhou et al. (2014); Hunt and Blake (2015b); Salhi and Loisel (2016)]. Given

Figure 2 :

 2 Figure 2: Simulated trajectory of y 1 (in black full line), y 2 (in blue dashed line), y 3 (in red dotted line), when α 2 = α 3 = 1, and β 2 = β 3 = 0, and the residuals are independent.

Figure 3 :

 3 Figure 3: Simulated trajectory of y 1 (in black full line), y 2 (in blue dashed line), y 3 (in red dotted line), when α 2 = α 3 = 1, and β 2 = β 3 = 0, and the residuals are identical.

Figure 4 :

 4 Figure 4: Simulated trajectory of y 1 (in black full line), y 2 (in blue dashed line), y 3 (in red dotted line), when α 2 = α 3 = β 2 = β 3 = 0.05, and the residuals are independent.

Figure 5 :

 5 Figure 5: Simulated trajectory of y 1 (in black full line), y 2 (in blue dashed line), y 3 (in red dotted line), when α 2 = α 3 = β 2 = β 3 = 0.05, and the residuals are identical.

Figure 6 :Figure 7 :

 67 Figure 6: The RMSFE of the Lee-Carter and STAR model applied to two populations separately.

Figure 8 :

 8 Figure 8: The out-of-sample forecasts of the average of y i,t from the STAR model and the Lee-Carter model.

Figure 9 :

 9 Figure9: The projected life expectancies at birth of the US and the UK population. The solid lines are the historical life expectancies and the mean forecasts from the twopopulation STAR model; the dashed lines are the mean forecasts of the single-population STAR model; the dotted lines are the 2.5% and the 97.5% quantiles of the two-population STAR model; dotted-dashed lines are the 2.5% and the 97.5% quantiles of the singlepopulation STAR model.

Figure 10 :

 10 Figure 10: The projected logarithm of crude death rates from the Lee-Carter model and the STAR model.

Figure 12 :Figure 13 :Figure 14 :

 121314 Figure 12: The smoothed and the original estimates for the UK population from the single-population model.

  Propagation of the shock i,t during the subsequent periods. At time t + 1, ages i, i + 1 and i + 2 are impacted; at time t + 2, ages i to i + 4 are impacted.

	ments on the principal diagonal, as well as the first two lower sub-diagonals.
	Such a constrained VAR can be regarded as a spatial-temporal autoregressive model
	[see e.g. Pace et al. (1998)], in the sense that different component processes (y i,t ) t
	are smoothed not only on the time domain, but also on the "space", that is the age
	domain. As a consequence, we expect that for a given t, components (y i,t ) and (y j,t ) are
	closed so long as i and j are closed. Indeed, the shock i,t enters, at time t + 1, in the
	expressions of y i,t+1 , y i+1,t+1 , y i+2,t+1 ; and enters, at time t + 2, in the expressions of
	y i,t+1 , y i+1,t+1 , y i+2,t+2 , y i+3,t+2 and y i+4,t+2 . More generally, it takes h years for a shock
	to reach out progressively to 2h + 1 different cohorts. As a consequence, although our
	VAR model is of first-order, it allows for the effect of the shocks to persist over several
	periods. This property is not satisfied by standard factor models, and explains the larger
	flexibility of our model with respect to the latter. Figure 1 illustrates the propagation
	of a mortality shock i,t between different cohorts.
	i
	i + 4			(i + 4, t + 2)
	i + 3			(i + 3, t + 2)
	i + 2			(i + 2, t + 2)
	i + 1			(i + 1, t + 2)
	i	(i, t)		(i, t + 2)
		t	t + 1	t + 2	t
	Figure 1:
					.6)
	where the (Granger causality) matrix R is lower triangular, and only has non null ele-

Table 2 :

 2

		06 0.07 2.39 0.19 0.18 0.07
	UK	0.47 0.34 5.21 0.50 0.86 0.32

  1t , 2t ) is an i.i.d. sequence of residuals, the Granger causality is equivalent to the following:Definition 2 (Granger's causality in a VAR model). Process (F t ) (resp. (Z t )) Granger causes (Z t ) (resp. (F t )) if and only if the coefficient r 21 (resp. r 12 ) is non zero.Thus this definition is directional. We can have situations where (F t ) Grange causes (Z t ), but not the converse. Moreover, this causality measures the cross-dependence between different periods. The next concept of causality is not directional, and concerns the instantaneous dependence between different time series.

							
	r 11 r 12 r 21 r 22	 	  Z t-1 F t-1	  +	 	2t 1t	  ,	(D.2)
	where (

Source: International Monetary Fund Global Financial Stability Report, April

Source: Population Trends N.145, Autumn 2011, Office of National Statistics.

See also[START_REF] Cairns | A Quantitative Comparison of Stochastic Mortality Models Using Data from England and Wales and the United States[END_REF],Plat (2009), etc.

The terminology "simultaneous equation" is in the sense that t's appear on both sides of the equations.

7 However, in general, even when |i -j| = 1 or 0, the correlation between i,t and j,t is non zero. Indeed, intuitively, i,t and its second-order neighbouring i+2,t are dependent since both terms are dependent with i+1,t

In the Lee-Carter model, in the first stage, parameters can be estimated via a single value decomposition.

We thank an anonymous referee for this insightful suggestion.

The estimate can be computed in closed form, in the same way as the estimate from the whole sample.

See: http://www.mortality.org/

The case of analysis of square matrix with nonnegative entries is similar to the case of positive matrix, see e.g.[START_REF] Seneta | Non-Negative Matrices and Markov Chains[END_REF]

This is because the mortality rates at higher ages i ≥ 2 depend not only on past mortality rates of other populations at the same age i, but also past mortality rates at an inferior age i -1.

In the sense of Lebesgue measure.

Similarly as the single population case, we have the following property concerning the co-movement of mortality of the two populations: Proposition 2. The difference of any two components of the (2×I)-dimensional process (y 1,i,t , y 2,i,t , i = 1, 2, ..., I) is stationary.

Proof. See Appendix 1.1. Therefore, our specification ensures both within-, and inter-population co-integration. At the same time, the dependence of parameters ρ 1,i , ρ 2,i on age i provides flexibility when it comes to capturing potential short-term mortality deviation at certain ages between the two populations. This necessity of reconciling the desirable long-term property and the empirical data is an important criteria for mortality models [see e.g. Hunt and Blake (2015a) for a discussion]. Our model provides such a trade-off.

The specification of the residuals is also extended to allow for inter-population dependence. In particular, let 1,i,t and 2,i,t be the residuals from population 1 and 2 for age i and time t, the dynamics of the residual is given by: j,I,t = θ j,I c j,I j,I-1,t + (1 -θ j,I ) -j,I,t + η j,I,t j,I-1,t = θ j,I-1 a j,I-1 j,I,t + c j,I-1 j,I-2,t + (1 -θ j,I-1 ) -j,I-1,t + η j,I-1,t j,I-2,t = θ j,I-2 a j,I-2 j,I-1,t + c j,I-2 j,I-3,t + (1 -θ j,I-2 ) -j,I-2,t + η j,I-2,t (3.3)

A J-population extension with J ≥ 3. Let us denote by y j,i,t the mortality of population j, j = 1, ..., J, at age x and time t. Then a natural idea is to generalize the two-population model [see equation (3.1)] into:

(3.4) populations simultaneously.

In the empirical study, we have fitted both the single-population and the twopopulation STAR model to the US and the UK population. We found that the STAR models are able to generate reasonable mortality forecasts in the long-run. Moreover, in the out-of-sample forecast, the STAR models outperforms the benchmark models (the [START_REF] Lee | Modeling and Forecasting US Mortality[END_REF] model in the single-population case, and the [START_REF] Li | Coherent Mortality Forecasts for a Group of Populations: An Extension of the Lee-Carter Method[END_REF] model in the two-population case).

A Appendix 1: Proofs

A.1 Proof of Proposition 2

First, from Equation (3.2) we get:

Next, by equations (3.1) and (3.2), we have y j,2,t+1 -y j,1,t+1 = ρ j,2 (1 -α j,2 )(y j,2,t -y j,1,t ) + (1 -ρ j,2 )(y -j,2,t -y -j,1,t ) + (ρ j,2 -ρ j,1 )(y j,1,t -y -j,1,t ) + m j,2 -m j,1 + j,2,tj,1,t stationary , (A.1) for j = 1, 2. From (A.1), we get a VAR model for the joint process:

Thus the stationary condition is that polynomial det(A -zId) has no roots outside the Definition 3 (Instantaneous causality in a VAR model). We say that there is instantaneous causality between processes (F t ) and (Z t ) if and only if the residuals 1t , 2t are not jointly independent.

For instance, if ( 1t , 2t ) is jointly normal, then there is instantaneous causality if and only if the variance-covariance matrix of the residuals is not diagonal.