Hong Li

Yang Lu
email: yang.lu@univ-paris13.fr

A Bayesian non-parametric model for small population mortality

Keywords: Autoregressive gamma, credibility, parameter uncertainty, two-population mortality model

This paper proposes a Bayesian non-parametric mortality model for a small population, when a benchmark mortality table is also available and serves as part of the prior information. In particular, we extend the Poisson-gamma model of Hardy and Panjer to incorporate correlated and age-specific mortality coefficients. These coefficients, which measure the difference in mortality levels between the small and the benchmark population, follow an age-indexed autoregressive gamma process, and can be stochastically extrapolated to ages where the small population has no historical exposure. Our model substantially improves the computation efficiency of existing two-population Bayesian mortality models by allowing for closed form posterior mean and variance of the future number of deaths, and an efficient sampling algorithm for the entire posterior distribution. We illustrate the proposed model with a life insurance portfolio from a French insurance company.

Introduction

Projecting mortality is of great concern for many disciplines, such as demographic analysis, social planning, policy-making, and actuarial practices. In recent decades, various methodologies for mortality modeling and projections have been proposed in the literature, such as the Lee-Carter model (1992), the Cairns-Blake- [START_REF] Cairns | A two-factor model for stochastic mortality with parameter uncertainty: theory and calibration[END_REF] model, and their various extensions (see e.g. [START_REF] Booth | Mortality modelling and forecasting: a review of methods[END_REF]Tickle 2008, Cairns et al. 2011a). Moreover, government agencies and societies of actuaries publish national or industry-level regulatory mortality tables on a regular basis. However, in many situations, the population of interest, such as the portfolios of pension funds or insurance companies, have limited size and/or short observation window. For such portfolios, the observed mortality experience is typically erratic, with hard-to-identify, or even unreliable longevity trend.

However, such small populations are often sub-populations of a larger population, such as a national population or a larger collective portfolio. The mortality pattern of the latter is usually much more regular, and easier to forecast. In this case, mortality table of the larger population may serve as a benchmark when modeling and forecasting mortality experience of the small population.

