Paul Doukhan
email: paul.doukhan@u-cergy.fr

Joseph Rynkiewicz
email: joseph.rynkiewicz@univ-paris1.fr

Yahia Salhi
email: yahia.salhi@univ-lyon1.fr

Optimal Neighborhoods Selection for AR-ARCH Random Fields with Application to Mortality *

Keywords: Mortality Rates, AR-ARCH Random Field, Estimation, Model Selection, BIC. Mathematics Subject Classification (2000) Primary: 60G70 Secondary: 60G10, 60F99

) investigate and illustrate a simple variant, called the three-level memory model. However, it is not clear which is the best parametrization to use for specific mortality uses. In this paper, we investigate the optimal model choice and parameter selection among potential and candidate models. More formally, we propose a methodology well-suited to such a random field able to select the best model in the sense that the model is not only correct but also most economical among all the correct models. Formally, we show that a criterion based on a penalization of the log-likelihood, e.g. the using the Bayesian Information Criterion, is consistent. Finally, we investigate the methodology based on Monte-Carlo experiments as well as real-world datasets.

Introduction

Recent advances in mortality modeling offer a wide variety of models to capture most of the stylized facts in human mortality dynamics. Most notably, there are the so-called factor-based models widely in use by practitioners, which have known an increasing recognition from the actuarial community, see [START_REF] Barrieu | Understanding, modelling and managing longevity risk: key issues and main challenges[END_REF] and the reference therein. Starting with the well-celebrated [START_REF] Lee | Modeling and forecasting US mortality[END_REF]'s model, most of the subsequent models rely on a factor-based decomposition of mortality surface. These factors are intended to handle the complex patterns of mortality evolution over time and capture essentially the age, period and cohort effects. Although these models are quite intuitive, their statistical properties are, however, not accurately understood. Moreover, the accuracy of such a framework has been widely discussed and it has been shown that the latter fails, in many cases, to cope with some essential stylized facts inherent in mortality dynamics. Indeed, the cohort effect is of paramount importance when dealing with mortality in some countries, e.g. England & Wales and United States. For instance, the factor-based models attempted to incorporate a specific factor dependent on the year of birth as to enhance the goodness-of-fit of models with only age and period effect. Even if this undesired remaining diagonal effect is, generally, accommodated, it is still unclear how such a cohorteffect can be interpreted and identified, e.g. [START_REF] Hunt | Robustness and convergence in the Lee-Carter model with cohort effects[END_REF]. On the other hand, these classic models are based on the assumption of the homoskedasticity of their residuals. However, the assumption of constant variance is always violated in practice as it is varying and does depend on the age and time.

A way to address these issues is to consider the surface of mortality improvements as a sole random field without any further assumption on the particular dependence structure neither the factors driving its evolution. This is done in [START_REF] Doukhan | A Class of Random Field Memory Models for Mortality Forecasting[END_REF] using a parsimonious approach involving a generalization of the well-known AR-ARCH univariate processes, see [START_REF] Doukhan | A fixed point approach to model random fields[END_REF]. The AR-ARCH for mortality surface is a random field where the improvement rates dynamics are described as an autoregressive random field over a given adjacent neighborhood. The errors of this improvement rates also are described as an autoregressive conditional heteroskedastic model. Each part of the model is characterized by a neighborhood giving the number of lags in years and ages to be taken into account. The model is introduced in its general form which is intended to accommodate to every population own characteristics, i.e. cohort effect, cross-cohorts dependence as well as conditional heteroskedasticity. For more insights, in Section 3, we fully describe the model and give some intuitions on its construction.

To estimate the parameters, [START_REF] Doukhan | A Class of Random Field Memory Models for Mortality Forecasting[END_REF] use a quasi maximum likelihood estimator (QMLE) and show its consistency and asymptotic normality. For numerical illustration, they use a simple variant without any empirical basis. Therefore, in the current paper, we investigate the optimal model choice and neighborhoods selection among potential and candidate models. More formally, we propose a methodology well-suited to such a random field able to select the best model in the sense that the model is not only correct but also most economical among all the correct models. Formally, in Section 3, after recalling some inference results in [START_REF] Doukhan | A Class of Random Field Memory Models for Mortality Forecasting[END_REF], we show that a criterion based on a penalization of the log-likelihood, e.g. the one using the Bayesian Information Criterion, is consistent.

The rest of the paper is organized as follows. In Section 2, after introducing the AR-ARCH random field model specification, we recall some theoretical results related to the existence of a stationary solution. In Section 3 we introduce the robust estimation procedure for the parameters based on the QMLE and recall some results on its consistency and asymptotic behavior. These will be crucial for defining a consistent procedure for model selection based a penalization of the QMLE. In this section, we will show that a specific penalization of the QMLE yields to a consistent criterion for selecting the optimal neighborhoods. Finally, in Section 4, we assess the effectiveness of the procedure using Monte-Carlo experiments. Next, we use real-world data-sets to illustrate the procedure and compare the performance of the selected model to some benchmark models.

2 Random Fields Memory Models for Mortality Improvement Rates 2.1 Classic Mortality Models. Denote by m (a,t) the crude death rate at age a and date t. Time is assumed to be measured in years, so that calendar year t has the meaning of the time interval [t, t + 1). For expository purpose and since we will be working with only a subset of historical data, we will henceforth re-index the observable ages by a = 0, 1, . . . , I -1 and the observable dates by t = 0, 1, . . . , J -1; where I and J are, respectively the number of ages and years. The crude mortality rate m (a,t) , among other actuarial quantities used to assess mortality, is of paramount importance when dealing life insurance and pension liabilities. Its future evolution is generally modeled using the so-called factor-based models. Among others, there is the so-called Lee and Carter. In their seminal paper, [START_REF] Lee | Modeling and forecasting US mortality[END_REF] postulated that the (log) mortality rates at different ages are captured by a common factor, and an age-specific coefficient with respect to this common trend. More precisely, we have for any a and t

log m (a,t) = α a + β a κ t + (a,t) , with (a,t) ∼ N (0, σ) (2.1)
where α a is the time average level of log m (a,t) at age a, κ t is the common factor also known as the period mortality effect and β a is the age-specific sensitivity coefficient with respect to κ t . The development of stochastic modeling of mortality was on the basis of the above model. Indeed, the variety of models that have been introduced in the literature, see [START_REF] Cairns | A two-factor model for stochastic mortality with parameter uncertainty: Theory and calibration[END_REF], [START_REF] Renshaw | A cohort-based extension to the Lee-Carter model for mortality reduction factors[END_REF], [START_REF] Cairns | A quantitative comparison of stochastic mortality models using data from England and Wales and the United States[END_REF] and [START_REF] Dowd | Evaluating the goodness of fit of stochastic mortality models[END_REF] among others, were a generalization of the Lee and Carter in the sense that they generally include additional parameters whether dependent on the age, the period or even on the age of birth. These factors are intended to capture the complex patterns of mortality evolution over time, age and cohort. As noted by [START_REF] Mavros | Stochastic Mortality Modelling: Key Drivers and Dependent Residuals[END_REF], "the number and form of these types of effects is usually what distinguishes one model from another ". However, some recent works show the limit of this mainstream approach, e.g. [START_REF] Giacometti | A comparison of the lee-carter model and ar-arch model for forecasting mortality rates[END_REF], [START_REF] Chai | A double-exponential garch model for stochastic mortality[END_REF], [START_REF] Hunt | Robustness and convergence in the Lee-Carter model with cohort effects[END_REF] and [START_REF] Mavros | Stochastic Mortality Modelling: Key Drivers and Dependent Residuals[END_REF] among others. A major drawback of these models relate in particular to the assumption of homoskedastic error terms (a,t) . In fact, the assumption of constant variance is always violated: the observed logarithm of central death rates is much more variable and the volatility is time varying, see e.g. [START_REF] Lee | Evaluating the performance of the Lee-Carter method for forecasting mortality[END_REF] and [START_REF] Gao | Dynamic mortality factor model with conditional heteroskedasticity[END_REF]. Furthermore, the mortality evolution is known to be related to the age of birth, see [START_REF] Willets | The cohort effect: insights and explanations[END_REF]. This is referred to as the cohort effect and translates the persistent of some shocks on mortality among cohorts. This is generally observed when plotting the residuals (a,t) of models such as (2.1) as an apparent diagonal structure which requires additional univariate cohort processes in some countries. This phenomenon has lead to various extensions of the initial Lee-Carter model, e.g. [START_REF] Renshaw | A cohort-based extension to the Lee-Carter model for mortality reduction factors[END_REF] or [START_REF] Cairns | A quantitative comparison of stochastic mortality models using data from England and Wales and the United States[END_REF] and the reference therein, by the inclusion of factors depending on the age of birth. The inclusion of additional univariate processes enhances the goodness-of-fit of the model but over-fit clearly the model and thus produces less reliable forecasts. Moreover, age-period-cohort models lack of identifiability and suffer from the limited interpretability of its structure, see [START_REF] Hunt | Robustness and convergence in the Lee-Carter model with cohort effects[END_REF]. On the other hand, the heteroskedasticity of mortality rates is not taken into account in such models as suggested by some recent empirical works, see [START_REF] Gao | Dynamic mortality factor model with conditional heteroskedasticity[END_REF] and [START_REF] Giacometti | A comparison of the lee-carter model and ar-arch model for forecasting mortality rates[END_REF] among others.

