Frédéric Dehais
email: frederic.dehais@isae.fr

Kevin Mandrick

Zarrin Chua

Mickaël Causse

Stéphane Perrey

Why a Comprehensive Understanding of Mental Workload through the Measurement of Neurovascular Coupling Is a Key Issue for Neuroergonomics?

mental workload, mental resources, neurovascular coupling, neuroergonomics, electroencephalography, near-infrared spectroscopy Raja Parasuraman, the father of Neuroergonomics (the crossroads of Ergonomics and Neuroscience, Figure 1) has opened the doors to new discoveries and techniques for advancing understanding of human behavior with the underlying brain mechanisms [START_REF] Parasuraman | Neuroergonomics: The Study of Brain and Behavior at Work[END_REF]. As of his death in 2015, a precise and objective definition of the concept of mental workload (MWL) had still not yet been formulated. In this opinion piece, we posit that MWL is associated through the measurement of neurovascular coupling (NVC); innovative neuroimaging methods is now capable of measuring such a phenomenon; all while highlighting Parasuraman's many contributions to this field.

BEYOND THE CONCEPT OF MENTAL WORKLOAD AND TOWARD MENTAL RESOURCES IN NEUROERGONOMICS

MWL measurement is an important issue in the Human Factors field, as seen through its ubiquitous presence in the literature. It is well acknowledged that an accurate assessment of MWL could help to reduce human error while improving human performance. The recently founded field of Neuroergonomics may help to reduce the ambiguity surrounding the MWL concept by providing data on its underlying neural processes. Neuroergonomics allows for the study of the human brain structure and function with respect to behavior during physical or cognitive performances in the workplace [START_REF] Mehta | Neuroergonomics: a review of applications to physical and cognitive work[END_REF]. The main goal of this interdisciplinary field is to integrate our understanding of the neural basis of cognition in relation to technologies and settings in complex daily life tasks.

However, Neuroergonomics does not yet provide a consensual and comprehensive explanation of the MWL. Despite being a roughly defined concept, there have been some formal attempts. Generally, MWL reflects how hard one's mind is working (under-over-loaded or occupied) at any given moment or how much mental effort it will cost for brain to meet given task demands [START_REF] Parasuraman | Neuroergonomics: research and practice[END_REF]. Furthermore, [START_REF] Parasuraman | Mental workload[END_REF] and Kramer and [START_REF] Karanasiou | Functional brain imaging using non-invasive non-ionizing methods: towards multimodal and multiscale imaging[END_REF] defined MWL as a set of mental and composite brain states that modulate human performance in different perceptual, cognitive, and/or sensorimotor skills. It is also considered as a construct used to reflect the relation between the demands of the environment (input load), the human characteristics (capacities), and the task performances (output performance). However, the notion of MWL is dissociated from performance as suggested [START_REF] Parasuraman | Neuroergonomics: The Study of Brain and Behavior at Work[END_REF][START_REF] Parasuraman | Neuroergonomics: research and practice[END_REF] has therefore coined the term "neuroergonomics" to refer to the inclusion of neuroscience in Human Factors/Ergonomics (HF/E). Neuroergonomics can be defined as the study of brain and behavior at work. Traditionally, ergonomics has not paid much attention to neuroscience or to the results of studies of the brain mechanisms underlying human perceptual, cognitive, affective and motor processes. To paraphrase the philosopher Mario Bunge (1980), until recently psychology (and HF/E) has been "brainless," whereas neuroscience has been "mindless." At the same time, neuroscience and its more recent offshoot, cognitive neuroscience, have been only partially concerned with whether their findings bear any relation to human functioning in real (as opposed to laboratory) settings, with the exception of applications to clinical disorders. Neuroergonomics is a response to this twin disregard.

by [START_REF] Ayaz | Optical brain monitoring for operator training and mental workload assessment[END_REF]. MWL presupposes that the consumption of true brain resources supports brain activity during work, suggesting a possible link between MWL and the key concept of mental resources. These two concepts can be treated by the intensity of the mental costs and be measured by the mental effort of performing tasks to predict operator performance. As stated by [START_REF] Cain | RTO-TR-HFM-121-Part-II-Human Factors Considerations in the Design, Use, and Evaluation of AMVE-Technology[END_REF] "As such, [MWL] is an interim measure and one that should provide insight into where increased task demands." Therefore, it is not possible to define MWL without also clearly characterizing mental resources.

Though it is generally admitted that mental resources are appreciable, multiple, independent, and limited [START_REF] Wickens | Multiple resources and mental workload[END_REF], most studies remain vague on their exact nature. One perspective is to think of mental resources as neural pathways. However, this oversimplification ignores the fact that mental resources exists in other forms. As a metaphor, an army may have efficient firepower, but without ammunition, a supply corps, and roads, it is useless. Similarly, the army of the brain has mental resources composed of neural pathways, energy supply, and irrigation (communication channel) to fuel mental effort, implemented by the mobilization of neurophysiological cellular processes in the operator's brain.

ENERGY MOBILIZATION OF NEUROVASCULAR COUPLING FOR THE OPERATOR'S BRAIN MACHINERY

The absence of consideration of the neurophysiological mechanisms in Neuroergonomics is certainly due to the difficulty in investigating them. Yet, there are real energy mobilizations that occur within the operator's brain machinery across several cellular levels to meet task demand. As previously compared to a super calculator or a computer, the brain machinery supports mental processors that need substantial and constant energy requirements. But the human brain is devoid of intra-cellular capacity for energy storage in oxygen, lactate, and glucose (even if small parts of glycogen exist). Fortunately, the demand for high-metabolic energy of the brain tissue is mainly regulated by complex but adequate energetic substrate delivery via a dense and redundant network of microvessels. Hence, metabolic demands are orchestrated by the blood supply hemodynamic response.

Since the first discoveries by [START_REF] Roy | On the regulation of the blood supply of the brain[END_REF], it has been possible to better understand the close spatiotemporal dynamics between the electrical activity of neuronal cells and the hemodynamic phenomenon that boost the local bloodstream circulation in localized arterioles and capillaries. The intimate neurofunctional relationship that concomitantly links the metabolically active neurons with the increasing oxygenation of the blood flow near of these cells reflects the functional hyperemia and is more widely known as neurovascular coupling (NVC). Simply, NVC is a tight temporal association of the neuronal activity with regional cerebral blood flow delivery. Understanding the fundamental cellular mechanisms underlying NVC is necessary to measure a dimension of the local brain machinery expenditure at work. The appraisal of the energetic costs required by NVC implies the assessment of mental resources. For instance, when an operator is engaged in a task, the mobilization of the neural pathways needs a synergistic support of massive astrocyte glial cells to fuel neurons and interneurons with oxygen and nutriments furnished by close capillaries.

NVC is observable due to changes in neuronal-astroglial and microvasculature activities, which occur in several steps. First, the measurable electrical neuronal activity (spiking and postsynaptic potential activity) is accompanied by synaptic neurotransmitter release (glutamate, GABA) with a neuronal-astroglial regional cerebral metabolic rate of oxygen consumption, mainly for regional cerebral metabolic rate of glucose demand. Second, this activity induces a cascading pathway involving the production and the release of powerful vasodilator metabolites by neurons and astrocytes and drives a chemical signal up to the vascular smooth muscle and pericytes cells along the microvessels which dilate the microvasculature. Third, the microvessels dilatation significantly modulates the regional cerebral blood activity (flow, volume, and oxygenation) which greatly exceeds the neuronalastroglial oxygen requirements, and results in a measurable overabundance of blood flow, hence, a local hyperoxygenation. Yet, the role of NVC as it contributes to the comprehension of the energy mobilization in response to mental resources is not common knowledge. The cellular measures of energy production, delivery, and utilization are crucial to understanding and interpreting NVC activity. How to clearly establish the role of NVC into the operator's brain machinery? One possible way would be to associate the level of correlates of NVC while interpreting the degree of task demand. It seems thus fairly possible that an accurate measurement of NVC, spatially and temporally and in terms of amplitude, would be a valuable neurophysiological marker for quantifying changes in brain activation. Although this statement is still reductionist (that NVC activity is proportional to operator's brain activity), this approach links the concept of human MWL and mental resources to objective neurophysiological measures for Neuroergonomics purposes.

