Alexandre Boumezoued
email: alexandre.boumezoued@milliman.com

Laurent Devineau
email: laurent.devineau@milliman.com

Individual claims reserving: a survey

Keywords: Non-life/General insurance, Stochastic reserving, Individual claims reserving, Poisson point processes, Process and estimation errors, Thinning algorithm

, and provides a consistent presentation of the modelling, inference, and forecasting (with simulation and closed-forms) of individual claims histories as well as aggregate quantities as the overall reserve for both RBNS and IBNR claims. Numerical illustrations are given based on real portfolio datasets, as well as comparisons with classical triangle-based methods.

Introduction

Since the work of [START_REF] Mack | Distribution-free calculation of the standard error of chain ladder reserve estimates[END_REF], stochastic reserving models based on runo triangles dominate the reserving practice, as surveyed in the book of [START_REF] Wüthrich | Stochastic claims reserving methods in insurance[END_REF] -the triangles are organized by origin period (e.g. accident, underwriting) and development period, and are built as an aggregation of input data made of individual claims paths. At around the same time, other papers proposed to account for a more precise description of the claims development. To our knowledge, [START_REF] Arjas | The claims reserving problem in non-life insurance: Some structural ideas[END_REF], [START_REF] Jewell | Predicting IBNYR events and delays: I. continuous time[END_REF], [START_REF] Norberg | Prediction of outstanding liabilities in non-life insurance[END_REF] and [START_REF] Hesselager | A Markov model for loss reserving[END_REF] are among the oldest references which introduced a proper probabilistic setting for individual reserving, in the context of increasing awareness of the power of a "point process" description of time patterns in many areas of application. In this setting, the time pattern of claims occurrences has systematically been investigated in the form of the Poisson model with time-varying intensity (a stochastic intensity for few exceptions with so-called Cox process), which will be the framework of interest in this paper. In their spirit, individual claims models provide a structural modelling of aggregate 'macro' quantities as a result of the 'micro' components of the individual modelfrom [START_REF] Norberg | Prediction of outstanding liabilities II. model variations and extensions[END_REF], a parametric specication of the continuous time model, which may be supported by physical reasoning, will automatically induce a parametrization of the discrete time model. More generally, the virtue of individual-based modelling emerges in a variety of problems from its structural nature, as it is the case of the separate evaluation of IBNR (Incurred But Not Reported) and RBNS (Reported But Not Settled) claims. Although some contributions focused on the development a coherent modelling framework based on runo triangles providing separate IBNR and RBNS reserves, see e.g. [START_REF] Verrall | Prediction of RBNS and IBNR claims using claim amounts and claim counts[END_REF], this appears in a more natural (and straightforward) way in individual claims modelling, see e.g. the discussion in [START_REF] Antonio | Micro-level stochastic loss reserving for general insurance[END_REF] and references therein, as it is also naturally the case for the valuation of non-proportional reinsurance treaties.

Since then, the individual based modelling of claims has gained recent interest in the literature. Following [START_REF] Norberg | Prediction of outstanding liabilities in non-life insurance[END_REF], [START_REF] Larsen | An individual claims reserving model[END_REF] revisited the Marked Point Process approach by including claims features to specify the model components. [START_REF] Antonio | Micro-level stochastic loss reserving for general insurance[END_REF] proposed to apply the modelling framework developed in [START_REF] Norberg | Prediction of outstanding liabilities in non-life insurance[END_REF] and [START_REF] Norberg | Prediction of outstanding liabilities II. model variations and extensions[END_REF] to a general liability insurance portfolio of a European insurance company, based on a set of distributional assumptions and a maximum likelihood estimation procedure. In their empirical investigation, and in the context of this specic case study, it is shown that individual reserving provides a better accuracy compared to selected aggregate models. Badescu et al. (2016b) and Badescu et al. (2016a) investigate the modelling of individual claims by means of a Cox process and the related inference; in this work, stochastic intensity is crucial to account for temporal dependence in claims arrival and settlement; the model focuses on number of reported and IBNR claims, without allowing for payments. Among the several authors who focused their attention into inferring the distribution of reporting lags as a rst step in the computation of a proper IBNR reserve, one can also refer to [START_REF] Guiahi | A probabilistic model for IBNR claims[END_REF], [START_REF] Jewell | Predicting IBNYR events and delays: I. continuous time[END_REF], [START_REF] Zhao | Applying copula models to individual claim loss reserving methods[END_REF], and more recently Verrall and Wüthrich (2016) -this latter paper provided additional interesting insights on the analysis and calibration of arrivals and delays distributions based on real data; indeed, it is observed a gain in using individual methods in environments characterized by non-stationarity.

It has also to be noted that other contributions chose another direction for individual models, away from a continuous time point process description: see in particular [START_REF] References Pigeon | Individual loss reserving with the multivariate skew normal framework[END_REF] for a discrete time formulation using stochastic indi-vidual development factors (see also [START_REF] Taylor | Individual claim loss reserving conditioned by case estimates[END_REF]), as well as [START_REF] Godecharle | Reserving by conditioning on markers of individual claims: a case study using historical simulation[END_REF] who rely on the Reserve by Detailed Conditioning (RDC) method in discrete time from [START_REF] Rosenlund | Bootstrapping individual claim histories[END_REF] based on detailed information on individual claims, and also [START_REF] Antonio | A multi-state approach and exible payment distributions for micro-level reserving in general insurance[END_REF] for a multi-state framework in discrete time to model the claims development process while allowing for exible payment distributions.

In this survey, we revisit the original probabilistic formulations of [START_REF] Norberg | Prediction of outstanding liabilities in non-life insurance[END_REF] and [START_REF] Hesselager | A Markov model for loss reserving[END_REF] as they provide a continuous-time and structural description of claims arrival and histories, and as we consider additionally that a higher potential emerges from these methods from a pure 'risk analysis' perspective.

Most of the papers investigated the benet accuracy of individual reserving from a pure empirical point of view, in the form of back-testing eciency assessment, leading to conclusions that were in each case limited to the data set used. Although this is crucial to get insights on the concrete order of magnitude of the gain of using individual method based on a given dataset, from our point of view, it has not been clearly investigated why do stochastic individual reserving methods outperform aggregate ones from a mathematical perspective. A lot of research work remains to be done in this direction. A rst attempt of a theoretical description of the gain in process error has been carried out in [START_REF] Huang | An individual loss reserving model with independent reporting and settlement[END_REF], leading to a comparison of individual and aggregate models under a common formalism, although simplied compared to the generic descriptions as provided in [START_REF] Norberg | Prediction of outstanding liabilities in non-life insurance[END_REF]. In this context, we argue that there is a need to revisit the mathematical foundations of individual claims modelling and provide professionals with a proper toolbox allowing for an improved risk assessment.

It is interesting to note that, although for triangle-based models the treatment of estimation error has early become a standard area of interest, dealing with parameter uncertainty in the context of individual reserving has been far less investigated. From an empirical perspective, references as [START_REF] Haastrup | Claims reserving in continuous time; a nonparametric bayesian approach[END_REF] (bayesian non-parametric) and Antonio and Plat (2014) (frequentist parametric) are among the few to address this question, although parameter uncertainty is a great area of improvement of accuracy of individual models as data available is "larger" compared to that of an aggregate triangle. Indeed, as noted by [START_REF] Antonio | Micro-level stochastic loss reserving for general insurance[END_REF], due to our large sample size, condence intervals for parameters are narrow. This is in contrast with runo triangles where sample sizes are typically very small and parameter uncertainty is an important point of concern.

As for simulation issues, Antonio and Plat (2014) rely on the specic form of a piecewise constant intensity by simulating occurrences as Poisson numbers in each (small) interval, or inverting cumulative distribution functions. In this paper, it is used a so-called thinning procedure to simulate the (past) occurrences of IBNR claims, the (future) development of RBNS and IBNR claims, as well as possibly claims which did not occur yet. We argue that such exact simulation method provides a simple and exible framework which is well suited for the simulation of claims occurrences, development and settlement.

The paper is structured in two parts. In Section 2, we develop the modelling framework for occurrences and reporting delays at the individual claims level, with a special focus on inferring and forecasting IBNR claims (past) occurrences. Section 3 enriches the previous model for claims payments development over time and settlement, allowing to precisely forecast the payment streams of both IBNR and RBNS claims. Each section develops the following components: the modelling framework, the simulation algorithm, the likelihood derivation and inference strategy, a set of macro results and closed-form formulas on the reserve and associated variance, as well as numerical examples. The paper ends with some concluding remarks in Section 4.

2 Occurrence and reporting

Stochastic individual claims model

Consider claims occurring at times (T n ) n≥1 given as the jump times of a (nonhomogeneous) Poisson process with intensity λ(t). Each claim which occurs at time T n is associated with a delay (waiting time) U n assumed to have distribution p U |Tn (du). In the modelling, this waiting time represents a reporting delay (the dierence between time of notication to the insurer and time of occurrence): this is a crucial quantity as it determines the set of claims which are observed in the dataset. Let us introduce an observation time τ -the set of claims observed is exactly the set of reported claims at time τ , denoted I R (τ ) and dened as

I R (τ ) = {(T n , U n ) such that T n + U n ≤ τ } . ( 1
)
The set of all claims which occurred before time τ , introduced as

I O (τ ) = {(T n , U n ) such that T n ≤ τ } ,
can then be split between reported and so-called incurred but not reported (IBNR) claims as

I O (τ ) = I R (τ ) ∪ I IBN R (τ ),
where

I IBN R (τ ) = {(T n , U n ) such that T n ≤ τ and T n + U n > τ } .
Milliman R&D 4/39 (i) for all non-overlapping measurable sets

B 1 , ..., B k of Ē = R + ×E, r.v. Q(B 1 ), ..., Q(B k )
are independent (innite r.v. being, as a constant, independent of any other),

(ii) for all measurable set B ⊂ Ē such that q(B) < +∞, Q(B) ∼ Poisson (q(B)) .

(iii) Q({0} × E) = 0 (this ensures no jump at time 0). Under this denition, the Poisson point measure is a counting measure on R * + , that is a.s. for each t > 0, Q({t} × E) ∈ {0, 1} (see e.g. Çnlar (2011), Theorem 2.17).

For the purpose of modelling occurrence and reporting in the present section, let us consider Q(dt, du) a Poisson point measure with intensity measure q(dt, du) = λ(t)dt ⊗ p U |t (du).

(2)

Then the following representation holds:

Q(dt, du) = n≥1 δ (Tn,Un) (dt, du).
The number of claims which occurred before τ , whatever their reporting delay, can be simply represented as

Q([0, τ ], [0, ∞)) = τ 0 ∞ 0 Q(dt,
du); this is the total sum of numbers in the completed triangle. This provides a straightforward split between reported claims and IBNR claims as

Q([0, τ ], [0, ∞)) = I R (τ ) Q(dt, du) Reported + I IBN R (τ ) Q(dt, du) Incurred but not Reported .
A key result which follows from i) and ii) of Denition 1 is given below. Proposition 1. The number of reported claims Q I R (τ ) and the number of IBNR claims Q I IBN R (τ ) are independent, and follow Poisson distributions with parameters q I R (τ ) and q I IBN R (τ ) respectively.

In particular, this result states that the expected IBNR number writes

q I IBN R (τ ) = τ 0 ∞ τ -t λ(t)p U |t (du)dt. (3) Milliman R&D 5/39
The frequency of reported claims Proposition 1 provides the distribution of the overall number of claims, both reported and IBNR claims, as a Poisson random variable which can be simulated in a straightforward manner. For our purpose, we are also interested in the frequency at which those claims occurred, that is on the distribution of the occurrences (instead of that of the overall number). Let us rst compute the probability distribution of reported claims and their associated marks.

To do this, we introduce the Poisson point measure of reported claims as

