Olivier Lopez

Xavier Milhaud

Individual reserving and nonparametric estimation of claim amounts subject to large reporting delays

Keywords: reserving, reporting delay, truncation, censoring, CART

Thanks to nonparametric estimators coming from machine learning, microlevel reserving has become more and more popular for actuaries. Recent research focused on how to integrate the whole information one can have on claims to predict individual reserves, with varying success due to incomplete observations. Using the CART algorithm, we develop new results that allow us to deal with large reporting delays and partially observed explanatory variables. Statistically speaking, we extend CART to take into account truncation of the data, and introduce plug-in estimators. Our applications are based on real-life insurance portfolios embedding Income Protection and Third-Party Liability guarantees. The full knowledge of the claim lifetime is shown to be crucial to predict the individual reserves efficiently.

Introduction

In non-life insurance, the final cost is rarely known when the claim occurs. In some cases, many years may pass before the amounts are settled. Standard procedures to predict the reserve, like Chain-Ladder (CL) approaches, rely on aggregated cumulative claim amounts. These techniques have the disadvantage that they do not exploit additional information one may have on the claims, and that may help to improve the prediction.

Our paper aims to use this information to better catch the heterogeneity of the data, and 1 get accurate predictions of the individual reserves. In this view, we propose a CARTbased approach [START_REF] Breiman | Classification and Regression Trees[END_REF]) that compensates for the presence of censoring and truncation, naturally present in this kind of data. CART allows us to deal with nonlinearities between the response and the explanatory factors while preserving simplicity of interpretation, and has therefore gained much popularity in actuarial sciences in the last decade (see Wüthrich [2018a], [START_REF] Baudry | A machine learning approach for individual claims reserving in insurance[END_REF] for recent papers related to CART and reserving). Note that other machine learning techniques were recently used for reserving applications, see for instance [START_REF] Duval | Individual loss reserving using a gradient boosting-based approach[END_REF] and Wüthrich [2018b].

This paper is based on two main ideas. First, the time before settlement gives meaningful information that helps to explain the final cost of claims. Second, the reporting delays can strongly affect the prediction of claim amounts. Indeed, right-censoring and left-truncation are known to have a significant impact when there is an underlying time phenomenon (Fleming and Harrington [2011a]). In reserving, most claims handling experts agree to say that the claim lifetime plays a crucial role to explain its final amount.

Typically, claims that take a long time to be settled are more likely to be expensive, as pointed by [START_REF] Maegebier | Valuation and risk assessment of disability insurance using a discrete time trivariate markov renewal reward process[END_REF], [START_REF] Spierdijk | Calculating loss reserves in a multistate model for income insurance[END_REF], [START_REF] Pitt | Modelling the claim duration of income protection insurance policyholders using parametric mixture models[END_REF], or [START_REF] Bluhm | Duration-based policy reserves[END_REF].

Calibrating a model only based on closed claims would thus lead to underestimations of the average claim amounts, since such a procedure would rely on observations with an overrepresentation of short lifetimes. Trying to better understand the relationship between models based on aggregate data (e.g. Chain Ladder) and those using individual data in the context of survival analysis has been done in numerous recent works; such as [START_REF] Bischofberger | Continuous chain-ladder with paid data[END_REF], [START_REF] Hiabu | On the relationship between classical chain ladder and granular reserving[END_REF], and [START_REF] Dolores | Continuous chain ladder: Reformulating and generalizing a classical insurance problem[END_REF]. They discuss the advantages and drawbacks of the most famous reserving techniques, and suggest to link the development of claims to the hazard rate of the cost distribution. This approach bridges the gap between both families of methods, and allow to clarify some assumptions underlying aggregate models. Nevertheless, no explanatory variables were considered in these works, contrary to the paper by [START_REF] Lopez | A tree-based algorithm adapted to microlevel reserving and long development claims[END_REF] who propose a tree-based estimator to predict the residual lifetimes of still open claims (Reported But Not Settled, or RBNS claims). We extend the latter work in three ways. Firstly, as we cannot estimate the final claim amount efficiently using partially observed risk factors (such as the claim lifetime), we introduce plug-in strategies to overcome this issue. Secondly, as large reporting delays lead to biased predictions, we must modify the CART estimation procedure. We thus determine consistent weights given to each observation to cancel (asymptotically) the bias caused by censoring and truncation. Finally, we use bootstrap resampling to estimate the predictive uncertainty of the estimators introduced thereafter.

The remainder of the paper is organized as follows. Section 2 describes the general framework and gives some data processing details. Our main contributions are presented in Section 3, including the theoretical extension that allows us to deal with both partially observed risk factors and large reporting delays. Finally, Section 5 is devoted to real data analyses. We compare our results to two competing approaches used in the insurance industry, namely the Collective Reserving Model (close to Chain Ladder for predicted reserves, see [START_REF] Wahl | The collective reserving model[END_REF]), and the semiparametric Cox model [START_REF] Cox | Regression models and life tables (with discussion)[END_REF]).

General framework and data management

Consider n claims, with amounts (M i ) 1≤i≤n . For some of them, the final amount M i is not observed, since the claim is still open. Indeed, the (random) time before the claim i is fully settled (the so-called claim lifetime), denoted by T i , is censored. Introducing some

censoring variables (C i ) 1≤i≤n , we define (Y i , δ i , N i , X i ) as i.i.d. replications of            Y = min(T, C), δ = 1 T ≤C , N = δM, X = (X (1) , ..., X (d) ) ∈ R d ,
where X are the claim covariates. Note that defining M and N this way amounts to consider one single payment per claim, which is unrealistic in reality. However, the information on payments already made for unsettled claims is most of time unavailable to actuaries, who are usually not in charge of collecting these information.

Remark 2.1. When one has additional information on M (e.g. partial payments already made for RBNS claims), taking it into account would be straightforward in our setting.

There basically exists two possibilities: integrate it as a covariate (N would be the cumulated amount already paid, and we know that M ≥ N ), or consider bivariate censoring.

These variables are easily built from original data. Let's say that d i is the date at which the ith claim occurs and s i the date at which it is fully settled. Moreover, introduce f i the date at which the claim stops being observed. We have T i = s i -d i , and