It is of course very important to tackle these limitations when considering a new modeling approach, but it is also essential to take into account the dependence structure between adjacent cohorts. Indeed, some recent works, and even common intuition, point out the importance of cross-cohorts correlation, see e.g. [START_REF] Loisel | core of longevity risk: hidden dependence in stochastic mortality models and cut-offs in prices of longevity swaps[END_REF] and [START_REF] Jevtić | Mortality surface by means of continuous time cohort models[END_REF]. In their empirical work, [START_REF] Loisel | core of longevity risk: hidden dependence in stochastic mortality models and cut-offs in prices of longevity swaps[END_REF] show that correlation among close generations is higher enough to be omitted. The same conclusions were drawn in the very recent work of [START_REF] Mavros | Stochastic Mortality Modelling: Key Drivers and Dependent Residuals[END_REF]. It is worth mentioning that the Lee and Carter's model in Equation (2.1) implicitly imposes a perfect dependence between mortality rates at the single age level. However, as noted by [START_REF] Loisel | core of longevity risk: hidden dependence in stochastic mortality models and cut-offs in prices of longevity swaps[END_REF], although common intuition suggests that correlation among close generations is high but should not necessary be perfect. A first attempt to tackle this last stylized fact was using vector-based formulation of the mortality rates for a given range of ages. More specifically, it does consists on the modeling of the time-dependent vector m (a,t) , for a = 1, . . . , I -1. For instance, [START_REF] Loisel | core of longevity risk: hidden dependence in stochastic mortality models and cut-offs in prices of longevity swaps[END_REF] and [START_REF] Lazar | A multivariate time series approach to projected life tables[END_REF] were among the first to consider a model based on the study of the multivariate time series of the differences of logits of yearly mortality rates taking into account inter-age correlations. Subsequent works focused on vector-autoregressive (VAR) models whether on a general framework without any assumption on the dependency between adjacent cohorts or with a given rigid dependency. For instance, in [START_REF] Mavros | Stochastic Mortality Modelling: Key Drivers and Dependent Residuals[END_REF] used a simple first order vector autoregressive process allowing for interaction between specific adjacent ages. This was applied to the residuals (a,t) for a factor-based model to cope with the remaining non-stationary effect mentioned above. However, the dependency was limited to immediate generations. In other words, the residuals of population aged a at time t are only explained using the ones aged a -1 and a + 1 at the same period. In the same line, [START_REF] Li | Coherent Forecasting of Mortality Rates: A Sparse AutoRegression Approach[END_REF] used the same idea but applied directly to the (transformed) mortality rates. Similar approaches were also considered when it comes to jointly modeling multiple populations, e.g. [START_REF] Gaille | Modelling mortality with common stochastic long-run trends[END_REF] and [START_REF] Salhi | Basis risk modelling: a co-integration based approach[END_REF].

2.2 Through A Random Field Framework. In contrast to this univariate factor-based framework, [START_REF] Doukhan | A Class of Random Field Memory Models for Mortality Forecasting[END_REF] approach the problem of modeling mortality rates by considering the whole surface as a sole random field without any further assumptions on the particular dependence structure neither the factors driving its evolution. Thus, unlike the mainstream models, such a framework is intended to accommodate not only the cross-cohorts dependence but also the conditional heteroskedasticity. This is amended by modeling the whole mortality (improvements) surface as a random field with a Markov property. The use of a random field for mortality modeling was, for example, considered in [START_REF] Millossovich | A bidimensional approach to mortality risk[END_REF] who propose a spatial random field model in which the covariance between adjacent cohorts is only a function of the distance between the ages, but not of ages. However, the introduced framework is developed in continuous time and it is inspired from the assumed similarities between the force of mortality and interest rates. It thus mainly aims at solving pricing problems for life annuities and longevity-linked securities, e.g. [START_REF] Olivieri | Life tables in actuarial models: from the deterministic setting to a Bayesian approach[END_REF] for a discussion. (i) Mortality Surface. In the sequel, we consider the process X s parameterized by the lattice points s = (a, t) and defined as the centered mortality improvement rates IR (a,t) = log m (a,t) /m (a,t-1) . Formally, we let

X s = IR s -IR, (2.2)
where IR is the average improvement rate over I × J. As noted before, empirical studies show that the differentiation of mortality rates (in the logarithmic scale) removes the non-stationarity, see for example [START_REF] Chai | A double-exponential garch model for stochastic mortality[END_REF]. This is also advocated by the models introduced above as the time-dependent factors are described by random walks with a constant drift. However, these models assume that IR s only depend on the observed age. Moreover, the conditional average of the improvement rates should not only depend on the age but also on the cohort, i.e. t -a, as well as the rates experienced in adjacent generations. First, this is to allow for capturing the cohort effect, i.e. the persistent effect of mortality shocks in the same cohort. Second, including experience from neighbor generations allows improving the assessment of the global mortality. Formally, this will account for the diffusion phenomenon well known in demographic theory. Indeed, some changes in health risk behaviors are adopted first among some cohorts, and then diffuse through the population. We can refer to this as learning or diffusive effects. Therefore, in order to account for correlations, in a natural way, across generations, a Markov property for the random field X s is needed. Formally, we assume that X s has interaction with a finite number of neighbors arranged in any manner around the point s = (a, t). This neighborhood is denoted by V ⊂ N 2 \ {0}.

Its shape is of paramount importance as it directly conditions the causality of the random field. In Figure 1, at the lattice points s = (a, t), here referred to as the point (0, 0), the neighborhood is depicted and corresponds to the grayed area, which obviously excludes the point s. This subset is causal in the sense of Doukhan and -5.

-4.

-3.

-2.

-1. 1. 2. 3. 4. 5.

-4.

-3.

-2.

-1.

1.

2.

3.

0

Age lag

Time lag

Figure 1: Causal random field. A bi-dimensional representation of the random field X s = IR s -IR, with s = (a, t). The grayed area repsents the causal neighbor V needed to characterize the evolution of X s .

Truquet (2007), that is, stable by addition.

(ii) Random Field Memory Models. The starting point of the approach introduced in Doukhan et al. ( 2017) is a formulation of the random field X in the sense of [START_REF] Doukhan | A fixed point approach to model random fields[END_REF]. More formally, let V ⊂ N 2 \ {0} be a given neighborhood, let Θ ⊂ R d be a set of possible parameters, for d a positive integer, and consider ξ = (ξ s ) s∈N 2 an independent and identically distributed (i.i.d.) random field, and F a given parametric function taking values in R: F : R V × Θ × R → R (endowed with a suitable norm • ), and consider solutions of the equation:

X s = F (X s-v ) v∈V , θ, ξ s , with s ∈ N 2 . (2.3) For Z ∈ R k (k ≥ 1) denote Z p = (E Z p ) 1/p
for an integer p ≥ 1. In the above equation the set V , as discussed above, refers to a neighborhood of X s used to characterize its evolution. Recall that the set V does not contain the point (0, 0), that is, (s -v) = s for any v ∈ V . This is abstract formulation in which X s is expressed in terms of its own past values and the present and past values of a sequence of i.i.d. random variables. In the case where the set V is causal, see [START_REF] Doukhan | A fixed point approach to model random fields[END_REF], the existence and uniqueness of a stationary solution of Equation (2.3) in L p (p ≥ 1) rely on the contraction principle, which can be summarized in the following two properties:

A-1 F (x 0 , θ, ξ) p < ∞ for some x 0 ∈ R V , A-2 F (x , θ, ξ) -F (x, θ, ξ) p < v∈V α v x v -x v for all x = (x v ) v∈V , x = (x v ) v∈V ∈ R V , where the coefficients α v are such that v∈V α v < 1.
The model in Equation (2.3) is a general formulation of random fields models with infinite interactions, which can be extended to the case where X takes values in R k and F takes values in R q , for integers k, q > 1. These models are not necessarily Markov, neither linear nor homoskedastic. Moreover the inputs do not need additional distributional assumptions. It is thus an extension of ARMA random fields which are special linear random fields, see e.g. [START_REF] Loubaton | Champs stationnaires au sens large sur Z 2 : Propriétés structurelles et modèles paramétriques[END_REF] and [START_REF] Guyon | Random fields on a network: modeling, statistics, and applications[END_REF]. Such an extension yields a novel random field which is capable of taking into account the nonlinearity and spatial dependency well adapted in our context. In other words, we can fit a particular spatial model of the form (2.3) to such data to give an appealing tools for investigating both spatiality and non-Gaussianity patterns in the evolution of mortality surface. Based on this abstract formulation, [START_REF] Doukhan | A Class of Random Field Memory Models for Mortality Forecasting[END_REF] proposed a specific form for the function F that is intended to tackle the various limitations encountered when using the factor-based approaches.

2.3 AR-ARCH-type random fields. We let s = (a, t) ∈ N 2 and consider some subsets V 1 , V 2 ⊂ N 2 \{0}. In view of different stylized facts discussed above and the aim of improving the mainstream approaches, [START_REF] Doukhan | A Class of Random Field Memory Models for Mortality Forecasting[END_REF] considered a Markov property characterizing both the conditional behavior of X s . First, they propose to decompose the specific features of the random field X s into a conditional mean m s and conditional variance σ 2 s defined as

m s = E X s | X (u,v) ; u < a, v ≤ t , and σ 2 s = Var X s | X (u,v) ; u < a, v ≤ t .
(2.4)

These conditional moments are then assumed to evolve as spatial autoregressive process, which has a structure very similar to the structure of the univariate AR model. More specifically, the conditional mean and variance, i.e. m s and σ 2 s in (2.4), are of the following form

     m s = v∈V1 β v X s-v , σ 2 s = α 0 + v∈V2 α v X 2 s-v , (2.5)