Q R (dt, du) = Q(dt, du)1 t+u≤τ ,
whose intensity measure is in fact q(dt, du)1 t+u≤τ .

Recall that the intensity measure of Q is q(dt, du) = λ(t)p U |t (du)dt; then the Poisson point measure Q R , represented as

Q R (dt, du) = n≥1 δ (T R n ,U R n ) (dt, du), is described as follows: Proposition 2. (i) The intensity of reported claims occurrences (T R n ) n≥1 is λ R (t) = λ(t)p U |t ([0, τ -t]), (4)
and at each time of reported claim occurrence T R n , the delay is drawn on

[0, τ -T R n ] with conditional distribution p U |T R n (du)/p U |T R n ([0, τ -T R n ]).
Indeed, the pure jump process N R t = Q R ([0, t), R + ) which counts the number of reported claims up to time t can be compensated by the associated intensity measure, so that the following process is a martingale:

N R t - (0,t)×R + q(ds, du)1 s+u≤τ = N R t - (0,t) λ(s)p U |s ([0, τ -s])ds,
which proves that the intensity is that given in Equation (4). The associated delay distribution then follows from the basic denition of the Poisson random measure Q R .

Simulation algorithm

The previous Proposition 2 characterizes the way the observed data sample can be resimulated: in particular, it states that occurrences of reported claims took place at frequency given in Equation (4). This formula will be useful to derive the likelihood in the next Subsection 2.3 dedicated to parameter inference.

Note that reciprocally, as the overall intensity is λ(t), Proposition 2 provides the occurrence frequency of IBNR claims which can be drawn independently from reported claims, as a non-homogeneous Poisson process with intensity

λ IBN R (t) = λ(t)p U |t ((τ -t, ∞)).
(5)

Once the parameters of the model λ(t) and p U |t (du) have been computed (see next Subsection 2.3), one is then interested in completing the observed data set based on occurence times of IBNR claims. Due to the stochastic nature of the model, an algorithm has to be set to practically 'draw' the times at which IBNR claims should have occurred in the observation period [0, τ ], which we denote T IBN R n . This can be achieved according to the so-called thinning algorithm which provides a very simple way of computing paths of any counting process with time-dependent intensity.

The thinning algorithm, which can be seen as a generalization of the acceptance/rejection method, works as follows:

1. First consider a bound λ on the intensity -take λ = sup t∈[0,τ ] λ IBN R (t).

2. Then draw the sequence (S k ) on [0, τ ] as a Poisson process with intensity λ -for example by drawing S k+1 -S k as iid increments that are exponentially distributed with parameter λ.

3.

For each time S k , draw a Bernoulli random variable with 1-value probability λ IBN R (S k )/ λ -the 1-value corresponds to the acceptance of the time as part of the sequence T IBN R n .

From a mathematical point of view, the thinning algorithm comes from a key martingale result which, again, can be conveniently formulated by means of Poisson random measures -this is detailed at the end of this subsection. Note that this general mathematical formalism for such procedure goes back to [START_REF] Kerstan | Teilprozesse Poissonscher prozesse[END_REF] and [START_REF] Grigelionis | The representation of integer-valued random measures as stochastic integrals over the Poisson measure[END_REF], although the thinning algorithm is most often referred to [START_REF] Lewis | Simulation of nonhomogeneous Poisson processes by thinning[END_REF] and [START_REF] Ogata | On Lewis' simulation method for point processes[END_REF] who proposed a more operational description of the process. A precise formulation is given in the context of Hawkes and population processes in Boumezoued (2016b). Note also that more advanced thinning procedures can be used, in order e.g. to improve its eciency, see e.g. [START_REF] Giesecke | Monte Carlo algorithms for default timing problems[END_REF]. In this paper, we restrict our attention to the original thinning approach which is, in its form, sucient to address eciently a wide variety of IBNR simulation contexts.

The thinning representation of counting processes Consider a Poisson point measure Ñ (ds, dθ) with intensity measure λds ⊗ dθ on the space R + × [0, 1], and denote by (F t ) the canonical ltration generated by Ñ . Let (λ t ) be a (F t )-predictable process such that a.s. for each t > 0, λ t ≤ λ. Then the following process (N t ) is a counting process with (F t )-predictable intensity λ t :

N t = (0,t] R + 1 [0,λs/ λ] (θ) Ñ (ds, dθ).
Indeed, N is clearly a counting process because each atom of Ñ is weighted 1 or 0. Also, since a.s. t 0 λ s ds < +∞, the martingale property for the Poisson point measure Ñ ensures that N t -t

0 1 0 1 [0,λs/ λ] (θ)dθ λds = N t - t 0 λ s ds is a (F t )-local
martingale, which shows that the counting process (N t ) has intensity (λ t ).

2.3

The likelihood Since claim occurrence is known only if the claim has been reported, it is crucial to be able to write the likelihood of the data sample made of the set I R (τ ) reported claims, see Equation ( 1). The partial observation scheme we face here lies in the family of classical censoring and truncation. More precisely, the occurrence T is observed if and only if T ≤ τ -U , therefore the occurrence time is right-truncated by the random variable τ -U .

Having established the occurrence intensity for reported claims, see Equation (4), we are now ready to derive the likelihood for the data sample observed between time 0 and time τ . Denote n R τ the number of claims reported before τ . In the individual dataset, one observes the occurrences

(t R n ) n R τ n=1
and the associated reporting delays

(u R n ) n R τ n=1 . The likelihood writes, with convention T R 0 = t R 0 , U R 0 = u R 0 = 0, P ∀1 ≤ n ≤ n R τ , T R n = t R n , U R n = u R n , and T R n R τ +1 > τ = P T R n R τ +1 > τ | T n R τ = t n R τ n R τ n=1 P T R n = t R n | T R n-1 = t R n-1 P U R n = u R n | T R n = t R n = exp - τ t n R τ λ R (s)ds n R τ n=1 λ R t R n exp - t R n t R n-1 λ R (s)ds p U |t R n (u R n ) p U |t R n ([0, τ -t R n ]) = exp - τ 0 λ R (s)ds n R τ n=1 λ R t R n p U |t R n (u R n ) p U |t R n ([0, τ -t R n ])
.

From the formula for λ R given in Equation (4), the likelihood nally writes

exp - τ 0 λ(s)p U |s ([0, τ -s])ds n R τ n=1 λ t R n p U |t R n (u R n ),
which leads to the following result:

Proposition 3. The log-likelihood of the observed sample writes

- τ 0 λ(s)p U |s ([0, τ -s])ds + n R τ n=1 ln λ t R n + ln p U |t R n (u R n ) .
Remark 1. (Dependency between occurrence and reporting) An interesting feature lies in the analysis of the case where delays are independent of occurrences, that is

p U |t (du) ≡ p U (du).
In this case, the likelihood writes

exp - τ 0 λ(s)p U ([0, τ -s])ds    n R τ n=1 λ t R n       n R τ n=1 p U (u R n )    ,
and it is interesting to note that due to the bias of observing reported claims only, a dependence between the occurrence and delay distribution still exists, and is captured in the rst term exp -

τ 0 λ(s)p U ([0, τ -s])ds .
That is why in practice the calibration of claims occurrence frequencies and delay distribution has to be performed simultaneously.

2.4

Some macro results

In the claims reserving practice, reported claims are grouped into years of occurrence i and year of reporting j; in this paper, both indices are integers starting at value 1, and we consider τ to be an integer (the beginning or end of a given year, or month, etc. depending on the time scale considered). The aim of such upper triangle is, again, to forecast IBNR counts (in the lower triangle). Formally, the aggregate incremental claim amount for claims which occurred in year i and are reported in their development year j is dened as follows, (see [START_REF] Hesselager | Modelling of discretized claim numbers in loss reserving[END_REF] for an extended study on this segmentation), as it is the case in the actuarial practice: for

1 ≤ i ≤ τ and 1 ≤ j ≤ τ , X i,j = n≥1 1 Tn∈[i-1,i) 1 Tn+Un∈[i+j-2,i+j-1) = n≥1 1 (Tn,Un)∈K i,j , (6)
where

K i,j = {(t, u) : t ∈ [i -1, i), t + u ∈ [i + j -2, i + j -1)}.
The aggregate cumulative claim amount is then dened as

C i,j = j k=1 X i,k .
From Equation (6) the following representation for incremental numbers holds:

X i,j = Q(K i,j ). (7)
The micro-macro result for claims counts is stated below:

Proposition 4. Milliman R&D 9/39

• The random numbers (X i,j ) are independent,

• For 1 ≤ i ≤ τ and 1 ≤ j ≤ τ , X i,j follows a Poisson distribution with parameter

         α i,1 = i i-1 λ(t)p U |t ([0, i -t)) dt if j = 1, α i,j = i i-1 λ(t)p U |t ([i + j -t -2, i + j -t -1)) dt for j ≥ 2. (8)
Proof of Proposition 4 As the sets (K i,j ) are disjoint, then by Denition 1 (i), the random variables (X i,j ) derived in Equation ( 7) are independent. By Denition 1 (ii), the random variable Q(K i,j ) has a Poisson distribution with parameter q(K i,j ), where the intensity measure q is dened in Equation ( 2). This concludes the proof.

With the previous result, it is clear that for counts, the aggregate model is still Poisson, which largerly justies why this model is mainly used for claims counts in practical applications. It can be emphasized here that the individual arrival process has a general time-dependent intensity, and that the (time-dependent) distribution of reporting delays is not specied, which shows the remarkable micro-macro consistency feature recovered here. In the following, we detail a particular case of interest.

Let us rst recall that the original Poisson model, which leads to the same reserves as the Chain-Ladder method, is in fact specied in a separable version of Equation ( 8) such that α i,j = β i γ j with some parameters β i and γ j , see [START_REF] Wüthrich | Stochastic claims reserving methods in insurance[END_REF], Section 2.3 -the separable form is the key assumption to recover the same reserves as Chain-Ladder. A question of practical implications is whether this separable assumption is valid -in our setting, this turns out to look for 'micro' specications which lead to the separable assumption. Interestingly, assuming a time-independent delay distribution p U |t ≡ p U is not sucient to recover a separable form. This shows that in general, a micro model with an homogeneous delay distribution does not lead to the separable Poisson assumption. Indeed, in this case the Poisson parameter in Equation ( 8) can be rewritten as, by a simple change of variable (for j ≥ 2 here),

α i,j = 1 0 λ(i -1 + s)p U ([j -1 -s, j -s)) ds,
which shows that the dependency between occurrence year i and development year j lies in the integral of a product of a i-dependent and a j-dependent function. A general sucient condition for the separable form to appear (see again Hesselager (1995)), is that the occurrence intensity can be written as an occurrence-year Milliman R&D 10/39 magnitude times a function on [0, 1] driving e.g. the intra-year seasonality, which formally covers the case where there exists some reference intensity (λ 0 (s)) 0≤s≤1 and constants (c i ) 1≤i≤τ such that for any

1 ≤ i ≤ τ and 0 ≤ s ≤ 1, λ(i -1 + s) = c i λ 0 (s).
Indeed, the Poisson parameter can then be written in the separable form (for j ≥ 2),

α i,j = c i 1 0 λ 0 (s)p U ([j -1 -s, j -s)) ds.
In other words, the separable form of the aggregate Poisson model is valid if the intra-year seasonality is similar from one occurrence year to the next.

Parameter uncertainty As stated in the previous result, the number of IBNR claims can be derived, in aggregate, as a Poisson random variable with given parameter. This is the so-called process error induced by the individual model. To measure parameter uncertainty, or equivalently estimation error, one has to address the impact of estimator randomness in the resulting output.

Let us denote generally X the total future claims number (or cash ows) -aim is to estimate the reserve which corresponds to the expectation of X -the estimator of the reserve is denoted X. To properly decompose the so-called mean square error of prediction dened in the following, one needs to specify the probability measures at stake:

• P the probability measure related to the uncertainty in the parameters (estimation error) λ(t) and p U |t (du) -these are estimated based on information on reported claims provided by the random measure Q R (dt, du),

• Q the probability measure related to the random measure Q IBN R (dt, du), which drives in particular the random number of IBNR (process variance).

Since by construction of the individual model the number of reported and IBNR claims are independent, this allows us to write the mean square error of prediction under the product measure

P ⊗ Q as M SEP ( X) = E P⊗Q ( X -X) 2 , = E P⊗Q ( X -E Q [X] + E Q [X] -X) 2 , = E Q (X -E Q [X]) 2 + E P X -E Q [X] 2 + 2E P⊗Q X -E Q [X] (E Q [X] -X) ,
where the last term reduces to zero as, again by independence,

E P⊗Q X -E Q [X] (E Q [X] -X) = E P X -E Q [X] E Q [E Q [X] -X] with E Q [E Q [X] -X] = 0.
Milliman R&D 11/39

2.4

Some macro results

This shows that the mean square error of prediction writes

M SEP ( X) = V ar Q (X) + E P X -E Q [X] 2 ,
and is decomposed into a pure process variance (rst term) and a squared bias which measures parameter risk (second term).

As the likelihood has been derived in this model, a convenient way of accounting for parameter uncertainty is to rely on asymptotic properties of maximum likelihood estimators, whose variance-covariance matrix can be computed as the opposite of the Hessian inverse. Such approach is discussed here in the standard parametrization of the individual model, although it can be extended to any piecewise-constant or continuous parametric specication of each model component.

In the standard setting, we consider homogenous Poisson claims arrivals, that is λ(t) = λ, as well as exponential and time-invariant reporting delay distribution, that is p U (du) = θ exp(-θu)du. In this framework, the log-likelihood writes (see Proposition 3):

ln L(λ, θ) = n R τ ln(λ) -λτ + λ θ (1 -e -θτ ) + n R τ ln(θ) -θ n R τ n=1 u n . (9)
As stated in Remark 1, the dependency between the occurrence frequency λ and the delay parameter θ remains, which represents the observation bias, here through the middle term of the form λ θ (1 -e -θτ ). The result about the mean-square error of prediction is stated below -it follows from a direct computation of the variancecovariance matrix Σ (opposite of the Hessian inverse) and the application of the Delta method; note that the rst term refers to process error whereas the second term relates to estimation error. Proposition 5. In the standard parametrization of the individual claims model, the total mean-square error of prediction can be approximated as

M SEP ( X) = φ( λ, θ) + ∇φ( λ, θ)Σ( λ, θ) ∇φ( λ, θ) T , where φ(λ, θ) = λ θ 1 -e -θτ , then ∇φ(λ, θ) = 1 λ φ(λ, θ), -1 θ φ(λ, θ) + λτ θ e -θτ
, and the variance-covariance matrix writes

1 D(λ, θ) n R τ θ 2 + λ θ 3 e -θτ (θ 2 τ 2 + 2θτ + 2) -2 1 θ 2 e -θτ (1 + θτ ) -1 1 θ 2 e -θτ (1 + θτ ) -1 n R τ λ 2
, with the determinant D(λ, θ) given by

D(λ, θ) = n R τ λ 2 n R τ θ 2 + λ θ 3 e -θτ (θ 2 τ 2 + 2θτ + 2) -2 - 1 θ 4 1 -e -θτ (1 + θτ ) 2 .
Remark 2. The result can be adapted for any piecewise constant parameters, in which framework the same structural form can be identied.

Milliman R&D 12/39 2.5 Numerical illustration 2.5 Numerical illustration

A numerical experiment

Let us start a numerical experiment with parameters λ = 500 (in average, 500 claims occur each year) and θ = 1/3 (the reporting delay is 3 years in average). This arbitrary setting allows us to easily illustrate key facts about the toy model.

In Figure 1, the empirical distribution of observed occurrences and delays is represented. At this stage, the sample bias appears: although the initial occurrences are distributed as an homogeneous Poisson process, the occurrences of reported claims are not Poisson anymore. Therefore, it is not possible to directly t any distribution on that of occurrences or delays. Instead, by maximizing the log-likelihood (9), we are able to compute MLEs for the frequency parameter λ and the delay parameter θ, see again Figure 1. This example illustrates how the true underlying parameters relate to the observed empirical distributions. In particular, it is shown in Figure 1 that the distance to the true frequency increases with time (due to the existence of reporting delays), this distance representing IBNR claims. In addition, it is shown that the true reporting delay distribution results in transportation of mass from short to longer reporting delays compared to the empirical distribution, due to the fact that larger delays are less likely to be observed -this will be illustrated on real data in the next paragraph.

Note that based on these parameters, the number of IBNR follows a Poisson distribution with a given parameter which takes a simple form in this numerical experiment; it is indeed equal to

τ 0 λ(s) 1 -p U |s ([0, τ -s]) ds = λ θ 1 -e -θτ .
In order to give a view on estimation error, we repeat the simulation of articial data and the maximization step, producing a sample of 1000 MLEs; these are depicted in Figure 2.

A practical example on personal injury data

The real dataset 'ausautoBI8999' is obtained from the R package CASdatasets [START_REF] Dutang | CASdatasets: Insurance datasets[END_REF]). It contains information on 22,036 settled personal injury insurance claims in Australia. These claims arose from accidents occurring from July 1989 through to January 1999. Note that claims settled with zero payment are not included. For the purpose of illustration, the model is calibrated on the 3-years observation period from month 50 to month 86, while considering only claims which occurred and were reported in the period. A piecewise constant arrival intensity λ(t) on bands of 3 months is chosen to account for non-stationarity, and an empirical specication of the delay distribution is set, augmented with an exponential tail as

p U (du) = n 0 k=1 p k 1 [k-1,k) (u) + 1 [n 0 ,∞) (u) 1 - n 0 k=1 p k α exp (-α(u -n 0 )) ,
where we set n 0 = 5 the number of rst months for which an empirical distribution is taken. The parameters obtained by the maximization of the log-likelihood established in Proposition 3 are represented in Figure 3. This shows that the model predicts a relatively stable claim occurrence for the last months, although a naïve observation of the reported claims subset would not lead to the same conclusion.

As a consequence, it is shown in Figure 3 that the observed delay distribution does not match the tted one, as the likelihood accounts for unobserved (higher) delays -this is observed in the gure as the density of the smallest delays is reduced in the nal distribution and transported in the tail. In this setting, the joint parameter estimation routine leads to 1,501 IBNR claims in average as an application of formula (3) (whereas the separate t of the observed occurrence and delay distributions would have led to an expected number of 565). The dispersion of IBNR due to process error can be easily drawn based on the macro consistency, see Proposition 1 -this is shown in Figure 4. This provides an input of the same 'nature' than an aggregate model would do. More original (and very useful in several practical applications) is the possibility to simulate the times at which those IBNR would have occurred, that is to complete the observed dataset (of reported claims) by that related to IBNR. One stochastic path of IBNR occurrence times performed by the thinning algorithm described in Subsection 2.2 is illustrated in Figure 5. As expected, this scenario shows that most IBNR occurred in the last months, in coherence with the dispersion shape of the IBNR intensity in Equation [START_REF] Dutang | CASdatasets: Insurance datasets[END_REF]. In addition, it is possible to add the uncertainty on the parameters to the simulation to produce a forecast comparable to aggregate alternatives such as that of [START_REF] Mack | Distribution-free calculation of the standard error of chain ladder reserve estimates[END_REF] -this has been discussed in Subsection 2.4. The nal stochastic forecast including both process and estimation errors is represented in Figure 6.

Finally, we summarize in Table 1 some comparison with Chain-Ladder estimates and the prediction, process and estimation errors by [START_REF] Mack | Distribution-free calculation of the standard error of chain ladder reserve estimates[END_REF] on triangles with varying granularity -that is we change the aggregation step form 6 to 1 months in the triangle. The numerical computations are performed based on the R package ChainLadder [START_REF] Gesmann | ChainLadder: Statistical Methods and Models for Claims Reserving in General Insurance[END_REF]. A rst result is that the Best Estimate values remain in the same orders of magnitude, although the Chain Ladder method always leads to an overestimation of the expected IBNR compared to the individual model, and moreover that the value depends on the aggregation step (although the underlying data is similar), as it remains to formulate dierent assumptions on the claims development in each case. In addition, one notices that process error is higher for the Mack method ; the individual model takes here advantage of its pure Poisson formulation. Finally, the estimation error for the Mack model is, as expected, very high for high aggregation step as it leads to fewer data points in the triangle. While decreasing the aggregation step, one increases the number of points which makes estimation error reduces (although the number of parameters increases as well, but not as much). In all cases, the coecient of variation for the number of IBNR stays at a higher level compared to that of the individual model. 3 Payments and settlement

IBNR distribution (process and estimation errors)

IBNR

Modelling the claims development

In this part, we develop the modelling framework allowing for claims payments development over time, which enriches the previous occurrence and reporting model developed in Section 2. In the literature, several studies based on triangles proposed a model in order to account for incurred amounts (in addition to paid) which include the reserve case estimates set by expert judgement, in order to guarantee a practical constraint to match ultimate paid and incurred forecasts -see [START_REF] Mack | Munich chain ladder: A reserving method that reduces the gap between IBNR projections based on paid losses and IBNR projections based on incurred losses, variance, volume 2, number 2[END_REF] proposed an extension to paid-incurred data. In this paper, we rather focus on pure paid claims and leave the extension of the framework surveyed here for further research.

In this section, the time origin 0 corresponds to a given claim's reporting. Our aim is to describe the claims development process from reporting to closing; we rely here on the convenient specication introduced in Antonio and Plat (2014). The claims development process is made of random times (V k ); at each random time, a mark E k is generated. The mark E k can take its values in the space {1, 2, 3}; let us introduce these three states:

• E k = 1 indicates that the claim is settled at time V k , without any payment at V k • E k = 2
indicates that the claim is settled at V k with a payment,

• E k = 3 indicates that a payment occurs at V k , without settlement of the claim.

Three intensity functions h 1 , h 2 and h 3 are introduced to govern the frequency at which events of type 1, 2 or 3 occur since reporting of the claim. If an event of type 2 or 3 occurs, a payment P k is generated with distribution which may depend on the time V k and the type of event, that is on the fact that the claim is closed or not.

Denition with a marked (non-homogenous and absorbed) Poisson process

Inter-arrival times can be seen as a competing risk model; we construct recursively the sequence (V k ) with V 0 = 0 and V k+1 = k +V k , dened as follows. First introduce

h i (V k + k ) 3 j=1 h j (V k + k ) .
This then proves that for a given claim, the sequence of the (V k , E k ) occurs as a marked Poisson process with occurrence rate 3 i=1 h i (t) and mark probability

p e|t (de) = 1 3 k=1 h k (t) 3 i=1 h i (t)δ i (de),
where δ i (de) is the Dirac mass at state i.

Finally, the true process is stopped if the generated mark is 1 or 2; one can dene S the settlement delay as the stopping time S = inf{V k : E k ∈ {1, 2}}, and consider the nal development process

{(V k , E k ), V k ≤ S}.
This construction lies in the framework of competing risk models. Although it provides an intuitive representation of inter-arrival times, this framework is not well suited to derive the closed-form formulas at stake. In the next subsection, we propose a continuous time Markov representation of the claims development, which will then allow us to rely on general results by Hesselager (1994).

Denition with a Markov process

The diculty in dening the claims development through a simple Markov process lies in the fact that 'state' 3 (payments without settlement) can occur at consecutive times. To overcome this issue, let us consider (X t ) a Markov process taking values in state space N * , and introduce the following new specication: if a jump occurs at time V k , then