C i = f i -d i .
In practice f i is often the same date for all the claims, due to data collection. In claim reserving applications, reporting delays must be carefully dealt with. Indeed, the database only reports claims that have been communicated, meaning that "Incurred But Not yet Reported" (the so-called IBNyR) claims are absent. This phenomenon relates to lefttruncation, but is different from the classical left-truncation in the literature. Introducing the reporting delay

τ i = r i -d i
where r i is the reporting date of claim i, only claims such that C i ≥ τ i are observed.

Figure 1 illustrates the case of an unknown claim, and reports the typical information one has on claim history. This truncation phenomenon is not standard, compared to the classical left-truncation model in survival analysis where it is usually assumed that observation only occurs when T ≥ τ (or when Y ≥ τ , see respectively [START_REF] Tsai | A note on the product-limit estimator under right censoring and left truncation[END_REF] and [START_REF] Sánchez Sellero | Uniform representation of product-limit integrals with applications[END_REF]). Here, the situation is slightly different: T < τ means that the claim has been settled before reporting. This can happen for instance if the indemnity has been fixed in advance, and the claim is thus stored in the database. To sum up, the observations are made of (Y i , δ i , N i , τ i , X i ) 1≤i≤n i.i.d. The distribution of this random vector is the conditional distribution of (Y, δ, N, τ, X)|C ≥ τ (indeed, by construction, all observations are such that C i ≥ τ i ).

As already mentioned, both censoring C and truncation τ induce bias in the analysis if they are not taken into account. The former leads to under-represent claims with high amounts, due to the positive correlation between the claim lifetime T and the final claim amount M . The latter impacts claims that occur just before the extraction of the data.

It is thus necessary to introduce a statistical framework in which such phenomenons can be considered, see Section 3. On top of that, our goal is twofold. First, we would like to understand the impact of T i and X i on M i , and typically estimate the final cost of a claim with characteristics X i and lifetime

T i , that is E[M i | T i , X i ]. Second, we aim to predict the final cost of RBNS claims, i.e. E[M i |T i ≥ Y i , δ i = 0, X i ],

Extensions to predict RBNS claims

We introduce in this section the statistical developments for the CART algorithm to be suited to reserving applications.

Managing reporting delays: new "Kaplan-Meier weights"

From a statistical viewpoint, the main difficulty when dealing with censored and truncated data is the fact that the classical empirical means become unadapted. If one wishes to estimate q φ = E[|φ(T, M, X)|] < ∞ (for a given measurable function φ), one may write q φ as an integral, that is

q φ = φ(t, m, x)dF (t, m, x),
where F (t, m, x) = P(T ≤ t, M ≤ m, X ≤ x), and plug a consistent estimator of F . In absence of censoring and truncation, a natural estimator of F is the empirical distribution function, leading to an estimation of q φ through an empirical mean. Now, consider that some observations are either partially observed or totally unobserved, due respectively to censoring and truncation. For simplicity matters, consider that we wish to estimate the marginal cumulative distribution function of some lifetime T , i.e. F T (t) = P(T ≤ t) = 1 -S(t), where S(t) is the survival function. In our case, the fact that truncation occurs when C ≤ τ impacts the way to estimate F T , and makes the estimator differ from the classical one under left-truncation [START_REF] Tsai | A note on the product-limit estimator under right censoring and left truncation[END_REF]). To introduce this new estimator, one needs to state the following assumption.

Assumption 1. τ is independent from (C, M, T, X) and C is independent of (M, T, X).

By analogy to standard techniques like Chain Ladder, this assumption means that we assume homogeneity in terms of time in all dimensions (occurence, development, calendar).

Proposition 1. Provided that Assumption 1 is satisfied, the cumulative distribution function of the lifetime T can be estimated by

FT (t) = 1 - Y i ≤t 1 - δ i 1 τ i <Y i n j=1 1 τ j <Y i ≤Y j
, where we assume that there is no ties among the variables (Y i ) 1≤i≤n .

Proof. For the sake of conciseness, the proof is moved to Appendix B. An alternative formula is also proposed in case of ties.

Note that FT is piecewise constant, and can therefore be written through the following additive formula:

FT (t) = n i=1 w i,n 1 Y i ≤t , with w i,n = δ i 1 τ i <Y i n j=1 1 τ j <Y i ≤Y j Y k <Y i 1 - δ k 1 τ k <Y k n j=1 1 τ j <Y k ≤Y j . (3.1)
Since they take truncation into account, these weights allow to deal with the presence of reporting delays. These new weights slightly differ from the ones of the classical productlimit estimator in presence of right-censoring and left-truncation, see [START_REF] Shulamith | Nonparametric estimation and regression analysis with left-truncated and right-censored data[END_REF].

Practically speaking, w i,n equals 0 when the observation is censored or truncated. The quantity w i,n can be thought as the weight to put on the ith observation to correct the bias caused by censoring and truncation. More weight is given to smallest and largest observations of T . Based on the idea of [START_REF] Sánchez Sellero | Uniform representation of product-limit integrals with applications[END_REF], one could use the same weight when it comes to estimate the joint cumulative distribution function F . Hence, we define

F (t, m, x) = n i=1 w i,n 1 Y i ≤t,N i ≤m,X i ≤x .
A natural estimator of q φ is thus given by

qφ = φ(t, m, x) d F (t, m, x) = n i=1 w i,n φ(Y i , N i , X i ). (3.2)
That being said, we can now introduce the weighting procedure used within the CART algorithm to deal with truncation and censoring.

Weighted regression-tree procedure

To extend the CART algorithm to predict the individual final claim amount M , we adopt the same framework as in [START_REF] Lopez | A tree-based algorithm adapted to microlevel reserving and long development claims[END_REF]. However, the weighting scheme differs because of the truncation due to reporting delays. Moreover, the use itself of the algorithm (denoted further wCART) for prediction purpose has to be adapted since T , considered as an explanatory variable of M thereafter, is partially observed for RBNS claims. For the paper to be self-contained, the complete version of wCART is given in Appendix C. Here, for brevity, we only remind to the reader the main principles.

Suppose that we want to estimate π(z

) = E[φ(M ) | Z = z],
where Z = (T, X). At each step of the algorithm, one determines a rule z = (t, x (1) , ..., x (d) ) → R j (z) to split the data and create two partitions of the covariate space (leading to the binary tree structure). That is, for each value of z, R j (z) equals either 0 or 1 depending on whether some conditions are satisfied by z. Of course R j (z)R j (z) equals 0 when j = j , and j R j (z) equals 1 to ensure that created subsets are exhaustive and disjoint. To select the best rule, we use a weighted quadratic loss so as to minimize within-node variances of the response φ(M ) in created subsets (instead of a quadratic loss in the classical CART algorithm, well-suited to the estimation of an expectation with fully observed variables).