Recent Neuroergonomics research has progressed in neurocognitive or neuroimaging-sensing instrumentation for determining operator states through the measurement of NVC activity associated with the degree of mental processes [START_REF] Parasuraman | Putting the brain to work: neuroergonomics past, present, and future[END_REF]. Tremendous advances have been made toward establishing approaches for portable neuroimaging equipment and brain activation measurements to assess sensitivity to NVC in human operators acting in realistic work environments. This development is especially the case in ambulatory functional neuroimaging methods such as functional near-infrared spectroscopy (fNIRS) and electroencephalography (EEG). To date, the aforementioned non-invasive brain imaging techniques are beginning to be well-established in the Neuroergonomics community. These advantages will be even more beneficial in the future as the coupling between these methods becomes more widespread.

ASSESSING NEUROVASCULAR COUPLING WITH FNIRS-EEG METHODS: AN OBJECTIVE NEUROERGONOMICS APPROACH FOR EVALUATION OF THE OPERATOR'S BRAIN ACTIVITY

Technological advances in opto-and electronic miniaturization have improved the portability and operational flexibility in brain imaging sensors, allowing for greater comprehension of the brain at work in real-world applications (aeronautics, automotive, robotics). fNIRS provides a continuous monitoring of the hemodynamic activity using near-infrared light transmitted between optodes. It infers the changes in the concentrations of oxygenated and deoxygenated hemoglobin in the cortical regions from scattering and absorption properties of light probing beneath the surface of the skull [START_REF] Perrey | Non-invasive NIR spectroscopy of human brain function during exercise[END_REF]. These two fNIRS signals have their origins in the metabolic response corresponding to a shift of oxygen consumed and the vascular response linked to a modulation of the microvasculature activity (dilatation). This hemodynamic response disrupts the regional cerebral blood flow and volume which exceeds oxygen intake (functional hyperemia) consumed by the recruited neuronal population. fNIRS responses characterize the operator's brain activity related to cerebral blood flow and cerebral tissue oxygenation changes over time [START_REF] Mandrick | Prefrontal cortex activity during motor tasks with additional mental load requiring attentional demand: a near-infrared spectroscopy study[END_REF][START_REF] Durantin | Using near infrared spectroscopy and heart rate variability to detect mental overload[END_REF][START_REF] Fishburn | Sensitivity of fNIRS to cognitive state and load[END_REF]. Good spatial localization can be derived if a high number of optodes are used in an array, but temporal resolution is coarse by the delayed nature of the hemodynamic response to cortical activity (few seconds).

On the other hand, EEG offers a fine temporal resolution (milliseconds) thus enabling detection of brief neuronal processes, but is limited in its capacity for spatial resolution, at least in real time even though dense array EEG permits source propagation localization. EEG uses scalp electrodes to capture weak electrical current fluctuations generated by inhibitory or excitatory postsynaptic potentials of a pool of neurons firing simultaneously in response to a stimulus. The electrophysiological roots of these signals correspond to the summation of the spontaneously and synchronously recruited neuronal population that contributes to the neuronal activity of the superficial layers of the cortex. EEG waves and eventrelated potentials signals are particularly strong candidates for objective measures of operator's brain activity at the workplace [START_REF] Parasuraman | Neuroergonomics: The Brain at Work[END_REF]. In general, fNIRS and EEG are complementary as they improve on each other's measurement weaknesses in terms of information content [START_REF] Fazli | Enhanced performance by a hybrid NIRS-EEG brain computer interface[END_REF]. Additionally, there is no noise cross-interference between fNIRS and EEG (light and electrical, respectively; [START_REF] Karanasiou | Functional brain imaging using non-invasive non-ionizing methods: towards multimodal and multiscale imaging[END_REF]. Therefore, simultaneous fNIRS-EEG signal acquisition would be suitable for assessing NVC in order to evaluate the operator's brain activity in ecological contexts [START_REF] Hirshfield | Combining electroencephalograph and functional near infrared spectroscopy to explore users' mental workload[END_REF][START_REF] Safaie | Toward a fully integrated wireless wearable EEG-NIRS bimodal acquisition system[END_REF]. However, it not should limit our understanding of the brain activity to only one perspective; looking at the brain at work with new tools and new eyes we could have new NVC comprehension during ecological context. Readers must note that the multimodality using fNIRS-EEG methods is a very promising approach in the investigation of where, when, and how much NVC exhibits energy mobilization during work. The spatiotemporal evolution of the functional neural connectivity and blood flow regulation through the scalp is permitted due to the recording of temporal electrical activity and spatial hemodynamic activity. Consequently, the evaluation of NVC distribution throughout the head becomes accessible. This measurement makes it possible to dynamically map the brain activity and identify the brain areas with the activated main NVC. Additionally, the assessment of the power of the electrical signal by EEG coupled with the amplitude of the hemodynamic signal by fNIRS will enable a better depiction of the intensity of the NVC, thus extrapolating the effectiveness of the metabolic effort of performing tasks. This view of the degree of extrapolated metabolic correlates as an indicator of the level of mental resources seems straightforward at first glance. However, the metabolic expenditure that fuels cognitive processes is the prerequisite for any mental resources and the assessment of operator's brain activity. The challenge now is to enhance the reliability of NVC measurement in situ with fNIRS-EEG methods.

THE FUTURE FOR NEUROERGONOMICS

It is clear that the extensive work of Parasuraman has left the scientific community in an excellent position to objectively define MWL and subsequently, mental resources, through the measurement of NVC activity. It is our opinion that NVC measurement could be achieved through the use of an efficient fNIRS-EEG coupling. In particular, there needs to be greater characterization of the energy mobilization of NVC with respect to neurophysiological mechanisms (neuronal-astroglial, metabolic and hemodynamic activity) and methods for its assessment in work settings [START_REF] Parasuraman | Neuroergonomics brain, cognition, and performance at work[END_REF]. There rests a great deal of work in Neuroergonomics before the development of a standard assessment approach of NVC with innovative neuroimaging technology for the evaluation of the operator's brain activity at work. In other words, there are still opportunities for the technological deployment of coupled hybrid devices (dry-electrodes EEG within a high density headset of fNIRS optodes). From a broader perspective, emerging devices must meet several criteria: discriminate different levels of workload; not interfere with the subject's work and environment; be accepted by the individual; be low cost with high portability; be easy to implement and to evaluate; be reproducible and reliable; and dissociate the mental workload from emotional processes (sensitivity and specificity). Theoretically, a multimodal fNIRS-EEG approach should help to investigate the interactions between different mental states and user behavior while taking into account the physiological processes. Further investigations are warranted to address newer assessments of the neurophysiological events of the operator's brain at work.

FIGURE 1 |

 1 FIGURE 1 | Illustration of Raja Parasuraman as he himself wrote in an article: "One of us[START_REF] Parasuraman | Neuroergonomics: The Study of Brain and Behavior at Work[END_REF][START_REF] Parasuraman | Neuroergonomics: research and practice[END_REF] has therefore coined the term "neuroergonomics" to refer to the inclusion of neuroscience in Human Factors/Ergonomics (HF/E). Neuroergonomics can be defined as the study of brain and behavior at work. Traditionally, ergonomics has not paid much attention to neuroscience or to the results of studies of the brain mechanisms underlying human perceptual, cognitive, affective and motor processes. To paraphrase the philosopher Mario Bunge (1980), until recently psychology (and HF/E) has been "brainless," whereas neuroscience has been "mindless." At the same time, neuroscience and its more recent offshoot, cognitive neuroscience, have been only partially concerned with whether their findings bear any relation to human functioning in real (as opposed to laboratory) settings, with the exception of applications to clinical disorders. Neuroergonomics is a response to this twin disregard.

FUNDING

This work was funded by the French Research National Agency, the French Defence Procurement Agency (ASTRID), and the AXA Research Fund.

AUTHOR CONTRIBUTIONS

Each of the authors has read and concurs with the content in the final manuscript. The first author (KM) wrote the majority of the manuscript. The other authors (ZC, MC, SP, and FD) have extensively reviewed and revised the manuscript from the first draft before giving final approval of the version to be submitted.

Conflict of Interest Statement:

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.


Yahia Salhi

Pierre-E Thérond

Age-Specific Adjustment of Graduated Mortality *

Keywords: Credibility, Mortality, Life Insurance, Graduation, Smoothing, Local Likelihood, Prediction

Recently, there has been an increasing interest of life insurers to assess their portfolios own mortality risk. The new European prudential regulation, namely Solvency II, emphasized the need to use mortality and life tables that best capture and reflect the experienced mortality, and thus policyholders' proper risk profile, in order to adequately quantify the underlying risk. Therefore, building a mortality table based on the experience from the portfolio is highly recommended and, for this purpose, various approaches have been introduced in the literature. Although, such approaches succeed in capturing the main feature, it remains difficult to assess the mortality when the underlying portfolio lacks of sufficient exposure.

In this paper, we propose to graduate the mortality curve using an adaptive procedure based on the local likelihood, which has the ability to model the mortality patterns even in presence of complex structures and avoid to rely on experts opinion. However, such a technique fails at proposing a consistent yet regular structure when for portfolios with limited deaths. Although the technique borrows the information from the adjacent ages, it is sometimes not sufficient to produce a robust life tables. In presence of such a bias, we propose to adjust the corresponding curve, at the age level, based on a credibility approach. This consists on reviewing, as new observations arrive, the assumption on the mortality curve.

We derive the updating procedure and investigate the benefits of using the latter instead of a sole graduation based on real datasets. Moreover, we look at the divergences in the mortality forecasts generated by the classical credibility approaches including Hardy-Panjer, the Poisson-Gamma model and Makeham framework on portfolios originating from various French insurance companies.

Introduction

For insurers, the assessment of the experienced mortality is of paramount importance. The new regulation and norms, Solvency II, shed light on the need of life tables that best reflect the experience of insured portfolios own mortality so as to reliably quantify the underlying risk. Insurers, in France for example, are used to rely on regulatory life tables for pricing and reserving purposes, which are sometimes too conservative. The use of inadequate life tables, being too conservative would considerably affect the financial profitability of the life insurance businesses as well as the insurer competitiveness. Among others, from a Solvency II perspective, using too conservative tables leads to two major impacts: (i) an increase of Best Estimate technical provisions (and thus a decrease of Basic Own-Funds); (ii) an increase of the base figure used for calculating the capital charge for mortality risk.

A natural and straightforward approach to handle this issue is to use the available data at the portfolio level and build an entity-specific mortality table. However, practitioners may face technical difficulties related to the size of the portfolios and the heterogeneity of the guarantees (for the same underlying risk). For instance, an insurer may have a fairly big portfolio but with policyholders holding different insurance contracts: pure endowment contracts, unit-linked contracts with minimum death guarantees, loan insurance and so on. In such a case, it is difficult to build a mortality table based on the sole experience of each product or guarantee. More precisely, the constructed table may not be able to represent the mortality profile of the policyholders thus failing in capturing the underlying risk profile. This should also be the case even if the mortality table is periodically updated with the incoming new data. Especially, the latter may induce significant impacts on the technical reserves if the table has to be updated more frequently over time. Also, if one draws a mortality table only based on the experience stemming from one product or guarantee, she shall have the face a problem of heterogeneity. The latter arises not only at the portfolio level but also for individual ages. In fact, the mortality profile is highly dependent on the age of the individuals and some age groups being poorly represented may alter the quantification of the mortality risk at each individual age.

In this paper, we consider an insurer with exposures to different policies and aiming at establishing an experience-based mortality table for each policy and age level, as individuals may differ substantially in their endowment for biometric risk (as showed by some empirical mortality studies, e.g. see [START_REF] Vaupel | The impact of heterogeneity in individual frailty on the dynamics of mortality[END_REF] and [START_REF] Hougaard | Life table methods for heterogeneous populations: distributions describing the heterogeneity[END_REF] among others). As a first step, we consider a graduation principle to build mortality rates at the insured portfolio level. There are usually two sorts of methods: non-parametric and parametric, see Forfar et al. (1988a) and [START_REF] Debón | A comparison of nonparametric methods in the graduation of mortality: Application to Data from the Valencia Region (Spain)[END_REF] for a comprehensive introduction to the use of both graduation techniques. The nonparametric framework is very useful in practice especially when there is sufficient data. This method relies on the use of kernel estimation techniques which were first used for graduation by [START_REF] Copas | Non-parametric graduation using kernel methods[END_REF] and [START_REF] Ramlau-Hansen | The choice of a kernel function in the graduation of counting process intensities[END_REF]. In the current work, we use local kernel-weighted log-likelihood techniques to graduate the observed mortality with attained age. It was considered in [START_REF] Tomas | A local likelihood approach to univariate graduation of mortality[END_REF]Tomas ( , 2012)), which introduced adaptive smoothing parameters choice with an application in life insurance. Local likelihood techniques have the ability to model the mortality patterns even in presence of complex structures and avoid to rely on experts opinion.