• X V k = 1 indicates that the claim is settled at time V k , without any payment at V k , • X V k = 2
indicates that the claim is settled at V k with a payment,

• X V k = j, for j ≥ 3, indicates that a payment occurs at V k , without settlement of the claim.

We consider that X 0 = 3 (the construction holds if we choose another integer higher than 3). Let us introduce the following transition intensities: for j ≥ 3:

• λ j,1 (t) = h 1 (t), • λ j,2 (t) = h 2 (t), • λ j,j+1 (t) = h 3 (t).
All other transition intensities are zero. Then one can prove that the construction leads to a similar development process as the one described in 3.1.1.

Milliman R&D 20/39 Of course, some constraints may be applied to h 4 to guarantee convergence and avoid innite payment cycles in the calibration and projection. A basic suggestion is to set h 4 (t) = 0 for t greater than some time t 0 since reporting, after which the claim can not be reopened.

The likelihood

We identify each claim by the index n, given by its observed occurrence T n ; for a given claim n we denote

• V (n) k
the time since reporting of the kth event in the payments development,

• S (n) the settlement delay for claim n,

• τ (n) the time during which the claims development is observed, which can be expressed as

τ (n) = min S (n) , τ -T n -U n , • E (n) k
the associated event type, taking values in {1, 2, 3},

• δ (n)
k (i) the indicator that the kth event for claim n is of type i ∈ {1, 2, 3},

• X (n) k
the payments (if any) related to the kth event,

• P (.) the probability density function of the payments distribution.

Proposition 6. The likelihood for the claims development component writes

n≥1 exp - τ (n) 0 (h 1 + h 2 + h 3 )(u)du k≥1 h 1 (V (n) k ) δ (n) k (1) h 2 (V (n) k ) δ (n) k (2) h 3 (V (n) k ) δ (n) k (3) × n≥1 k≥1 δ (n) k (1) + P (X (n) k ) δ (n) k (2) + δ (n) k (3) . ( 10
)
Remark 4. The use of the asymptotic normality of maximum likelihood estimators is still possible in this framework, in a parametric or non-parametric (piecewise constant) specication of the intensity parameters. The structure of the log-likelihood shows however that the parameter estimates for each type of event are decorrelated (in the precise present specication)-this contrasts with the occurrence and reporting parameters, as discussed in Subsection 2.3, for which we proposed a more detailed discussion.

Simulation algorithm

The simulation algorithm of each claims path lies, again, in the eld of general thinning procedures for marked population processes -we refer the reader to Boumezoued (2016a) for more details on the formalism and the related thinning representation. The algorithm is described here in the spirit of the Marked Poisson process representation of the claims path presented in 3.1.1, although it can also be used to generate general Markov models closer to the construction in 3.1.2. The algorithm works as follows for generating paths from a given time (set at zero) after reporting to any arbitrary time ν:

1. First consider bounds hi for each event-specic intensity function h i , i ∈ {1, 2, 3}, such as hi = sup t∈[0,ν] h i (t).

2. Then draw a sequence (U k ) on [0, ν] as a Poisson process with intensity h1 + h2 + h3 -for example by drawing U k+1 -U k as iid increments that are exponentially distributed with parameter h1 + h2 + h3 .

3. For each time U k , draw a Bernouilli random variable with 1-value probability (h 1 (U k ) + h 2 (U k ) + h 3 (U k )) / h1 + h2 + h3 -the 1-value corresponds to the acceptance of the U k as a time of event (1, 2 or 3) for the individual claims path.

4. To determine the type of event, draw another random variable with values in {1, 2, 3} with probability for the i th value given by h

i (U k )/ (h 1 (U k ) + h 2 (U k ) + h 3 (U k ))
5. In case of event related to payment (2 or 3), draw the payment with the appropriate distribution P (.).

Some macro results

3.4.1 Closed-form formulas for a single claim development [START_REF] Hesselager | A Markov model for loss reserving[END_REF] described the path of a given claim as a general continuous time Markov process with state space N, and derived semi-explicit rst and second order moments for X (t) (u, v), which denotes the total payments in (t + u, t + v] for a claim which occurred at time t. In the model proposed by [START_REF] Hesselager | A Markov model for loss reserving[END_REF], the claims path (after reporting) from one state to another is described by a process S (t) (u) with values in N * (we removed 0 which represents the IBNR state), with transition Milliman R&D 22/39 probabilities p (t) mn (u, v) and associated transition intensities λ (t) mn (u). Payments are made at the transition times of this dynamics; for a transition m → n at a time u since reporting, a payment Y (t) mn (u) is drawn, payments being independent from each other, with mean y (t) mn (u) and standard deviation σ (t) mn (u). In the following analysis of a single claim, we omit without any restriction the index t. We will rely on the following results by [START_REF] Hesselager | A Markov model for loss reserving[END_REF].

Proposition 7. For all j ∈ N * , the mean and variance of the payments related to the claim in state j after u times units of development can be expressed as

E [X(u, ∞) | S(u) = j] = m∈N * n∈N * \{m} ∞ u p jm (u, v)λ mn (v)y mn (v)dv, (11)
Var (X(u, ∞) | S(u) = j) = m∈N * n∈N * \{m} ∞ u p jm (u, v)λ mn (v) σ 2 mn (v) + r 2 mn (v) dv, (12)
where

r mn (v) = y mn (v) + E [X(v, ∞) | S(v) = n] -E [X(v, ∞) | S(v) = m] . (13)
In the next paragraphs, we derive the formulas for the mean and variance (process error) under our particular (non-homogeneous) Markovian setting as described in Subsection 3.1.2. Let us recall that in this setting, intensities can be set at zero, except for the following with j ≥ 3: λ j,1 (t) = h 1 (t) (settlement without payment), λ j,2 (t) = h 2 (t) (settlement with payment) and λ j,j+1 (t) = h 3 (t) (payment without settlement).

Expected payments along a claim path In our setting, it is clear that X(u, ∞) is zero given that the claim is settled, that is given that S(u) = 1 or S(u) = 2. Moreover, without any change in the path distribution, it is sucient to study X(u, ∞) given that S(u) = 3. Based on Proposition 7, we then get after some computation:

E [X(u, ∞) | S(u) = 3] = ∞ u h 2 (v)y 2 (v) + h 3 (v)y 3 (v) m≥3 p 3m (u, v) dv.
Note that m≥3 p 3m (u, v) is the probability that the Markov process stays in the set {3, 4, ....}, whose transition to states 1 and 2 are h 1 and h 2 respectively ; therefore

m≥3 p 3m (u, v) = exp - u s (h 1 (v) + h 2 (v)
)dv -this leads to the following result: Proposition 8. The expected payments after some time s for a single claims path can be written as

µ(s) = E [X(s, ∞) | S(s) = 3] = ∞ s y 2 (u)h 2 (u) + y 3 (u)h 3 (u) exp - u s (h 1 (v) + h 2 (v))dv du. (14)
Milliman R&D 23/39

3.4 Some macro results

Remark 5. An alternative way of computing the quantities µ

k (u) := E [X(u, ∞) | S(u) = k],
as suggested in Hesselager (1994) (see Remark 3.1), is to rely on Thiele's dierential equation given as

µ k (s) = - j∈N * \{k} λ kj (s)(y kj (s) + µ j (s) -µ k (s)).
This equation, stated here in a general Markov framework, provides a convenient numerical scheme for computing transition probabilities in more general individual claims models. This can be proved based on the representation given in Equation (11) -as the proof does not appear in [START_REF] Hesselager | A Markov model for loss reserving[END_REF], we propose to detail it in Appendix to guarantee self-consistency of the results included in the present survey.

In our case, denoting µ 3 ≡ µ, Thiele's dierential equation reduces to

µ (s) = -y 2 (s)h 2 (s) + y 3 (s)h 3 (s) + h 1 (s) + h 2 (s) µ(s), (15)
whose solution satises for s ≤ t

µ(s) = µ(t) exp - t s (h 1 (u) + h 2 (u))du + t s y 2 (u)h 2 (u) + y 3 (u)h 3 (u) exp - u s (h 1 (v) + h 2 (v))dv du. (16)
Letting t → ∞, we recover the closed form in Equation ( 14) for the expected payments on a single claim path.

Remark 6. In the homogeneous Markovian setting, under which h i (s) ≡ h i and y i (s) ≡ y i , the expected payments along a claims path write

µ(s) = µ(0) = y 2 h 2 + y 3 h 3 h 1 + h 2 .
Variance of claims payments From Equation (12), and by analogous reasoning as in the previous part on expected payments, one can state the following result:

Proposition 9. The variance of total payments after some time s along a single claims path can be written as

γ(s) = Var (X(s, ∞) | S(s) = 3) = ∞ s H(u) exp - u s (h 1 (v) + h 2 (v))dv du. ( 17
)
where the function H, which depends on time since reporting, is given by

H = h 1 µ 2 + h 2 σ 2 2 + (y 2 -µ) 2 + h 3 σ 2 3 + y 2 3 . (18)
Milliman R&D 24/39

(u) = Var (X(u, ∞) | S(u) = k), then γ k (s) = - j∈N * \{k} λ kj (s) σ 2 kj (s) + r 2 kj (s) + γ j (s) -γ k (s) ,
where r kj (s) is dened in Equation (13). Denote γ(s) ≡ γ 3 (s); in our case, Thiele's dierential equation reduces to

γ (u) = -H(u) + γ(u) (h 1 (u) + h 2 (u)) ,
which allows to recover the result in Proposition 9.

Remark 8. In the homogeneous Markovian setting, under which h i (s) ≡ h i and y i (s) y i , the variance of the total payments along a claims path writes

γ(s) ≡ γ = h 1 µ 2 + h 2 (σ 2 2 + (y 2 -µ) 2 ) + h 3 (σ 2 3 + y 2 3 ) h 1 + h 2 ,
where

µ = y 2 h 2 + y 3 h 3 h 1 + h 2 .
3.4.2 Closed-form formulas for the claims population IBNR reserve Recall that times and reporting delays are given as a Poisson point measure Q(dt, du) with intensity measure q(dt, du) = λ(t)dt ⊗ p U |t (du). Therefore the ultimate payments for IBNR claims at observation time τ can be expressed as

X IBN R τ := τ 0 ∞ τ -s X (s) (0, ∞)Q(ds, du). (19)
The expected payments are therefore given as

E X IBN R τ = τ 0 ∞ τ -s µ (s) (0)λ(s)p U |s (du)ds,
where µ (s) (0) ≡ µ(0) is expressed in Equation ( 14). This leads to the following result:

Proposition 10. The total IBNR reserve writes

E X IBN R τ = τ 0 dsλ(s) ∞ τ -s p U |s (du) × ∞ 0 y 2 (u)h 2 (u) + y 3 (u)h 3 (u) exp - u 0 (h 1 (v) + h 2 (v))dv du. (20)
Based on the representation (19), we can also derive the variance of the IBNR payments as Remark 9. Note that this is not a trivial result; it holds according to the fact that the family of processes X (s) (u, ∞) for occurrence times s is independent from the Poisson point measure Q. To prove it, let us use this independence and construct an extended Poisson measure M (ds, du, d X) with intensity measure m(ds, du, d X) = λ(s)ds ⊗ p U |s (du) ⊗ p X|s,u (d X), where we use a general notation p X|s,u (d X) for the density distribution of the R-valued random variable X (s) (0, ∞). Based on this Poisson point measure, the following representation holds:

Var X IBN R τ = τ 0 ∞ τ -s E X (s) (0, ∞) 2 λ(s)p U |s (du)ds.
X IBN R τ = τ 0 ∞ τ -s R XM (ds, du, d X), (22)
and from the classical properties of Poisson point measure, one gets

Var X IBN R τ = τ 0 ∞ τ -s R X2 m(ds, du, d X),
which then reduces to Equation (21).