where (β v ) v∈V1 , α 0 and (α v ) v∈V2 are real positive coefficients. Here, V 1 and V 2 are the neighborhoods needed to characterize the evolution of conditional mean and variance of each X s in terms of its own past values and the present and past values of the adjacent cohorts. By doing so, we incorporate the observed heteroskedasticity in real-world datasets as well as the persistence of mortality shocks over adjacent cohorts. The main idea is that σ 2 s , the conditional variance of X s , has an spatial autoregressive structure and is positively correlated to its own recent past and to recent values of some adjacent cohorts defined by the set V 2 . Therefore, the error terms are not equal across ages, in extent that they could be larger for some ages or ranges of ages than for others. Given these first two conditional moments, [START_REF] Doukhan | A Class of Random Field Memory Models for Mortality Forecasting[END_REF] introduced a random field model of the form:

X s = m s + σ s • ξ s , (2.6)
where (ξ s ) s∈I×J is an independent and identically distributed (centered) random field. Here, the dynamics of m s and σ s are given in Equation (2.5). More specifically, the latter can be written in a more tractable form as follows

X s = ξ s • α 0 + v∈V2 α v X 2 s-v + v∈V1 β v X s-v . (2.7)
This is a class of conditionally heteroskedastic random field with an autoregressive (AR) component where the innovation is defined from a standardized random field ξ s given as an autoregressive conditional heteroskedasticity (ARCH), see Equation (2.5). The combined model is referred to as the AR-ARCH random field and it is a generalization of the now widely used AR-ARCH models for random processes, see [START_REF] Engle | Autoregressive conditional heteroscedasticity with estimates of the variance of united kingdom inflation[END_REF]. As already mentioned, such a specification is in accordance with some recent results obtained for mortality series that exhibit conditional heteroskedasticity and conditional non-normality features, which are significant empirical long-term mortality structures, see e.g. [START_REF] Giacometti | A comparison of the lee-carter model and ar-arch model for forecasting mortality rates[END_REF], [START_REF] Chai | A double-exponential garch model for stochastic mortality[END_REF] and [START_REF] Chen | Multi-population mortality models: A factor copula approach[END_REF].

Using the notation from [START_REF] Doukhan | A fixed point approach to model random fields[END_REF], recalled in the previous subsection, the parameters of the model are then denoted θ = α 0 , (α v ) v∈V2 , (β v ) v∈V1 where V 2 is the neighborhood for the conditional variance and V 1 is the the neighborhood characterizing the conditional mean. Note that the model in Equation (2.7) satisfies the assumptions A-1 and A-2 needed to ensure the existence and the uniqueness of a stationary solution.

Indeed, note that the function F is given by

F (x, θ, z) = z α 0 + v∈V2 α v x 2 v + v∈V1 β v x v , for x = (x s ) s∈N 2 , θ = (α 0 , (α v ) v∈V2 , (β s ) v∈V1
), and z ∈ R.

Then, letting µ be the law of the centered i.i.d. random field (ξ t ) t∈N 2 , Equation (2.6) can be expressed in the following form

X s = H (ξ s-t ) t∈N 2 , with H ∈ L p R N 2 , µ ,
so that the conditions for the existence and uniqueness of a stationary solution can be written as

ξ 0 p < ∞, κ p ≡ ξ 0 p s∈V1 |α s | + s∈V2 |v s | < 1, p ≥ 1. ( C0
)
For p = 2, a weaker condition than the previous one is given by

ξ 0 2 < ∞, Eξ 0 = 0, κ 2 2 ≡ ξ 0 2 2 s∈V1 |α s | 2 + s∈V2 |v s | 2 < 1, p = 2. (C1)
3 Estimation, Asymptotic Inference and Model Selection

In this section, we study the maximum likelihood estimator (MLE) for the parameter θ in the AR-ARCH random field model in Equation (2.7). More precisely, we consider an approximation of the MLE called Quasi-Maximum Likelihood Estimator (QMLE) and prove its consistency, together with the asymptotic normality. To this end, we assume throughout the sequel that the errors ξ s is an i.i.d. zero-mean and centered Gaussian random field.

3.1 Quasi Maximum Likelihood Estimator (QMLE). For the sake of simplicity we consider that the observation are real numbers (i.e. taking values in R k , with k = 1), the results exposed here can be extended to more general systems, but this remains a topic for future research. For a real function g and an integer p, let us write g(X) p := (E |g(X)| p ) 1 p , the L p norm. Suppose that the random field is observed for

I = {s-v, v ∈ V 1 ∪ V 2 } (the initial states) and s ∈ O ⊂ N 2 .
Let us write T the number of observations in O, we assume that the following equation can always be written using s ∈ O and s -v ∈ O ∪ I. Conditionally to I the quasi log-likelihood of the model can be written as:

L T (x s , s ∈ O; θ) = 1 T s∈O - 1 2 ln α 0 + v∈V2 α v x 2 s-v - x s -v∈V1 β v x s-v 2 2 α 0 + v∈V2 α v x 2 s-v . (3.1)
We will consider the estimator based on maximizing the above function (QMLE) over the set Θ, which will be denoted θ T , i.e.

θ T = arg max θ∈Θ L T (x s , s ∈ O; θ) . (3.2)
We say that θ T is the QMLE estimator. In order to study the consistency of this estimator we assume the following properties:

H-1 Finite second order moment, i.e. E(X 2 s ) < ∞.
H-2 The model is identifiable, in the sense that for each set of parameters

(β v ) ∈V1 , (α v ) ∈V2 and (β v ) ∈V1 , (α v ) ∈V2 we have v∈V1 β v X s-v - v∈V1 β v X s-v 2 = 0 ⇔ (β v ) v∈V1 = (β v ) v∈V1 ,
and

α 0 + v∈V2 α v X 2 s-v α 0 + v∈V2 α v X 2 s-v -1 a.s. = 1 ⇔ (α v ) v∈V2 = (α v ) v∈V2 .

H-3

The set of possible parameters Θ is compact and the true parameter θ 0 of the model (2.7) belongs to the interior of Θ.

H-4 E X 4 s < ∞.
These are classic assumptions required for consistency and asymptotic normality of the QMLE used in [START_REF] Doukhan | A Class of Random Field Memory Models for Mortality Forecasting[END_REF]. The hypothesis H-1 assumes that the variance of the random field X is finite, which is in line with the object of interest, i.e. mortality improvements. On the other hand, H-2 ensures the identifiability of the model, which is a critical condition for the consistency as we may see later on this section. This will impose that the quasi likelihood L T (x s , s ∈ O; θ) has a unique maximum at the true parameter value θ 0 over the compact set Θ from assumption H-3. Based on these assumptions, [START_REF] Doukhan | A Class of Random Field Memory Models for Mortality Forecasting[END_REF] derive a first result which deals with the asymptotic behavior of the (quasi)-likelihood recalled in the following theorem.

Theorem 1. For the model in Equation (2.7), under the assumptions H-1, H-2 and H-3

L(θ) = lim T →∞ L T (x s , s ∈ O; θ)
exists a.s. for all θ ∈ Θ and is uniquely maximized at θ 0 . This is the first step toward showing that the QMLE is consistent for the parameters of the jointly parametrized conditional mean and conditional variance. To this end, notice that the set of continuous functions G defined as