The graduated mortality can be then used as such to project future insurance liabilities related to the underlying population. However, the evolution of the flow of data related to latest available information is not taken into account. This should be, for example, used to update the graduated mortality. However, if one decides to redo a graduation procedure including the new data, the forecasts are likely to be unstable; adding potential volatility to the underlying reserves and capital charges.

Therefore, the primary contribution of this paper is the incorporation of heterogeneity into the graduated mortality table model by introducing an unobserved variate for individual differences in each attained age. Such an approach has been considered in [START_REF] Salhi | A credibility approach of the Makeham mortality law[END_REF] but with different graduated curve. The latter uses a parametric model, i.e. Makeham law, to first built the mortality curve and then applies a credibility to the portfolio-sensitive parameter. Other approach have been also introduced in the literature but work directly on the death counts, e.g. [START_REF] Hardy | A credibility approach to mortality risk[END_REF] and [START_REF] Bühlmann | A course in credibility theory and its applications[END_REF] among others. Unlike the classical approaches that focus on the update of the aggregate deaths recorded over the whole portfolio, the proposed adjustment approach is intended to enhance the predictive ability of the graduated mortality using a credibility-based revision at the age-level and not on the aggregate portfolio level. By doing so, we implicitly account for the heterogeneity observed at the age-level and has beneficial outputs in pricing and reserving fields when considering the insurance guarantees according to the age.

This note is organized as follows. Section 2 has still an introductory purpose. It specifies the notation and assumptions used in the following. We also describe the graduation procedure based on the local-likelihood approach and recall some statistical inference results used in the sequel. Section 3 introduces the credibility approach to the graduated mortality mortality. We specify the model and make connection with the recent literature and derive the main tools needed to fully characterize the next period prediction of mortality rates when a (multiplicative) credibility factor is taken into account. Finally, Section 4 presents an application with experience data originating from French insurance companies. Finally, some remarks in Section 5 conclude the paper.

Notation, Assumptions and Preliminaries

2.1 Notation and Assumptions. We suppose that we have at our disposal age-specific mortality statistics originating from K portfolios (or companies) over the period [1, T i ], with i ∈ {1, • • • , K}. For each portfolio i, we let L i x,t denote the number of individuals at attained age x during calendar year t and D i x,t represents the number of deaths recorded from an exposure-to-risk E i x,t measuring the time during which individuals are exposed to the risk of dying. To each of the observations j, we associate the dummy variable δ j indicating if the individual j dies or not over one year. That is, for each j = 1, . . . , L i x,t , and age-band [x, x + 1], δ j = 1 if individual j dies and null otherwise. Next, we define the time lived by individual j before her (x + 1) th birthday by τ j and assume that we have at our disposal i.i.d. observations (δ j , τ j ) for each of the L i t,x individuals. Then, it is straightforward that:

L i x,t j=1 τ j = E i x,t and L i x,t j=1 δ j = D i x,t , for i = 1, . . . , K .
Furthermore, we assume that the age-specific forces of mortality are constant within bands of time and age, but allowed to vary from one band to the next, ϕ i x+τ ,t+ξ = ϕ i x,t for 0 ≤ τ < 1 and 0 ≤ ξ < 1. We also denote by p i x,t the probability that an individual aged x in calendar year t originating from company i reaches age x + 1, and by q i

x,t = 1 -p i x,t the corresponding probability of death. Under the assumption of piecewise constant forces of mortality, we have for integer ages x and calendar years t, p i x,t = exp -ϕ i x,t and ϕ i x,t = -log p i x,t Moreover, the likelihood can be expressed in the following form:

L ϕ i x,t = exp -E i x,t ϕ i x,t ϕ i x,t D i x,t .
The associated log-likelihood is then given by

ϕ i x,t = log L ϕ i x,t = -E i x,t ϕ i x,t + D i x,t log ϕ i x,t .
Maximizing the log-likelihood ϕ i x,t gives ϕ i x,t = D i x,t /E i x,t which coincides with the central death rates m i x,t . Then it is apparent that the likelihood ϕ i x,t is proportional to the Poisson likelihood based on the following assumption