Not surprisingly, the weights w i,n given by (3.1) are considered to compensate for censoring and truncation. Finally, each set of rules R = (R 1 , ..., R K ) is associated with a tree-based estimator of the regression function, that is

πR (Z) = K j=1 πj R j (Z), where πj = n i=1 w i,n φ(N i ) R j (Z i ) n i=1 w i,n R j (Z i )
.

(3.3) Remark 3.1. Without any stopping rule, this procedure ends up with a complex tree which is very likely to overfit the data. This estimator is thus not satisfactory to estimate π. To get an estimator with lower dimension, a pruning step is required. Let K(R) denote the number of leaves of a subtree: the pruning approach consists of minimizing the following penalized loss,

n i=1 w i,n (φ(N i ) -πR (Z i )) 2 + α K(R) n ,
where α > 0 is a tuning parameter, usually chosen through cross-validation.

These building and pruning strategies lead to consistent estimators which are asymptotically unbiased, thanks to the properties of Kaplan-Meier estimators (see Lopez et al.

[2016] and references therein). Section 4 illustrates such consistency via simulations.

Nonetheless, for reserving applications, this estimator cannot be applied directly. Indeed, the explanatory vector Z is not fully observed, because T can be censored or truncated.

We therefore propose some possibilities to overcome this issue.

Strategies to compute the reserves of RBNS claims

Take φ(M ) = M . From the selected tree estimator, it is possible to deduce a predictor of M for claims where Z is fully observed. However, in the context of RBNS claims, T is censored. It means that Z is partially observed. Denote by y the observed duration.

To provide a reasonable estimator of M when δ = 0, one thus needs to take into account that T ≥ y. Consider the two following approaches.

(A) Bayes -the best prediction given the available data is

E[M | Y = y, δ = 0, X = x],
which can also be written as

M = E[M | T ≥ y, X = x] = E[M 1 T ≥y | X = x] E[1 T ≥y | X = x] . (3.4)
Two trees are built using wCART, to estimate the numerator and denominator.

(B) Plug-in -build one tree π by wCART to estimate π(t,

x) = E[M | T = t, X = x].
Then fit a model for T | T ≥ y, X = x, from which a prediction T (y, x) can be computed. Finally, predict M using the plug-in principle: M = π( T (y, x), x).

In the latter approach, several prediction models T can be used, among which any machine learning prediction model adapted to censoring and truncation. In the sequel, we consider the three following cases:

(B1) : build two different wCART estimators to estimate r 1,y (x) = E[T 1 T ≥y |X = x] and r 2,y (x) = E[1 T ≥y |X = x],
and compute T (y, x) = r1,y (x)r 2,y (x) -1 ;

(B2) : use a simplified prediction of T , assuming that T does not depend on X:

T (y, x) = T (y) = t1 t>y d FT (t) 1 t>y d FT (t) = n i=1 w i,n Y i 1 Y i ≥y n j=1 w j,n 1 Y j ≥y .
Although inconsistent in full generality, the simplicity of this approximation is expected to favor good behavior of this strategy.

(B3) : get T (y, x) from the semiparametric Cox model, see [START_REF] Cox | Regression models and life tables (with discussion)[END_REF].

For comparison purposes, we consider two other naive competitors, where the input information about T is not realistic:

(B4) : predict M by M = π(y, x), where y is the observed censored duration;

(B5) : predict M by M = π(r(x), x), where r(x) is an estimator of r(x) = E[T | X = x]
obtained by wCART.

Given M , the individual reserves as well as the global reserve can easily be obtained.

In our setting, this amounts to keep M i as the estimator of the ith individual reserve, since we do not have information on partial payments already made. Strategies (B4) and

(B5) are expected to underestimate the individual reserves, since they do not consider the full information about T (i.e. T > y,

knowing that E[T | X] ≤ k + E[T -k | T ≥ k, X]
for all k ≥ 0). On the contrary, strategies (A), (B1), (B2) and (B3) take it into account, and are thus expected to provide well-estimated reserves on average.

8

First we investigate the practical behaviour of our tree-based estimator with truncated and censored data, where truncation is defined as in Section 2. For the sake of simplicity, we consider the case where we are interested in the distribution of the lifetime T , thus

focusing on estimating π 0 (x) = E[T | X = x].
Consider the following simulation scheme:

1. draw n i.i.d. replications (X 1 , . . . , X n ) of the covariate, with X i ∼ U(0, 1); 2. draw n i.i.d. lifetimes (T 1 , . . . , T n ) following an exponential distribution such that