G = g (x s , (x s-v ) v∈V1 , (x s-v ) v∈V2 ; θ = - 1 2 ln α 0 + v∈V2 α v x 2 s-v - x s -v∈V1 β v x s-v 2 2 α 0 + v∈V2 α v x 2 s-v , θ ∈ Θ , is Glivenko-Cantelli for the L 1 -norm
√ T θ T -θ 0 L -→ N 0, A -1 0 B 0 A -1 0 , where          A 0 = -E ∂ 2 g(X s , (X s-v ) v∈V1 , (X s-v ) v∈V2 ; θ 0 ) ∂θ 2 , B 0 = E ∂g(X s , (X s-v ) v∈V1 , (X s-v ) v∈V2 ; θ 0 ) ∂θ ∂g(X s , (X s-v ) v∈V1 , (X s-v ) v∈V2 ; θ 0 ) ∂θ T . (3.3)
Let us give some interpretations of the above results. (i) Asymptotic normal distribution result in Theorem 2 allows us to approximate the distribution of the estimated parameters θ. This can be used, for instance, to obtain confidence intervals, to conduct hypothesis testing. (ii) Results recalled in Theorem 2 are also of paramount importance when it comes to mortality forecasting. As the amount of data at our disposal is, generally, limited, the parameter estimates using the QMLE in (3.2) most inevitably be subject to some degree of uncertainty. As demonstrated by [START_REF] Cairns | A two-factor model for stochastic mortality with parameter uncertainty: Theory and calibration[END_REF], the parameter uncertainty forms a significant element of the uncertainty in forecasts of future mortality. Thus, Theorem 2 allows us to quantify this uncertainty on parameters based on their asymptotic distribution. These different results deal mainly with the estimation of the parameters. [START_REF] Doukhan | A Class of Random Field Memory Models for Mortality Forecasting[END_REF] also introduce an estimation method based on moments matching. However, the latter is less effective than the one based on the likelihood maximization. Moreover, there is no statistical inference results for the moment method which makes it less interesting for our application. On the other hand, the question of model choice remained open. Indeed, the authors introduce an illustrative model considering two distinct neighborhoods for the conditional mean and the conditional variance. Formally, the two neighborhoods considered in [START_REF] Doukhan | A Class of Random Field Memory Models for Mortality Forecasting[END_REF] are given by V 1 = {(1, 1)} and V 2 = {(1, 0), (0, 1)}. To some extent, numerical applications show that this simple variant is able to outperform classic models in predicting both mortality rates and residual life expectancies. However, the particular form of the neighborhoods does not have any empirical basis as the spatial dependence between adjacent cohort may change from one country to another. It is quite obvious that enlarging the sets V 1 and V 2 will fit the data better, in particular, if one compares the QMLE attained by each model. Moreover, there is no quantitative tool developed for such a class of models to optimally select the best model, in the sense that full feature of a given (mortality) population is represented. In general, when various models are plausible, the best model should be fitted within each of the competing types and compared using an information criterion. The latter encompasses the (quasi) likelihood attained by each model but penalizes models that are over-parametrized. The choice of the penalty is of critical importance to ensure the consistency of the selected model. In other words, the convergence of the selected neighborhoods V 1 and V 2 to the true ones. In our case, in the literature, there is no theoretical results on the appropriate criterion to be used to ensure this consistency. Therefore, in the next section, we will show that, provided some technical assumptions on the asymptotic behavior of the penalty term, the consistency of the model selection is fulfilled by some of the classic procedures.

Optimal Model Selection.

The previous subsection recalls some results of [START_REF] Doukhan | A Class of Random Field Memory Models for Mortality Forecasting[END_REF] dealing with inference about the model (2.6). The model being general, it is not clear which is the best parametrization to use as there are a number of potential, or candidate, models. Here, being best is in the sense that the model is not only correct but also most economical, meaning that it is simplest among all correct models. In this section we deal with the problem of model selection and propose a methodology well-suited to the random field specification in Equation (2.6). The results recalled in Theorem 2 allows to prove the tightness of the likelihood ratio test statistic (LRTS) and so to prove consistency of an estimator of the true neighborhoods for the conditional mean and the conditional variance using an information criterion like the Bayesian Information Criterion (BIC). First we define the LRTS and give its asymptotic behavior.

(i) Likelihood ratio test statistic (LRTS). Suppose we wish to test the null hypothesis H 0 : θ ∈ Θ 0 versus the alternative H 1 : θ ∈ Θ 1 where θ 0 is (a segment of) an affine linear subspace of (a segment of) a bigger affine linear space Θ 1 . Write θ T,0 and θ T,1 , the (quasi)-maximum likelihood estimators if the parameter set is Θ 0 and Θ 1 respectively. Here, we assume that the true value of the parameter θ 0 is an inner point of Θ 0 (so of Θ 1 ).

According to the definition of the quasi log-likelihood (3.1), the LRTS denoted by ∆ T is given as follows:

∆ T = -2T sup θ∈Θ1 L T (x s , s ∈ O; θ) -sup θ∈Θ0 L T (x s , s ∈ O; θ) .
(3.4)

Let us recall that a family of random sequences {Y n , n = 1, 2, . . .} will be said to be o P (1) if for every δ > 0 and ε > 0 there exists a constant N (δ, ε)

such that P (|Y n | < ε) ≥ 1 -δ for all n ≥ N (δ, ε).
A direct application of Theorem 2 and classic results on QMLEs (see [START_REF] Vuong | Likelihood ratio tests for model selection and non-nested hypotheses[END_REF]) yields:

∆ T = √ T θ T,1 -θ T,0 T A 0 √ T θ T,1 -θ T,0 + o P (1) T →∞ -→ dim(Θ1)-dim(Θ0) i=1 λ i χ 2 i (1). (3.5)
Where χ 2 i (1) 1≤i≤dim(Θ1)-dim(Θ0) are independent χ 2 variables with one degree of freedom and (λ i ) 1≤i≤dim(Θ1)-dim(Θ0) are strictly positive real numbers. Note that there is no analytical computation of the distribution of a weighted sum of χ 2 variables. However, this result shows that the LRTS is tight (bounded in probability) for our model and we can assess the true model thanks to a suitable information criterion.

(ii) Penalized QMLE. Let us denote V 1 and V 2 the maximal possible neighborhoods for the conditional mean and the conditional variance of model (2.7). A sub-model of

X s = ξ s α 0 + v∈V2 α v X 2 s-v + v∈V1 β v X s-v ,
will be specified by setting some parameters α v and β v to 0. Let us write V 0 = V 0 1 ∪V 0 2 the smallest neighborhood such that the model (2.7) realizes the true regression function. For a fixed neighborhood V = V 1 ∪ V 2 , let |V | be the number of elements of the neighborhood associated to non-zero parameters α v or β v and write Θ V the set of these non-zero parameters in model (2.7). In order to optimally select the adequate model we need to introduce another quantity that aims at penalizing the less economical models. Indeed, penalization is a classic approach to select a model. In short, penalization chooses the model minimizing the sum of the empirical risk, i.e. how well the model fits data, and of some measure of complexity of the model (called penalty). In our case, let us define the penalized QMLE as follows:

U T (V ) = T × sup θ∈Θ V L T (x s , s ∈ V; θ) -a T (|V |), (3.6)
and write V the neighborhood maximizing (3.6). Here, a T (|V |) is a penalty which is a function of |V |, i.e. the model complexity. In the sequel, the selection approach consists in assuming the existence of a true model of minimal size and aims at finding it. More formally, this corresponds to the existence of the set

V 0 = V 0 1 ∪ V 0 2
such that the model in Equation (2.7) realizes the true regression function. Then, we are aiming at localizing these using a selection criterion. To this end, the following assumptions are needed to ensure that the selected neighborhoods tend almost surely (with probability one) to the true model when the dimension of the data goes the infinity:

I-1 V 1 , V 2 (maximal possible neighborhoods) are finite set and V 0 1 ⊂ V 1 , V 0 2 ⊂ V 2 .
I-2 a T (.) is increasing, a T (k 1 ) -a T (k 2 ) goes to infinity as T tends to infinity as soon as k 1 > k 2 and a T (k)/T tends to 0 as T tends to infinity for any k ∈ N.

Theorem 3. Under the assumptions of Theorem 2-(i), I-1 and I-2, V converges in probability to V 0 .

Proof. Let us write V = V 1 ∪ V 2 . Applying (3.5), yields to

P(| V | > |V 0 |) ≤ |V| k=|V 0 |+1 P U T (k) ≥ U T (k 0 ) = |V| k=|V 0 |+1 V, |V |=k P T × sup θ∈Θ V L T (x s , s ∈ V; θ) , -sup θ∈Θ 0 L T (x s , s ∈ V; θ) ≥ a T (|V |) -a T (|V 0 |) T →∞ -→ 0, since T × sup θ∈Θ V L T (x s , s ∈ V; θ) -sup θ∈Θ 0 L T (x s , s ∈ V; θ) is bounded in probability and a T (|V |)-a T (|V 0 |)
goes to infinity as T tends infinity (Assumption I-2). Now, we can write

P(V 0 ⊂ V ) ≤ V,V 0 ⊂V P sup θ∈Θ V L T (x s , s ∈ V; θ) -sup θ∈Θ 0 L T (x s , s ∈ V; θ) ≥ a T (|V |) -a T (|V 0 |) T ,
and, if

V 0 ⊂ V , sup θ∈Θ V L T (x s , s ∈ V; θ) -sup θ∈Θ 0 L T (x s , s ∈ V; θ) converges in probability to sup θ∈Θ V L(θ) - L(θ 0 ) < 0, so P(V 0 ⊂ V ) T →∞ -→ 0 because (a T (|V |) -a T (|V 0 |))/T T →∞ -→ 0. Finally, V converges in probability to a neighborhood such that | V | ≤ |V 0 | and V 0 ⊂ V , hence V T →∞ -→ V 0 in
probability, which is the desired the result.

Note that it is easy to see that assumption I-2 is fulfilled by the Bayesian Information Criterion (BIC). Therefore, according to the above theorem, the BIC is a consistent selection criterion for our model.

(iii) Bayesian Information Criterion. As noted before, model selection is the choice of an element V ⊂ V 1 ∪ V 2 (neighborhood) using the data records X s with s ∈ O. To do this, we use the Bayesian Information Criterion (BIC) as a decision function that penalizes the (quasi) log-likelihood at rate log T by the model dimension, i.e. a T (|V |) = |V | log T . More formally, we choose the model (neighborhoods) as

V = arg min V ⊂V1∪V2 U T (V ).
This choice represents a compromise between goodness-of-fit and a simple, interpretable model. Also, thanks to Theorem 3, such a criteria selects the true model as V → V 0 in probability. In addition, the use of the BIC allows for comparative analyses to the existing models. Indeed, most of the recent work on mortality modeling rely on this criteria to rank the models, 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] and [START_REF] Dowd | Evaluating the goodness of fit of stochastic mortality models[END_REF] among others. However, we have to adjust the (quasi)-likelihood in order to include the information at the edges of the space I. More specifically, one should let the parameters α ν and β ν equal to zero for all data that are not contributing to the QMLE.

(iv) Search Algorithm. The selection of the optimal model needs to browse over the maximal sets V 1 and V 2 looking into all the possible combination of points s = (a, t) ∈ V 1 ∪ V 2 separately for the conditional mean and the conditional variance. For instance, if the neighborhoods V 1 and V 2 consist, respectively, of n and m elements; then the implementation of the selection procedure for this random field model has a 2 n+m complexity. In other terms, when n = m = 3, e.g. V 1 = V 2 = {(1, 0), (1, 1), (0, 1)}, we need to test 64 models. Even if this is computationally very simple, the time complexity can make the model choice rapidly infeasible. In fact, when -4.