D i x,t ∼ P E i x,t ϕ i x,t . (H0)
Thus, it is equivalent to work on the basis of the true likelihood or on the basis of the Poisson likelihood as recalled in [START_REF] Delwarde | Construction de tables de mortalité périodiques et prospectives[END_REF]. Accordingly, under the assumption of constant forces of mortality between integer values of x and t, we consider, henceforth, (H0) to take advantage of the Generalized Linear Models (GLMs) framework.

In the sequel we frequently use the following notation.

Let D i x • = Ti t=1 D i x,t and E i x • = Ti t=1 E i x,t
be respectively the company specific time aggregated death counts and exposures. For convenience, we write

D i = (D i x1 • , • • • , D i xn • ) and E i = (E i x1 • , • • • , E i xn •
) the vectors of aggregated deaths and exposures for company i and also denote

ϕ i = (ϕ i x1• , • • • , ϕ i xn•
) . Thus, the bullet indicates the summation over the corresponding index in the notation introduced above.

Local Likelihood

Smoothing of Mortality Tables. The crude estimates of mortality for each company i might be subject to random fluctuations. Typically, the latter may arise when the underlying exposure to risk is small. Therefore, the common practice in life insurance is whether to graduate mortality using a parametric approach or use some non-parametric smoothing techniques, see Forfar et al. (1988a) and [START_REF] Debón | A comparison of nonparametric methods in the graduation of mortality: Application to Data from the Valencia Region (Spain)[END_REF] for a comprehensive introduction to the use of graduation techniques. Here, we consider a non-parametric graduation method of experienced data originating from each company following the approach introduced in Tomas and Planchet (2013a,b). We use the hypothesis of Poisson distributed deaths as in (H0) and assume that the intensity ϕ i x,t of the portfolio i is related to a reference ϕ ref x,t as follows:

ϕ i x,t = exp(f (x))ϕ ref x,t , (1)
where f is an unspecified, smooth and deterministic function of the age x. The latter allows to link the mortality of the company i to the baseline at the attained age level x. We define an estimator ϕ i x,t of ϕ i x,t based on the data from the i'th subgroup. To do so, a non-parametric estimation of f is needed as developed in Tomas and Planchet (2013a,b). More formally, we consider the local kernel-weighted log-likelihood method to estimate the smooth function f (x). Such an approach was used to graduate life tables with attained age context in [START_REF] Delwarde | Modèles linéaires et additifs géneralisés, maximum de vraisemblance local et méthodes relationelles en assurance sur la vie[END_REF], [START_REF] Debón | A comparison of nonparametric methods in the graduation of mortality: Application to Data from the Valencia Region (Spain)[END_REF] and [START_REF] Tomas | A local likelihood approach to univariate graduation of mortality[END_REF]. One of the key advantages of such an approach is its ability to model the mortality patterns even in the presence of complex structures and avoid to rely on experts opinion in graduating mortality, thanks to the flexible functional form of the regression curve. The starting point of our estimation is the specification of the GLM under the assumption (H0) that allows to write the first moment of D i x,t as follows:

log E[D i x,t ] = log E i x,t + log ϕ i x,t = log E i x,t + log ϕ ref x,t + f (x), (2)
where the term log E i x,t is an offset. Such a model has been widely considered in the literature. For example, in [START_REF] Currie | Fitting models of mortality with generalized linear and non-linear models[END_REF], the function f has the parametric form f (x) = β 0 + β 1 x for some unknown parameters β 0 and β 1 . In our case, there is no model assumption about f . The natural extension of this idea is to suppose that for each age x k , with k = 1, • • • , n, f (x k ) is a p th degree polynomial in x j 's, where x j is an element in the neighborhood of x k . Formally, denoting

x j = 1, x k -x j , • • • , (x k -x j ) p and β i = β 0 , • • • , β p
we can write f (x k ) in the following form f (x k ) = x j β i . Next, to estimate β i , we maximize the likelihood of the model. To do this, we add a weight on the observations at each age. Following Tomas (2012), we use a non-negative weight function, i.e. w j = w j (x k ), that depends on the distance between the observations and the fitting point x j , i.e.

ω j =    ω j (|x j -x k |/h), if |x j -x k |/h ≤ 1, 0, otherwise, (3)
where h is a smoothing parameter determining the radius of the neighborhood of x k used in the smoothing. It gives the bandwidth of the neighbor used in the kernel. For instance, the smallest h the thiner is the voisinage that contributes to the likelihood at each attained age. Examples of the weighing function ω j are given in Tomas (2011, Table 1). In our case we consider a Gaussian kernel and use these weights on the characterization of the likelihood. Accordingly, it follows that the log-likelihood, for portfolio i, under the assumption (H0) writes:

L β i |D i , E i , (ω j ) j = Ti t=1 n j=1 ω j D i xj ,t x j β i -E i xj ,t ϕ ref xj e x j β i -log D i xj ,t ! + D i xj ,t log E i xj ,t ϕ ref xj , = n j=1 ω j D i xj ,• x j β i -E i xj ,• ϕ ref xj e x j β i + C i ,
where C i is a constant offset which does not depend on the parameters β i . Hence, taking the derivative with respect to each β i (p), for p = 1, • • • , l, leads to the following system of equations,

n j=1 ω j (x j -x i ) p D i xj ,• -E i xj ,• ϕ ref xj e x j β i = 0,
which can be rewritten in a matrix notation as follows:

(X i ) W i (D i -ϕ i ) = 0, (4)
with ϕ i being non-linearly dependent on β i , X i is the matrix

X i =       1 x 1 -x i (x 1 -x i ) 2 • • • (x 1 -x i ) n 1 x 2 -x i (x 2 -x i ) 2 • • • (x 2 -x i ) n . . . . . . . . . . . . . . . 1 x n -x i (x n -x i ) 2 • • • (x n -x i ) n       , (5)
and W i is the weight matrix with elements w k,j given as w k,j = ω j for k = 1, • • • , n. The solution of the above equation, i.e. estimations, must be obtained numerically using an iterative algorithm as Nelder-Wedderburn, Newton-Raphson algorithms or the fisher scoring methodology, see [START_REF] Tomas | A local likelihood approach to univariate graduation of mortality[END_REF] for further development. From these, we can get the estimation of β i and ϕ i denoted, henceforth, by β i and ϕ i 2.3 Inference of the Graduated Mortality. The aim of this subsection is to characterize the statistical feature of the estimators considered above. We recall some results from [START_REF] Tomas | A local likelihood approach to univariate graduation of mortality[END_REF] well known in the nonparametric smoothing literature, see e.g. [START_REF] Tibshirani | Local likelihood estimation[END_REF] and [START_REF] Wand | Kernel smoothing[END_REF], regarding particularly the variance of the graduated mortality and the expected behavior of these estimations. Accordingly, using theoretical results concerning bias and variance the estimator ϕ i is shown to be asymptotically design-unbiased and consistent. It is, for instance, recalled in [START_REF] Tomas | A local likelihood approach to univariate graduation of mortality[END_REF] that the smoothed mortality rates ϕ i are unbiased estimators of ϕ i in the sense that:

E[ ϕ i ] ≈ ϕ i . ( 6
)
The above approximation is based on the inspection of the mean squared errors, which are commonly used to assess the bias of the estimation in such a framework. Expressions of the latter are available in the classical textbooks and the readers is referred to (Tomas 2011, Sec. 3.2.) which provides an approximation to the bias of the estimator ϕ i . Unlike the linear model fitting, there is no exact expression for moments of ϕ i x due to the non-linearity in Equation ( 4). Using a multivariate version of Taylor series expansion around β i allows to use classical results on the inference of the estimated parameter β i . Note that this approximation depends on the bandwidth of the neighborhood h used in the kernel. More precisely, the bias is decreasing with the bandwidth. This is, particularly, reasonable in practice, because a large bandwidth induces a miss-fitting of the local polynomials and hence also the sum of squared residuals. To derive the second order moment of ϕ i , a variance approximation based on Taylor linearization is also generally suggested and shown to be consistent, see [START_REF] Loader | Local regression and likelihood[END_REF]. More precisely, we have the following expression for the variance:

Var( ϕ i ) = ϕ i ⊗ S i S i , (7)
with " ⊗ " being the pointwise division operator. Here, the matrix S i is given as follows:

S i =       s i 1 (x 1 ) s i 2 (x 1 ) s i 3 (x 1 ) • • • s i n (x 1 ) s i 1 (x 2 ) s i 2 (x 2 ) s i 3 (x 2 ) • • • s i n (x 2 ) . . . . . . . . . . . . . . . s i 1 (x n ) s i 2 (x n ) s i 3 (x n ) • • • s i n (x n )       , (8)
with rows

s i (x j ) = (s i 1 (x j ), s i 2 (x j ), • • • , s i n (x j )) = (X i W i X i ) -1 X i W i
, where W i is the weight matrix and X i is given in Equation ( 5) 2.4 Small Sized Porfolios and Heterogeneity. It is worth mentioning that the relational (proportional) model considered in Equation ( 1) implicitly accounts for differential mortality that may arise due to portfolio specific features, e.g. particular socioeconomic groups involved, income level, etc. Indeed, insured portfolios show a typical behavior compared to a national mortality. The mortality of insured population is significantly lower than the national population from which it is drawn. However, when it comes to the study of the mortality at a single portfolio level, some stylized facts arise, which might compromise the efficiency of the graduation procedure. For instance, insured population are generally of small size, so none or very few deaths are observable at some ages. Therefore, the use of the model ( 1) and the local-likelihood based estimation procedure advocates using the information stemming from the adjacent ages to construct the mortality curve. This learning procedure will enhance the determination of the mortality at a given age. However, when successive ages lack of information, the approach exposed above will need a large bandwidth h for the estimator to access distant ages with sufficient and reliable information. By doing so, we increase the bias surrounding the smoothed curve. In fact, as noted above, the mean of squared errors measuring the bias due to the local regression increase with the bandwidth h.

Due to these different sources of uncertainty we suppose that the true mortality curve ϕ i x , for x = x 1 • • • , x n , is known up to an unobservable multiplicative factor Θ i

x . In other words, the portfolios examined should be regarded as a sample of the reference. Estimates based on the data will be subject to sampling errors and the smaller the group is, the greater will be the relative random errors in the number of deaths and the less reliable will be the resulting estimates. This argument is extended to include the bias stemming at the attained age level due the consideration exposed above. Thus, if one has estimated the curve using the non-parametric approach, the true curve is an adjustment of the latter as multiplied by the random and non-observable parameter Θ i

x . This reasoning is closely related to the multiplicative frailty model that was introduced by [START_REF] Vaupel | The impact of heterogeneity in individual frailty on the dynamics of mortality[END_REF]. This concept is well established in biostatistics and survival analysis, see e.g. [START_REF] Hougaard | Life table methods for heterogeneous populations: distributions describing the heterogeneity[END_REF] and [START_REF] Aalen | Heterogeneity in survival analysis[END_REF] among others. In the actuarial literature frailty models have also been used. An overview of the use of frailty models to handle the heterogeneity phenomenon in life insurance is given in [START_REF] Olivieri | Heterogeneity in survival models. Applications to pensions and life annuities[END_REF]. In our case, we address this problem using the credibility theory, see [START_REF] Bühlmann | A course in credibility theory and its applications[END_REF] for more details.

Company-Specific Relative Risk Level

We suppose that there are K companies, with a portfolio of observed individuals of ages ranging from x 1 to x ni . Here, the n i 's might be all different to be in line with the insurance practices. This kind of information structure is similar to the so-called unbalanced framework used in actuarial science. For the sake of readability, without loss of generality, we will henceforth assume similar observed age groups for all companies, i.e. n = n 1 = • • • = n K . With a slight variation to the model, however, it can be easily extended to the unbalanced case.

3.1 The credibility model. Given the specific parameterization of the problem, one may think of the K portfolios as a subset of the reference population and thus each population is characterized by a risk profile. The latter is due to the heterogeneous sizes of the portfolios as well as the underlying guarantees (for the same underlying risk). These sources of heterogeneity might also induce an age varying risk profile within the same portfolio. Therefore, for a company i, we let the vector

Θ i = diag(Θ i x1 , • • • , Θ i xn ) be her relative risk level. For x ∈ {x 1 , • • • , x n }, each Θ i
x characterizes the age-specific risk level, which are unobservable random variables.

The primal objective is to characterize the force of mortality of each company i at a specific age x through the proportional relationship introduced in Subsection 2.2, i.e.

ϕ i = Θ i α, (1)
where

α = (α x1 , • • • , α xn ) such that for j = 1, • • • , n, we have α xj = f (x j )ϕ ref xj .
This model suggests that for each company i, the age-specific experienced force of mortality varies around the baseline α x , which can be seen as a reference or best-estimate mortality. This fluctuation is modeled by a heterogeneity parameter Θ i

x capturing the individual properties (heterogeneity) of each company at attained age x. Thus, using new incoming data should allow to update the the next period mortality ϕ i

x by adjustment following model in Equation ( 1). The approach is to first find an estimator ϕ i x of ϕ i x = Θ i x α x for each company i using the likelihood-based approach introduced in Section 2.

Henceforth, the notation ϕ i x |Θ i x refers to the estimation of the quantity α x = f (x)ϕ ref x , which, by abuse of language, referred to as the estimated mortality conditional on the risk profile. [START_REF] Nielsen | Credibility weighted hazard estimation[END_REF] and [START_REF] Gustafsson | Quantifying operational risk guided by kernel smoothing and continuous credibility[END_REF] proposed a similar framework to quantify the operational risk, where an adequate quantity related to losses from different lines of business share a common baseline. In our setting, we generalize their approach to the case where the mortality is structured by age (line of business) and intensities of mortality are dependent within the same group. Indeed, the random variables ϕ i x1 , • • • , ϕ i xn are assumed to be dependent, namely, the force of mortality of one age does directly impact those of other ages. This is mainly due to the graduation of mortality at a given age, which weights up over the adjacent age groups, see Section 2. This dependency will be explored later on this section. Finally, in order to characterize the next period mortality level, we make use of the credibility theory. For this purpose and using the usual credibility setting, we shall make the following assumptions:

(A1) The random vectors Θ i are independent across companies and ages. Moreover,

for i = 1, • • • , K, Θ i x 's are identically distributed with E[Θ i ] = I n and Var(Θ i ) = σ,
where σ is a diagonal matrix with elements σ x and I n is the identity matrix.

(A2) The random vectors

(ϕ i , Θ i ), i = 1, • • • , K, are independent across companies (A3) ϕ i x1 , • • • , ϕ i xn are conditionally independent given Θ i .
The first assumption (A1) ensures that the baseline mortality produces the a priori expected number of deaths under the model assumption (H0), in the sense that

E[D i x,• ] = E[Θ i x α x ] = α x .
The assumptions (A2) means that the risk profiles are independent over portfolios. In other words, the successive realizations of the mortality intensity (so as the death counts) for any portfolio are independent of each other except through the risk parameter. Finally, assumption (A3) translates the dependency of the mortality over ages. The latter is only captured by the vector Θ i . Conditionally on the latter the forces of mortality at the age level are independent.

As noted before, ϕ i |Θ i is the conditional local-likelihood estimator of the intensity in Equation (1) based on the data from the i th portfolio as developed in Section 2. In view of the the assumptions (A1)-(A3), it is important to recall that conditional on the knowledge of the risk profile Θ i the theoretical properties of ϕ i are identical to those the of local-likelihood estimator considered in Subsection 2.3. This will be used, among others, in the following lemma, in order to state some fundamental features of the dependence structure.

[A revoir la notation pour le produit] Lemma 1 Under assumptions (A1), (A2) and (A3) and the notation above, we have (i) The first order moment of ϕ i is given by

E[ϕ i ] = α.
(

) (ii) The variance matrix of ϕ i |Θ i , denoted Σ i (Θ i ) = Var( ϕ i |Θ i ), is given by Σ i (Θ i ) = ϕ i S i S i . 2
Hence, the variance Σ i = Var( ϕ i ) can be written as

Σ i = S i S i + σ α (4)
(iii) The covariance of ϕ i x with ϕ i x is given by

Cov ϕ i x , ϕ i = σ x α i x ) 2 e δx , (5)
with δ x = j if x = x j and e j is the vector with all 0's except for a 1 in the j th coordinate.

Proof. To show these results, we make an intensive use of the law of total variance. (i) Equation ( 2) is a direct consequence of assumption (A1) which gives E[ϕ i |Θ i ] = Θ i α.

(ii) The conditional variance Σ i (Θ i ) is directly derived from the calculus in Subsection 2.3. Hence, to check (4), the law of total variance gives