To this end, different solutions have been proposed in the literature. The first branch of literature is based on deterministic approaches, such as relational models (see e.g. [START_REF] Himes | A relational model of mortality at older ages in low mortality countries[END_REF][START_REF] Brouhns | A Poisson log-bilinear regression approach to the construction of projected life tables[END_REF][START_REF] Tomas | Prospective mortality tables: taking heterogeneity into account[END_REF], as well as methods based on the distortion of probability distributions (see e.g. Denuit et al. 2007, Bienvenüe and[START_REF] Bienvenüe | Iterative adjustment of survival functions by composed probability distortions[END_REF]. These models specify a deterministic, linear, or non-linear functional relationship between the mortality curves of the two populations. Their drawbacks are that, first, the functional relationship typically involves a small number of parameters,1 which limits the flexibility of the model; second, and more importantly, the deterministic relationship does not capture the impact of parameter uncertainty on pricing/reserving/capital requirement calculation, which could be substantial for small populations.

The second class of models are Bayesian (see e.g. [START_REF] Dowd | A gravity model of mortality rates for two related populations[END_REF], Cairns et al. 2011b[START_REF] Antonio | Bayesian Poisson log-bilinear models for mortality projections with multiple populations[END_REF][START_REF] Van Berkum | A Bayesian joint model for population and portfolio-specific mortality[END_REF]. While they are flexible, and account for parameter uncertainty in a natural way, these models are currently rather computationally intensive. While credibility theory has been proposed to simplify such models (see e.g. [START_REF] Hardy | A credibility approach to mortality risk[END_REF][START_REF] Olivieri | Life tables in actuarial models: from the deterministic setting to a Bayesian approach[END_REF][START_REF] Salhi | A credibility approach of the Makeham mortality law[END_REF], the computational gain usually relies on strong parametric restrictions.

This paper proposes a Bayesian non-parametric model for the mortality of a small population, when mortality table of a larger population is exogenously given as the benchmark. The proposed model has a flexible structure, while at the same time being less computationally intensive than the existing Bayesian multi-population models. In the proposed model, the small population is assumed to follow the same longevity trend as the benchmark population, but has different, age-specific mortality levels. These levels are captured by a set of correlated and age-specific coefficients. A priori, the mortality level of the small population is the same as the benchmark population for all ages, and its forecasts depend solely on that of the latter. As mortality data accumulate, the estimation of the mortality coefficients are updated by the Bayesian rule, and the small population has increasing weights when forecasting its own mortality.

From the computational side, the model gives closed form expressions of the posterior moments, such as mean and variance, of the next period's total death count of the small population. Moreover, in an insurance context, where the observations are claims with heterogeneous payments, closed form expressions are available for moments of the total claim amount as well. These closed form expressions significantly reduce the computational requirement of our model compared to existing Bayesian models. Therefore, we are able to investigate, with limited computational cost, the impact of the prior distribution on the expected future death counts and claim amount and the corresponding uncertainty. These parameters include, among others, the size of the population of interest, the length of the observation window, and the choice of the prior distribution. From the flexibility side, our approach is non-parametric, in the sense that the mortality coefficients, which measure the difference of mortality level, are age-specific. Moreover, the mortality coefficients are correlated across different ages, with stronger correlation when two ages are closer. Such correlation leads to smoother projected mortality curves. Compared to conventional non-parametric models such as [START_REF] Olivieri | Life tables in actuarial models: from the deterministic setting to a Bayesian approach[END_REF], the Bayesian aspect of the model allows us to conduct stochastic extrapolation of the mortality coefficients to ages with no historical risk exposure.

The proposed model has also implications for regulation and risk management. Firstly, the closed form expressions allow us to conveniently decompose the total variance into two parts: the diversifiable part due to sampling variation, and the non-diversifiable part due to parameter uncertainty. In the empirical analysis, it is shown that, for a mildly large life insurance portfolio with around 15,000 policyholders and a fairly long observation period (11 years), the contribution of the parameter uncertainty to the total variance can still be rather substantial (about 12%). This result illustrates the danger of under-estimating the risk if a deterministic relational model is used. Secondly, our model may be suitable for building an internal model under the new insurance regulations, such as Solvency II. These regulations encourage the use of proprietary mortality tables that better reflect company-specific risks. Nevertheless, their implementation has to be monitored carefully, in order to prevent insurers from regulatory arbitrage. When a regulatory life table (often built with some risk margin on top of the 'best estimate') serves as the benchmark, our model leads to projections that are closer to the insurer's own mortality experience when its mortality experience is more reliable (e.g. a larger portfolio size with a longer observation window).

In a numerical illustration, we apply the non-parametric Bayesian model to a French life insurance portfolio. In particular, we compute the one-year-ahead predictive expectation and variance of the total death counts, and illustrate the impact of the portfolio size and the choice of prior distributions on the risk management of the portfolio. We illustrate, for a given prior distribution of the mortality coefficients, how a larger portfolio gains more weight on its own mortality experience when projecting future mortality, whereas a smaller portfolio relies more heavily on the benchmark mortality table. Besides, we also discuss the joint impact of the portfolio size and the prior distribution, which measures the degree of belief on the benchmark mortality table.

The paper is organized as follows. Section 2 reviews the Poisson-gamma model of [START_REF] Hardy | A credibility approach to mortality risk[END_REF]. Section 3 discusses the drawbacks of this benchmark approach, and introduces the Bayesian non-parametric model. Section 4 provides the Bayesian updating formula for the mortality forecasts, in particular the closed form expression of the conditional expectation and variance of the total death count (claim amount). Section 5 illustrates the model with a real life insurance portfolio. Section 6 concludes. Technical proofs are gathered in Appendices.

The benchmark Poisson-gamma model

The Poisson-gamma credibility model of [START_REF] Hardy | A credibility approach to mortality risk[END_REF] is based on the conjugacy property of the gamma distribution with respect to the Poisson likelihood function, inspired by studies of nonlife insurance portolios (see e.g. [START_REF] Dionne | A generalization of automobile insurance rating models: the negative binomial distribution with a regression component[END_REF]. Assume that an insurance portfolio (the small population) contains mortality experience of I ages and T years. Ages are indexed by 1, 2, …, I, with indices '1' and 'I' referring to the smallest and the largest ages in the portfolio, rather than the biological age 1 and I. The number of death at each age x and year t is conditionally Poisson distributed (see e.g. [START_REF] Czado | Bayesian Poisson log-bilinear mortality projections[END_REF]:

d x,t ∼ P(θ x μ 0,x,t e x,t ), ∀x = 1, . . . , I, t = 1, . . . , T, (2.1)
where

• e x,t is the exposure of the portfolio at age x and year t.

• μ 0,x,t is the force of mortality of a benchmark population at the same age and year. This can be, for instance, mortality of the corresponding national population, a collective portfolio, or a regulatory mortality table. In general, mortality of this population can be stochastic with a longevity trend, i.e. mortality rates at each age are non-stationary, and decrease over time. μ 0,x,t can be estimated from mortality data of the benchmark population, or exogenously given, e.g. by regulators. • the age-specific mortality coefficient, θ = (θ 1 , . . . , θ I ) , are time-invariant. They capture the relative risk of the portfolio with respect to the benchmark population. For instance, if θ x is smaller than 1, then the insurance portfolio has lower expected mortality at age x than the benchmark. The mortality coefficients are unobserved and stochastic, inducing uncertainty on the portfolio's mortality.

The assumption that θ is constant means that the mortality ratios between the two populations are stable. In other words, the small population shares the same longevity trend with the benchmark population. This assumption is also used in the subsequent literature (see e.g. [START_REF] Olivieri | Life tables in actuarial models: from the deterministic setting to a Bayesian approach[END_REF][START_REF] Salhi | A credibility approach of the Makeham mortality law[END_REF][START_REF] Van Berkum | A Bayesian joint model for population and portfolio-specific mortality[END_REF]. This assumption is motivated by the following arguments:

• It allows the mortality improvements in these two populations to be coherent in the long-run, which is an important criterion for modeling mortality of closely related populations. 2 • When the small population has a reduced observation window, which is often the case in insurance application, it is extremely difficult to reliably derive the mortality trend of the small population based on its own data. Therefore, it is more reasonable to assume that the two populations share the same mortality trends.

Let us now specify the prior distribution of θ. Model (2.1) is first proposed by [START_REF] Hardy | A credibility approach to mortality risk[END_REF], who assume that θ x is common for all ages. Specifically, θ x = θ for all x's, where θ has the prior distribution γ (ν, ν). In this way, the prior mean of θ equals to 1, and thus E[d x,t ] = μ 0,x,t e x,t . In other words, the portfolio is expected to have a priori the same mortality experience as the benchmark population. As is shown by [START_REF] Salhi | A credibility approach of the Makeham mortality law[END_REF], this condition is often not satisfied in practice. In particular, the common θ assumption may lead to over-estimation of d x,t at certain ages, and under-estimation at some others.

As an alternative to the common θ assumption, [START_REF] Olivieri | Life tables in actuarial models: from the deterministic setting to a Bayesian approach[END_REF] and [START_REF] Van Berkum | A Bayesian joint model for population and portfolio-specific mortality[END_REF] assume that θ is a vector of i.i.d. random variables, each with the prior distribution γ (ν, ν). This assumption is Bayesian non-parametric, in the sense that the number of parameters vary with the number of ages included in the analysis. In the sequel, we refer to the model with i.i.d. θ assumption as the benchmark Poisson-gamma model.

As observations accumulate over time, the conditional distribution of each θ x is regularly updated using Bayes' rule. For each T, denote by D x = (d x,1 , . . . , d x,T , e x,1 , . . . , e x,T , μ 0,x,1 , . . . , μ 0,x,T ) the set of death counts and exposure of the portfolio and the force of mortality of the benchmark population at age x, the conditional density of D x given θ x is:

l(D x |θ x ) = T t=1 l(d x,t |θ x ) = T t=1 (θ x e x,t μ 0,x,t ) d x,t exp ( -θ x e x,t μ 0,x,t ) (d x,t + 1) ∝ θ T t=1 d x,t x exp -θ x T t=1 e x,t μ 0,x,t .
(2.2)

The following closed form expression of the posterior distribution is the key contribution to the popularity of the benchmark Poisson-gamma model. Using Equation (2.2), the conditional distribution of θ x given the death counts is:

l(θ x | D x ) ∝ θ ν-1+ T t=1 d x,t x exp -θ x ν + T t=1 e x,t μ 0,x,t .
(2.3) That is, l(θ x | D x ) follows the gamma distribution γ (ν + T t=1 d x,t , ν + T t=1 μ 0,x,t e x,t ). Thus, the predicted mean of next year's death count is: (2.4) Therefore, the posterior expected number of deaths can be interpreted as the prior expected number, e x,T+1 μ 0,x,T+1 , multiplied by an average of 1 and the ratio between the observed and the expected number of deaths,

E[d x,T+1 |D x ] = μ 0,x,T+1 e x,T+1 ν + T t=1 d x,t ν + T t=1 μ 0,x,
T t=1 d x,t T t=1 μ 0,x,t e x,t
. In this model, parameter ν captures the belief of the benchmark mortality table: if ν is much larger than the cumulative sum, T t=1 μ 0,x,t e x,t , then we have:

E[d x,T+1 |D x ] ≈ e x,T+1 μ 0,x,T+1 .
In other words, we have strong belief that mortality of the small population would be close to that of the benchmark population. Contrarily, if ν is small, we have: .

E[d x,T+1 |D x ] ≈ e x,
In this case, the posterior expected death count depends on the prior expectation, e x,T+1 μ 0,x,T+1 , multiplied by an adjustment factor,

T t=1 d x,t T t=1 μ 0,x,t e x,t
, which reflects the difference between the observed and the expected number of deaths. In other words, predictive death counts in the small population are heavily driven by its own historical pattern.

Let us now compare the variance of the posterior distribution with that of the prior distribution. We have:

V[θ x |D x ] = ν + T t=1 d x,t (ν + T t=1 μ 0,x,t e x,t ) 2 , V[θ x ] = 1 ν . (2.5) Therefore, we have V[θ x ] > V[θ x |D x ] if and only if: ν T t=1 d x,t -2 T t=1 μ 0,x,t e x,t < T t=1 μ 0,x,t e x,t 2 .
A sufficient condition of this latter inequality is:

T t=1 d x,t < 2 T t=1
μ 0,x,t e x,t .

(2.6)

When the mortality of the small population is lower, or close to that of the benchmark population, Equation (2.6) holds with a large probability. Hence, the posterior variance V[θ x |D x ] is generically smaller than the prior variance V[θ x ]. In other words, information updating lowers the uncertainty of the conditional distribution of l(θ x |D x ).

Let us now discuss the difference of the Poisson-gamma model with the relational models, such as the Brass model (see e.g. [START_REF] Brouhns | A Poisson log-bilinear regression approach to the construction of projected life tables[END_REF]. The latter model assumes that the force of mortality of a population of interest, μ 1 , is related to that of a benchmark population, μ 2 , via an affine relationship:

f (μ 1 (x, t)) = θ 1 f (μ 2 (x, t)) + θ 2 ,
where f is some link function, and θ 1 and θ 2 are typically deterministic and constant across different ages. The Poisson-gamma model can be viewed as having a similar spirit, in the sense that f (x) = x, θ 2 = 0. However, θ x is age-specific and stochastic, rather than deterministic and common across ages.

The autoregressive gamma-Poisson model

The benchmark Poisson-gamma model assumes independence between the mortality coefficients θ . As a consequence, the posterior distribution of θ x depends only on the small population's mortality experience at age x, and the resulting curve of E[θ x |D x ] can be rather erratic. In particular, at ages with small (or no) exposure, E[θ x |D x ] is close to the prior mean, 1; whereas at ages with a large exposure, E[θ x |D x ] could be rather different than 1.

In this paper, we relax the independence assumption on the mortality coefficients by introducing the autoregressive gamma (ARG) process. This allows the posterior distribution of each θ x to depend not only on mortality experience at age x, but the whole sample. Moreover, the ARG process leads to stochastic and smooth extrapolation of the portfolio mortality table to ages with no historical experience, rather than a point estimate given by the deterministic extrapolation methods (see e.g. [START_REF] Salhi | A credibility approach of the Makeham mortality law[END_REF].

The ARG process is the exact time discretization of the continuous time Cox-Ingersoll-Ross process (see [START_REF] Nieto-Barajas | Markov beta and gamma processes for modelling hazard rates[END_REF]Walker 2002, Gouriéroux and[START_REF] Gouriéroux | Autoregressive gamma processes[END_REF]. It has gamma marginal distribution, and nests the benchmark Poisson-gamma credibility model. To define the ARG process, we can start from either end of the age group, i.e. θ 1 or θ I , and construct the process iteratively with an auxiliary counting process N = (N 1 , . . . , N I ) . For example, the construction procedure starting from θ 1 is as follows:

• θ 1 follows a gamma distribution, γ (ν, ν).

• Given θ 1 , introduce a latent count variable, N 1 , which follows a Poisson conditional distribution with parameter βθ 1 . Given (N 1 , θ 1 ), θ 2 then follows the conditional gamma distribution γ (ν + N 1 , β + ν). • The next latent count variable, N 2 , follows the conditional Poisson distribution P(βθ 2 ), and θ 3 follows γ (ν + N 2 , ν + β), etc.

To summarize, we have the following causal chain:

θ 1 → N 1 → θ 2 → N 2 → θ 3 • • • (3.1)
This is a Markov chain, in the sense that the conditional distribution of each variable given variables on its left hand side depends only on its nearest left neighbour. For instance,

l(N 2 |θ 2 , N 1 , θ 1 ) = l(N 2 |θ 2 ), l(θ 2 |N 1 , θ 1 ) = l(θ 2 |N 1 ).
Moreover, N and θ are both Markov processes with respect to their own history. For instance, the conditional distribution of N 3 given N 1 , N 2 only depends on N 2 . Alternatively, we can define the joint distribution of θ from the oldest age x = I, and construct the process in the reverse direction:

• • • , θ I-2 ← N I-2 ← θ I-1 ← N I-1 ← θ I ,
where N I-1 follows P(βθ I ), and θ I-1 follows γ (ν + N I-1 , ν + β), etc. It will become clear in Proposition 1 that process θ has the time reversibility property. Intuitively, this property means that these two methods of construction lead to the same joint distribution of θ, and hence the definition of its dynamics does not depend on the choice of the direction.

This specification has a similar spirit to the lognormal autoregressive specification proposed by van [START_REF] Van Berkum | A Bayesian joint model for population and portfolio-specific mortality[END_REF], in the sense that both models allow correlation between different θ x 's. However, under the log-normal model, moments of the posterior distribution need to be computed using simulations.

It is also worth mentioning that the counting variables N x 's are latent and have no direct link with the observed mortality counts d x,t 's. They can be interpreted as (age-dependent) regimes, since the larger N x , the larger θ x in expectation. Although the ARG process can also be introduced without these count variables (see e.g. Lu Forthcoming) for an introduction via the conditional Laplace transform], they are crucial to the closed form posterior moment expressions. Now we summarize some properties of the ARG process that are useful for our application.

Proposition 1 (see Gouriéroux and Jasiak 2006):

(1) (Marginal distribution) For each age x, variable θ x is marginally γ (ν, ν) distributed.

(2) (Serial correlation) For each h, we have:

E[θ x+h |θ x ] = 1 -ρ h + ρ h θ x , (3.2) V[θ x+h |θ x ] = (1 -ρ h ) 2 ν + 2 ρ h (1 -ρ h ) ν θ x . (3.3)
As a consequence, we have

corr[θ x , θ x+h ] = ρ h ,
where ρ := β β+ν . When β is equal to zero, θ 1 , θ 2 ,…,θ x are independent, and we get the benchmark Poisson-gamma model in Section 2 as a special case.

(3) (Reversibility) The dynamics of the process θ is reversible, that is, in the reverse direction, the process θ I , θ I-1 , …, θ 1 has the same Markov, ARG dynamics. As a consequence, we have also:

E[θ x-h |θ x ] = 1 -ρ h + ρ h θ x , ∀h ∈ N. (3.4)
The proof of these properties can be found in [START_REF] Gouriéroux | Autoregressive gamma processes[END_REF], and is thus omitted. The first property on marginal distribution is due to the fact that the gamma distribution is a conjugate prior to the Poisson likelihood function. This conjugacy property greatly simplifies the updating formula and explains why gamma distribution is the invariant distribution of the process. It also implies that the proposed model is a direct generalization of the benchmark model. Indeed, if we condition only on the death counts at age x, we still have the linear prediction formula:

E[d x,T+1 |D x ] = e x,T+1 μ 0,x,T+1 ν + T t=1 d x,t ν + T t=1 μ 0,x,t e x,t
.

The second property shows that the auto-correlation of the process θ is the same as an AR(1) process. This explains the terminology ARG. The larger ρ, the smoother the (prior) trajectory of θ , and as a consequence, the smoother its posterior trajectory. Moreover, the dependence structure between the different risk parameters θ x 's is spatial, in the sense that the closer the two ages x 1 and x 2 , the stronger the correlation between θ x 1 and θ x 2 .

The conditional expectation formula (3.2) allows us to extrapolate the mortality table to ages older than I. Indeed, the expected value of θ I+h given observations at all ages is:

E[θ I+h |D I ] = E E[θ I+h |θ I , D I ] | D I = E E[θ I+h |θ I ] | D I = 1 -ρ h + ρ h E[θ I |D I ],
(3.5)

where

D x = (D 1 , D 2 , . . . , D x ).
In particular, D I is the whole dataset of death counts and exposure of the portfolio, and the corresponding death probabilities of the benchmark population. The extrapolation method introduced above is smooth, while being able to take into account uncertainties in θ x . First, for a fixed h, if ρ is close to 0, then the RHS of Equation (3.5) is close to 1, i.e. the projected mortality at the advanced ages are nearly not affected by observed mortality of the portfolio. On the other hand, if ρ is close to 1, then this term is close to E[θ I |D I ]. Therefore, with appropriate choice of ρ, the correlation structure of θ ensures that the projected mortality in each year vary smoothly across ages, rather than having a jump at age I + 1. Second, compared to existing deterministic extrapolation methods, this method takes into account the uncertainty of the extrapolated mortality coefficients. For instance, by Equation (3.3), the variance of the extrapolation at age I + h is:

V[θ I+h

|D I ] = V E[θ I+h |θ I ]|D I + E V[θ I+h |θ I ]|D I = ρ 2h V[θ I |D I ] + (1 -ρ h ) 2 ν + 2 ρ h (1 -ρ h ) ν E[θ I |D I ]. (3.6)
Hence, when h = 0, the RHS in Equation (3.6) is equal to V[θ I |D I ]; whereas when the extrapolation horizon h is large, the RHS converges to

1 ν = V[θ I ].
Finally, the reversibility property of the ARG process allows us to extrapolate the mortality of the portfolio to lower ages without historical data.

The Bayesian updating algorithm

After specifying the prior distribution of θ , it remains to compute their joint posterior distribution given all historical data, l(θ|D I ). The introduction of the serial correlation between θ x 's makes this task more difficult than the benchmark case. In particular, the posterior distribution is no longer gamma for each x, as opposed to Equation (2.3). This problem is usually tackled using computationally intensive MCMC techniques (see e.g. [START_REF] Pitt | Filtering via simulation: auxiliary particle filters[END_REF][START_REF] Czado | Bayesian Poisson log-bilinear mortality projections[END_REF][START_REF] Pedroza | A Bayesian forecasting model: predicting US male mortality[END_REF][START_REF] Antonio | Bayesian Poisson log-bilinear models for mortality projections with multiple populations[END_REF][START_REF] Koopman | Numerically accelerated importance sampling for nonlinear non-Gaussian state-space models[END_REF]. In the rest of this section, we formulate our model in a state-space form, and derive closed form expressions of the posterior predictive mean and variance using the chain structure of the ARG process (Equation (3.1)). A simple sampling algorithm for the joint distribution l(θ|D I ) is also proposed.

The state-space representation

A state-space model is a time series model in which there is a latent process in parallel to the observed processes. In our application, this 'time series' refers to the age dimension, and the latent state variables are the mortality coefficients, θ , whereas the observed variables are the observed death counts d x,1 , . . . , d x,T . The system has the following state-space representation 3 :

In Figure 1, the vectors (d x,1 , . . . , d x,T )'s are i.i.d. across x conditional on the benchmark force of mortality, (μ 0,x,t ) x,t , the portfolio exposure (e x,t ) x,t , and θ . The arrow between θ x and d x,t , t = 1, . . . , T indicates that each d x,t depends on θ x via Equation (2.1).

3 The state space based on the reverse order, i.e. from θ I to θ 1 , can be constructed in the similar manner.

In this state-space representation, transition density of each θ x can be derived by conditioning on the counting variable:

l(θ x |θ x-1 ) = ∞ n=0 P[N x-1 = n]l(θ x |n, θ x-1 ) = ∞ n=0 e -βθ x-1 (βθ x-1 ) n (n + 1) (β + ν) ν+n (ν + n) e -(β+ν)θ x θ ν+n-1
x , for x = 2, . . . , I. This expansion has infinitely many terms, but can be approximated by the sum of the K + 1 first terms. In other words, we neglect the possibility of the count variable N x-1 taking a value larger than K. Thus we get:

l(θ x |θ x-1 ) = ∞ n=0 P[N x-1 = n]l(θ x |n, θ x-1 ) = ∞ n=0 e -βθ x-1 (βθ x-1 ) n (n + 1) (β + ν) ν+n (ν + n) e -(β+ν)θ x θ ν+n-1 x ≈ K n=0 e -βθ x-1 (βθ x-1 ) n (n + 1) (β + ν) ν+n (ν + n) e -(β+ν)θ x θ ν+n-1 x ≡ F (θ x-1 )G(θ x ), ( 4.1)
where F and G are both (K + 1) dimensional vectors given by

F(x) = e -βθ x-1 (βθ x-1 ) 0 (0 + 1) , . . . , e -βθ x-1 (βθ x-1 ) K (K + 1) , G(x) = (β + ν) ν+0 (ν + 0) e -(β+ν)θ x θ ν+0-1 x , . . . , (β + ν) ν+K (ν + K) e -(β+ν)θ x θ ν+K-1 x .
Clearly, when K goes to infinity, the approximation formula becomes exact. In practice, the value of K is fixed to be mildly large to ensure a good quality of approximation. The choice of K will be discussed in Appendix 3. Given a fixed K, the conditional distribution, l(θ x |D I ), can be computed using an algorithm similar to the standard Hamilton filter for Markov switching model. The posterior densities are summarized in the following proposition. Proposition 2:

l(θ x |D I ) = P (I -1)l(D I |θ x )G(θ x ) P (I -1)Q(I) , if x = I, (4.2) and l(θ x |D I ) = P (x -1)G(θ x )F (θ x )l(D x |θ x )Q(x + 1) P (I -1)Q(I) , if x < I, (4.3)
where:

• P (0) is the row vector of the marginal distribution of N x , i.e. a negative binomial distribution truncated at K. P j (0) ∝ C j+ν-1 j p ν (1p) j , j = 0, . . . , K, (4.4)

where p = 1 -ρ. • The subsequent value of P (x) is obtained recursively and forward by:

P (x) = P (x -1)M x ,
• The sequence Q(x) is obtained recursively and backward by:

Q(I) = ∞ 0 l(D I |θ I )G(θ I )dθ I Q(x -1) = M x Q(x).
• the age-indexed matrix M x is given by:

M x = ∞ 0 G(θ x )F (θ x )l(D x |θ x )dθ x .
(4.5)

This approach has recently been suggested by [START_REF] Creal | A class of non-Gaussian state space models with exact likelihood inference[END_REF]. For the sake of completeness, we provide the proof of this proposition in Appendix 1. The interpretation of this result is that given D I , the process N is still a (time-inhomogeneous) Markov chain, with transition matrix from N x-1 to N x given by, up to a multiple constant, the matrix M x . Moreover, the vectors P(x) and Q(x) are, up to appropriate normalization, the vector of probabilities of N x given D x = D 1 , . . . , D x and D x = D x , D x+1 , . . . , D I , respectively.

Thus, based on Proposition 2, the computation of the smoothing density necessitates only the computation of the matrices M x , whose (i, j)th entry has the following explicit expression:

M x (i, j) = ∞ 0 (β + ν) ν+i (ν + i) e -(β+ν)θ x θ ν+i-1 x e -βθ x (βθ x ) j (j + 1) θ T t=1 d x,t x exp -θ x T t=1 e x,t μ 0,x,t dθ x = (β + ν) ν+i β j (ν + i) (j + 1) (ν + i + j + T t=1 d x,t ) β + ν + β + T t=1 e x,t μ 0,x,t ν+i+j+ T t=1 d x,t
, ∀i, j = 0, . . . , K.

(4.6)

We can also remark that the smoothing density l(θ x |D I ) is a linear combination of gamma densities γ (δ + k), where k varies between 0 and 2K. As a consequence, the lth moment of this smoothing density E[θ x |D I ] is:

E[θ l x |D I ] = P (x -1)M l,x Q(x + 1) P (I -1)Q(I) , with M l,x = ∞ 0 θ l x G(θ x )F (θ x )l(D x |θ x )dθ x . (4.7)

Predictive mean and variance of the one-year-ahead death count

From Equation (4.7), we obtain the posterior prediction of the one-year-ahead total death count:

E x d x,T+1 |D I = x μ 0,x,T+1 e x,T+1 E θ x |D I . (4.8)
For the conditional variance, we first decompose it into the sum of the conditional variance of individual d x,T+1 's and the covariance of d x,T+1 for different x:

V x d x,T+1 |D I = x V[d x,T+1 |D I ] + 2 x<x CoV[d x,T+1 , d x ,T+1 |D I ].
(4.9)

By the law of total variance, the variance term V[d x,T+1 |D I ] can be further decomposed into:

V d x,T+1 |D I = E V(d x,T+1 |θ x , D I )|D I + V E(d x,T+1 |θ x , D I )|D I = E μ 0,x,T+1 e x,T+1 θ x |D I + V[μ 0,x,T+1 e x,T+1 θ x |D I ]. (4.10)
The first term on the RHS is the contribution of the idiosyncratic variation, and can be derived from (4.7), whereas the second term is the contribution of the uncertainty around the unobserved mortality coefficient θ x , and is equal to:

V[μ 0,x,T+1 e x,T+1 θ x |D I ] = μ 0,x,T+1 2 e 2 x,T+1 E[θ x |D I ] 2 -E[θ 2 x |D I ] .
Note that this formula has been derived under the implicit assumption that the future value of μ 0,x,T+1 is deterministically given. This can be motivated, for instance, when μ 0,x,T+1 is given by the regulator, which usually does not provide uncertainties around the mortality table. In many situations, however, uncertainties in longevity improvements of the benchmark population should be taken into account when forecasting death counts in the insurance portfolio. Accounting for uncertainties of the benchmark mortality rates is straightforward, under the assumption that μ 0,x,T+1 is independent of θ x for each x. In this case, we have:

V[μ 0,x,T+1 e x,T+1 θ x |D I , F T ] = e 2 x,T+1 E[μ 0,x,T+1 |F T ] 2 E[θ x |d I ] 2 -E[μ 0,x,T+1 2 |F T ]E[θ 2 x |D x ] , (4.11)
where F T denotes the mortality history of the benchmark population, i.e. (μ 0,x,t ) x=1,...,I,t=1,...,T .The independence assumption between each μ 0,x,T+1 and θ x indicates that differences in mortality level between the two populations are stable, regardless the mortality improvement in the benchmark population. In other words, the insurance portfolio has the same mortality trend as the benchmark population, i.e. mortality trends in the two populations are coherent [START_REF] Li | Coherent mortality forecasts for a group of populations: an extension of the Lee-Carter method[END_REF]. Under the independence assumption, computing the variance of the total death counts in the portfolio amounts to computing the first two conditional moments of μ 0,x,T+1 for all x. This can be conducted by simulation when a stochastic mortality model, such as the [START_REF] Lee | Modeling and forecasting US mortality[END_REF] model, is specified for the benchmark population.

Similarly, the covariance term in (4.9) is given by:

Cov[d x,T+1 , d x ,T+1 |d I ] = μ 0,x,T+1 e x,T+1 μ 0,x ,T+1 e x ,T+1 Cov[θ x , θ x |D I ] = μ 0,x,T+1 e x,T+1 μ 0,x ,T+1 e x ,T+1 E[θ x θ x |D I ] -E[θ x |d I ]E[θ x |D I ] .
(4.12)

Thus, the computation of the covariance term involves solving for E[θ x θ x |D I ], which is given in the following proposition. Proposition 3: For two ages x < x , E[θ x θ x |D I ] is given by:

E[θ x θ x |D I ] = P (0) M 1 • • • M x-1 M 1,x M x+1 • • • M x -1 M 1,x M x +1 • • • M I-1 ∞ 0 l(d I |θ I )G(θ I )dθ I P (0) M 1 • • • M I-1 ∞ 0 l(d I |θ I )G(θ I )dθ I ,
with M x and M 1,x given in Equations (4.5) and (4.7), respectively.

Proof: See Appendix 2.

Finally, when uncertainties of μ 0,x,T+1 are incorporated, the covariance term becomes:

Cov[d x,T+1 , d x ,T+1 |d I ] = e x,T+1 e x ,T+1 E[μ 0,x,T+1 μ 0,x ,T+1 |F T ]E[θ x θ x |D I ] -e x,T+1 e x ,T+1 E[μ 0,x,T+1 |F T ]E[μ 0,x ,T+1 |F T ]E[θ x |D I ]E[θ x |D I ].
(4.13)

Accounting for heterogeneous policy face values

Besides the number of death, we can also calculate the posterior moments of the total claim amount.

For illustration, we consider a portfolio of policies only with death benefits. In the simplest scenario where all policies have the same face value, multiplying the forecast of the total death count by this face value (resp. the squared face value) yields the next year's expected total claim amount (resp. the variance of the claim).

In practice, however, insurance policies have different face values. Suppose that there are e x,t policies with age x at year t. Their face values are chosen independently by their policyholders, and we denote by V x,t the set of these e x,t observed face values. Then the sum of the face values of the d x,t death events at age x and time t, A x,t , has the compound Poisson representation:

A x,t = d x,t i=1 z i,x,t , (4.14)
where z i,x,t are mutually exclusive draws from the set of face values V x,t . When e x,t is large, these draws can be approximately regarded as i.i.d., and the mean (resp. variance) of each z i,x,t is well approximated by the empirical mean m x,t (resp. empirical variance σ 2 x,t ) of these observed face values, which are easily computable. Thus, we have:

E[A x,t+1 |D I ] = m x,t E[d x,t+1 |D I ] and V[A x,t+1 |D I ] = V E[A x,t+1 |d x,t+1 , D I ]|D I + E V[A x,t+1 |d x,t+1 , D I ]|D I = m 2 x,t V[d x,t+1 |D I ] + σ 2 x,t E[d x,t+1 |D I ].
Similarly, the covariance is given by:

Cov[A x,t+1 , A x ,t+1 |D I ] = m x,t m x ,t Cov[d x,t+1 , d x ,t+1 |D I ].
In practice, m x,t and σ 2 x,t are unknown, but can be estimated from historical claim experience. Finally, the third and fourth conditional moments of I x=1 A x,t given D I can also be obtained in a similar manner, using the law of total cumulance. These computations involve the third and fourth order (co-)moments of d x 's (and hence those of θ x 's), which have closed form expressions as well (see Section 4.2). They can be used to approximate the conditional distribution of the total claim amount using the moment matching technique (see e.g. Ghysels and Wang 2014). Besides, one can also obtain the posterior distribution of the total claim amount via simulation, which is discussed in the next subsection.

Simulating trajectories of θ

Our previous closed form predictive formulae concern the mean and variance of the next year's death count. However, other risk metrics, such as Value-at-Risk (VaR) and Conditional-Tail-Expectation (CTE), may also be interesting in solvency calculation. These risk metrics can be calculated using simulation. This can be done by first simulating trajectories from the posterior distributions of θ , then, for each simulated θ, simulating the number of death for each age. Furthermore, the simulation exercise is useful when the forecast horizon is longer than one year, which is particularly relevant for annuity or pension businesses.

Simulation of the posterior trajectory in standard state space models is based on computationally intensive particle filter, or MCMC (see e.g. [START_REF] Pedroza | A Bayesian forecasting model: predicting US male mortality[END_REF]). In our model, the auxiliary counting variable can be used to simplify the sampling procedure (see [START_REF] Nieto-Barajas | Markov beta and gamma processes for modelling hazard rates[END_REF]Walker 2002, Creal 2017). The basic idea is that, in order to simulate trajectories of θ , we first simulate a trajectory of N from its posterior distribution. Simulating trajectories of N is simpler, since each N x is discrete and has an (approximated) finite state space {0, . . . , K}. Then, in the second step, we simulate each θ x conditional on d x , N x , and N x+1 . The detailed sampling algorithm is provided in Appendix 4.

Application

In this section, we fit the ARG model to a portfolio of endowment policies (that is, policies with death benefits) from an anonymous French life insurance company. The data contain about 15,000 male policyholders over an observation window of 11 years (1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007). The portfolio is not closed, that is, new policyholders can enter the portfolio after the inception date, while existing policyholders can choose to lapse their contract and leave the portfolio at any time during the observation period. In this paper, the French national male mortality data during the same period is used as the benchmark mortality table. The benchmark mortality data are downloaded from the Human Mortality Database. 4For illustration purpose, mortality development of the benchmark population is assumed to be deterministic.

Data

Figure 2(a) displays the age-specific exposure of the portfolio in year 1997. The exposure is adjusted for lapsation.5 Policyholders' ages range between 10 and 102, and has a two-modal distribution. We see that the majority of policyholders are active lives. However, the second peak at around 70 indicates that recently retired people tend to buy endowment policies as well, possibly because of heritage motive. In particular, the proportion of policyholders above age 65 in the portfolio is around 27% in year 1997, while the analogous number is 11% in the French male population. Therefore, uncertainties regarding the total number of deaths in the portfolio is likely to be substantial, due to the higher elderly proportion. Figure 2(b) displays the age-specific mortality rates of the portfolio from year 1997 to 2007. Due to higher longevity risks, the mortality rates are more erratic at older ages. Moreover, the small exposures at these ages could also contribute to the uncertainties. Finally, no clear mortality improvement is observed in the insurance portfolio, probably due to the short observation window. This motivates the use of a benchmark population with reliable mortality trend.

Impact of the age dependence

Since θ are correlated, the posterior distribution of each θ x depends not only on D x , but also on other ages. To illustrate the impact of the auto-correlation coefficient ρ, we plot in Figure 3 the posterior mean of the process (θ x ) x=1,...,110 for ρ = 0.1, 0.9, and ν = 10. For comparison, the posterior means Several observations can be made from Figure 3. First, the posterior means of θ x are smaller than 1 for most ages and all three model specifications. This is reasonable, since mortality rates of life insurance policyholders are in general lower than those of the national average. Second, the posterior mean of θ x converges to 1 when x reaches either the minimal age, 1, or the maximal age, 110. This is due to the fact that exposure in the portfolio decreases drastically when approaching these limiting ages, and thus the posterior means of θ x 's are converging to their prior means (see Equation (3.5)). Third, the convergence speed is the fastest in the independent case, and the slowest with ρ = 0.9. The reason is that the correlation structure of θ leads to a smoothing effect on their posterior means. Moreover, this effect is more substantial at ages with a small exposure and when ρ is larger. In particular, for ρ = 0.9, the spikes on the independent trajectory (dashed line) near age 30 and 105 have almost been eliminated.

Table 1. The posterior predictive mean, variance, and 99.5% value-at-risk of the total number of deaths next year in the portfolio using different values of ρ and ν. The predictive mean (row 1) and variance (row 2) are obtained from closed-form expressions derived in Section 4, while the VaR (row 3) is based on simulation. ν 1 5 1 0 ρ 0.1 0.5 0.9 0.1 0.5 0.9 0.1 0.5 0.9

E x d x,

Projected total number of deaths

After forecasting θ, we forecast the one-year-ahead total death count. In Table 1, we show, for different values of ρ and ν, the posterior mean, variance, and 95% VaR of next year's total number of deaths.

The first two are calculated using the closed form formulae in Sections 4.1 and 4.2. The 95% VaR is obtained by 100,000 simulated trajectories of θ based on the sampling algorithm in Section 4.4. Generically speaking, when ρ is fixed, the predictive mean and variance are increasing in ν. The reason is that ν is a belief parameter and measures the weight assigned to the benchmark mortality, which has in general higher age-specific mortality rates (see Equation (2.4)). Therefore, the larger ν, the larger the weight of the benchmark mortality, and so are the forecasted θ. Moreover, with fixed ν, we see that the predictive mean and variance are decreasing in ρ. Indeed, larger ρ implies stronger correlation between θ x 's, and thus lower posterior means for the θ x 's with small exposure. As a result, the projected mortality rates in the insurance portfolio become smaller when ρ increases.

Variance decomposition of the total death count

Besides projecting the mean and variance of the insurance portfolio's death count in the next period, the ARG model allows us to evaluate the extent to which parameter uncertainty on θ affects the projection. Specifically, the conditional variance of the total number of deaths is given by:

V x d x,T+1 |D I = E x d x,T+1 |D I + V x d x,T+1 |D I -E x d x,T+1 |D I .
(5.1)

The first part on the right hand side of Equation (5.1) corresponds to the (diversifiable) contribution from the sampling variation, whereas the second part is the (non-diversifiable) contribution from the uncertainty on θ x , i.e. the parameter uncertainty. By law of large numbers, the diversifiable part is decreasing (relative to the conditional mean) in the size of the portfolio. As for the second, nondiversifiable part, when the portfolio size and/or the observation window increases, this uncertainty of the underlying mortality rates becomes smaller as well. 6Table 1 shows that, although the sampling variation is larger, the contribution of the uncertainty on θ x accounts for, roughly speaking, 12% of the estimated number of death, E

x d x,T+1 |D I . This number is far from negligible, even when the data covers an observation window of 11 years, which is relatively long for one insurance product. Therefore, by the arguments of Section 2 (see e.g. Equation (2.3)), we can expect even larger contribution from the uncertainty of θ x when the observation window is shorter.

To further illustrate the realized economy due to the decrease of the parameter uncertainty, we estimate the model using different observation windows, and compare the results. In particular, for different T's smaller than the maximal possible value T max = 11, we use the first T years of observations to estimate E[θ x |D I ], then compute the predictive mean and variance of the total death count at year T + 1. Figure 4 displays the posterior mean of θ with T = 1, 5, and 11. We see that the posterior mean of θ becomes smaller when T increases. This result is intuitive, since as the observation accumulates, portfolio's own mortality rate, which is lower than the benchmark population, has larger impact on the estimation of θ .

Let us now examine how the predictive mean and variance of the total number of deaths evolves when the observation window increases from T = 1 to T = 11. Figure 5 shows the predictive mean and variance of next year's death count as a function of T. First, we observe a hump shape for both the mean and the variance. This is the result of two opposite effects. On the one hand, total exposure in the portfolio decreased by around 15% over the 11 years, which would reduce the projected death count. On the other hand, the average age of policyholders had been increasing over time, which lead to higher average mortality rates. More importantly, we see that increasing the observation window clearly reduces the predictive variance relative to the predictive mean, especially when T is small. In other words, the contribution of parameter uncertainty to the variance of total death count becomes less substantial as observation accumulates.

Finally, we illustrate the simulation of the trajectories of θ . Figure 6 displays 100 simulated trajectories of θ drawn from the posterior distribution l(θ|D I ), with ν = 10 and ρ = 0.5. We observe that the average of these trajectories correspond roughly to Figure 3(b), i.e. E[θ x |D I ]. Moreover, as expected, the trajectories are more erratic at lower and higher ages, due to the lack of exposure.

Implications on regulation and reinsurance pricing

Let us now discuss the implications of the ARG model on regulations and pricing, as well as how to determine the values of ν and ρ. In this paper, we assume that ν and ρ are exogenously given. In a frequentist framework, however, it is possible to estimate ν and ρ from the insurance portfolio. Roughly speaking, at age x, θ x is approximately the ratio between the historical death counts of the insurance portfolio and the benchmark population. Therefore, a reasonable proxy of θ x would be t d x,t t μ x,t e x,t

. Using this approximation, we can calibrate the values of ν and ρ. Nevertheless, from the regulatory point of view, this frequentist approach is not recommended. Indeed, the benchmark mortality table for life insurance policies (resp. life annuities ) is often a 'conservative' table which overestimates (resp. underestimates) the mortality rates of the portfolio. Thus if the previous estimation approach is directly applied, we would get raw estimate of θ x which are usually smaller than 1 for a life insurance portfolio. This is contrary to our prior assumption that

E[θ x ] = 1.
So why do we make, at the beginning, this unitary mean assumption? Indeed the possibility for the insurer to set its own proprietary life table has to be monitored very carefully by the regulator, in order to prevent them from regulatory arbitrage. Without the prior assumption E[θ x ] = 1, it would be very tempting for insurers to apply the ARG model, or other similar models, to obtain a proprietary life table, so long as the portfolio mortality level seems to be lower than that implied by the regulatory table. regardless of the reliability of such a conclusion. Using such proprietary life tables is very risky if the observed mortality difference does not hold in the future. In fact, the use of proprietary table should only be allowed when there is strong evidence that the risks of the insurer's portfolio is significantly lower than implied by the benchmark life table. This strong evidence might be the fact that the portfolio is sufficiently large, or that a lower mortality level has been observed over a sufficiently long period.

The Bayesian approach proposed in this paper is more adapted, as it allows to automatically evaluate the reliability of a portfolio's own mortality experience. In our approach, the prior mean, E[θ x ], is fixed rather than estimated. As a result, even when the portfolio's mortality rates seem to be lower, we would still get the posterior means of θ x 's close to 1 if i) the portfolio size is not sufficiently large; ii) or the observation window is not long enough. Contrarily, the predicted mortality would be close to the insurer's own mortality experience when the portfolio size is large, or when the observation window is long. In other words, under the Bayesian approach, the more reliable the proprietary mortality experience, to a larger extent the insurer can take it into account when building the internal model. Therefore, to ensure coherence, the values of ρ and ν, which characterize the prior distribution of θ , should be fixed by the regulator at the industry level, rather than be estimated for each portfolio, and their values will reflect the regulator's risk aversion (via parameter ν), as well as its emphasis on smooth mortality forecasts (via parameter ρ). This approach is also in line with the existing literature (see e.g. [START_REF] Olivieri | Life tables in actuarial models: from the deterministic setting to a Bayesian approach[END_REF]Pitacco 2012, van Berkum et al. 2017).

Besides solvency calculation, our model is also useful for reinsurance pricing, in which case the role of the regulator is replaced by the reinsurer, and the regulatory mortality table is replaced by the reference pricing table of the reinsurer. In this case, the values of ρ and ν could be estimated by the reinsurer, as it has access to a wide range of insurance portfolios. This is the approach adopted by [START_REF] Hardy | A credibility approach to mortality risk[END_REF] and [START_REF] Salhi | A credibility approach of the Makeham mortality law[END_REF]. As a result, the Bayesian approach allows to derive reinsurance quotes that take into account both the portfolio's own mortality risk and the industry experience.

In the context of a reinsurance quote, the willingness of the reinsurer to price the contract using the insurer's own mortality experience depends on its knowledge of the latter. If the insurer's portfolio is small, then the reinsurer will put more credibility on a conservative reference mortality table. When the size of the insurer's portfolio becomes larger, the degree of the asymmetric information becomes smaller, and the reinsurer becomes more willing to take into account the insurer's own mortality experience. Note, also, that when an industry-level data is used, the conditions E[θ x ] = 1 for each x ensures that this pricing approach is fair for the reinsurer. That is to say, when the reinsurer underwrites lots of portfolios, it can expect that the average cost across all the portfolios is equal to that predicted by the benchmark mortality table.

To conclude this section, we provide an example to illustrate the impact of parameter ν and the portfolio size on the credibility of its own mortality experience. Specifically, we create a small portfolio by randomly drawing policyholders from the whole portfolio. The small portfolio consists of around 2,800 policyholders in year 1997, thus its size is about 20% of the whole portfolio. The age composition of the small portfolio is similar to the whole portfolio, and is given in Figure 7.

Figure 8 displays the posterior means of θ for the whole portfolio and the small portfolio with different values of ν and ρ = 0.5.7 First, we see that the posterior means from both portfolios are more volatile with a smaller ν. The reason is that, when the prior belief is weaker (with a smaller ν), the portfolio's own mortality experience has heavier influence on the posterior distribution of θ . Second, with ν = 5 and 10, i.e. when the posterior means are less volatile, the whole portfolio has in general lower posterior means of θ x . Therefore, given a fixed set of ρ and ν, a larger portfolio could gain more credibility on its own mortality experience in mortality forecasting.

Conclusion

This paper proposes a Bayesian non-parametric mortality model for a small population, where mortality of a larger population is exogenously given as the benchmark. The proposed model captures differences in mortality level between the small and the benchmark population by an autoregressivegamma (ARG) process. Compared with existing deterministic models, the proposed model is more flexible whereas at the same time ensures smoothness of the mortality projection across different ages. Moreover, the model is associated with closed form expressions of the (posterior) moments, and a computationally efficient sampling algorithm of the whole posterior predictive distribution, of future death counts of the small population. This makes the ARG model more accessible to practitioners.

In a numerical application, we apply the ARG model on a French life insurance portfolio. In particular, we examine how the portfolio size and the length of observation window would affect the weight assigned to the portfolio's own mortality experience in mortality forecasting. Further, we discuss the impact of the prior distribution, which measures the degree of belief on the benchmark mortality table, on the risk management of the insurance portfolio.

Finally, while the primary application in this paper relates to mortality forecasting of a small population, the proposed model can also be used as building block of a two-population Bayesian mortality model (see e.g. [START_REF] Antonio | Bayesian Poisson log-bilinear models for mortality projections with multiple populations[END_REF]. Indeed, the estimation of Bayesian mortality models typically require substantial computational efforts (see also [START_REF] Czado | Bayesian Poisson log-bilinear mortality projections[END_REF], and the computational costs are even higher in the two-population case. The ARG model provides a tractable distribution for the conditional mortality rate of one population, when the mortality of the other population is assumed to be exogenously given. This could substantially reduce computational costs resulting from simultaneously modeling of two populations.

Disclosure statement

No potential conflict of interest was reported by the authors.

ORCID

Yang Lu

http://orcid.org/0000-0002-7141-1048

where the normalization constant to divide is to ensure that the components of P (D x ) sum up to unity.

In other words, P (D x ) is equal to, up to a multiplication constant, the vector P (x) introduced in Proposition 2: P (x) ∝ P (D x ), ∀x.

Equation ( A3) is just a matrix formulation of the following elementary identity, obtained via the Bayes' rule:

P[N x = j|D x ] = K i=0 P[N x = j|N x-1 = i, D x ]P[N x-1 = i|D x-1 ] ∝ K i=0 G i (θ x )F j (θ x )l(D x |θ x )dx P N x-1 = i|D x-1 .
It is then easily checked that this formula leads to the filtering formula (4.2).

A.2. The smoothing formula (4.2)

Let us consider the joint distribution: that is the ratio between the RHS of Equation (A5) and its integral with respect to θ x .

Appendix 2. Proof of Proposition 3

It suffices to show that the joint smoothing distribution l(θ x , θ x |d I ) is equal to: Algorithm 1 (see Creal 2017, Section 3.3):

l(θ x , θ x |D I ) ∝ P (0) M 1 • • • M x-1 G(θ x )F (θ x )l(d x |θ x ) M x+1 • • • M x -1 × G(
• Draw N I-1 using the vector of elementary probabilities p j ∝ P j (I -1)Q j (I), ∀j = 0, . . . , K, where vectors P and Q are introduced in Section 4.2. In other words, N I-1 is drawn from the distribution of N I-1 |D I . • For each x ∈ [1, I -2] and given N x+1 , . . . , N I , draw N x using the vector of elementary probabilities: Once we have a trajectory of N, we can simulate a trajectories of θ from the distribution θ|N I , D I . Using the Markov chain structure, we can remark that θ x depends on N I and D I only via its left and right neighbours, N x and N x+1 (except for the initial and the final value), as well as the death counts D x . In particular, the sampling algorithm of θ is summarized as follows. , that is the gamma distribution with shape parameter N 1 + T t=1 d 1,t and rate parameter β + T t=1 e 1,t μ 0,1,t . In terms of complexity, the sampling of a trajectory of θ involves only the simulation of 2I variables, all of which have closed form (approximated) density functions. This means that the complexity is significantly lower than standard MCMC methods, which rely on the (often slow) convergence of some auxiliary Markov chains to sample variables with untractable distributions.

Figure 1 .

 1 Figure 1. The age-indexed state-space model.

Figure 2 .

 2 Figure 2. The exposure and mortality rates of the portfolio at different ages and different years.

Figure 3 .

 3 Figure 3.Posterior mean of θ x with age dependence (full line, corresponding to ρ = 0) and without age dependence (dashed line, that is, ρ = 0) with ν = 10.

Figure 4 .

 4 Figure 4. Posterior mean of θ x for different lengths of observation window: T = 1 in full line, T = 5 in dashed line, and T = 11 in dot-dash line.

Figure 5 .Figure 6 .

 56 Figure 5. Predictive mean (in full line) and variance (in dashed line) of total death count in the next year, for different lengths of observation window.

Figure 7 .Figure 8 .

 78 Figure 7. The age-specific exposure in the small portfolio in year 1997.

  l(D I , θ I , D I-1 , θ I-1 , . . . D x+1 , θ x+1 , θ x |D x ) = l(θ x |D x )l(θ x+1 |θ x )l(D x+1 |θ x+1 ) • • • l(D I-1 |θ I-1 )l(θ I |θ I-1 )l(D I |θ I ) = l(θ x |D x )F (θ x )G(θ x+1 )l(D x+1 |θ x+1 )F (θ x+1 )G(θ x+2 )l(D x+2 |θ x+2 ) • • • F (θ I-1 )G(θ I )l(D I |θ I ) (A4)By integrating out θ x+1 , . . . , θ I , we obtain:l(D I , D I-1 , . . . , D x+1 , θ x |D x ) = l(θ x |D x )F (θ x ) ∞ 0 G(θ x+1 )l(d x+1 |θ x+1 )F (θ x+1 )dθ x+1 × • • • ∞ 0 G(θ I-1 )F (θ I-1 )l(d I-1 |θ I-1 )dθ I-1 ∞ 0 l(D I |θ I )G(θ I )dθ I , (A 5 )where the term in the first pair of brackets is a row matrix, the term in the last pair of brackets is a column matrix, and all the intermediate terms are square matrices that have already been computed during the filtering stage. Finally, by Bayes' formula, the smoothing density is equal to: l(θ x |D I ) = l(D I , D I-1 , . . . , D x+1 , θ x |D x ) l(D I , D I-1 , . . . , D x+1 |D x ) ,

  θ x )F (θ x )l(D x |θ x ) M x +1 • • • M I-1 ∞ 0 l(D I |θ I )G(θ I )dθ I , (B 1 )where the normalization constant is equal toP (0) M 1 • • • M I-1 ∞ 0 l(D I |θ I )G(θ I )dθ I .The proof is similar as Proposition 2. Indeed, from Equation (A4), we integrate outθ x+1 , • • • , θ x -1 , θ x +1 , • • • , θ I ,to obtain the counterpart of Equation (A5). Then by Bayes' formula we have: l(θ x , θ x |D I ) = l(D I , D I-1 , . . . , D x+1 , θ x , θ x |D x ) l(D I , D I-1 , . . . , D x+1 |D x ) .

  p j ∝ P j (x) l(D x+1 |θ x+1 )G j (θ x+1 )dθ x+1 × l(D x+1 |θ x+1 )G j (θ x+1 )F Nx+1 (θ x+1 )dθ x+1 , j = 0, . . . , K.In this equation, the termP j (x) l(D x+1 |θ x+1 )G j (θ x+1 )dθ x+1 is, up to normalization constant, the probability P[N x = j|D x+1 ], whereas the second integral l(D x+1 |θ x+1 )G j (θ x+1 )F Nx+1 (θ x+1 )dθ x+1is the (N x+1 , j)th element of matrix M x+1 , that is (up to normalization constant), the transition probability of N x+1 given N x and D x+1 .

  is sampled from the distribution: l(θ I |N I , D I ) ∝ G NI-1 (θ I )l(D I |θ I ) ∝ exp -(β + ν)θ I -θ I T t=1 e I,t μ 0,I,t θ ν+NI-1+ T t=1 dI,t I, that is θ I depends only on D I and N I-1 , and has the gamma conditional distribution with shape parameter ν + N I-1 + T t=1 d I,t and rate parameter β + ν + T t=1 e I,t μ 0,I,t . • for x = 2, . . . , I -1, θ x is sampled from the distribution:l(θ x |N I , D I ) ∝ F Nx-1 (θ x )G Nx (θ x )l(D x |θ x ) ∝ exp -(2β + ν)θ x -θ x T t=1 e x,t μ 0,x,t θ Nx +Nx-1+ν+ T t=1 dx,t x, that is the gamma distribution with shape parameter N x + N x-1 + ν + T t=1 d x,t and rate parameter 2β + ν + T t=1 e x,t μ 0,x,t . • θ 1 is sampled from the distribution: l(θ 1 |N I , D I ) ∝ F N1 (θ 1 )l(D 1 |θ 1 ) ∝ exp -βθ 1 -θ 1

See however,[START_REF] Tomas | Multidimensional smoothing by adaptive local kernel-weighted log-likelihood: application to long-term care insurance[END_REF],[START_REF] Cadena | Semi-parametric accelerated hazard relational models with applications to mortality projections[END_REF] for recent advances in flexible relational models.

Coherence assumption is typically imposed when the modeled populations have similar socioeconomic characteristics[START_REF] Li | Coherent mortality forecasts for a group of populations: an extension of the Lee-Carter method[END_REF]), or when one population is a sub-population of another[START_REF] Dowd | A gravity model of mortality rates for two related populations[END_REF]).

See: http://www.mortality.org/. The application of the ARG model is not restricted to the national mortality table. Regulatory or industrial aggregate mortality tables can also be used as the benchmark mortality table.

For example, if a policyholder lapses his/her policy at sometime between 1 January and 31 December in year t, the year t exposure only counts the time elapsed between 1 January and the time of lapse, but not one entire year.

A similar conclusion has been reached in[START_REF] Hardy | A credibility approach to mortality risk[END_REF], who report that a large company (say, with 1 million policies) can release up to 50 % of the regulatory capital compared to a small company (say, with 50,000 policies). See also equation (2.5) for a discussion.

The comparison is qualitatively similar for different values of ρ.

Appendix 1. Proof of Proposition 2

The derivation of the smooth distribution is provided in Section 3.3 in [START_REF] Gouriéroux | State-space models with finite dimensional dependence[END_REF] for a general class of models with finite dimensional dependence. For the sake of completeness, we provide in this appendix a simplified proof using the latent switching interpretation of the ARG process. The proof is decomposed into two parts. We first derive the filtering formula of θ in Appendix A.1, then derive the smoothing formula in Appendix A.2.

A.1. The filtering formula

Let us compute recursively the filtering distribution l(θ x |D x ), i.e. the conditional distribution of θ x given deaths occurred at ages non larger than x. As discussed in the paper, it is easier to work with the discrete state variable N x in the derivation. Therefore, we introduce the filtering distribution of N x given death count observations up to age x:

(A1)

Then we can remark that the link between the filtering distribution of N x and that of θ x is the following:

This formula is a direct consequence of the Bayes formula and thus its proof is omitted. Hence, in order to obtain the filtering density, it suffices to obtain the vector of probabilities P(D x ). Let us start by providing its initial value P(0), when the conditioning set is degenerated. That is to say, P(0) is the stationary distribution of N x (truncated by upper bound K). Since θ x has a gamma γ (ν, ν) marginal distribution and N x is conditionally Poisson P(βθ x ), the marginal distribution of N x is Negative Binomial with size parameter ν and probability of success p = 1 -ρ. Hence we have formula (4.4).

Let us now provide a simple recursive updating formula to update P(D x ). We have the following property: Lemma 1: The recursive updating formula of the sequence (P (D x )) is given by:

Appendix 3. The choice of K

In this paper, we let K, the truncated upper bound of the stationary distribution of N x , to be 200 for each x. Figure C1 displays the stationary distribution of N x with ν = 1, 5, 10, and ρ = 0.9. We see that K = 200 seems sufficient, as probability of N x exceeding 200 is negligible for all ν's. From the relation β = νρ 1-ρ and Equation (4.4), we have that the prior distribution of N x has a smaller mean and variance when either ρ or ν decreases. Therefore, K = 200 would be sufficient for the prior distribution of N x with any pair of (ρ, ν) ≤ (0.9, 10). Figure C2 displays the posterior distribution of N x for both the whole portfolio and the small portfolio with ρ = 0.9 and ν = 10. We see that the posterior distribution of N x becomes more concentrated for ages with nonzero observed number of deaths. For the rest of the choice of ρ and ν in the paper, the posterior distributions have smaller means and variances. Therefore k = 200 is also sufficient for the posterior distribution of N x in our analysis.

Appendix 4. Sampling of θ

In order to simulate trajectories of θ from the posterior distribution θ |D I , let us first simulate trajectories of N. This is conducted backwardly as follows.