-3.

-2.

-1. 1.

2.

3.

-3.

-2.

-1.

1.

2.

0 (0, 0) (-1, 0) (-1, -1) (0, -1) (0, -2) (-1, -2) (-2, -2) (-2, -1) (-2, 0) Figure 2: The initial (maximal) candidate neighborhoods V 1 ≡ V 2
given as a collection of points (i, j) ∈ {(1, 0), (1, 1), (0, 1), (1, 2), (2, 1), (2, 2), (0, 2), (2, 0)}.

considering larger sets the number of models grows drastically. Indeed, when one considers maximal sets V 1 and V 2 as depicted in Figure 2 for instance, we should test 2 16 = 65536 models, making this particular step very time-consuming. In fact, the selection procedure involves the optimization of the (quasi)-likelihood in order to estimate the parameters of each candidate model. For instance, if one relies on some statistical softwares, such as R, the (non-linear) optimization routines in most packages will be sufficient to find the values of the parameters θ that maximize the (quasi)-likelihood. This can be, for instance, implemented based on some derivative-free algorithms as the Nelder-Mead using heuristics to search for the maximizers under the constraint of stationarity (C0) or (C1). However, these methods have excessive run-time and the characterization of the optimal model can take some days or even months. Therefore, in order to enhance this step we can supply the optimization algorithm with the closed form formulas of the gradient and Hessian of the likelihood function g, see Appendix A1 .

(v) Spatial Autocorrelation Function. In order to enhance the computation time, we can choose to reduce the dimension of the maximal neighborhoods. To this end, we will build an empirical procedure, which does not aim at robustness, on the basis of the well celebrated Box & Jenkins procedure. The latter advocates the inspection of the autocorrelation function both for the time-series as well as their residuals. The pre-selection of the maximal lags relies on the confidence on the level of the autocorrelation. Therefore, in the sequel, we consider the generalization of the autocorrelation for random field, called the spatial autocorrelation function. Formally, denote the covariance function of the random field X as C(s, s -h) = Cov(X s , X s-h ), where s, s -h ∈ I. Since X is stationary, we have C(s, s + h) = C(0, h) = C 0 (h) for all s ∈ I, h, and a stationary covariance function C 0 . Now, let S(h; I) = {s : s ∈ I, s + h ∈ I} be the collection of pairs that are distanced by h and denote |S(h; I)| its cardinality, i.e. the number of distinct pairs in S(h; I). Under the constant (zero)-mean assumption, the following statistic

A(h; I) = 1 |S(h; I)| s∈S(h;I) X s X s-h , (3.7)
is an estimator of the covariance C 0 (h) = C(0, h) for a spatial lag h ∈ R2 + and ρ(h) = C 0 (h)/C 0 (0) is the estimator of the autocorrelation for the random field X. This functional will be used later on to determine, empirically, the candidate maximal sets V 1 and V 2 . A more concise inspection of the asymptotic behavior of the autocorrelation in Equation (3.7) for AR-ARCH random field is, however, beyond the scope of the current paper and will be explored in future work. records would induce a misspecification of the future mortality decline. In fact, the mortality improvements were systematically lower at the beginning of the last century till the end of the second world war so as the cohort effect. Therefore, in the following empirical analysis we restrict the mortality data historical data and do include the recent period. More precisely, we use the period 1970 -2016. The age band will be fixed to 55 -89 to be in line with recent literature in mortality modeling and thus enables readers to compare the estimation results to other papers more readily.

4.2 Real-World Data Application. In order to illustrate our model selection methodology, we consider the mortality data from the Human Mortality Database (HMD) 2 for the male population of France (FRA), England & Wales (GBP) and the United States (USA). We have at our disposal, the exposures and deaths sample ranging from 1970 to 2016, and for 35 ages, 55 to 89 inclusive.

(i) Model Selection. Before proceeding to the selection of the optimal neighborhoods using the criterion in Subsection 3.2, we explore the autocorrelation function of the random field X using the definition in Equation (3.7). The aim is to reduce the dimensionality of the maximal neighborhoods V 1 and V 2 characterizing the two conditional moments. Using larger neighborhoods will not only require heavy machine calculation and thus take significant time to test all the candidate models but it is more likely to overfit the data. Therefore, we inspect the autocorrelation function to determine potential maximal lags in years and ages. This is motivated by the univariate time-series practice which advocates the examination of the so-called ACF and PACF plots3 in order to identify the maximum order (lag) of the models, based on the characterization of the asymptotic behavior of the autocorrelation. In most univariate models, the autocorrelation decreases with a given rate suggesting the nullity beyond a certain lag. The latter is used as the maximum lag. Unfortunately, in our case, there is no available asymptotic results for the behavior of the autocorrelation in 3.7, in particular when the underlying random field is described using the AR-ARCH structure in Equation (2.7). Developing a specific test for this particular model is challenging but goes beyond the scope of the current paper. Here, we will only use a visual inspection of the autocorrelation function to determine the maximal neighborhoods V 1 and V 2 without any theoretical basis. If a test or the asymptotic behavior of the autocorrelation were available, one has to plot the confidence regions with a given probability and thus select the candidate models within these regions. Figure 3 shows the autocorrelation function (ACF) of the random field X for the three populations (bottom panel) and the ACF of the squares of the random field (top panel). For each population, the squares are much more autocorrelated than the initial improvement which is in line with the intuition behind ARCH structure. The first orders autocorrelation, in year and ages, lie in the interval [0.4, 0.6] in the immediate neighborhoods, and they gradually decline after some lags. These autocorrelations are sufficiently large to be taken into account locally. In view of the amount of data at our disposal, we may think that the computation the autocorrelation for large lags are not consistent as few pairs (X s , X s-h ) are available. To recap, these figures are in supporting the initial hint on the autoregressive behavior of the mortality improvement rates together with the inter-generational dependence. On the other hand, the three populations show substantial evidence of ARCH effects given the structure of the conditional variance, as judged by the autocorrelations of the squared process in the top panel of Figure 3. Indeed, the ARCH behavior is generally characterized by a little correlation in the process itself and yet profound correlation in the squared process, see Figure 3.

For the purposes of the present article, and based on the considerations above, it suffices to consider a narrow characterization of each maximal neighborhood. Thus, for convenience and sake of applicability, the maximal neighborhoods used throughout the sequel are those displayed in Figure 2. Although this does not capture the whole dependence as described in Figure 3, we will check later on this section if the selected model captures conveniently the stylized behavior in each population's mortality. The grayed area corresponds to parameters that were not selected by the procedure and are excluded from the neighborhoods.

Age lags

Year lags 0 -1 -2 0 -1 -2 - 0.014 1 - - - - - - Cond. Variance -0.005 0.001 1 - - - 0.051 - - 0 -1 -2 0.036 0.008 1 - - - - - - - -0.
In Figure 4, we report the results of the selection procedure. A set of 2 16 models were investigated and the one maximizing the BIC criterion is selected. The parameters of the neighborhood with "-" do not belong to the optimal neighborhood. The remaining parameters are estimated and their values are reported.

(ii) Diagnostic Checks. In order to assess the validity of the fitted models for these populations, we consider some standard diagnostics checks. Formally, we investigate the distributional assumption on ξ s under the correct model. Therefore, given the estimation of the parameters ( α 0 , ( α v ) v∈V2 , ( β v ) v∈V1 ) for optimal neighborhoods V 1 and V 2 in Figure 4, we derive the corresponding implied residuals given by ξ s = m s / σ s , where the conditional mean and variance m s and σ s are given in Equation (2.5). These are assumed to behave as i.i.d. standard normal random field. Figure 10 reports the QQ-plots of ξ s depicting the empirical quantiles against those of a standard normal distribution. We notice that the plots are nearly straight lines except some outliers. However, the sample size is 2030, so the outliers are a very small fraction of the data. Thus, a Gaussian distribution would be suitable for the residuals ξ s . Moreover, the histograms of the residuals, shown in Figure 11, appear to be normal with slight deviations at the tails for all populations. Further evidence may be obtained by applying normality tests such as those in [START_REF] Razali | Power comparisons of Shapiro-Wilk, Kolmogorov-Smirnov, Lilliefors and Anderson-Darling tests[END_REF]. These are some of the frequently used statistical test for residual analysis. Failing these tests indicates that the residuals are not normal, which is an indication that the selected model is not consistent with the data. Results are displayed in Table 2, which shows that for most of the tests the residuals do validate the hypothesis of normality. Furthermore, we want to ensure that the residuals do not exhibit any serial correlation which may violate the underlying assumption of independence. In Figure 9, we plot the spatial autocorrelation of the standardized residuals, as implied by the selected model using the definition in Equation (3.7). None of the plots provide any evidence of a persistent autocorrelation. Indeed, regardless of the lag in the age and year, the autocorrelation does not exceed 0.18 and is equal or close to 0 in most cases. This being said, we can conclude that the residuals for the three populations do satisfy the assumptions. Hence, the whole information embedded in each mortality surface is perfectly captured by the selected models and the remaining information consists only of white noise. On the other hand, we further examine the fitting of the selected models using residual heatmaps. The right panel of Figure 5 displays the residuals against the calendar year and age. In the left panel we also depict the residuals of the Lee-Carter model fitted to the same data. These figures allow to compare the goodness-offit behavior of the models. The plots of the Lee-Carter residuals indicate an obvious remaining structure not captured by the model. In particular, we observe clusters of high and low residuals, which is a specific feature that characterizes heteroscedastic fields. Indeed, this means that "large changes tend to be followed by large changes, of either sign, and small changes tend to be followed by small changes" a well known effect in financial literature also called volatility clustering. In actuarial literature, such an effect is referred to as a cohort effect and represents the tendency of some group of cohorts to experience more rapid improvements mortality rates than generations born either side of this period. However, this effect is captured by our model as we can see from Figure 5, which is also intended to inspect if any pattern is still present and did not perfectly captured by the selected models. Indeed, we observe that these residuals are reasonably randomly distributed along both year and age dimensions. The heatmaps only show small traces of diagonal patterns around the cohorts born in 1920 and 1918 for the French population, but these are barely noticeable. These remaining isolated effects may be related to anomalies in the datasets. In fact, [START_REF] Cairns | Phantoms never die: living with unreliable population data[END_REF] highlighted the presence of such anomalies and errors in the data in period mortality tables in the form of isolated cohort effects for years of birth around the two world wars, that is 1915, 1920, 1940 and 1945. We can, thus, conclude that all the selected models succeed to capture all the main features present in the initial mortality improvement rates. (iii) Predictive Performance. The estimated models are used to predict future evolution of mortality. We look at the mortality rates predicted by our model and the Lee-Carter model in Equation (2.1). Here, we only validate the performance based on an out-of-sample analysis instead of an in-sample inspection. This is because the benchmark model is over-parametrized and might undoubtedly outperform over the AR-ARCH.