We therefore get the following result:

Proposition 11. The variance of the IBNR reserve (process error) writes

Var X IBN R τ = (γ(0) + µ(0) 2 ) τ 0 ∞ τ -s λ(s)p U |s (du)ds, (23)
where µ(0) and γ(0) can be computed from Equations ( 14) and ( 17) respectively.

Based on the representation (22), it is also possible to compute the Laplace transform of the total IBNR payments, which involves the distribution p X|s,u (d X) as a whole. Indeed, from Poisson point measures properties as well, one gets Proposition 12. The Laplace transform of the IBNR reserve writes

E exp -θX IBN R τ = τ 0 ds ∞ τ -s p U |s (du) R p X|s,u (d X) exp 1 -exp -θ X .

RBNS reserve

The future payments for the RBNS claims write

X RBN S τ = τ 0 τ -s 0 X (s) (τ -s -u, ∞)1 S (s) (τ -s-u) ∈{1,2} Q(ds, du).
Let us denote by F τ the information about reported claims (including RBNS) up to time τ ; then the RBNS reserve writes

E X RBN S τ | F τ = τ 0 τ -s 0 E X (s) (τ -s -u, ∞)1 S (s) (τ -s-u) ∈{1,2} | F τ Q(ds, du).
Using the independence between claims paths and the Markov property of each single path, we get the following result Milliman R&D 26/39

3.5

Numerical illustration Proposition 13. The RBNS reserve writes

E X RBN S τ | F τ = τ 0 τ -s 0 µ(τ -s -u)1 S (s) (τ -s-u) ∈{1,2} Q(ds, du), (24)
where µ(.) is given in Equation ( 14).

Finally, the variance of RBNS future payments is given in the following proposition, as that of a sum (integral) of independent trajectories: Proposition 14.

Var X RBN S τ | F τ = τ 0 τ -s 0 γ(τ -s -u)1 S (s) (τ -s-u) ∈{1,2} Q(ds, du).

Numerical illustration

The numerical illustration is based on a real portfolio, in which 8 groups of claims can naturally be identied by means of an observed covariate. For purpose of condentiality, the covariates possible values are labeled from A to H. We represent in Figure 7 the observed occurrence times and reporting delays used for calibration (all groups merged in the illustration). The sample is made of 19,870 observed claims, split into 6,121 settled claims and 13,749 RBNS claims. The occurrence and reporting component of the model is taken as the standard parametrization in each group -that is, a constant occurrence frequency is chosen (which depends on claims covariate) as well as a constant reporting parameter (which again depends on claims group) for the exponentially distributed delay. The estimated occurrence and reporting parameters are illustrated in Figure 8 and compared to the 'naive' estimates based on a separate t of frequency and reporting parameters. Note that this is the average reporting delay (in years) which is represented here (the inverse of each claims group reporting parameter). As expected, the joint estimation of occurrence and reporting parameters leads to a higher occurrence frequency and a higher average delay compared to a separate estimation.

Based on these parameters, it is possible to derive the distribution of the number of IBNR -this is given in Figure 9 (with the associated prediction error) and compared with the Mack model (with Gaussian distribution assumption). Again here, similar conclusions hold as in 2.5.2: the Mack model provides a higher IBNR (unitary) reserve (4,485) compared to the individual claims model (4,334), while providing a larger total condence interval. Note that here the Mack method is applied to the full triangle (all 8 groups merged) in order to guarantee reasonable forecasts -note also that the joint application of the Mack method to both 8 triangles would require an extended framework to account for correlations, see [START_REF] Braun | The prediction error of the chain ladder method applied to correlated run-o triangles[END_REF] and [START_REF] Merz | Prediction error of the chain ladder reserving method applied to correlated run-o triangles[END_REF].

Let us now focus on the claims payment streams, and recall that the overall events frequency in the claims path (whatever their type) is driven by h 1 + h 2 + h 3 and that the weight of each parameter relates to the event probability, as described in 3.1.1. We depict the frequency parameters h 1 , h 2 and h 3 in Figure 10, which are assumed to be constant in a rst step of the analysis. This rst indicates that both settlement frequencies and payments without settlement are of the same order of magnitude, although one can identify in more details that group C is characterized by both a high settlement frequency (shorter claims history) and a low payment frequency (more time between two payments). This contrasts with group B, with a larger time to settlement, and more frequent payments during the claims paths. In between, claims such as that of group F are characterized by a medium average time to settlement and times between payments, compared to the other groups. Note that as almost no payments at settlement occurred, the estimates for parameter h 2 are close or equal to zero, depending on the group considered.

To rene the analysis, we propose to t yearly piecewise constant parameters h 1 , h 2 and h 3 , to better account for the timing of the kind of events which are observed here. Note that maximum likelihood estimators in this piecewise constant setting can be derived in closed form based on the likelihood representation (10) -the calibration results are depicted in Figure 11. This setting appears quite informative (except for h 2 whose values are negligible) as a proper structure can be exhibited among the several claims groups considered, namely:

• The settlement frequency h 1 is low for the two rst years, then increases for years 3 and 4, then a slight decrease is observed and nally an increase or decrease in the latest years up to 10, depending on the group considered. In particular, the overall increase of the settlement intensity from earliest years of development to the latest ones means that the probability for an IBNR claim to be settled in its rst year of development is lower than that of an RBNS claim with already 2 years of development, which in comparison is more likely to be settled. Therefore such conclusions contradict a simplied homogeneous Markov framework, which would lead in particular to consider that expected future payments are the same for RBNS and IBNR claims, see Remarks 6 and 8. Note however that the illustrative computations for the Markov case can serve as a building block for the piecewise constant intensity model (worth mentioning that it is not homogeneous Markov).

• The payment (without settlement) frequency h 3 shows a dierent pattern, as it is high for the rst year, then deceases to the second year, and nally increases (overall) up to the latest year. This means in particular that, the settlement intensity being given (let us think about it as constant for the interpretation), payments are more frequent in the rst and the last year, and increases from the second year to the last year. Also, this means that e.g. payments are more likely to occur (still in this interpretation framework) in the rst year of an IBNR claim compared to the second year of a one-year old RBNS claim.

As for the observed payments, we also account for the covariate information and t a separate distribution for each claim group. Moreover, due to the nature of the payments distribution, we use a mixture of a log-normal distribution and an exponential in the tail -the density is written as

P (x) = p 1 xσ √ 2π exp - (ln x -µ) 2 2σ 2 1 x≤x 0 + (1 -p)γ exp (-γ(x -x 0 )) 1 x>x 0 ,
with x 0 a xed threshold separating attritional and large claims, which depends on the group considered -the results are depicted in Figure 12. The frequency and payments parameters being tted, this allows us to compute by closed-form µ(s) and γ(s), the expected payments and related variance along a single claim path conditional on development length, see Equations ( 14) and (17) in Subsection 3.4. The results are depicted in Figure 13 -they can be described as follows: as a function of claim development, the expected payments decrease in the 3 rst years, due to the combination of a decrease of the payment frequency and an increase (spike) in the settlement frequency h 1 (see again Figure 11); then the expectation increases again, due to the reduction in the settlement frequency after the third year barrier is attained. This illustrates the benet in analyzing the timing of claims development, and in particular the bias when considering simple frequency × cost models.

On this topic, the computation of the IBNR and RBNS reserves is informative, which can be performed in closed-form based on Equations (20) and (24). Indeed, although the expected number of IBNR represents around 20% of the claims numbers to be paid (both RBNS and IBNR), which would be the ratio in terms of reserve in a homogeneous Markov setting, in terms of amounts this represents here around 50%, which illustrates the importance of quantifying IBNR claims accurately for lines of business with extended reporting delays -this is depicted in Figure 14.

Also, the comparison with the Chain-Ladder prediction is given in Figure 15. As for the IBNR unitary reserve (see 2.5.2), the Chain-Ladder method provides higher reserves for all (except one) groups.

Finally, we propose in Figure 16 a detail of the coecients of variation (CV) for estimation and process errors for each claim group, and we compare the results of the individual model with that of the Mack model. First, the graph shows a reduction in process variance in the individual model, whose stochastic formulation is dierent (and apparently less overdispersed) compared to the original discrete time Markov model of [START_REF] Mack | Distribution-free calculation of the standard error of chain ladder reserve estimates[END_REF]. More importantly, it appears that the use of the individual claims model provides a great reduction in estimation error when Milliman R&D 29/39 payments are modelled, which strengthens the results obtained in 2.5.2 for IBNR claim number. This key conclusion, which has also been observed in other studies as referenced in this survey, can be interpreted based on two facts: rst, the more natural formulation seems to allow for a better t of the several individual model building block compared to the Mack model, and second, the larger amount of data available for each part of claims history combined with this smart parametrization makes estimators have a low variance. On this basis, we argue that such modelling framework represents a promising (and complementary) method to assess future claims development with improved accuracy.

Concluding remarks

In this survey, we proposed a unied framework to address the stochastic modelling of individual claims for reserving purposes, based on the earliest developments of [START_REF] Norberg | Prediction of outstanding liabilities in non-life insurance[END_REF] and [START_REF] Hesselager | A Markov model for loss reserving[END_REF] in this direction. The model is built upon two major components, driving the occurrence and reporting process, as well as the claims payments path until settlement. It is shown in this paper how such framework can be used for calibration, risk analysis and simulation in a variety of contexts based on a coherent modeling framework. We argue that individual-based methods, which allow to analyze, understand and forecast risk in a more precise way compared to triangle-based approaches while taking advantage of detailed individual information, is worth implementing in current risk management practices. Further work remains to be done on the exploration of the numerous potentialities of individual models, in particular by calibrating more exible models involving advanced dependence structures between occurrence, reporting and payments (timing and amount). Also, the comparison with classical methods based on run-o triangles suggests mathematical perspectives which are out of scope of the present survey, but needed to continue improving our understanding of the link with more classical models.

As for the second term, let us focus on that in Equation (26) which can be rewritten as This concludes the proof.

Milliman R&D 39/39

Figure 1 :Figure 2 :

 12 Figure 1: Distribution of observed occurrences (left) and delays (right) -each graph includes the underlying theoretical distribution in red.

Figure 3 :Figure 4 :

 34 Figure 3: Observed and tted distribution related to the subsample extracted from the dataset 'ausautoBI8999' is obtained from the R package CASdatasets.

Figure 5 :

 5 Figure 5: One stochastic path of IBNR occurrence times.

  , Merz and Wüthrich (2010), Happ and Wüthrich (2013), as well as Miranda et al. (2012) and Martínez Miranda et al. (2015). From an individual perspective, in the spirit of the modelling framework developed in Pigeon et al. (2013), Pigeon et al. (2014)

Figure 7 : 39 Figure 9 : 39 Figure 10 :

 73993910 Figure 7: Occurrence times and reporting delays used for calibration

Figure 11 :Figure 12 :

 1112 Figure 11: Computation of the frequency parameters h 1 and h 3 driving settlement and payments respectively, in a time dependent setting.

Figure 13 :

 13 Figure 13: Expectation (left) and variance (right) of future payments on a single claim path, conditional on claim group and development length.

Figure 14 :

 14 Figure 14: Proportion of IBNR claims in all claims to be payed (both IBNR and RBNS): numbers and amounts.

Figure 15 :

 15 Figure 15: Comparison with the Chain-Ladder method

-

  k∈N * \{j} λ jk (s)(R(s | k) -R(s | j)), = -k∈N * \{j} λ jk (s)

  reserve has to be set to cover the future development of both reported (but not settled) claims, as well as IBNR claims. The prediction of IBNR counts and occurrences is a rst step to achieve the latter objective. The overall reserve computation is adressed in the next Section 3 which enriches the individual model for claims payments development.Poisson point measure representation The Poisson measure formalism is a powerful tool to analyze quantities based on Poisson times and associated marks. We recall below the denition of Poisson point measures. Denition 1. A Poisson point measure Q(ds, dy) with intensity measure q(ds, dy) = dsµ s (dy) on R + × E (with µ s some sigma-nite measure) is a random measure taking values in N ∪ {+∞} satisfying:

	2.1	Stochastic individual claims model