Σ i = E Var( ϕ i |Θ i ) + Var E[ ϕ i |Θ i ] , = E[ ϕ i ] ⊗ S i S i + α ⊗ Var(Θ i ) = α ⊗ S i S i + σ .
(iii) Finally, to prove (5), notice that Cov(ϕ i x , ϕ i |Θ i ) = 0. Thus,

Cov ϕ i x , ϕ i = Cov E[ϕ i (x)|Θ i ], E[ ϕ i |Θ i ] + E Cov(ϕ i x , ϕ i |Θ i ) , = Cov(Θ i x α i x , Θ i ⊗ α) = σ x α x ) 2 e δx ,
where the last equality follows from the independence assumption in (A1).

3.2

The next-period linear per-age mortality estimator. The goal is to predict the future force of mortality for each company i at the age level. Therefore, we will be looking for the inhomogeneous credibility predictor corresponding to the linear estimators of ϕ i x . We thus solve the following optimization problems:

min c i 0,x ,c i x E ϕ i x -c i 0,x -c i x ϕ i 2 , ( 6
)
where c i 0,x ∈ R and c i x ∈ R n . This formulation suggests adjusting the next period force of mortality at a given age using the information stemming from the other age groups. This should enhance the prediction for ages with low or sparse information using the credibility in ages of high information. Based on Proposition 1, we can easily derive the inhomogeneous credibility estimators of ϕ i . Indeed, we state the following proposition.

Proposition 1 The point estimate of the linear factors in (6) can be written as follows

c i 0,x = 1 n -c i x α and c i x = σ x α i x 2 (Σ i ) -1 e δx . ( 7
)
The next period predicted mortality (estimator) ϕ i of ϕ i is given by