In Figure 7, we display the 95% confidence intervals as well as the median prediction of the death rates from the two models for selected ages 65, 75 and 85. First of all, we observe that the projected mortality rates continue the historical downward trend for the considered populations. In comparison with the forecasts of the Lee-Carter, we notice that the mortality rates predicted by our model show substantial differences at higher ages, i.e. 75 and 85 for instance. This phenomenon is also observed by [START_REF] Doukhan | A Class of Random Field Memory Models for Mortality Forecasting[END_REF] using a simple variant of the model (2.6). For lower ages, at first sight, we see that the AR-ARCH model provides forecasts that are very close to the benchmark, which suggests that there is no tendency to a general out-performance at these age levels. This is due to the boundary condition which makes the bias on the projection of the mortality rates made at lower ages spread over the diagonals. However, the impact of such a bias has to be minimized as the cohort effect and cross-generation dependency is limited at lower ages. This is why the model produces significantly different forecast at high ages. Moreover, we observe that the evolution of mortality rates exhibits downward drifts comparable to the benchmark for all populations with the exception of the UK population where our model anticipates a more drastic fall of mortality compared to the benchmark. Regarding the confidence intervals, we can observe that those related to the benchmark are quite narrow. This may be explained by the rigid structure of the Lee-Carter as it was widely discussed in the literature. For instance, in [START_REF] Lee | The Lee-Carter method for forecasting mortality, with various extensions and applications[END_REF], the Lee-Carter model is shown to generate overly narrow confidence intervals and generally result in an underestimation of the risk of more extreme outcomes. On the other hand, the forecasting intervals produced by the AR-ARCH are substantially wider than those of the Lee-Carter. The difference is even more pronounced for the higher ages, i.e. 75 and 85. This is mainly due, as discussed earlier, to the increasing uncertainty at higher ages, which is not captured by the benchmark that assumes a constant variance regardless of the age. The AR-ARCH allows to cope with this drawback as the variance of the mortality improvements is also modeled as an autoregressive random field. For a closer examination, we use the forecasts provided by the out-of-sample analysis and derive the corresponding remaining period life expectancies. Figure 7 shows 95% confidence intervals for life expectancies forecasts for ages 65, 75 and 85. We also depict with solid lines the median life expectancies of these projections. First, we can see that the AR-ARCH random field model produces forecasts at least equivalent to the LC at lower ages. However, for higher ages, our model suggests a more increasing life expectancy compared to the Lee-Carter. As soon as the confidence intervals are concerned, we remark that those related to the benchmark are quite narrow and underestimate the evolution of life expectancy. Therefore, the AR-ARCH forecasts should encompass a broader range of probable outcomes. The most striking feature of the AR-ARCH model is its ability to predict very high life expectancies at higher ages. It is, then, important to consider biological reasonableness of such a behavior, see [START_REF] Cairns | A two-factor model for stochastic mortality with parameter uncertainty: Theory and calibration[END_REF]. As noted by [START_REF] Liu | The Locally Linear Cairns-Blake-Dowd Model: A Note On Delta-Nuga Hedging Of Longevity Risk[END_REF], this may be interpreted in terms of a limiting behavior of life expectancies. In other terms, the future life expectancy should not exceed a certain range. However, there is no common census of experts on such a limit. Therefore, in order to validate the behavior of the two models we fit the models to data on the period ranging from 1970 to 1999 for each population, and forecast the future improvement rates from year 2000 to 2016. Next, we recover the corresponding mortality rates as well as the life expectancies.

Besides the visual inspection of the forecast, we introduce some quantitative metrics to assess the forecast accuracy. In fact, as noted before, looking at Figure 7 may not reveal the actual performance over the whole mortality surface. Therefore, we consider an aggregate representative metric that summarizes the performance of each model. More precisely, we considered the mean absolute forecast error (MAFE) and root mean squared forecast error (RMSFE). These criteria measured the closeness of the forecasts in comparison with the actual values of the variable being forecast, regardless of the direction of forecast errors. These are defined as follows:

MAFE = 1 |J | s∈J m s -m s , and RMSFE = 1 |J | s∈J (m s -m s ) 2 ,
where J the set used to assess the forecasts and |J | its size. Here, m s represents the actual holdout sample and m s the forecasts, see Figure 12 in Subsection B.2. In addition, we are also considering a measurement for the confidence intervals accuracy. As discussed above, we wanted to check if the outputs of the considered model encompasses sufficiently large amount of observed quantities. Also, the model produced wide confidence intervals and one could question the coherence of such a behavior. In fact, this can be a drawback if the forecast interval is too wide compared to the realized scenarios. To do this, we introduce a scoring rule for evaluating the pointwise interval forecast accuracy defined by [START_REF] Gneiting | Strictly proper scoring rules, prediction, and estimation[END_REF]. Therefore, as the above analysis, we consider the symmetric 95% prediction interval and we let m u s and m l s be, respectively, its lower and upper bounds. Then, letting s = (a, t), we define the scoring metric for a given forecast horizon t and age a as follows

IS (a,t) (95%) = m u s -m l s + 2 5% m l s -m s 1 ms< m l s + 2 5% ( m s -m u s ) 1 ms> m u s .
As noted by [START_REF] Shang | Point and interval forecasts of age-specific life expectancies: A model averaging approach[END_REF], the optimal interval score is achieved when the actual realized mortality rate lies between predicted bounds at the 95% level, with the distance between the upper and lower bounds being minimal. In order to summarize the above metric, defined per age and forecast horizon, we consider the average score over the considered ages for each horizon. First, on Table 3, we compare the point forecast errors between the benchmarks and the optimal AR-ARCH model discussed so far, both for mortality rates q a and life expectancy forecasts. We can observe that the AR-ARCH model outperforms the other models in terms of point and forecast accuracy as measured by the MAFE and RMSFE. This corroborates the the above discussion and confirms that the model can behave in a similar fashion, for some ages, but outperform in average both for life expectancy and probability rates. When comparing the interval forecast accuracy, see Figure 8, we observe that our model outperforms the others. In particular, the longer the forecast horizon, the larger the deviation between the AR-ARCH and the benchmarks, see also Figures 13 and 12 in Subsection B.2. Moreover, we can remark that the score remains stable over the forecast horizon unlike the others whose score diverge over the horizon. However, for the French population, our model is very close to the Lee and Carter. This could mainly due to the fact that the mortality dynamics for this particular population do not exhibit a clear cohort effect.

Concluding Remarks

In this paper, we addressed the problem of optimal and robust selection of neighborhoods for the random field models introduced in [START_REF] Doukhan | A Class of Random Field Memory Models for Mortality Forecasting[END_REF]. The latter is a parsimonious alternative approach for modeling the stochastic dynamics of mortality improvements. We have shown that a procedure based on the penalization of the QMLE is asymptotically consistent, in the sense that it selects the true model (neighborhoods) as far as the sample size goes infinity. This particular behavior depends on the form of the penalty and we show that a classic BIC criterion fulfills the conditions that amend to the consistency. This being shown, we explored numerically the procedure robustness and consistency for selecting the true model. For this, we draw conclusions on the efficiency of the BIC-based penalization using Monte-Carlo experiments. We implement the procedure in a simple case, where the maximal neighborhoods are predefined. Next, we investigate the procedure using real-world datasets of mortality rates. As shown in [START_REF] Doukhan | A Class of Random Field Memory Models for Mortality Forecasting[END_REF] using even a simple model specification of the AR-ARCH structure, the selected model for the considered population corroborate the conclusions in [START_REF] Doukhan | A Class of Random Field Memory Models for Mortality Forecasting[END_REF]. In fact, our model outperform the classic factor-based approaches even those with a cohort component.

Finally, we should note that the AR-ARCH structure can be implemented to tackle the problems encountered when using factor-based approaches in order to characterize the evolution of the residuals. Indeed, based on the work of [START_REF] Mavros | Stochastic Mortality Modelling: Key Drivers and Dependent Residuals[END_REF] which proposes a simple vector auto-regressive model (VAR) to model the residuals, one can hence consider the AR-ARCH as a parsimonious and coherent alternative for residuals modeling. Indeed, as noted throughout the paper, the latter does not impose any predefined structure on the dependency of the underlying random field and does tackle the problem of the heteroskedasticity. These were not considered by the VAR model used in [START_REF] Mavros | Stochastic Mortality Modelling: Key Drivers and Dependent Residuals[END_REF].