T i ∼ E(β = α 1 1 1 X i ∈[a,b[ + α 2 1 1 X i ∈[b,c[ + α 3 1 1 X i ∈[c,d[ + α 4 1 1 X i ∈[d,e] ).
3. draw n i.i.d. censoring times (Pareto-distributed:

C i ∼ Pareto(λ, µ)) and n i.i.d.
truncation times, uniformly distributed (τ i ∼ U(γ, δ));

4. from the simulated lifetimes, censoring and truncation times, get for all i the actual observed lifetime Y i and the indicator δ i = 1 T i ≤C i ;

5. compute the weights w i,n from the entire generated sample

(Y i , C i , τ i , δ i ) 1≤i≤n .
Parameter values are stored in Table 1, and descriptive statistics corresponding to the various simulated datasets are available in Table 4 of Appendix D. To each simulated sample, we fit a regression tree with the algorithm of Section 3.2. Then, we compute the weighted squared errors given by W SE i = w i,n (π l(i) -π 0 (X i )) 2 , where πl(i) is deduced from (3.3) for the ith observation that belongs to leaf l(i), and where we know that π 0 (X i ) = 1/β. To gain some robustness, we repeated 100 times the simulation scheme above to compute empirical means of W SE i , leading to the M W SE. We also considered different values for (λ, µ) and (γ, δ) in the censoring and truncation processes so as to measure the impact of both phenomenons on the procedure's performance (see Table 1).

Figure 2 and Table 2 report the results. While the sample size remains low, the censoring and truncation phenomenons seem to significantly impact the performance. This impact seems to be complex when both truncation and censoring are present. However, as soon as there are enough observations, the weighted CART estimator with weights w i,n provides satisfactory results.

Group-specific means Component probabilities

Censorship rate Truncation rate

α 1 α 2 α 3 α 4 [a, b[ [b, c[ [c, d[ [d,

Applications

To be in line with regulatory constraints and current practices in insurance companies, the assessment of reserves is considered on a quarterly basis (note that we changed this time step to check the stability of our conclusions, which was confirmed). Building the database and running the estimations is thus performed every three months, with updated policyholders' features and claim characteristics. Each time, we split our n-sized data into two independent subsamples to assess the prediction power of the methods: the learning set with n l observations (representing two third of the dataset), and the test set containing the remaining n t observations (from which the predictive power can be assessed). In practice, we backtest our reserve predictions. To do so, we only consider the closed claims (among the n t observations) at the last observed date t 1 . This way, the true outcomes M and the true global reserve P are available and can be compared to the predictions.

10

The predictive uncertainty of our estimators M and P , where P stands for the global estimated reserve, is approximated using bootstrap resampling (see [START_REF] Björkwall | Non-parametric and parametric bootstrap techniques for age-to-age development factor methods in stochastic claims reserving[END_REF],

Pinheiro et al. [2003], [START_REF] England | Analytic and bootstrap estimates of prediction errors in claims reserving[END_REF]). Knowing that P = n t i=1 (1 -δ i ) M i , we will consider:

• the overall standardized error, defined by = ( P -P )/P ,

• and the root mean squared deviation RM SD, defined by

RM SD = √ M SE = 1/ n t i=1 (1 -δ i ) n t i=1 (1 -δ i )( M i -M i ) 2 .
The first quantity indicates the overall quality of the reserving method, whereas the second one gives insights about the accuracy of the prediction of individual reserves. We will also report the bootstrap estimate of the relative standard deviation (RSD), defined by

RSD = 1/ P 1 B -1 B b=1 ( Pb -P ) 2 , where P = (1/B) B
b=1 Pb , with Pb the estimated reserve on the bth bootstrap sample. In the sequel, the number of bootstrap samples is set to B = 1000. For each of the above indicators, the superscript (.) will be used to precise which strategy was used to make the predictions.

Income Protection

Short-term disability insurance was designed to protect the policyholders against the loss of some revenue. In our context, the coverage can last up to three years, meaning that the duration of payments T is capped. Here, predicting the final claim amount is similar to predicting the residual lifetime in the disability state. Without loss of generality, say that the insurer has to pay 1A C for each insured day (i.e. M = T ). We wish to predict the global reserve at various settlement dates, by summing predicted remaining claim lifetimes using strategy (A) (with M replaced by T , see Section 3.3 and Remark C.1 in Appendix C).

Other strategies presented in Section 3.3 are useless, since T is the response in this case.

We compare our results to two famous approaches: i) the Cox model (recommended by regulators for such insurance risk), and ii) the Chain Ladder (CL) model applied on RBNS claims only (since there is no reporting delays here, there are no IBNyR claims).

About the database

Our database reports 65 670 claims related to income protection guarantees over six years, from 01/01/2006 to 12/31/2011. For each claim, we know the gender of the policyholder (14 455 males, 51 215 females), her socio-professional category (2 406 managers, 62 799 employees and 465 others), her age when the claim occured, the duration in the disability state, the commercial network (three kinds of brokers: 28 662 "Net-A", 4 890 "Net-B " and 32 118 "Net-C "), and the cause (57 131 sicknesses, and 8 539 accidents) that triggered the coverage. The censoring rate equals 7.2% at the end of the observation period. The mean observed duration in the disability state is 100 days (beyond a deductible of 30 days), with a median of 42 days and a standard deviation of 162 days. If necessary, more details about the data can be found in [START_REF] Lopez | A tree-based algorithm adapted to microlevel reserving and long development claims[END_REF], Section 3.

Evolution of the predicted reserve

We predict the global reserve P every quarter, from 01/01/2008 to 10/01/2009. The overall results are stored in Table 3, with additional details about the datasets given in the top part of Table 8 in Appendix E.4. As an illustration, Figure 3 shows the evolution of the prediction error for each strategy.

First, notice that P is strongly underestimated in 2008, with both (CL) and (A). On 01/01/2008, the error (CL) reaches 60%, whereas it roughly equals 44% using (A). Clearly, the (CL) model does not take into account the censoring phenomenon adequately, which causes large underestimations of the final claim amounts. Although (A) is supposed to appropriately deal with censoring, largest observed lifetimes equal two years at the beginning of 2008. Given that asymptotic properties of our tree estimator are guaranteed once the observations almost entirely map the domain of possible values for T (which is not the case here, since some of the claims will last up to three years in practice), it is not surprising to underestimate the individual reserves M (and thus P ). On the contrary, the Cox model [START_REF] Cox | Regression models and life tables (with discussion)[END_REF]) seems to provide satisfactory results for all the settlement dates. Its semiparametric specification allows us to anticipate longer lifetimes since the beginning, thanks to the baseline hazard component. Still, these good results should be moderated since they originate from favourable circumstances: low censoring rate (see Table 8), no indication that the proportional hazards (PH) assumption may not be reasonable (looking at scaled Schoenfeld residuals), and bounded lifetimes that cover a narrow interval ([0, 12] quarters), ensuring that the response T has low variance.

Second, the errors of (CL) and (A) decrease as time passes. This is in line with expec- 8). The indicators RSD (A) and RSD (Cox) show that the variance of the prediction of P are similar for both strategies, whereas there is a significant difference between them when looking at the RM SD indicator. In particular, the Cox model seems to provide more accurate individual predictions here.

Third, for most recent settlement dates, (A) and Cox tend to overestimate P whereas (CL) still underestimates P ( (CL) = -24%). These statements make sense, since the portfolio is observed until 12/31/2011. Getting closer to this date, the percentage of fully observed claims increases and the actual global reserve decreases. This explains why the (CL) reserve gets closer to the reality. However, such improvements would not be experienced in practice, since P would have no reason to decrease so much (except in case of run-off business). Concerning the results for strategy (A), our backtesting approach implies that considered claims for the prediction of the reserve P are settled on 12/31/2011. The estimations thus face a selection bias, due to an overrepresentation of claims with short developments (whatever X). Since our estimator is based on past information, it anticipates longer developments on average for those claims. Hopefully, this is therefore a non-issue, and the same applies to the Cox model. To make sure about that, we ran the estimations on the whole database (without only selecting closed claims on 12/31/2011), and confirmed that the prediction error remained low and stable.

To sum up, the predictions by strategy (A) and the Cox model give similar results on this example, provided that the main underlying characteristics of the insured risk have been observed (nearly three years of historical information here). In addition, although intermediary results (regression coefficients for Cox, and tree estimators for (A)) were not presented here for conciseness, both models accord with designating the policyholder's age as the most discriminant risk factor to explain T . This is good news since this is in line with what risk experts do observe in short-term disability.

Third Party Liability (TPL) insurance

In this application, we aim to estimate individual reserves M based on T and X, where

T is censored and where claims are subject to large reporting delays. We compare the methodologies (A) and (B1-5) proposed in Section 3.3 to the Collective Reserving Model (CRM), recently introduced by [START_REF] Wahl | The collective reserving model[END_REF]. The CRM enables to split the reserves dedicated to RBNS and IBNyR claims, allowing for a more fair comparison of RBNS reserves than when using Chain Ladder (CL). We consider the ausautoBI8999 dataset providing claims in motor insurance1 . A lot of claims have long development times, causing specific claim management processes and atypical loss development triangles (an example of such triangle is provided in Appendix E.2).

Brief description of the database under study

The dataset is made of 22 036 settled automobile bodily injury claims in Australia. These Knowing that the mean claim lifetime equals 558 days and that accident dates start in mid-1989, it looks necessary to omit the data before 01/01/1993 in our study to avoid issues related to data quality. Finally, our database reports 16 822 claims.

Prediction of individual reserves for RBNS claims

First, claim amounts have to be inflated to the most recent date (01/01/1999). The annual inflation rate, estimated by the method of [START_REF] Lopez | A censored copula model for micro-level claim reserving[END_REF] (Appendix A), equals 0.39%.

We wish to give quarterly predictions of the cost of RBNS claims, i.e. E[M | T > y, X],

between 09/30/1996 and 06/30/1997. We use the strategies (A) and (B1-5) to do so; which obviously lead to different estimators of the individual reserves. More precisely, strategies of type (B) give the same tree estimator π(t, x), but reserve predictions will differ for one simple reason: rebuilding the information on T is not considered analogously. Anyway, for all the settlement dates studied, looking at the tree π indicates that the claim lifetime has been detected as the risk factor with the strongest impact on the settled claim amounts.

For instance, on 03/31/1997, Figure 4 shows that the most discriminant threshold for claim duration is 838 days, meaning that claims that last more than this threshold before being closed are expected to cause significantly higher final claim amounts (on average 75,000$, as compared to 19,000$ otherwise). Note that this threshold is somewhat stable: it varies by less than 5% depending on the settlement date in practice. Finally, in this example, the population was divided using seven segmentation rules; only based on the claim lifetime, the reporting delay, and the legal representation. Other characteristics such as the number of injured people were not selected, justifying our initial beliefs.

Concerning the reserve predictions, they are in line with our expectations. Figures 5 and6 give the whole picture of the results. For detailed numerical results and additional information on the datasets, please refer respectively to Table 7 in Appendix E.3 and bottom part of Table 8 in Appendix E.4. The best prediction method seems to be (B1), since the prediction error remains low and stable as compared to others. The error (B1) is about ±7% (except for two settlement dates, but more on this later). Unless this may seem substantial, the explanation is twofold : learning samples are of limited size, and the censoring rate is high (between 30% and 55%, depending on the settlement date).

Globally, our reserve estimates are more accurate than insurance practice (CRM model).

Indeed, the CRM (closely related to Chain Ladder) systematically underestimates the reserve, which is not surprising in case of claims with long developments. Improvements related to the latest settlement dates are once again fictive, as explained in Section 5.1.2.

The results from other strategies can be further analyzed. Recall that strategy (B4)

considers
CRM (A) (B1) (B2) (B3) (B4) (B5)
Figure 5: Evolution of the prediction errors for each strategy, for all the settlement dates. not integrate the information on elapsed time. Given that T

= E[T | X] underestimates y + E[T -y | T ≥ y, X],
these results were clearly expected. This strategy is even worse than (B4) since it systematically makes the same mistake and does not really benefit from newer information coming from latest experience (at most recent settlement dates).

Strategy (B3), estimating T (given that T > Y ) by the semiparametric Cox model before plugging it into the model that predicts M , reveals very unstable results. As can be seen in Figure 5, the associated prediction error (B3) sometimes explode (see for instance the predicted global reserve on 09/30/96). In fact, a deeper analysis on intermediate results about Cox modelling shows that the crucial underlying PH assumption is strongly violated. In this context, it is very unlikely that predictions on T can be trusted, meaning that our final estimator can not rely on such predictions. Now focusing on strategies (A), (B1) and (B2), the corresponding prediction errors seem to behave the same. In terms of computation power, (B1) is more demanding than (A). However, it is easier to estimate

E[T |T > y, X = x] than E[M |T > y, X = x],
since T is a much lower dispersed random variable than M . As expected, those methodologies tend to overestimate the global reserve because they suffer from the lack of data related to high values of y, which makes the denominator of (3.4) tend to zero. This effect is smoothed when using (B1), thanks to the plug-in step. In addition to being much simpler, strategy (B2) shows good performance, which makes it attractive. However, because it does not integrate the information on X, it should not be recommended if the portfolio composition (in terms of the distribution of X) is subject to significant changes. For latest settlement dates, all strategies taking into account that T ≥ y tend to overestimate the reality, because of the selection bias due to backtests (see also the discussion in Section 5.1.2). The bump of the prediction errors for 1997-09-30 1997-12-31 1998-03-31 1998-06-30 1996-09-30 1996-12-31 1997-03-31 1997-06-30 all these strategies on 12/31/97 is caused by significantly different distributions of claim durations in the learning and test sets. To mention an example of such difference, the third quartiles respectively equal 760 and 730 days. This leads to overestimations when predicting T (given that T ≥ y) on the test set, hence on M .

(A) (B1) (B2) (B3) (B4) (B5) (A) (B1) (B2) (B3) (B4) (B5) (A) (B1) (B2) (B3) (B4) (B5) (A) (B1) (B2) (B3) (B4) (B5)
Concerning the prediction uncertainty of the strategies considered, Figure 6 highlights that one should have little confidence in predictions coming from (B3). Strategies (A) and

(B4) globally seem to provide predicted reserves with little variance, which makes sense since they are one-step prediction methods (without plug-in). Finally, (B1) and (B2) lead to similar prediction uncertainties, higher than when using other strategies (except (B3)).

In a nutshell, the strategy (B1) seems to outperform all other methods in terms of reserve prediction quality, and shows stable results on our data. This estimator is supposed to be asymptotically unbiased, with acceptable variance. We do not pretend it to be the best choice whatever the case, but proved through situations implying different sample sizes and censoring rates that it remains interesting.

Conclusion

We proposed different methodologies to perform individual claim reserving, based on regression trees. Our contribution is twofold: our modelling enables to take into account the reporting delays appropriately, and shows the good way to implement plug-in estimators to get reasonable estimations of the reserves. The main features of these approaches are the possibility to use all available information on a claim to predict its final state. In other words, the information on the time since occurrence of the claim is appropriately and fully integrated in the model in our framework. To go further, this work could be improved in several ways. Our applications are mainly a picture of the reserve at some point of time. In particular, no dynamic readjustment of the reserve -due to new information or events that affect the claim -is considered. Nevertheless, our technique may be easily modified to incorporate this, as long as the required information is available.

Among other possible improvements, let us mention the possibility to use random forests (i.e. aggregations of regression trees) to stabilize the results, since the CART algorithm is known to be sometimes sensitive to the introduction of new data. The drawback would be a loss of intelligibility of the obtained model. Moreover, some assumptions could be relaxed. In particular, a key assumption in our work is the independence between τ and the other variables of the model. It is possible to easily relax this assumption by making τ depend on the covariates X, and then compute a stratified version of the estimator (where the computation of the weights is done separately on groups of observations belonging to the same stratum, as in [START_REF] Galimberti | A weighted kaplan-meier estimator for matched data with application to the comparison of chemotherapy and bone-marrow transplant in leukaemia[END_REF]).

Let M i denote the observed claim amount (before removing inflation). In general, the data reports M i instead of M i . We assume that log

M i = βd i + log M i , (A.1)
where d i is the fiscal year at which the ith claim is observed, and β is an inflation factor that is going to be estimated using our data. We assume that (M i ) 1≤i≤n are i.i.d. and independent of (d i ) 1≤i≤n . The dates d i take their value in {0, ..., D} (where D + 1 periods, e.g. years, are observed). Then, proceed as follows:

• Compute m i,j , defined as the average of the fully observed claims that occurred on the fiscal year d i , and which are settled after j years (for (i, j) such that

d i + j ≤ t 1 ,
where t 1 is the last observed date). Let n i,j denote the number of such claims.

• Under (A.1), we know that log m i,j ≈ β j d i + α j , where

α j = E[log M i |T i = j].
For each j, we compute βj the weighted least-square estimator of the slope β based on the points (m i,j , d i ) i:d i +j≤t 1 . More precisely, one solves (α j , βj ) = arg min

α j ,β j i:d i +j≤t 1 n i,j (log m i,j -α j -β j d i ) 2 .
• Finally, our estimator of β is given by β = j n 1/2 j βj j n

1/2 j , with n j = i n i,j .

For each claim i such that δ i = 1 (settled claims), we thus consider Mi = M i e -βd i as an estimator of the final claim amount M i , once removed the inflation effect. In practice, one will use this amount in applications.

B Estimator of S(t) in our censoring-truncation model

Let us recall that, for a discrete variable A taking value at point {a 1 , ..., a k }, its survival function S A (t) = P(A ≥ a) can be written as

S A (t) = k j=1 (1 -λ A (a j )) , (B.1) with λ A (t) = - dS A (t) S A (t) .
sistent estimator obtained from the data. If the variable A is not discrete, it can still be approximated by a discrete distribution where the (a j ) 1≤j≤k are replaced by the value of the complete observations (in our case, the uncensored observations). This is the basis of the construction of Kaplan-Meier and other product-limit based estimator, see for example Fleming and Harrington [2011b].

Hence, our aim is to determine a consistent estimator of λ T (t) = -dS T (t) S T (t) .

Let L(t) = P(τ ≤ t), S C (t) = P(C ≥ t), α = P (τ < C), and

S 1 (t) = E [δ1 τ ≤Y 1 Y ≤t |τ < C] .
By basic computations, we obtain dS

1 (t) = -α -1 S C (t)L(t)dS T (t). Indeed, S 1 (t) = α -1 E 1 τ ≤t 1 T ≤T E[1 max(τ,T )<C |τ, T ] = α -1 E [1 τ ≤T 1 T ≤t S C (max(τ, T ))] = α -1 E [1 τ ≤T 1 T ≤t S C (T )] = -α -1 t 0 S C (y)L(y)dS T (y),
where we used that (T, τ ) is independent of C for the second line, and the independence between T and C for the last line. On the other hand, let

S 2 (t) = E [1 τ <t<Y |τ < C] . We have S 2 (t) = α -1 S T (t)S C (t)L(t). Since S 2 (t) = α -1 E 1 τ <t 1 t<T 1 max(t,τ )<C = α -1 E 1 τ <t E 1 t<T 1 max(t,τ )<C |τ = α -1 E [1 τ <t S T (t)S C (max(t, τ ))] = α -1 E [1 τ <t S T (t)S C (t)] , hence - dS T (t) S T (t) = dS 1 (t) S 2 (t) .
The quantities S 1 and S 2 can be estimated consistently by

Ŝ1 (t) = 1 n n i=1 δ i 1 τ i <Y i 1 Y i ≤t and Ŝ2 (t) = 1 n n i=1 1 τ i <t≤Y i .

E.3 Detailed results on the calibration of reserves

We provide here the detailed results for the application on TPL guarantees. When using the CRM model, we only report the part of the reserve dedicated to RBNS claims.

  which would enable us to determine appropriate individual reserves for still open claims. Remark 2.2. As in other reserving techniques, one has to remove the inflation affecting historical information on M to ensure that the amounts are i.i.d. Various methodologies can be used, among which using external data when available. Otherwise, one could follow the procedure by Lopez [2019] (see Appendix A for corresponding technical details).

Figure 1 :

 1 Figure 1: Case of an unknown claim. The observation period begins at t 0 and stops at t 1 .

Figure 2 :

 2 Figure 2: MWSE as a function of the logarithm of the sample size (linear interpolation for n = 100, 1000, 5000).

Figure 3 :

 3 Figure 3: Evolution of the prediction errors (in %) depending on the methodology; for all the settlement dates under consideration.

  claims arose from accidents occurring from July 1989 to January 1999. The database contains event dates (accident, reporting, closing), operational time (indicator of claim management difficulties), type of injury, number of injured people, potential legal representation of the policyholder, and aggregated settled claim amount. In Appendix E.1, Table 6 summarizes descriptive statistics on these variables, as well as other created by-hand variables useful for our study (e.g. reporting delay, or claim duration). Additionally, Fig-ure7shows that reporting started very late as compared to accident dates, which is very likely due to data collection. In particular, there are no settled claims before 01/01/1993.

Figure 4 :

 4 Figure 4: Optimal tree π following (B)-type strategies, on 03/31/1997. The lifetime ('En-dObsW') appears as the most important explanatory variable to predict claim amounts.

Figure 6 :

 6 Figure 6: Boxplot of predictions of P at all settlement dates (B=1000 bootstrap samples).

Table 1 :

 1 Parameters involved in the simulation scheme.

	e]	10%	30%	50%	5%	15%

Table 2 :

 2 Mean weighted squared errors w.r.t. the censoring rate and sample size, in the case where the truncation rate equals 15%.

	Truncation	Censoring	Sample		Group-specific MWSE		MWSE
	rate (%)	rate (%)	size (n)	Group 1	Group 2	Group 3 Group 4	Global
			100	2.4480	20.8320000 1.61175 2.12600	2.23375
		10%	1000	0.4705	11.7646667 0.03925 0.00975	1.10795
			5000	0.0034	0.6306333 0.00120 0.00020	0.05720
			100	5.9335000 53.97467	2.1965	2.32850	5.24670
	15%	30%	1000	4.4461667 15.79833	0.0375	0.00225	1.84215
			5000	0.1077333	1.21450	0.0006	0.00005	0.11877
			100	5.639833 13.9091667 2.64475 4.27975	2.03135
		50%	1000	2.069667 10.3565000 0.13850 0.15400	1.12775
			5000	0.098900 0.9647667 0.00100 0.00000	0.09480

of the MWSE in function of the number of observations

	Censoring 50% -Truncation 5%
	Censoring 30% -Truncation 5%
	Censoring 10% -Truncation 5%
	Censoring 50% -Truncation 15%
	Censoring 30% -Truncation 15%
	Censoring 10% -Truncation 15%

Table 3 :

 3 Actual reserve P , prediction P , and corresponding predictive uncertainties. tations, as more and more information become available. One year later (on 01/01/2009),

		01/01/08 04/01/08 07/01/08 10/01/08 01/01/09 04/01/09 07/01/09 10/01/09
	P	378 817	363 899	384 703	382 289	400 365	391 806	380 500	342 298
	P (CL)	151 017	166 614	193 593	207 677	243 701	242 688	254 947	259 834
	(CL)	-60.1%	-54.2%	-50%	-45%	-39%	-38%	-33%	-24%
	P (Cox)	406 559	359 710	386 701	366 381	414 068	388 272	389 268	378 820
	(Cox)	7.3%	-1.2%	0.5%	-4.2%	3.4%	-0.9%	2.3%	9.5%
	RSD (Cox)	2.1%	5.6%	1.9%	1.8%	1.6%	2.4%	0.9%	0.8%
	RM SD (Cox)	243	234	241	230	224	210	201	187
	P (A)	211 357	227 088	263 030	312 400	402 398	384 361	387 525	374 133
	(A)	-44.2%	-42%	-31.6%	-18.3%	0.5%	-1.9%	1.8%	9.3%
	RSD (A)	3.2%	3%	3.7%	2.2%	2.3%	2.1%	0.9%	1%
	RM SD (A)	417	430	444	405	383	376	377	371

(A)

decreases steeply to reach 0.5%, whereas it still equals 39% with (CL). The prediction of P thus improves much faster with (A) than with (CL). Strategy (A) benefits from the fact that larger observed lifetimes are now included in the learning set, allowing to significantly improve the quality of our estimator. On the contrary, (CL) still suffers from the proportion of censored lifetimes (almost 11.9% on 01/01/2009, see Table

Table 6 :

 6 Example of loss triangle.

		jan-99						jan-99						jan-99
								jan-98
		jan-98												jan-98
								jan-97
		jan-97						jan-96						jan-97
	Closing date	jan-96					Reporting date	jan-95					Closing date	jan-96
								jan-94
		jan-95						jan-93						jan-95
								jan-92
		jan-94												jan-94
								jan-91
			jan-91	jan-93	jan-95	jan-97	jan-99		jan-91	jan-93	jan-95	jan-97	jan-99	jan-92	jan-94	jan-96	jan-98
				Accident date					Accident date				Reporting date
	Figure 7: Dates of events in ausautoBI8999 dataset. Reporting only started mid-1993.
		d1	d2		d3	d4	d5	d6	d7		d8	d9	d10 d11 d12 d13 d14	d15 d16
	01/01/93	0.1	27		57	237	343 199 184 524 563	381	428	250	449	651	606	19
	04/01/93	0	36		171	162 909 332 324 454 208	389	192	862	344	275	119
	07/01/93	5	61		228	395 210 370 353 501 805	188	330	235	849	358
	10/01/93	2	49		139	169	98	571 327 647 655	687	829	342	88
	01/01/94	5	112	101	210	209	299 223 729 264	855	370	120
	04/01/94	2	15		236	372 240 371 792 136 652	604	176
	07/01/94	8	152	688	350	437	364 580 327 507	228
	10/01/94	19	125 304 347	157 439 452 300 143
	01/01/95	61	102 254 166	154 414 492 686
	04/01/95	11	82		129 225 310 477 148
	07/01/95	15	79		67	245 361 719
	10/01/95	0.9	43		234 228 139
	01/01/96	0	112	263	42
	04/01/96	5	91		53
	07/01/96	7	16
	10/01/96	0

E.2 Example of loss triangle in the TPL insurance application

We give the (non cumulated) payments for RBNS claims on 12/31/1996 (settlement date), leading to a CL reserve of 26 165 560$. Numbers are given in thousands and rounded to make the reading easier.

Table 7 :

 7 Quarterly assessment of the global reserve at various settlement dates.

		09/30/96	12/31/96	03/31/97	06/30/97	09/30/97	12/31/97	03/31/98	06/30/98
	P	100 872 881 96 545 407 97 351 400 87 497 571 78 389 452 71 489 554 58 207 140 56 493
	P CRM	7 212 119	8 014 483	9 015 614	9 758 563	11 935 499 14 650 041 15 397 966 20 617
	(CRM )	-92.8%	-91.7%	-90.7%	-88.8%	-84.7%	-79.5%	-73.4%	-63.5%
	P (A)	115 776 538 94 284 076 103 701 446 108 631 498 96 881 249 114 245 317 92 722 102 71 456
	(A)	+14.8%	-2.3%	+6.5%	+24.2%	+23.6%	+59.8%	+59.3%	+26.5%
	RSD (A)	14.5%	5.3%	6.7%	6.1%	5.1%	11.8%	4.1%	3.8%
	RM SD (A)	131 477	118 117	103 363	74 249	61 914	82 596	67 906	149 120
	P (B1)	103 890 209 91 167 753 88 837 467 89 498 892 99 104 425 94 607 399 65 703 004 60 192
	(B1)	+3%	-5.6%	-8.7%	+2.2%	+26.4%	+32.3%	+12.9%	+6.5%
	RSD (B1)	43.5%	11.7%	16.7%	29%	16.5%	19.1%	26.5%	18.9%
	RM SD (B1)	142 831	119 204	106 326	83 161	65 001	95 649	76 673	155 189
	P (B2)	117 001 710 111 788 049 123 314 660 98 413 820 113 376 802 125 713 715 81 641 459 66 327
	(B2)	+16%	+15.8%	+26.7%	+12.5%	+44.6%	+75.8%	+40.3%	+17.4%
	RSD (B2)	33%	8.1%	17.7%	32.5%	14.3%	15.4%	40.2%	31.9%
	RM SD (B2)	140 472	119 905	107 727	84 696	67 241	102 941	95 601	158 678
	P (B3)	210 330 566 102 857 540 113 318 807 98 456 074 102 297 911 268 349 985 92 537 067 71 268
	(B3)	+108%	+6.5%	+16.4%	+12.5%	+30.5%	+275%	+59%	+26.2%
	RSD (B3)	70%	6.1%	20.3%	39.4%	12.6%	42.6%	58.7%	28.7%
	RM SD (B3)	176 961	119 356	107 170	84 024	67 402	200 421	108 349	159 456
	P (B4)	66 954 771 52 828 868 58 790 353 63 305 968 58 620 699 58 700 374 55 884 349 41 978
	(B4)	-33.6%	-45.3%	-39.6%	-27.6%	-25.2%	-17.9%	-4%	-25.7%
	RSD (B4)	9.5%	3.2%	4.6%	5.1%	2.8%	6.3%	20.5%	2.5%
	RM SD (B4)	132 730	120 251	106 166	78 661	62 152	81 923	93 140	154 126
	P (B5)	56 792 306 41 315 697 43 034 713 53 340 431 39 257 861 35 175 163 42 180 221 32 297
	(B5)	-43.7%	-57.2%	-55.8%	-39%	-49.9%	-50.8%	-27.5%	-42.8%
	RSD (B5)	27%	42.6%	21.3%	39.2%	62%	28.5%	29.5%	28.8%
	RM SD (B5)	132 114	122 838	110 157	81 709	68 191	81 638	69 368	157 515
	E.4 Details about subsamples at each settlement date

Table 8 :

 8 Actual global reserve at various settlement dates, and other statistics about the data under consideration.

Available in the R package CASdatasets.

Acknowledgments

The authors would like to thank the two anonymous referees for their constructive comments that contributed to significantly improve the paper. This work was conducted within the Research Chair DIALog under the aegis of the Risk Foundation, a joint initiative by CNP Assurances and ISFA, Université Claude Bernard Lyon 1 (UCBL).

This leads to the following estimator of the survival function,

when there is no ties.

In case of ties, let (t 1 , ..., t k ) denote the distinct values taken by (Y i ) 1≤i≤n , leading to

C Modification of the CART algorithm

We present here the weighted CART algorithm (wCART), used throughout the paper to take into account censoring and truncation phenomenons.

Step 1: R 1 (z) = 1 for all z = (y, x), and n 1 = 1 (corresponds to the root node).

Step k+1: Let (R 1 , ...R n k ) denote the rules obtained at step k. For j = 1, ..., n k ,

• if all observations such that δ i R j (Y i , X i ) = 1 have the same characteristics, then keep rule j as it is no longer possible to segment the population;

• else, rule R j is replaced by two rules R j1 and R j2 determined in the following way:

for each component Z (l) of Z = (Y, X) (l = 1, ..., d + 1), define the best threshold l) to split the data, such that z (l) = arg min z (l) s(R j , z (l) ), with

.

Then, select the best component to consider, that is l = arg min l s(R j , z (l) ).

Define the two new rules R j1 (z) = R j (z)1

• Let n k+1 denote the new number of rules.

Stopping rule: stop if n k+1 = n k .

Remark C.1. Sometimes, M may be a deterministic function of T . The algorithm is thus easily adapted by replacing N by Y , and Z = (Y, X) by X. In this simpler situation (where only one single censored/truncated variable has to be predicted), competing approaches include survival trees and forests, see [START_REF] Hemant Ishwaran | Random survival forests[END_REF] and [START_REF] Molinaro | Tree-based multivariate regression and density estimation with right-censored data[END_REF].

D Simulations

We give here the characteristics of the simulated samples used in the simulation study.

E Applications E.1 Descriptive statistics for the motor insurance dataset

Depending on the type of the variable, we give different indicators: for categorical variables, exposure for each category is provided. Concerning numerical variables, the minimum, the maximum, the median, the mean, and the standard deviation are given. Table 5: Descriptive statistics on available information for ausautoBI8999 dataset.

We also show that reporting started much later than accidents, which led us to remove oldest observations since it is not clear whether these information are reliable.


Adriana Rebelo

Ilse Eidhof

Vivian Cintra

Léna Guillot-Noel

Claudia Pereira

Dagmar Timmann

Andreas Traschütz

Ludger Schöls

Giulia Coarelli

Alexandra Durr

Biallelic loss-of-function variations in PRDX3 cause cerebellar ataxia