A

Table 1 :

 1

	number

Best estimate of IBNR number and related process, estimation and prediction errors (standard deviations). The Mack model is applied to a triangle with varying aggregation window, expressed in months in parenthesis.

  Consider the random times independent from each other for any kind of event and any k. Then construct k = min(S 1 ) and specifyE k such that E k = i if and only if k = S i k . Standard results ('clock lemma') show that the rate of k is 3 i=1 h i (V k + t), and that conditionally on V k and k , the probability that E k = i is

	k , S 2 k , S 3

• S 1 k a random time with rate h 1 (V k + t), associated with event 1,

• S 2 k a random time with rate h 2 (V k + t), associated with event 2,

• S 3 k a random time with rate h 3 (V k + t), associated with event 3. k

  Dealing with reopenings) Once stated in a general Markov setting with countable state space, it is possible to enrich the model to account for other events than occurrence, reporting, payments and settlement. A major issue in many datasets is driven by reopenings, which are observed on several claims path. A specic reopening frequency h 4 (t) can be introduced, which allows to go from settlement (state 1 or 2) to new payments (state 3) leading to an other cycle of new payments (states 4 and more) until settlement. In the model, this amounts to consider that transition intensities λ 1,3 (t) and λ 2,3 (t) are non zero, equal to h 4 (t).

	3.2	The likelihood
	Remark 3. (

  As in Remark 5, one can alternatively use Thiele's dierential equation satised by the variance of payments; denote γ k

	3.4	Some macro results
	Remark 7.

  km (s, t) -p jm (s, t)) λ mn (t)f mn (t)dt, * \{j} λ jk (s) (p km (s, t) -p jm (s, t)) dt. Now, let us compute the following dierential:

		∞
	m∈N * n∈N * \{m} ∞ (p = -s λ mn (t)f mn (t)
	m∈N * n∈N * \{m}	s
	∂p jm (s, t)
	∂s
	= lim

k∈N h→0 p jm (s + h, t) -p jm (s, t) h = lim h→0 p jm (s + h, t) -k∈N * \{j} p jk (s, s + h)p km (s + h, t) h = k∈N * \{j} lim h→0 p jk (s, s + h) h (p jm (s + h, t) -p km (s + h, t)) = k∈N * \{j} λ jk (s) (p jm (s, t) -p km (s, t)) .

Milliman R&D

Appendix

In this appendix, we propose the proof of Thiele's dierential equation (discussed in 3.4) satised by any prospective reserve on a (general) individual Markov model, as introduced by Hesselager (1994).

Proposition 15. For any set of continuous maps f mn (.) for all n = m, the prospective reserve conditioned on state j at time s, where the rst term reduces to -k∈N * \{j} λ jk (s)f jk (s), since p jm (s, s) = 1 m=j .

Milliman R&D


Boris Nectoux
email: boris.nectoux@enpc.fr

Sharp estimate of the mean exit time of a bounded domain in the zero white noise limit

We prove a sharp asymptotic formula for the mean exit time from an open bounded domain D ⊂ R d for the overdamped Langevin dynamics

.

Setting and main results

Let us consider (X t ) t≥0 the stochastic process solution to the overdamped Langevin dynamics in R d :

dX t = -∇f (X t )dt + √ 2ε dB t , (1)
where f ∈ C ∞ (R d , R) is the potential function, ε > 0 is the temperature and (B t ) t≥0 is a standard d-dimensional Brownian motion. The overdamped Langevin dynamics can be used for instance to describe the motion of the atoms of a molecule or the diffusion of impurities in a crystal (see for instance [START_REF] Matkowsky | The exit problem: a new approach to diffusion across potential barriers[END_REF]Sections 2 and 3] or [START_REF] Chandrasekhar | Stochastic problems in physics and astronomy[END_REF]). One of the major issues when trying to have access to the macroscopic evolution of the system from simulations made at the microscopic level is that the process ( 1) is metastable: it is trapped during long periods of time in some regions of the configuration space. This implies that it typically reaches a local equilibrium of these regions long before escaping from them. These regions are called metastable regions (see [START_REF] Bovier | Metastability: a potential-theoretic approach[END_REF]Chapter 8]) and the move from one metastable region to another is typically associated with a macroscopic change of configuration of the system. The average time it takes for the process [START_REF] Berglund | Kramers' law: validity, derivations and generalisations[END_REF] to leave a metastable region is given by the Eyring-Kramers formula (see [START_REF] Hänggi | Reaction-rate theory: fifty years after kramers[END_REF]). In this work, we would like to prove, in a typical geometric setting (see [H-D] below), that the average time it takes for the process [START_REF] Berglund | Kramers' law: validity, derivations and generalisations[END_REF] to leave a metastable region satisfies in the small temperature regime (ε → 0) a kind of Eyring-Kramers formula even in the degenerate case when arg min ∂D f does not consists of a finite number of non degenerate critical points of f | ∂D .

To this end, let us consider a C ∞ bounded open set D ⊂ R d and introduce

τ D c = inf{t ≥ 0|X t ∈ D c } (2)
where D c = R d \ D, the first exit time from D. The framework we consider in this work is the following:

Assumption [H-D]: D ⊂ R d is a C ∞ bounded open set and f ∈ C ∞ (R d , R).
The function f satisfies ∂ n f > 0 on ∂D (where n is the unit outward normal to ∂D). Moreover, f has a unique critical point x 0 in D which is non degenerate and which satisfies f (x 0 ) = min D f .

Under the assumption [H-D],

it is proved in [13, Theorem 4.1] that for any x ∈ D:

lim ε→0 ε log E x [τ D c ] = min ∂D f -f (x 0 ).
See also [START_REF] Mathieu | Spectra, exit times and long time asymptotics in the zero-whitenoise limit[END_REF] for the study of the asymptotic behaviour of the law of ε log τ D c in the limit ε → 0 and [START_REF] Matkowsky | The exit problem: a new approach to diffusion across potential barriers[END_REF] for formulas obtained with formal computations. Let us mention [START_REF] Berglund | Kramers' law: validity, derivations and generalisations[END_REF] for a review of the different techniques used to obtain asymptotic estimates on the mean exit time from a domain in the limit ε → 0 in various geometric settings and for an extension of the Eyring-Kramers formulas in some degenerate cases when D = R d . In this paper, we prove a sharp asymptotic formula on the mean exit time from D in the limit ε → 0. Our main result is the following.

Theorem 1. Let us assume that the assumption [H-D] holds. Then, for any compact set K ⊂ D, it holds in the limit ε → 0 and uniformly with respect to x ∈ K:

E x [τ D c ] = (2πε) d 2 det Hess f (x 0 ) ∂D ∂ n f (σ)e -1 ε f (σ) dσ e -1 ε f (x 0 ) 1 + O(ε) ,
where dσ is the Lebesgue measure on ∂D.

Remark 1. Under some assumption on f | ∂D , an asymptotic estimate of the term

∂D ∂ n f (σ)e -1 ε f (σ)
dσ in the limit ε → 0 can be obtained with Laplace's method.

As a consequence of Theorem 1, one obtains an estimate in the limit ε → 0 on the first eigenvalue of the infinitesimal generator of the diffusion ( 1)

L ε = -ε∆ + ∇f • ∇. (3)
with homogeneous Dirichlet boundary conditions on ∂D. Let us recall that since

D ⊂ R d is a C ∞ bounded open set and f ∈ C ∞ (R d , R), the operator L ε with domain H 2 (D) ∩ H 1 0 (D) on L 2 (D, e -f (x)
ε dx) is self-adjoint, positive and has compact resolvent, where L 2 (D, e -f (x) ε dx) is the completion of the space C ∞ (D) for the norm

φ ∈ C ∞ (D) → D |φ| 2 e -1 ε f .
Its smallest eigenvalue is denoted by λ ε > 0. Theorem 1 together with [8, Corollary 1] (which is recalled in Section 2.2 below) imply the following estimates on λ ε .

Corollary 1. Let us assume that the assumption [H-D] holds. Then, in the limit ε → 0:

λ ε = det Hess f (x 0 ) ∂D ∂ n f (σ)e -1 ε f (σ) dσ (2πε) d 2 e 1 ε f (x 0 ) 1 + O(ε) .
Let us mention that sharp estimates of the smallest eigenvalues of L ε have been obtained in [START_REF] Di Gesù | The exit from a metastable state: concentration of the exit point on the low energy saddle points[END_REF][START_REF] Di Gesù | Sharp asymptotics of the first exit point density[END_REF][START_REF] Helffer | Quantitative analysis of metastability in reversible diffusion processes via a Witten complex approach: the case with boundary[END_REF][START_REF] Lelièvre | Low temperature asymptotics for quasistationary distributions in a bounded domain[END_REF] in the Dirichlet case and in [START_REF] Peutrec | Small eigenvalues of the Neumann realization of the semiclassical Witten Laplacian[END_REF] in the Neumann case when f | ∂D is a Morse function (i.e. when all the critical points of f | ∂D are non degenerate). When D = R d , we refer to [START_REF] Bouchet | Generalisation of the eyring-kramers transition rate formula to irreversible diffusion processes[END_REF][START_REF] Bovier | Metastability in reversible diffusion processes. I. Sharp asymptotics for capacities and exit times[END_REF][START_REF] Bovier | Metastability in reversible diffusion processes. II. Precise asymptotics for small eigenvalues[END_REF][START_REF] Helffer | Quantitative analysis of metastability in reversible diffusion processes via a Witten complex approach[END_REF][START_REF] Landim | A Dirichlet and a Thomson principle for nonselfadjoint elliptic operators, metastability in non-reversible diffusion processes[END_REF][START_REF] Michel | About small eigenvalues of Witten laplacian[END_REF]. Corollary 1 gives a general formula on the asymptotic estimate of λ ε which allows in particular, under the assumption [H-D], to deal with the case when f | ∂D is not a Morse function. For example, direct consequences of Theorem 1 are the following:

• Let us assume that f is constant on ∂D: f (z) ≡ f 1 for all z ∈ ∂D. Then, for any compact set K ⊂ D, it holds:

E x [τ D c ] = (2πε) d 2 det Hess f (x 0 ) ∂D ∂ n f (σ)dσ e 1 ε (f 1 -f (x 0 )) 1 + O(ε) ,
in the limit ε → 0 and uniformly with respect to x ∈ K. Moreover, one has in the limit ε → 0

λ ε = det Hess f (x 0 ) ∂D ∂ n f (σ)dσ (2πε) d 2 e -1 ε (f 1 -f (x 0 )) 1 + O(ε) .
• Let us assume that there exists k ∈ N * such that arg min ∂D f = {z 1 , ..., z k } and for all j ∈ {1, ..., k}, z j is a non degenerate critical point of f | ∂D . Then, for any compact set K ⊂ D, it holds:

E x [τ D c ] = √ πε k j=1 det Hess f | ∂D (z j ) ∂ n f (z j ) det Hess f (x 0 ) e 1 ε (f (z 1 )-f (x 0 )) 1 + O(ε)
in the limit ε → 0 and uniformly with respect to x ∈ K. Moreover, one has in the limit ε → 0