Other types of models may be considered and a reference for random processes is [START_REF] Doukhan | Stochastic models for time series[END_REF]. First we will work out the case of integer valued models of the GLM type as worked in [START_REF] Ferland | Integer?valued GARCH processes[END_REF], since this naturally corresponds to count individuals of a certain age at a given date. and for v and v in V 2

∂ 2 ∂β v ∂β v g(X s , (X s-v ) v∈V1 , (X s-v ) v∈V2 ; θ 0 ) = - X s-v X s-v α 0 0 + v∈V1 α 0 v X 2 s-v
.

Figure 3 :

 3 Figure 3: Autocorrelation. Autocorrelation of the random field and its squared values given in absolute value following the definition in Equation (3.7).

Figure 4 :

 4 Figure4: Model Selection. Estimation of the parameters for the optimal model maximizing the penalized QMLE. The grayed area corresponds to parameters that were not selected by the procedure and are excluded from the neighborhoods.

Figure 5 :

 5 Figure 5: Standarized Residuals. Heatmaps for the standarized residuals produced by the Lee-Carter (left panel) and the selected AR-ARCH model (right panel) for the French, England & Wales and United States male populations.

Figure 6 :Figure 7 :

 67 Figure6: Mortality Rates. Forecast mortality rates at age 65, 75 and 85 for the France, England & Wales and US populations produced by both the selected AR-ARCH (black) model andLee-Carter (red). The 95% prediction intervals are also included.

Figure 8 :

 8 Figure 8: Confidence Interval Forecast Score. The mean interval score across ages for different forecasting years for the LC, AR-ARCH and CBD two-factor models.

  Then, applying Theorem 5.7 of van der Vaart (1998),[START_REF] Doukhan | A Class of Random Field Memory Models for Mortality Forecasting[END_REF] demonstrate that the QMLE is consistent and asymptotically normally distributed, provided that the random field has a finite fourth moment, i.e. assumption H-4. Consistency. If the assumptions H-1, H-2 and H-3 hold, then the QMLE estimator θ T in Equation (3.2) is consistent, in the sense that for each (β v ) v∈V1 , (β v ) v∈V1 , (α v ) v∈V2 and (α v ) v∈V2 we have

	Theorem 2 (Doukhan et al. (2017)).
	(i) θ T	P -→ θ 0 .

since the set of possible parameters Θ is a compact set, i.e. assumption H-3.

(ii) Asymptotic normality. Under assumptions H-1, H-2, H-3 and H-4,

Table 2 :

 2 Normality tests. We perform normality tests of the residuals ξ s using the Shapiro-Wilk (SW) test, Anderson-Darling statistics (AD), the Cramér-von Mises test (CVM), Pearson-χ 2 (PC) test and the Shapiro-Francia (CF).

	SW	AD	CVM	PC	SF
	USA 0.3428 0.7048 0.7991 0.8523 0.3947
	FRA 0.6171 0.5212 0.5438 0.6279 0.4322
	GBP 0.6114 0.7626 0.8548 0.9309 0.4827

Table 3 :

 3 Point Forecast Accuracy. Forecast errors between the LC, AR-ARCH and the CBD (see[START_REF] Cairns | A two-factor model for stochastic mortality with parameter uncertainty: Theory and calibration[END_REF]) two-factor model both for mortality rates and life expectancies over the period2000-2016.

	AR-ARCH	LC	CBD

In our case, the execution time for implementing the neighborhood selection procedure on the maximal neighborhoods in Figure

using R was expected to take too much time. Therefore, we have solved this by considering parallel calculations and making use of external compiled languages to speed up the selection routine.

The data was downloaded from www.mortality.org on December 2016.

These stand, respectively, for the autocorrelation function and the partial autocorrelation function.

* The work is supported by the Institut des Actuaires (French Institute of Actuaries). P. Doukhan's work has been developed within the MME-DII center of excellence (ANR-11-LABEX-0023-01). Y. Salhi's work is supported by the BNP Paribas Cardif Chair DAMI "Data Analytics & Models for Insurance". The views expressed herein are the author's owns and do not reflect those endorsed by BNP Paribas.

Numerical and Empirical Analyses

In this section, we investigate the performance of the BIC based model selection. First, we run a Monte-Carlo experiment to assess the statistical properties developed in the previous section. Next, we consider a real-world dataset of mortality rates and explore the optimal neighborhood characterizing the mortality improvement rates dynamics. Hence, we forecast future evolution of mortality and compare these outputs to classic models discussed in Subsection 2.1.

4.1 Simulation Study. In this section, we examine the finite sample properties of the model selection using Monte-Carlo simulations. With no loss of generality, we generate data from the model in Equation (2.6) with i.i.d. r.v. ξ s ∼ N (0, 1), given the specification of the neighborhoods in Figure 2 and an arbitrary set of parameters that fulfill the stationarity constraint. The Monte-Carlo analysis serves as an illustration of the robustness of the QMLE estimation as well as the model selection. Using the notation of the previous sections, we investigate the impacts of varying the sample size T on the performance of these procedures.

Building on the previous considerations, we set the sample size T = I × J = 3000 (I = 30 ages and J = 100 years), 30 × 40 and use 1000 replications, we then study the performance of the selection procedure. In Table 1, we report the outputs of this analysis. The two panels show how many times (in frequency) the BIC selected each model out of 1000 independent simulations for different sizes of the dataset. The true model is presented in bold where the 1's correspond to the parameters belonging the true neighborhood and 0's are those excluded. We reported the percentages of samples in which each competing models provide a better penalized fit than the other competing models. The model selection method performed reasonably well in identifying the true model, though its ability to recover the true model (in bold) increases with the sample size. For small sized samples (lower part of Table 1), the model selection methods recovered at most 42.3%. When the sample size is large (upper part of Table 1), the model selection recovered at most 64.8%. We notice also that the number of competing models decreases with the sample size. Note that the consistency of the BIC selection criteria in Theorem 3 is achieved (the model selection method identifies (almost surely) the true neighborhoods) for big enough datasets, i.e. I × J → ∞, which is in line with the result shown in Table 1.

Table 1: Model Selection. Monte-Carlo assessment of relative performance of the selection method across data sizes based on 1000 replications.

In the light of the above remarks, a particular attention should be paid to the implementation of this model for mortality dynamics. In fact, to be able to select the true model (neighborhoods) we need to have relatively large datasets to achieve the consistency requirement. However, we are then left with a challenging trade-off between enlarging the observed age and year periods or losing the consistency of the selection procedure. The former can undoubtedly be a good compromise to achieve the identification (almost surely) of the true neighborhoods but may alter the relevance of the fitted model. Indeed, including very old mortality improvement Appendix A The gradient and the Hessian of the (quasi)-likelihood function g

The derivatives of the function g needed for the optimization problem can be derived explicitly, after some algebra, and are given as follows. These appear also in Equation (3.3). First, the different elements of the gradient ∂g/∂θ are given by:

Similarly, the elements of the Hessian matrix ∂ 2 g/∂θ 2 are given by:

and for v and

Appendix B Additional Figures Lee-Carter (red). The 95% prediction intervals are also included. The crude mortality over the observed period is also presented (black circles).


Jean Martinez
email: jean.martinez@umontpellier.fr

Jean-Alain Fehrentz

From JMV-1843 to Macrilen TM

Growth hormone deficiency (GHD) is a severe pathology that greatly affects the quality of life, and increases morbidity and mortality of patients owing to the augmentation of cardiovascular events. Treatment of GHD is challenging, mainly because there is no specific characteristic sign or symptom that can be used to make a clear diagnosis. There is need for an unequivocal diagnosis of GHD to avoid unnecessary treatment with GH, because the available provocative tests (GH stimulation tests) are not specific and sensitive enough, and are contraindicated in some patients.

Ghrelin is an endogenous peptide that stimulates GH secretion by interacting with a G-proteincoupled receptor named ghrelin receptor (GH secretagogue receptor 1a, GHS-R1a). Given this, a GH stimulation test using ghrelin or its analogues appears to be attractive. In this paper, a modified tripeptide first named JMV-1843 in the laboratory is briefly presented. It is potent and selective in stimulating the release of GH and is orally active. It has been recently commercialized for the diagnosis of adult growth hormone deficiency under the tradename Macrilen TM . The test using this compound appears to be reliable, well tolerated and simple.

Introduction

Adult growth hormone deficiency (AGHD) is usually caused by childhood GHD or by various conditions of the hypothalamus and pituitary gland resulting from tumours and their treatments, including surgery and irradiation [1]. AGHD results in severe morbidity associated with reduced quality of life, including reduced exercise tolerance as a consequence of fat mass increase, osteopenia, glucose intolerance and insulin resistance, as well as altered cardiovascular functions, which all decrease survival [2]. In the United States, approximately 60 000 people have been diagnosed with AGHD, while 6000 new cases are reported each year [3]. In Europe, about 51 000 adults are estimated to have been diagnosed with GHD [4]. However, because the symptoms are non-specific and there is no well-defined clinical marker of GHD in adults, a method to diagnose AGHD is necessary.

Results and Discussion

The detection of GHD is challenging, because growth hormone (GH) secretion occurs episodically in a pulsatile pattern, and is influenced by age, gender and body mass index [5,6]. Consequently, the diagnosis should be based on an appropriate provocative GH stimulation test (GHST) using a stimulus that can rapidly induce release of GH, which can be detected in blood samples [7].