ϕ i = (I n -α(Σ i ) -1 σα α + α(Σ i ) -1 σα ϕ i . ( 8
)
Proof. Let us first derive the intercept c i 0,x . To do this, we develop the expectation in Equation ( 6) and take the derivative with respect to c i 0,x . This yields to the following equality:

c i,0 + (c i x ) E[ ϕ i ] = 1.
On the other hand, differentiating the expectation in Equation ( 6) with respect to the vector c i x gives rise to the following variance

Var ϕ i x -(c i x ) ϕ i ,
needed to fully characterize the solution. The latter can be computed using results in Lemma 1. Indeed, we can write:

Var ϕ i x -(c i x ) ϕ i = Var(ϕ i x ) -2(c i x ) Cov(ϕ i x , ϕ i ) + (c i x )Σ i (c i x ) .
Taking the derivative with respect to the vector

(c i x ) yields 2Cov(ϕ i x , ϕ i ) -2Σ i c i x = 0.
The terms Σ i and Cov(ϕ i x , ϕ i ) are given in Lemma 1, which concludes the proof. Note that we are able to estimate all the components needed to characterize the next-period intensity ϕ i , except for the variance σ i . Remarking that ϕ i is an estimator of ϕ i = α ⊗ Θ i , we can write Θ i = ϕ i α, with " " being the pointwise division operator. Therefore, a natural choice of the estimator of σ is

σ = Θ i -I n Θ i -I n . (9)
We can now derive the Lemma 2 The optimal credibility estimator of Θ i is given by

Θ i = (I n -α(Σ i ) -1 σα 1 n + α(Σ i ) -1 σα Θ i ,
and the next period prediction of ϕ i can be approximated by α ⊗ Θ i .

Remark 1 The adjustment procedure described in Proposition 1 and Lemma 2 can written for each individual age x in the following form:

ϕ i x = (1 -z i x )α x + z i x ϕ i x ,
where z i x is the credibility factor that can be written in the following form:

z i x = (α x ) 2 σ 2 i (α x ) 2 σ 2 i + ϕ i x s i (x) 2 .
Here, recall that s i (x) 2 = n j=1 s i j (x) 2 and measures the reduction in variance of the smoothed mortality curve ϕ i x .

Remark 2 All the ingredients required to implement the credibility approach in Lemma 2, in order to predict the next-period estimator, are already determined. However, we still need to characterize an estimation for α. To do this, we borrow the same procedure considered in [START_REF] Nielsen | Credibility weighted hazard estimation[END_REF], which amend to estimate α as linear weighted average over the portfolios.

Numerical Analysis

4.1 Source of Data. The data come from studies conducted by Institut des Actuaires. These studies include in total 14 portfolios covering the period 2007-2011 with each companies contributing data for at least 4 of a possible 5 years. Table 1 presents the observed characteristics of the male population of the portfolios. For this dataset, we are considering respectively T i = 3 and T i = 4 for all companies. The remaining years serve to test the predictive feature of the model through an in-sample analysis. The considered analysis follows similar lines as in [START_REF] Salhi | A credibility approach of the Makeham mortality law[END_REF], which also exploits the same dataset. Therefore, the age band for all companies ranges from x = 30 to 95 years old. Figure 1 shows the age distribution of the portfolios, i.e. the aggregate number of individuals exposed to risk at each attained age. It graphically depicts the heterogeneity observed between the portfolios with insured holding different policies. On the other hand, the baseline mortality ϕ ref x is a market table, denoted IA20131 , constructed for the French insurance market provided by Institute des Actuaires, see Tomas and Planchet (2013c). Before proceeding to the implementation of the methodology developed in the previous sections, we must should look deeper into the particular feature of our dataset.. Specifically, we must focus on those that may arise specific concerns when it comes to the graduation of a mortality table using the smoothing procedure considered in Section 2. As previously reported, the experienced mortality does not only suffer from a small sample size but also the under-representation of those within some age groups. This is typically the case of portfolio P2, see Table 1 and Figure 1. In fact, we have a small size of 7589 individuals with only 2% aged under 60. This is the case for portfolios P7, P8 and to some extent P10, but with a bigger exposure. For these portfolios, the use of the smoothing procedure in Subsection 2.2 has the advantage of borrowing the information in age bands where the exposure is substantially larger. However, this may alter the local properties generally required for a mortality curve. Indeed, enlarging the smoothing window h, giving access to far distant ages, may amend the increasing mortality over ages, which is not only a very much sought behavior but also a biologically reasonable quality.

4.2 Entity-Specific Graduated Mortality. In order to implement the likelihood based graduation approach, we need to identify the fitting variables. To this end, we first choose a weight function More precisely, we should characterize the smoothing parameters h and p referring, respectively, to (i) the window-width or bandwidth h which determines how fast the weights decrease and (ii) the degree of the polynomial used for smoothing. To do so, we have to balance the trade-off between bias and variance on the estimation of the mortality as they move on the opposite direction in terms of the smoothing parameters. Formally, we compute a number of candidate fits and use a criteria to select, among the fits, the one with the lowest score. For instance, we use a generalization of the Akaike information criterion (AIC) based on the deviance function taking into account the likelihood as well as the degrees of freedom of the fitting model. The smoothing parameters are selected based on graphical diagnostics. Figure 2 displays the AIC scores against the fitted degrees of freedom for the portfolios 1, 6 and 14. The interpretation of such a graphic is as follows. A two-degree of freedom represents a smooth model with very little flexibility while 10 degrees of freedom represents a noisy model showing many features. It also aids comparability as we can compute criteria scores for other polynomial degrees or for other smoothing methods and added them to the plot. We select the smoothing parameters at the point when the criterion reaches a minimum or a plateau after a steep descent. In the left panel of Figure 2, the graphic represents the selection procedure for the portfolio P1. This consists of scores for various polynomial degrees. The optimal choice of the smoothing parameters corresponds to the point on the plan where with the lowest AIC score. For portfolio P1, the lowest score corresponds to a cubic fit with 5 degrees of freedom leading to a smoothing window of 7 points. In Table 2, we report the optimal parameters for each portfolio.

In Table 2, we can see that for some portfolios the optimal choice of the input controls induce high level of degrees of freedom, i.e. portfolios P1 and P5. This is to say that the corresponding "smoothed" curves ϕ i , i = 1, 5, will be noisy showing many feature. Indeed, the degree of freedom is a qualitative proxy for the regularity of graduated mortality curve as the smoothness inversely to

1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3
Values of the AIC for various polynomial degrees q 2 the degree of freedom. This feature can already be deduced from the limited amount of information (exposures) that are at our disposable for these portfolios. For the remaining portfolio, the degrees of freedom are relatively small. In the following, we will implement the credibility approach described in the last section to assess the impact of the latter on the graduated mortality curves.

4.3 Next-Period Mortality Rate. Here, we consider the mortality experience over the period 2007 -2010 upon which we calibrate the smoothing procedure considered in the above subsection.

For each portfolio, we build a graduated mortality table ϕ i and aim at adjusting the latter for the next period projection. For each age x, the graduated mortality gives a candidate rate for the next period, i.e. ϕ i x . The insurer has the possibility of whether to rely on this rate or adjust the latter given the experience stemming from the other rates at other ages. In other words, the mortality used for the next period forecasts can be adjusted using the credibility formula in Equation (8). To do this, we estimate the different quantities needed to implement (8) as follows: (i) Following Remark 2, the expected mortality rate α can be estimated as follows:

α = xn x=x1 d i=1 E i x,• ϕ i x / d i=1 E i •,• .
(ii) The weight loading matrix S i is given as an output of the graduate step and can be estimated using Equation ( 8).

(iii) The diagonal matrix σ i relies on the variance of Θ i x 's which might be estimated given that

Θ i x = ϕ i
x / α, and thus we can write

σ i x = d i=1 ϕ i x / α 2 /d - d i=1 ϕ i x / α /d 2 .
Figure 3 depicts, respectively, the graduated mortality over the period 2007 -2010 as described above and the next-period (2011) mortality rates using the credibility formula in Equation ( 8). Similarly, Figure 4 represents the next-period predicted deaths using the two mortality rates. In these figures we grayed areas (ages) where the relative difference between the smoothed mortality and its adjusted counterpart exceeds a 10% level. More precisely, this corresponds to the ages x where

| ϕ i x -ϕ i x |/ ϕ i x > 0.1.
At a first glance, we remark that the credibility adjustment does change the mortality rate and overall propose a smoother curve compared to the initial one, and this is even evident when dealing with portfolios with small portfolios and high degrees of freedom. In fact, when we deal with portfolios such as P5, where the exposure-to-risk as well as the underlying deaths are very limited, the smoothing approach fails to capture the mortality structure and the output of the procedure proposed in Section 2 are very irregular and noisy. Indeed, as noted above such a procedure need information stemming from adjacent ages when a particular age lacks of sufficient exposure. The case of P5 is very appealing of the limit of the non-parametric smoothing as the limitation on the information is shared over the ages. This is why the corresponding degree of freedom is high and the AIC is low, see Table 2, and explains the irregular curve (dashed line) for the smoothed mortality.

On the other hand, the degrees of freedom given as tr(S i S i ) provide information on the credibility of the smoothed curve ϕ i . In fact, as we can see in Equation ( 8) or in a more tractable way as in Remark 1, the higher the degrees of freedom, i.e. tr(S i S i ) = xn x=x1 s i (x) 2 ; the smallest is the weight attached to the smoothed curve ϕ i (in aggregate). At the age level, one should look the component-wise of the degrees of freedom, i.e. the individual variance s i (x) 2 . That being said, we can conclude that the parameter driving the adjustment at the age levels is vanishing meaning that the adjusted mortality rate ϕ i

x is close to the reference α x . It comes as no surprise, then, to find that the adjusted curve tends to obliterate this undesired effect thanks also to the information coming from other ages but from different portfolios.

The visual inspection of the credibility based mortality curve shows that the regularity is preserved avoiding the limitation of the sole smoothing procedure discussed above. For some portfolios, such as portfolio P12, the regularization based on the credibility attached to each age level enhance the prediction of the future mortality. Indeed, the smoothed mortality based on past observations suggest a local distortion of the curve for ages ranging from 60 to 80. This particular feature is however not observed in the mortality curve for the year 2011 and thus the credibility based curve has a better fitting. This can also be observed in Figure 4, where the predicted deaths using ϕ i , for i = 12, are (visually) more in line with the observations. The same conclusions, in the grayed area, can be drawn for the other portfolios.

4.4 Proximity Between the Observations and the Model. Besides the visual inspection of the proposed adjustment, and in order to understand the impact of the the latter, we will use some known statistics to quantify the proximity between the observations and the outputs of the two curves considered in Figure 3 and Figure 4. We assess the overall deviation with the observed mortality by comparing criteria measuring the distance between the observations and the models with the χ 2 , i.e. Forfar et al. (1988b), the mean average percentage error (MAPE), see for instance Age Mortality Rate -8 -6 -4 -2 qq qq q q q q q qq q q qq q q q q q qq q q q q q q q q q qq q q q q q q qq q q q q q q q qqq qq q qq q q q q q q q Portfolio 1 40 50 60 70 80 90 q q q q qq q q q q q q q q q q q q q q q q q q q q q q q Portfolio 2 qq q q qq q qqq q q q q q q q qq q q q q q q q q q q qq q qq q q q Portfolio 4 40 50 60 70 80 90 q q q q q q q qq qq q q q q q q Portfolio 5 qqq qq q q qq q q q q q q qq q q q q q q q qqq qq qq qq qqq q q q q qq q q qq q q q qq Portfolio 7 q q qq q q q q q q q q q qq q qqq q qq qqq q qqq q q q q qq q Portfolio 8 qq q q q q q q qq qq q q q qq qq qq q qq q qqqq q q qqqq qqqq q qq qq qqq q q q q q q q qq q q q qq Portfolio 9

-8 -6 -4

-2 q q q q q q q q q qq q q q qq q q q q q qq q q q q q q q q qq q q q qqqq q q q q q q Portfolio 10 -8

-6 -4

-2 40 50 60 70 80 90 qq q qq q qq q q q q q qq q q q q q qqq q q qqq q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q Portfolio 11 q q q q q q qq q q q q q q q q qq q q q q q q q q q q q q qq qq q q q q qq qq q q q Portfolio 12 40 50 60 70 80 90 q q q q q q q q q q q q q q q qq q q q q qq qq q q q q q q q q qqq q q q qq q q qq q q q qq q q q qq qq qq qq q qq q q Portfolio 13 q q qq qq q q q q qq q q q q qq q q q q q q q q q q q q q qqq qq q q q q q q q q q q q q q q qq q q q Portfolio 14 q Raw Smoothed Smoothed+Adjusted

Age

Deaths 0 5 10 15 20 25 30 q qq q qq q q q q q qq q q qq q q q q q qq q q q q qq q qq q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q Portfolio 1 40 60 80 0 1 2 3 4 5 6 qqqqqqqqqqqqqqqqqqqqqqqqqqqqqqqq q q q qq q qq q q q q qq qq q q q q q q q q q q q q q qq q q q Portfolio 2 0 2 4 6 8 10 qq qqq q q q qq q qqq q q q q q q q q qq q q q qq q q qqq q q q q q qq qq q qqq qq q q q qqqqqqqqqqqqqq Portfolio 4 40 60 80 0.0 0.5 1.0 1.5 2.0 2.5 3.0 qqqqqqqq q qqqq q qq q qq qq q q q q q qq q q q q q q qqq q qqqqqqqqq q q q q qqqqqqqqqqqqqqq Portfolio 5 0 100 200 300 400 qqqqqqqqqqqqqqqqqqqqqqqqqqqq qqqq qqqqq q q q q q q q q q q q q q q q qqqqqq qq qqqqqq Portfolio 7 0 20 40 60 qqqqqqqqqqqq q qqqqqqq q q qq qq q qqq q qqqq qq q q q q q q qq q q q q q q q q q q q q q q q q q q q q q Portfolio 8 0 20 40 60 qqq qq q q q q q q q q qq qq q q q q q qq qq q q q q q q q q q q q qq q q q q q q q q qq q q q q q q q q q q q q q q qq q Portfolio 9 0 10 20 30 40 50 qqqqqqqqqqqq qq qq q q qqq q q q qq q q qq q q q qq q q q q q q qq q q q q q q q q q q q q q q q q q q q qq q q Portfolio 10 0 5 10 15 20 25 40 60 80 qq q qq q qqq q q q q q q q q q q q q q q q q q q q q q q q q q qq qq q q q q q q q q q q q q qq q q q qq q q q q q q qq q Portfolio 11 0 2 4 6 qq q q q qq q q q qq q qq q q q q q q q q q qq q q qq q q q q q qq q q q qq qq qq qq q qq q q q q q qq q q q q qq qq Portfolio 12 0 10 20 30 40 50 40 60 80 q q q q q q q q q q q q q q q qq q q q q q q qq q q q q q q q q q q q q q q q q q q qq q q q qq q q q q q q q q q q q q q q q q Portfolio 13 0 5 10 15 20 qqqqqqqqq q q qq q qq q q q q qq q q q q q qq q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q Portfolio 14 q Raw Smoothed Smoothed+Adjusted et al. (2002); as well as the standardized mortality ratio (SMR), i.e. the ratio of deaths observed to those predicted. In addition, we find useful to use the SMR test proposed in [START_REF] Liddell | Simple exact analysis of the standardised mortality ratio[END_REF] and the likelihood ratio test to asses of significance of the SMR and to its interval estimation.

The tests and quantities summarizing the proximity between the observations and the model, for each portfolio i at calendar year t = 2011, are described in the following.

(i) The χ 2 i allows to measure the quality of the fit of the model. It writes,

χ 2 i = xn x=x1 D i x,t -L i x,t q i x 2 L i x,t q i x (t) 1 -q i x (t)
.

(ii) The MAPE is the average of the absolute values of the deviations from the observations,

MAPE i = xn x=x1 D i x,t /L i x,t -q i x / D i x,t /L i x,t xn x=x1 D i x,t × 100.
(iii) The SMR is computed as the ratio between the observed and fitted number of deaths:

SMR i = xn x=x1 D i x,t xn x=x1 L i x,t q i x (t)
.

Hence, if SMR i > 1, the fitted deaths are under-estimated and vice-versa if SMR < 1. Note that we can consider the SMR as a global criterion which does not take the age structure into account, compared to the χ 2 and the MAPE for instance. Figure 5 summarizes the above mentioned tests and quantities giving the overall deviation between the observations and the adjustment analysis for the portfolios P1 to P14 (except 3 and 6 which does not contain observations of year 2011) obtained by the smoothing approach together with the credibility adjustment procedure. When looking at criteria and quantities which take the age structure of the error into account, the credibility approach has an important benefit compared to the sole graduated curve. The quality of the fit increases, sometimes drastically, i.e. portfolio P1, in terms of having the minimum χ 2 and MAPE values, see the last panels of Table 3. Also, the credibility adjustment exhibits the highest p-value for the likelihood ratio test. Even when we consider a global indicator of the quality of the fit such as the SMR which does not take the age structure into account, the proposed procedure seems to perform better than the graduated curve. However, notice that the impact of adjustment in smaller when the portfolios are quiet big. This is already remarkable when one operated visual checks as mentioned earlier.

Comparaison with Classical

Approaches. Here, we wish to compare our model to the [START_REF] Hardy | A credibility approach to mortality risk[END_REF] and Poisson-Gamma credibility analysis applied to our mortality datasets. Moreover, we intend to compare our results to a similar approach introduced in [START_REF] Salhi | A credibility approach of the Makeham mortality law[END_REF], where the graduation of mortality is based on a parametric method, i.e. Makeham law, and the credibility theory is used to adjust the latter with incoming new data experience. The first two approaches focus on the actual to expected mortality ratio, in aggregate level, as a key observation. The adjustment is directly applied to this quantity. Specifically, the a priori expected number of deaths for each portfolio is updated at each period given the credibility weight on the observations coming from this portfolio and the one computed on the basis of the other portfolios. these different approaches in comparison the one exposed in this paper. We first note that the Hardy-Panjer and Poisson-Gamma approaches produce relatively similar graduations as the tests suggest sensibly similar outputs. However, we notice some differences with the Makeham credibility model which displays more favorable results. This is already outlined in [START_REF] Salhi | A credibility approach of the Makeham mortality law[END_REF] and this may be explained by the age-specific adjustment but thanks to the structural feature added by the Makeham parametric model. When it comes to the tests that are sensitive to the age structure, we notice that the credibility-based Makeham model offers a outstanding adjustment as the tests are favorable compared to the Hardy-Panjer and Poisson-Gamma approaches.

Concluding Remarks

In this paper we proposed a methodology to adjust the graduated mortality table that uses an adaptive smoothing procedure based on the local likelihood. The adjustment is based on the credibility weighting technique of the smoothed curves and a reference. Out approach takes into account the age specific heterogeneity that may arise in real world datasets. We thus consider updating the mortality for each age based on the upcoming past information from the same age but also the neighboring ages. The inclusion of the neighboring ages is crucial as the particular smoothing procedure used in this paper add a dependency between the single ages. Based on classical results on the inference of the smoothing procedure we derived the closed form formulas needed to adjust the mortality.

The proposed methodology is shown to outperform compared to the classical credibility approaches that does not take into account the age structure of the portfolio. This is in line with the recent work in this field as mentioned by [START_REF] Salhi | A credibility approach of the Makeham mortality law[END_REF]. Even when the age structure is accounted for the methodology developed in this paper has an important benefit. This is mainly due to the underlying curve built using an adaptive procedure compared to the parametric model considered in [START_REF] Salhi | A credibility approach of the Makeham mortality law[END_REF].

We should note that the model proposed can be investigated in order to quantify mathematically the errors induced in the assessment of the next-period mortality curve. This amend to consider the uncertainty stemming from the estimation of the different variables used in the updating procedure. There are also several practical issues we do not address here such as the impact on the pricing of life insurance contracts. These are open questions that we openly acknowledge and leave for future research. Table 3: Tests and quantities summarizing the deviation between the observations and the model

A Additional Tables

Figure 1 :

 1 Figure 1: Distribution of age groups in the portfolios.

Figure 2 :

 2 Figure 2: AIC scores for various polynomial degrees for portfolios P1, P6 and P14

Figure 3 :

 3 Figure 3: The next-period mortality rate based on the graduation of mortality (solid line) as well as its credibility based adjustment (dashed line) based on the period 2007 -2010. Both predictions are compared to the observed mortality rate over the period 2011 (black dots).

Figure 4 :

 4 Figure 4: The next-period deaths prediction based on the graduation of mortality (solid line) as well as its credibility based adjustment (dashed line) based on the period 2007 -2010. Both predictions are compared to the observed deaths over the year 2011 (black dots).

Felipe

Figure 5 :

 5 Figure 5: Distance between the observations and the models (top panel) together with their mean average percentage error between (middle panel) and the standarized mortality ratio (bottom panel).

Table 1 :

 1 Observed characteristics of portfolios population.

		Period of observation Beginning End	Mean age In Out	Average exposure	Mean age at death	Size
	P1	01/01/07	12/31/11	36.96 39.74	2.77	68.78	616390
	P2	01/01/07	12/31/11	69.3 73.35	4.05	80.34	7589
	P3	01/01/07	12/31/10	40.16 43.1	2.94	71.77	80086
	P4	01/01/07	12/31/11	37.5 41.13	3.63	54.08	93165
	P5	01/01/07	12/31/11	36.9	39.1	2.2	59.31	21540
	P6	01/01/07	12/31/10	48.5 52.11	3.62	82.34	847469
	P7	01/01/07	12/31/11	66.65 71.29	4.64	73.68	89507
	P8	01/01/07	04/13/11	67.51 71.38	3.86	80.72	78650
	P9	01/01/07	06/30/11	45.97 49.6	3.62	73.17	1556150
	P10	01/01/07	12/31/11	62.97 67.64	4.67	79.77	132990
	P11	01/01/07	12/31/11	38.89	42	3.11	56.44	420405
	P12	01/01/07	12/31/11	37.05 39.2	2.15	57.41	904020
	P13	01/01/07	12/31/11	43.01 46.89	3.88	71.03	848757
	P14	01/01/07	12/31/11	50.12 54.16	4.04	72.37	233488

Table 2 :

 2 Local-likelihood smoothing parameters' optimal choice

		P1	P2	P3	P4	P5	P6	P7
	AIC	98.90 29.34 79.68 78.91 61.10 106.73 68.37
	Degree	2	1	1	1	3	3	1
	DF	18	3	3	6	16	10	8
		P8	P9	P10	P11	P12	P13	P14
	AIC	74.11 83.96 53.44 70.86 74.82 82.70 78.37
	Degree	2	1	2	1	1	3	2
	DF	4	7	5	6	6	7	5

  Table 3 presents the tests and quantities summarizing the overall deviation between the observations for

This table is derived from mortality trends originating from the INSEE table, French national bureau of statistics, and covers the period

2007-2060.

* This work has been supported by the BNP Paribas Cardif Chair "Data Analytics & Models in Insurance" and the ANR project "LoLitA" (ANR-13-BS01-0011). The views expressed in this document are the authors owns and do not necessarily reflect those endorsed by BNP Paribas Cardif.