λ ε = 1 √ πε k j=1 ∂ n f (z j ) det Hess f (x 0 ) det Hess f | ∂D (z j ) e -1 ε (f (z 1 )-f (x 0 )) 1 + O(ε) .
In particular, if f | ∂D is a Morse function, one recovers the results of [START_REF] Di Gesù | The exit from a metastable state: concentration of the exit point on the low energy saddle points[END_REF][START_REF] Di Gesù | Sharp asymptotics of the first exit point density[END_REF][START_REF] Helffer | Quantitative analysis of metastability in reversible diffusion processes via a Witten complex approach: the case with boundary[END_REF] on the first eigenvalue λ ε .

2 Change of coordinates in a neighborhood of ∂D

In this section, one constructs coordinates which will be useful for the computations in Section 4. The construction of these coordinates heavily depends on the assumption

∂ n f > 0 on ∂D.
In all this section, we assume that the assumption [H-D] is satisfied.

Eikonal solution near ∂D

Let us start with the following lemma.

Lemma 1. Let us assume that the assumption [H-D] holds. Then, there exists a neighborhood of ∂D in D, denoted by V ∂D , such that there exists

Φ ∈ C ∞ (V ∂D , R) satisfying      |∇Φ| 2 = |∇f | 2 in D ∩ V ∂D Φ = f on ∂D ∂ n Φ = -∂ n f on ∂D. (4)
Moreover, one has the following uniqueness results

: if Φ is a C ∞ real valued function defined on a neighborhood Ṽ of ∂D satisfying (4), then Φ = Φ on Ṽ ∩ V ∂D . Finally, V ∂D can be chosen such that Φ > f on V ∂D \ ∂D and ∇(Φ -f ) = 0 on V ∂D .
Proof. Let z ∈ ∂D. Using [11, Theorem 1.5] or [12, Section 3.2] and thanks to the fact that ∂ n f > 0 on ∂D, there exists a neighborhood of z in D, denoted by

V z , such that there exists Φ ∈ C ∞ (V z , R) satisfying      |∇Φ| 2 = |∇f | 2 in D ∩ V z Φ = f on ∂D ∩ V z ∂ n Φ = -∂ n f on ∂D ∩ V z .
Moreover, V z can be chosen such that the following uniqueness result holds: if a function Φ ∈ C ∞ (V z , R) satisfies the previous equalities, then Φ = Φ on V z . Now, one concludes using the fact that ∂D is compact and can thus it can be covered by a finite number of these neighborhoods (V z ) z∈∂D . Finally, since

∂ n (Φ -f ) = -2∂ n f < 0 on ∂D, V ∂D can be chosen such that Φ > f on V ∂D \ ∂D and ∇(Φ -f ) = 0 on V ∂D .

Definition of the coordinate x d

In this section, one defines coordinates near ∂D which will be convenient in the upcoming computations in Section 3. Let us now consider Φ the solution to (4) on the neighborhood V ∂D of ∂D as introduced in Lemma 1. Let us define on V ∂D :

f + = f + Φ 2 and f -= Φ -f 2 . (5)
Using Lemma 1, it holds on V ∂D \ ∂D: f -> 0 and one has on V ∂D :

∇f -• ∇f + = 0. ( 6
)
Let us now consider δ > 0 such that

V δ := {x ∈ D, 0 ≤ f -(x) ≤ δ} ⊂ V ∂D .
For any x ∈ V δ , the dynamics

   γ x (t) = - ∇f - |∇f -| 2 (γ x (t)) γ x (0) = x (7)
is well defined (since from Lemma 1, one has on V ∂D , ∇f -= 0) and is such that γ x (t x ) ∈ ∂D, where t x = inf{t, γ x (t) ∈ ∂D}. This is indeed a consequence of the fact

that d dt f -(γ x (t)) = -1 < 0 on [0, t x ).
Proposition 1. The application

Θ : V δ → ∂D × [0, δ] x → (γ x (t x ), t x ) defines a C ∞ diffeomorphism. The inverse application of Θ is Ψ : (z, x d ) ∈ ∂D × [0, δ] → γ z (-x d ).
Remark 2. Let us mention that the application Ψ has been introduced locally in [START_REF] Helffer | Quantitative analysis of metastability in reversible diffusion processes via a Witten complex approach: the case with boundary[END_REF] and have also been used in [START_REF] Di Gesù | Sharp asymptotics of the first exit point density[END_REF].

Let us now give some properties of the function Ψ which are used in the sequel. Using the fact that Ψ(z, x d ) = γ z (-x d ), one obtains that for all z ∈ ∂D and x d ∈ [0, δ]:

∇ x d Ψ(z, x d ) = d dx d γ z (-x d ) = ∇f -(z, x d ) |∇f -(z, x d )| 2 . (8)
Thus, one has for all z ∈ ∂D:

∇ x d Ψ(z, 0) = - 1 ∂ n f (z, 0) n, ( 9
)
where n is the unit outward normal to ∂D. Moreover, using the fact that Ψ(z, 0) = (z, 0) for all z ∈ ∂D and n = -∇x d |∇x d | together with [START_REF] Di Gesù | The exit from a metastable state: concentration of the exit point on the low energy saddle points[END_REF], it holds for all u ∈ T z ∂D and for all v ∈ R:

dΨ (z,0) (u + vn) = u + v ∂ n f (z, 0) n, (10)
and thus: jac

Ψ(z, 0) = 1 ∂ n f (z, 0) , (11)
where jac Ψ is the determinant of the jacobian matrix of Ψ. Finally, by construction (since ρ k (y) = 1 [START_REF] Freidlin | Random Perturbations of Dynamical Systems[END_REF] such that for all k ∈ {1, ..., N }, there exist smooth coordinates x ∈ R d-1 defined by a

d dt f -(γ x (t)) = -1) x d (x) = f -(x) and one has {x d = 0} = ∂D, {x d > 0} = D∩V δ and V δ = x = Ψ(z, x d ) ∈ D, 0 ≤ x d ≤ δ . (12
C ∞ mapping Γ k : supp ρ k → R d-1 z → x . ( 14
)
The coordinates x ∈ Γ k (supp ρ k ) are then extended in a neighborhood of supp ρ k in D, as constant along the integral curves of

γ (t) = ∇f - |∇f -| 2 (γ(t)), for t ∈ [0, δ]. The function x → (x , x d ) (where, we recall, x d (x) = f -(x)) thus defines a smooth system of coordinates in a neighborhood V k of supp ρ k in D. Let us define for all (x , x d ) ∈ Γ k (supp ρ k ) × [0, δ], Υ k (x , x d ) := Ψ Γ -1 k (x ), x d ( 15
)
where Ψ is introduced in Proposition 1. Notice that it holds for all (x

, x d ) ∈ Γ k (supp ρ k )× [0, δ], Jac Υ k (x , x d ) := Jac Ψ Γ -1 k (x ), x d Jac Γ -1 k (x ) 0 0 1 . ( 16
)
where Jac Υ k is the jacobian matrix of Υ k . In this system of coordinates, the metric tensor

G k (x , x d ) = t Jac Υ k (x , x d ) Jac Υ k (x , x d ) writes: G k : (x , x d ) ∈ Γ k (supp ρ k ) × [0, δ] → Gk (x , x d ) 0 0 (G k ) dd (x , x d ) ( 17
)
where Gk is a C ∞ square matrix of size d -1 and (G k ) dd is a C ∞ positive function.

Let us prove [START_REF] Karlin | A second course in stochastic processes[END_REF]. Let us denote by x = (x 1 , ..., x d-1 ). Since by construction, for all

(x , x d ) ∈ Γ k (supp ρ k ) × [0, δ], f -(Υ k (x , x d )) = x d , one has: ∀j ∈ {1, ..., d -1}, ∇ x j Υ k (x , x d ) • ∇f -(Υ k (x , x d )) = 0. ( 18
)
Moreover, from ( 16) and ( 8), one has for all (x , x d )

∈ Γ k (supp ρ k ) × [0, δ], ∇ x d Υ k (x , x d ) = ∇f -(Υ k (x , x d )) |∇f -(Υ k (x , x d ))| 2 . (19)
Then, from ( 18) and ( 19), it holds

∀j ∈ {1, ..., d -1}, (G k ) j,d = ∇ x j Υ k • ∇ x d Υ k = 0.
This proves [START_REF] Karlin | A second course in stochastic processes[END_REF]. Furthermore, from ( 10) and ( 16), one has:

for all x ∈ Γ k (supp ρ k ), (G k ) dd (x , 0) = 1 ∂ n f (x , 0) 2 . ( 20
)
Finally, a consequence of ( 6) is that d dt f + (γ x (t)) = 0, where γ x satisfies (7) and thus, in the system of coordinates (x , x d ), the functions f + and f write:

f + (x , x d ) = f + (x , 0) and f (x , x d ) = f + (x , 0) -x d , (21)
where with a slight abuse of notation, one denotes

f (Υ k (x , x d )) (resp. f + (Υ k (x , x d ))) by f (x , x d ) (resp. by f + (x , x d )).
3 Potential theory and mean exit time of D

Potential theory

Let us recall the main results from Potential theory which are used in this work. These results can be found for instance in [START_REF] Bovier | Metastability: a potential-theoretic approach[END_REF]. Let us denote by C = B(x 0 , r 0 ) ⊂ D a closed ball centred at x 0 and of radius r 0 > 0 chosen such that B(x 0 , r 0 ) ∩ V δ = ∅ where V δ is given by [START_REF] Evans | Partial differential equations[END_REF]. Let h C,D c be the unique weak solution in

H 1 (R d ) of the elliptic boundary value problem      L ε v = 0 on D \ C v = 0 on D c v = 1 on C,
The function h C,D c is called the equilibrium potential of the capacitor (C, D c ) (as denoted in [4, Section 2]). From elliptic regularity estimates (see for instance [12, Theorem 5, Section 6.3]), the function

h C,D c belongs to C ∞ (D \ C). Therefore, it holds h C,D c ∈ H 1 (D) ∩ C ∞ (D \ C).
Using the Dynkin's formula (see for instance [START_REF] Karlin | A second course in stochastic processes[END_REF]Theorem 11.2]), one has for all x ∈ D,

h C,D c (x) = P x [τ C < τ D c ], (22)
where From [4, Section 2] (see equation (2.27) there), one has the following relation:

τ C = inf{t ≥ 0|X t ∈ C}
∂C E z [τ D c ] e -1 ε f (z) e C,D c (dz) = D e -1 ε f (x) h C,D c (x) dx. ( 23
)
Let us now define, as in [4, Section 2] (see equation (2.13) there), the capacity associated with (C, D c ):

cap C (D c ) = ∂C e -1 ε f (z) e C,D c (dz). ( 24
)
3.2 A first asymptotic estimate on the mean exit time of D

The following results from [8, Corollary 1] and [7, Theorem 2] will be useful in the sequel.

Proposition 2. Let us assume that the assumption [H-D] holds. Let K ⊂ D be a compact set. Then, there exists c > 0 such that it holds in the limit ε → 0 and uniformly with respect to x ∈ K:

λ ε E x [τ D c ] = 1 + O(e -c ε ), and h C,D c (x) ≥ 1 -e -c ε ,
where, we recall, for all

x ∈ D, h C,D c (x) = P x [τ C < τ D c ] (see (22)).
Remark 3. In [8, Corollary 1], the result on

λ ε E x [τ D c
] is not stated with an error term. However, in view of the proof of [8, Corollary 1], the error term is O(e -c ε ) and is uniform with respect to x in a compact subset of D.

Proposition 2 implies that in the limit ε → 0 and uniformly with respect to x ∈ K:

E x [τ D c ] = E x 0 [τ D c ](1 + O(e -c ε )). ( 25
)
We are now in position to obtain a first estimate on the mean exit time of D. Using (25), ( 23) and ( 24), for any compact set K ⊂ D, there exists c > 0 such that in the limit ε → 0:

E x 0 [τ D c ] = D e -1 ε f (x) h C,D c (x)dx cap C (D c ) (1 + O(e -c ε )). Moreover, since h C,D c ≡ 1 on C, h C,D c ≤ 1 on D, f (x) ≥ max C f > f (x 0 ) for all
x ∈ D \ C and using Laplace's method (since x 0 is non degenerate), one obtains that there exits c > 0 such that in the limit ε → 0:

D e -1 ε f (x) h C,D c (x)dx = C e -1 ε f (x) dx + O(e -1 ε (f (x 0 )+c) ) = (2πε) d 2 det Hess f (x 0 ) e -1 ε f (x 0 ) 1 + O(ε) .
Thus, one has the following result.

Lemma 2. Let us assume that the assumption [H-D] is satisfied. Then, in the limit ε → 0:

E x 0 [τ D c ] = (2πε) d 2 det Hess f (x 0 ) cap C (D c ) 1 + O(ε) , (26)
where τ D c is defined by (2) and cap C (D c ) by (24).

To prove Theorem 1, it remains to give an estimate on cap C (D c ) in the limit ε → 0. This is the purpose of the next section.

Proofs of Theorem 1 and Corollary 1

In this section, one obtains sharp lower and upper bounds on the capacity cap C (D c ). From [4, Section 2], one has the following variational principle:

cap C (D c ) = ε D\C ∇h C,D c (x) 2 e -1 ε f (x) dx (27) = inf h∈H C,D c ε D\C ∇h(x) 2 e -1 ε f (x) dx,
where

H C,D c = h ∈ H 1 (R d ), h(x) = 1 for x ∈ C, h(x) = 0 for x ∈ D c }.
Formula (27) holds since the function h C,D c is a minimizer of the functional