There are six main tests for evaluating AGHD [8,9] including the insulin tolerance test (ITT), the combined arginine-GH-releasing hormone (ARG-GHRH) test and the glucagon test, which are used for adults. The ITT is considered to be the reference test, although it is limited to certain patients because of the risk of hypoglycaemia; patients must not have ischaemic heart disease or seizures [8,9]. However, all GHSTs have limitations in the diagnosis of AGHD concerning their sensitivity, specificity and/or reproducibility. In addition, they usually require parenteral administration by trained personnel, with a follow-up period in monitored facilities.

Several factors influence GH release. Among them, a potent 28-amino acid acylated peptide named ghrelin (Figure 1) produced mainly by the stomach, was isolated and identified by Kojima et al. in 1999 [10]. By binding to the 7-transmembrane G-protein-coupled receptor (GHS-R1a) [11], ghrelin induces GH secretion. The wide range of pharmacological and physiological properties of ghrelin have been described [12].

Since 1996, in collaboration with Dr R. Deghenghi, Europeptides (France) and several laboratories in Europe, including the Universities of Milano and Torino (Italy), we have been involved in a European programme (Peptido, Eureka n°1923) to design and synthesize peptides and mimetics able to stimulate GH release. The goal of this project was to identify an active compound as short as possible, selective, stable and orally effective, with possible therapeutic applications including the diagnostic test for detection of AGHD among GH-related disorders, and in some non-GH-related disorders associated with the wide range of biological activities of ghrelin [13,14].

In our laboratory, we have synthesized a series of analogues based on compound EP-51389, a potent GH-releasing peptide developed by R. Deghenghi at Europeptides (Figure 1) [15]. Based on our previous studies on gastrin [16] and bombesin [17], we have designed and synthesized a series of modified peptides, among them compound JMV-1843 (Figure 1), which, according to its pharmacological profile, was selected for further studies [18].

The chemical characteristics of JMV-1843 reside in its C-terminal formyl-amino function that mimics the C-terminal amide of EP-51389, which according to our previous studies on bombesin analogues should confer a certain degree of oral bioavailability [17], and on a gem-diamino tryptophan that improves enzymatic resistance [16].

Tested in competition binding assays using LLC PK-1 cells transiently overexpressing the cloned hGHS-R1a human receptor [18], human pituitary glands and hypothalamic membranes expressing the hGHS-R1a receptor [19], compound JMV-1843 displaced 125 I-labelled ghrelin in a dose-dependent manner and with high affinity. It exhibited approximately the same potency as ghrelin (in the 10 nM range) and the potent hexapeptide ghrelin agonist hexarelin [20], a wellcharacterized potent growth hormone secretagogue (GHS).

In vivo, JMV-1843 administered subcutaneously (s.c.) in rats was able to stimulate GH release dose dependently to the same extent as hexarelin. In dogs, JMV-1843 administered intravenously (i.v.) or given orally (1 mg/kg) stimulated significant GH release [18]. Compound JMV-1843 (then named EP1572) was then given i.v. (1.0 μg/kg as a bolus) and orally (0.06, 0.125, 0.25 and 0.5 mg/kg body weight) to two young male volunteers. At a dose of 0.06 mg/kg body weight given orally, JMV-1843 induced a remarkable and prompt increase in GH levels [19,21] that persisted for 3 hours. A further randomized, placebo-controlled study was conducted on 36 healthy volunteer male participants to evaluate both the pharmacological profile and GHreleasing activity of increasing single oral doses of JMV-1843 [22].

The results of this study indicated that oral administration of JMV-1843 (0.5 mg/kg) to healthy male volunteers resulted in a rapid and dose-dependent release of GH, with maximal GH release between 50 and 75 minutes, which lasted for >2 hours after administration. A significant GH release was detected even at a dose of 0.125 mg/kg. No adverse effects were observed and, in contrast to ghrelin, there was only a marginal effect on prolactin levels, while adrenocorticotropin (ACTH), cortisol, and ghrelin levels, as well as plasma concentrations of glucose and insulin, were not significantly affected.

Although far from being proven, this selectivity could be the consequence of the activation of different G proteins that modulate different intracellular signalling pathways and induce different biological responses and physiological processes in a given tissue. This specific G-protein activation could be the result of either homodimerization of the GHS-R1a or of its heterodimerization with another G-protein-coupled receptor expressed in the same tissue [23,[START_REF] Damian | Proc. Natl. Acad. 6 Sci[END_REF].

According to these preliminary results, Aetherna Zentaris selected compound JMV-1843 for further clinical investigation in 25 hospitals around the world: 377 participants were enrolled for the phase 3 clinical trials of oral administration of macimorelin (JMV-1843) as a diagnostic agent for AGHD.

Although the oral bioavailability of macimorelin is not high, the results of this multicentre study indicate that, at a dose of 0.5 mg/kg body weight, macimorelin was efficient in stimulating GH release in the blood to a level approximately 1.4-fold higher than with the ITT, and it was well tolerated. It was concluded that macimorelin was useful in the diagnosis of AGHD; usually after collection from the patients of four blood samples, at 30, 45, 60, and 90 minutes, to reduce the risk of false-positive results. There was no need of an extended follow-up period, and no need of parenteral administration [START_REF]Committee for Medicinal Products for Human Use (CHMP)[END_REF].

The most common adverse effects associated with macimorelin, which appeared in <5% of patients, were a change in the sense of taste, headache, diarrhoea, nausea and fatigue. During the initial phase 3 clinical trial, a drug-related alteration of cardiac rhythm (QT prolongation) was observed in one patient, which disappeared within 24 h. This participant had been taking citalopram, a drug known to be associated with QT prolongation. Accordingly, the use of macimorelin with drugs known to induce torsades de pointes, including antipsychotic products, should be contraindicated. This observation required a second phase 3 clinical investigation more particularly oriented towards this adverse effect. Finally, a single orally administered dose of macimorelin (usually 0.5 mg/kg) was recognized not to represent a particular risk and to be acceptable. In addition, because cytochrome P450 (CYP3A4) is involved in the metabolic degradation of macimorelin, one should be aware of the potential for false positives when diagnosing AGHD in patients treated with CYP3A4 inducers. In summary, these clinical trials demonstrated that orally available macimorelin was reliable in the diagnosis of AGHD, with minimal adverse effects when compared with the other tests available, and with an acceptable safety profile. Details of the phase 3 clinical trials are available at https://www.clinicaltrials.gov/ct2/show/NCT02558829?term=macimorelin. However, the use of macimorelin for young patients and for adults with extreme and morbid obesity remains to be demonstrated [START_REF] Gabreanu | [END_REF].

In December 2017, following the successful clinical trial results, the US Food and Drug Administration (FDA) granted their approval for commercialization of macimorelin in the USA and Canada, under the tradename Macrilen TM , for the diagnosis of AGHD [27]. In January 2019, the European Medicines Agency (EMA) also approved the marketing authorization for macimorelin for the diagnosis of AGHD.

A summary of the most significant events that occurred during the development of JMV-1843 to macimorelin and Macrilen TM are reported in Table 1.

Conclusion

Macimorelin (Macrilen TM ) is a novel active GH secretagogue that has an interesting pharmacological and clinical profile for diagnosis of AGHD, with several advantages compared with the established tests. Macimorelin appears to provide a simple and accurate diagnosis of AGHD, which has remained a challenge, and which is required for identification of patients that need GH treatment. The diagnosis of AGHD with macimorelin is highly competitive with the existing tests, which lack sensitivity and selectivity, and for which the parenteral administration route and the follow-up after administration of the drug are restrictive.

In addition, taking into account the range of biological activities of ghrelin, macimorelin can be envisioned in other clinical applications. These include treatment of growth retardation in children, and cachexia associated with chronic diseases such as AIDS and cancer. In this respect, a phase 1/2 study started recently, comprising an open-label clinical trial involving 24 children and adolescents under the age of 18 years. Macimorelin is also being used in a phase 2A clinical trial for the treatment of cancer-induced cachexia. However, whether or not macimorelin is beneficial in clinical applications other than the diagnosis of AGHD remains to be demonstrated.

October 31, 2018 :February 2019 :

 20182019 Novo Nordisk acquired the rights of Macrilen TM (Macimorelin) in the USA and Canada. November 2018: The Committee for Medicinal Products for Human Use (CHMP) gave a positive opinion of Macimorelin (Macrilen TM ). January 2019: The EMA granted marketing authorization for Macimorelin (Macrilen TM ). Phase 1/2 paediatric clinical trial initiated with Macrilen ® (Macimorelin) by Aeterna Zentaris GmbH and Novo Nordisk.

Figure 1 .

 1 Figure 1. Structure of ghrelin, hexarelin, EP-51389 and JMV-1843

Acknowledgments

Studies performed at the University of Montpellier were supported by CNRS, University of Montpellier, EU Eureka Project (Peptido Project 1923), Europeptides and Aetherna Zentaris. The authors would also like to thank all the individuals that were involved in the JMV-1843 project allowing the compound to reach the market.

Conflict of interests

The authors declare no conflicts of interest.

TOC

JMV-1843 (tradename Macrilen TM ) is a modified tripeptide, which stimulates the release of growth hormone (GH) by interacting with the ghrelin receptor. Macrilen TM is the first orally available compound commercialized for diagnosing adult growth hormone deficiency (EGHD). The test appears to be reliable, well tolerated and simple. The story of JMV-1843 in the laboratory to Macrilen TM in the market is summarized.