h ∈ H C,D c → ε D\C ∇h(x) 2 e -1 ε f (x) dx.
Using this variational principle, one can get a sharp upper bound on cap C (D c ) by choosing a suitable function h ∈ H C,D c .

Upper bound on cap

C (D c )
In this section, one gets a sharp upper bound on cap C (D c ). Let V δ be defined by [START_REF] Evans | Partial differential equations[END_REF] and let h ∈ H C,D c . From 27, one has

cap C (D c ) ≤ ε V δ ∇h(x) 2 e -1 ε f (x) dx + ε D\V δ ∇h(x) 2 e -1 ε f (x) dx. (28)
From ( 13), ( 14), ( 15) and ( 17), one has:

ε V δ ∇h(x) 2 e -1 ε f (x) dx = ε N k=1 x ∈Γ k (supp ρ k ) ρ k (Γ -1 k (x )) × δ 0 t ∇h(x , x d )G k (x , x d ) -1 ∇h(x , x d ) e -1 ε f (x ,x d ) jac Υ k (x , x d ) dx d dx (29)
where t ∇ = (∂ x , ∂ x d ), Υ k is defined by [START_REF] Helffer | Quantitative analysis of metastability in reversible diffusion processes via a Witten complex approach[END_REF], G k is the tensor metric associated with the change of variable x = Υ k (x , x d ) (see [START_REF] Karlin | A second course in stochastic processes[END_REF]) and jac Υ k = √ det G k is the jacobian of Υ k . Let us now consider the following function:

x d ∈ [0, δ] → g(x d ) =
∂ x h(x , x d ) = 0 and ∂ x d h(x , x d ) = d dx d g(x d ).
From ( 17), ( 21), (28), and (29) together with the fact that ∇h = 0 on

D \ V δ , one cap C (D c ) ≤ ε N k=1 x ∈Γ k (supp ρ k ) ρ k (Γ -1 k (x )) e -1 ε f + (x ,0) ε 2 1 -e -δ ε 2 × δ 0 e -x d ε (G k ) dd (x , x d ) -1 jac Υ k (x , x d )dx d dx .

Now let us notice that for any function

ϕ ∈ C ∞ (Γ k (supp ρ k ) × [0, δ], R * + ), one has in the limit ε → 0: δ 0 ϕ(x , x d )e -x d ε dx d = ε ϕ(x , 0) 1 + O(ε) , (30)
uniformly with respect to x ∈ Γ k (supp ρ k ). Thus, applying (30) with ϕ = (G k ) -1 dd jac Υ k , it holds in the limit ε → 0:

cap C (D c ) ≤ N k=1 x ∈Γ k (supp ρ k ) ρ k (Γ -1 k (x ))e -1 ε f + (x ,0) 1 -e -δ ε 2 × (G k ) dd (x , 0) -1 jac Υ k (x , 0)dx 1 + O(ε) .
Finally, using (10), ( 16) and ( 20), it holds in the limit ε → 0:

cap C (D c ) ≤ N k=1 x ∈Γ k (supp ρ k ) ρ k (Γ -1 k (x ))e -1 ε f (x ,0) ∂ n f (x , 0) jac Γ -1 k (x )dx 1 + O(ε) .
Therefore, since from (5) and Lemma 1, f (x , 0) = f (x) for all x = Υ k (x , 0) ∈ ∂D, one has following result.

Lemma 3. Let us assume that the assumption [H-D] is satisfied. Then, it holds in the limit ε → 0:

cap C (D c ) ≤ ∂D ∂ n f (σ) e -1 ε f (σ) dσ 1 + O(ε) , (31)
where, we recall, cap C (D c ) is defined by (24).

Let us now give a sharp lower bound on cap C (D c ).

Lower bound on cap C (D c )

In this section, one gets a sharp lower bound on cap C (D c ). Let V δ be defined by [START_REF] Evans | Partial differential equations[END_REF]. Using ( 27), ( 13), ( 14), ( 15) and ( 17), one has:

cap C (D c ) ≥ ε V δ ∇h C,D c (x) 2 e -1 ε f (x) dx ≥ ε N k=1 x ∈Γ k (supp ρ k ) ρ k (Γ -1 k (x )) δ 0 L k (x , x d ) dx d dx (32) with L k (x , x d ) := ∂ x d h(x , x d ) 2 (G k ) dd (x , x d ) -1 e -1 ε f (x ,x d ) jac Υ k (x , x d ).
Let us define for k ∈ {1, ..., N } and (

x , x d ) ∈ Γ k (supp ρ k ) × [0, δ]: χ k (x , x d ) := (G k ) dd (x , x d ) -1 jac Υ k (x , x d ). ( 33
)
The function χ k satisfies min

Γ k (supp ρ k )×[0,δ] χ k > 0. (34)
Let us consider k ∈ {1, ..., N } and x ∈ Γ k (supp ρ k ). Then, it holds: Then, using (30), ( 33), ( 20), ( 10) and ( 16), one has in the limit ε → 0:

δ 0 L k (x , x d ) dx d = δ 0 ∂ t h C,D c (x , t) 2 χ k (x , t) e t ε dt ≥ inf g∈H 1 (0,δ) g(0)=0 g(δ)=h C,D c (x ,δ) δ 0 d dt g (
where, we recall, cap C (D c ) is defined by (24).

Theorem 1 is then a consequence of (31) and (39) together with (26) and (25). Corollary 1 is a consequence of Theorem 1 and Proposition 2.

) 2 . 3

 23 Metric associated with the change of variable x = Ψ(z, x d ) Let us consider (ρ k ) k∈{1,...,N } ∈ C ∞ (∂D, [0, 1]) N a partition of unity of ∂D: for all y ∈ ∂D, N k=1

  and τ D c is defined by (2). Let us denote by G D be the Green function of L ε associated with homogeneous Dirichlet boundary conditions on ∂D. The equilibrium measure e C,D c associated with (C, D c ) (see [4, Section 2] and more precisely the equation (2.10) there) is defined as the unique measure on ∂C such that h C,D c (x) = ∂C G D (x, y)e C,D c (dy).

1 -e -x d ε 1 -

 11 e -δ ε , which satisfies g(0) = 0 and g(δ) = 1. Let h : V δ → R be such that h • Ψ(z, x d ) := g(x d ), for all (z, x d ) ∈ ∂D × [0, δ].The function h is then extended by 1 in D \ V δ and by 0 outside D. Thus, h belongs toH C,D c since C ⊂ D \ V δ . For all k ∈ {1, ..., N } and for all (x , x d ) ∈ Γ k (supp ρ k ) × [0, δ],denoting with a slight abuse of notation h • Υ k by h, one has h(x , x d ) = g(x d ) and then for any (x , x d ) ∈ Γ k (supp ρ k ) × [0, δ]:

t 0 χ 2 χ 0 ∂ 1 ε 0 χ 2 e - 1 ε 2 N

 02010212 k (x , s) -1 e -s ε ds δ 0 χ k (x , s) -1 e -s ε ds h C,D c (x , δ).The set K = {g ∈ H 1 (0, δ), g(0) = 0 and g(δ) = h C,D c (x , δ) is a closed convex subset of H 1 (0, δ) and the functionalF : θ ∈ H 1 (0, δ) → k (x , t) e t ε dt ε → 0: cap C (D c ) ≥ ε N k=1 x ∈Γ k (supp ρ k ) ρ k (Γ -1 k (x )) e -1 ε f + (x ,0) × δ x d g * x (x d ) 2 χ k (x , x d )e x d dx d dx = ε N k=1 x ∈Γ k (supp ρ k ) ρ k (Γ -1 k (x )) h 2 C,D c (x , δ) e -1 ε f + (x ,0) × δ k (x , x d ) -1 δ 0 χ k (x , s) -1 e -s ε ds x d dx d dx ≥ ε(1 -e -c ε k=1 x ∈Γ k (supp ρ k ) ρ k (Γ -1 k (x )) e -1 ε f + (x ,0) × δ 0 χ k (x , x d ) -1 δ 0 χ k (x , s) -1 e -s ε ds 2 e -1 ε x d dx d dx .

Lemma 4 .

 4 cap C (D c ) ≥ N k=1 x ∈Γ k (supp ρ k ) ρ k (Γ -1 k (x ))e -1 ε f (x ,0) ∂ n f (x , 0) jac Γ -1 k (x )dx 1 + O(ε) .Therefore, since from (5) and Lemma 1, f (x , 0) = f (x) for all x = Υ k (x , 0) ∈ ∂D, one has the following lower bound on cap C (D c ). Let us assume that the assumption [H-D] is satisfied. Then, it holds in the limit ε → 0:cap C (D c ) ≥ ∂D ∂ n f (σ) e -1 ε f (σ) dσ 1 + O(ε) .

Acknowledgements. I am very grateful to T. Lelièvre and D. Le Peutrec for their suggestions and help. This work is supported by the European Research Council under the European Union's Seventh Framework Programme (FP/20072013)/ERC Grant Agreement number 614492.

is continuous and from (34), it is strongly convex. Furthermore, since for all u ∈ K, u(0) = 0, there exists C > 0 such that for all g ∈ K,

Thus, using in addition (34), there exists c > 0 such that for all g ∈ K,

Let us consider a sequence (g n ) n≥0 ∈ K N such that lim n→∞ F (g n ) = inf K F . Then, from (37), (g n ) n≥0 is a bounded sequence in H 1 (0, δ) and thus converges for the weak topology of H 1 (0, δ) towards some g ∈ H 1 (0, δ). Since F is continuous and convex on K, it is a lower semi-continuous function for the weak topology in H 1 (0, δ). Therefore, inf K F ≤ F (g) and since g ∈ K, g is a minimizer of F on K. Finally, because F is strongly convex, g is the unique minimizer of F on K. Let α ∈ R and ϕ ∈ C ∞ c (0, δ). Then, it holds g + αϕ ∈ K and thus

Thus, g ∈ H 1 (0, δ) is a weak solution to the following one dimensional Dirichlet problem on (0, δ):

From (34), one can use the Lax-Milgram Theorem which implies that there exists a unique solution in H 1 (0, δ) of (38). Clearly, this solution is given by

and thus g = g * x . This concludes the proof of (36). Using (32), ( 35) and (36) together with the second statement in Proposition 2, there exists c > 0 such that in the limit