Robust bootstrap procedures for the chain-ladder method

TL;DR

This paper develops and compares robust bootstrap methods for the chain-ladder claims reserving technique, addressing outlier sensitivity and improving inference accuracy in insurance reserve estimation.

Contribution

It introduces and evaluates new robust bootstrap procedures tailored for the chain-ladder method, enhancing reliability in the presence of outliers.

Findings

Robust bootstrap methods outperform classical ones with outlier data.

The proposed procedures provide more accurate standard error estimates.

Application to real data demonstrates practical effectiveness.

Abstract

Insurers are faced with the challenge of estimating the future reserves needed to handle historic and outstanding claims that are not fully settled. A well-known and widely used technique is the chain-ladder method, which is a deterministic algorithm. To include a stochastic component one may apply generalized linear models to the run-off triangles based on past claims data. Analytical expressions for the standard deviation of the resulting reserve estimates are typically difficult to derive. A popular alternative approach to obtain inference is to use the bootstrap technique. However, the standard procedures are very sensitive to the possible presence of outliers. These atypical observations, deviating from the pattern of the majority of the data, may both inflate or deflate traditional reserve estimates and corresponding inference such as their standard errors. Even when paired with a robust chain-ladder method, classical bootstrap inference may break down. Therefore, we discuss and implement several robust bootstrap procedures in the claims reserving framework and we investigate and compare their performance on both simulated and real data. We also illustrate their use for obtaining the distribution of one year risk measures.

Tables

Table 1. Table 1: Run-off triangle (incremental claim amounts).

$\frac{\text{development year}}{\text{accident year}}$	1	2	…	$j$	…	$J-1$	$J$
1	$Y_{11}$	$Y_{12}$	$\dots$	$Y_{1j}$	$\dots$	$Y_{1,J-1}$	$Y_{1J}$
2	$Y_{21}$	$Y_{22}$	$\dots$	$Y_{2j}$	$\dots$	$Y_{2,J-1}$
$⋮$	$\dots$	$\dots$	$\dots$	$\dots$	$\dots$
$i$	$Y_{i1}$	$Y_{i2}$	$\dots$	$Y_{ij}$
$⋮$	$\dots$	$\dots$
$I$	$Y_{I1}$

Table 2. Table 2: claims data from Taylor and Ashe (1983).

	1	2	3	4	5	6	7	8	9	10
1	357,848	766,940	610,542	482,940	527,326	574,398	146,342	139,950	227,229	67,948
2	352,118	884,021	933,894	1,183,289	445,745	320,996	527,804	266,172	425,046
3	290,507	1,001,799	926,219	1,016,654	750,816	146,923	495,992	280,405
4	310,608	1,108,250	776,189	1,562,400	272,482	352,053	206,286
5	443,160	693,190	991,983	769,488	504,851	470,639
6	396,132	937,085	847,498	805,037	705,960
7	440,832	847,631	1,131,398	1,063,269
8	359,480	1,061,648	1,443,370
9	376,686	986,608
10	344,014

Table 3. Table 3: Simulated run-off triangle.

	1	2	3	4	5	6	7	8	9	10
1	9891	7914	6092	4842	3812	3012	2370	1842	1522	1139
2	10373	8210	6624	5304	4146	3218	2486	1939	1525
3	10890	8611	6702	5345	4119	3328	2502	2054
4	11463	9036	7091	5596	4316	3429	2783
5	12050	9385	7424	5828	4425	3533
6	12675	9793	7594	6026	4623
7	13116	10185	7928	6315
8	13473	10652	8234
9	14205	11101
10	14598

Table 4. Table 4: Pearson residuals obtained after fitting a Poisson GLM to the data in Table 3 .

	1	2	3	4	5	6	7	8	9	10
1	-62.0	8.4	10.3	9.9	-61.1	97.9	-11.6	-6.5	-7.8	-0.0
2	-12.0	53.3	55.0	83.4	-27.3	-82.3	12.9	15.4	9.8
3	-63.3	9.2	12.0	11.8	-63.7	96.7	-13.3	-5.6
4	-3.7	63.8	61.5	54.9	-26.3	-81.6	18.8
5	-0.6	66.7	65.2	57.6	-26.4	-82.2
6	120.9	-129.2	-109.2	-96.5	89.1
7	-46.3	28.1	29.1	27.0
8	-41.3	35.6	35.0
9	-29.3	46.7
10	0.0

Table 5. Table 5: Robust weights obtained after fitting a robust Poisson GLM to the data in Table 3 .

	1	2	3	4	5	6	7	8	9	10
1	1.00	1.00	1.00	1.00	1.00	0.00	1.00	1.00	1.00	1.00
2	0.70	1.00	0.81	0.05	0.73	1.00	1.00	1.00	1.00
3	1.00	1.00	1.00	1.00	1.00	0.00	0.90	1.00
4	1.00	1.00	1.00	1.00	1.00	1.00	1.00
5	1.00	1.00	1.00	1.00	0.75	0.76
6	0.00	1.00	1.00	1.00	0.00
7	1.00	1.00	1.00	1.00
8	1.00	1.00	1.00
9	1.00	1.00
10	1.00

Table 6. Table 6: Coverage percentages for the benchmark model using the Cordeiro adjusted residuals 𝒓 C superscript 𝒓 𝐶 \boldsymbol{r}^{C} in the bootstrap process.

		75%	90%	95%	99.5%
	TCL	71.50	86.60	92.00	98.50
Clean	CB	72.60	88.00	91.80	98.70
	WB	66.50	78.60	84.80	93.60
	FRB	69.60	83.50	89.20	96.70
	IFB	72.00	86.80	92.40	98.70
	TCL	100.00	100.00	100.00	100.00
Cont 1	CB	87.00	100.00	100.00	100.00
	WB	65.50	79.30	83.80	94.40
	FRB	70.00	85.70	91.30	97.90
	IFB	63.10	80.80	88.40	97.60
	TCL	100.00	100.00	100.00	100.00
Cont 2	CB	92.00	100.00	100.00	100.00
	WB	74.90	85.70	90.00	96.60
	FRB	77.40	89.40	93.50	98.50
	IFB	76.20	91.80	98.30	100.00

Table 7. Table 8: Computation times expressed in minutes.

	TCL	CB	WB	IFB	FRB
mean	0.98	6.75	19.47	1.03	0.016
st. dev.	0.06	2.90	14.44	0.05	0.002

Table 8. Table 9: Coverage results for Schiegl model.

		75%	90%	95%	99.5%
Clean	TCL	71.84	87.09	92.66	98.62
	FRB	71.02	85.75	91.34	97.90
	IFB	67.13	80.50	86.67	95.51
Cont 1	TCL	100.00	100.00	100.00	100.00
	FRB	65.37	83.18	90.31	98.14
	IFB	55.43	72.31	80.89	93.79
Cont 2	TCL	0.00	100.00	100.00	100.00
	FRB	59.23	79.12	87.47	97.35
	IFB	47.68	65.64	74.91	90.89

Table 9. Table 11: Real run-off triangle from Schedule P data.

	1	2	3	4	5	6	7	8	9	10
1988	794	569	497	409	221	55	62	9	0	0
1989	696	499	429	437	227	154	10	57	0	0
1990	799	690	391	529	326	34	4	37	45	0
1991	796	405	234	287	131	263	438	0	1	4
1992	742	498	368	199	346	49	1	0	0	1
1993	637	518	140	231	79	34	19	89	5	0
1994	703	427	364	194	98	40	15	6	61	3
1995	580	375	227	132	121	18	10	52	45	6
1996	563	266	106	203	71	30	23	4	3	19
1997	324	193	173	30	36	49	44	19	3	0

Table 10. Table 12: Weights from the robust chain-ladder method for the Schedule P data.

	1	2	3	4	5	6	7	8	9	10
1988	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00
1989	1.00	1.00	1.00	1.00	1.00	0.63	1.00	1.00	1.00
1990	1.00	1.00	1.00	0.93	1.00	1.00	1.00	1.00
1991	1.00	1.00	1.00	1.00	1.00	0.22	0.07
1992	1.00	1.00	1.00	0.83	0.57	1.00
1993	1.00	0.80	0.76	1.00	0.91
1994	1.00	1.00	1.00	1.00
1995	1.00	1.00	1.00
1996	1.00	1.00
1997	1.00

Table 11. Table 13: Pearson residuals obtained from the classical chain-ladder algorithm for the Schedule P data.

	1	2	3	4	5	6	7	8	9	10
1988	-0.94	0.70	6.35	2.76	-1.06	-5.44	-5.82	-4.28	-0.02	0.00
1989	-3.25	-1.34	3.71	5.13	-0.04	4.40	-10.17	4.28	0.02
1990	-2.78	3.96	-0.54	7.45	4.45	-7.92	-11.37	0.09
1991	-0.56	-6.05	-7.07	-3.51	-6.77	14.33	27.72
1992	0.20	0.10	1.77	-6.77	9.12	-5.19
1993	3.13	7.16	-7.27	-0.87	-6.56
1994	1.73	-0.77	3.69	-5.45
1995	1.05	0.07	-1.64
1996	3.03	-3.70
1997	-0.00

Table 12. (a)

	1	2	…	$J$
1	$Y_{11}^{\left(b\right)}$	$Y_{12}^{\left(b\right)}$	$\dots$	$Y_{1J}^{\left(b\right)}$
2	$Y_{21}^{\left(b\right)}$	$Y_{22}^{\left(b\right)}$	$\dots$
$\dots$	$\dots$	$\dots$
$I$	$Y_{I1}^{\left(b\right)}$

Table 13. (a)

	1	2	…	$J$
1	$Y_{11}^{\left(b\right)}$	$Y_{12}^{\left(b\right)}$	$\dots$	$Y_{1J}^{\left(b\right)}$
2	$Y_{21}^{\left(b\right)}$	$Y_{22}^{\left(b\right)}$	$\dots$
$\dots$	$\dots$	$\dots$
$I$	$Y_{I1}^{\left(b\right)}$

Table 14. (b)

	1	2	…	$J$
1	$Y_{11}^{\left(b\right)}$	$Y_{12}^{\left(b\right)}$	$\dots$	$Y_{1J}^{\left(b\right)}$
2	$Y_{21}^{\left(b\right)}$	$Y_{22}^{\left(b\right)}$	$\dots$	${\hat{Y}}_{2J}^{\left(b\right)}$
$\dots$	$\dots$	$\dots$	$\dots$	$\dots$
$I$	$Y_{I1}^{\left(b\right)}$	${\hat{Y}}_{I2}^{\left(b\right)}$	$\dots$	${\hat{Y}}_{IJ}^{\left(b\right)}$

Table 15. (c)

	1	2	…	$J$
1	$Y_{11}^{\left(b\right)}$	$Y_{12}^{\left(b\right)}$	$\dots$	$Y_{1J}^{\left(b\right)}$
2	$Y_{21}^{\left(b\right)}$	$Y_{22}^{\left(b\right)}$	$\dots$	${\hat{Y}}_{2J}^{\left(b\right)}$
$\dots$	$\dots$	$\dots$	$\dots$	$\dots$
$I$	$Y_{I1}^{\left(b\right)}$	${\hat{Y}}_{I2}^{\left(b\right)}$	$\dots$	${\tilde{Y}}_{IJ}^{\left(b\right)}$

Equations

$$R=\sum _{i=2}^n\sum _{j=n-i+2}^n{Y}_{ij}.$$

$${\hat{f}}_j=\frac{{\sum }_{i=1}^{n-j+1}{C}_{ij}}{{\sum }_{i=1}^{n-j+1}{C}_{i,j-1}}.$$

$${\hat{C}}_{i,n-i+2}$$

$${\hat{C}}_{ik}$$

$$Y_{ij}\stackrel{i.i.d.}{\sim }\text{mbox}Poisson({\zeta }_i{\xi }_j),$$

$$\log \left({\mu }_{ij}\right)=\tau +{\alpha }_i+{\beta }_j=x_{ij}^t\mathbf{\theta }.$$

$$r_{ij}=\frac{y_{ij}-{\mu }_{ij}}{\sqrt{{\mu }_{ij}}}\text{for }1⩽i,j⩽n\text{ and }i+j⩽n+1,$$

$${\mathbf{r}}^E=\sqrt{\frac{N}{N-p}}\mathbf{r},$$

$$r_{ij}^P=\frac{r_{ij}}{\sqrt{1-h_{kk}}}$$

$$E[\mathbf{r}]$$

$$\text{mbox}Var[\mathbf{r}]$$

$${\mathbf{r}}^C=\frac{\mathbf{r}-E[\mathbf{r}]}{\sqrt{\text{mbox}Var[\mathbf{r}]}}.$$

$${\mathbf{y}}^{(b)}={\mathbf{r}}^{(b)}\sqrt{\mathbf{\mu }}+\mathbf{\mu }.$$

$$\sum _{i=1}^n\sum _{j=1}^{n-i+1}\mathbf{\psi }(y_{ij},{\mu }_{ij})=\sum _{i=1}^n\sum _{j=1}^{n-i+1}\left[{\psi }_c(r_{ij})w(x_{ij})\frac{1}{V^{1/2}({\mu }_{ij})}{\mu }_{ij}^{\prime }-a(\mathbf{\theta })\right]=0.$$

$$\psi_{c}(r_{ij})=\left\{\begin{array}[]{ll}r_{ij}&|r_{ij}|\leqslant c\\ c\cdot\textrm{sign}(r_{ij})&|r_{ij}|\geqslant c,\end{array}\right.$$

$$\sum _{i=1}^n\sum _{j=1}^{n-i+1}\left[\left({\psi }_c(r_{ij})-E\left[{\psi }_c(r_{ij})\right]\right)w(x_{ij})\sqrt{{\mu }_{ij}}{x}_{ij}\right]=0$$

$$E\left[{\psi }_c(r_{ij})\right]=$$

$$+\sqrt{{\mu }_{ij}}\left[\text{mbox}P\left(Y_{ij}=j_1\right)-\text{mbox}P\left(Y_{ij}=j_2\right)\right]$$

$$\frac{{\psi }_c(r_{ij})}{r_{ij}}.$$

$${\hat{\mathbf{\theta }}}_{\text{FRB}}^{(b)}=\hat{\mathbf{\theta }}+{\left[I_p-\nabla g_N(\hat{\mathbf{\theta }})\right]}^{-1}\cdot \left(g_N^{(b)}(\hat{\mathbf{\theta }})-\hat{\mathbf{\theta }}\right),$$

$${\psi }_N(\hat{\mathbf{\theta }})=\sum _{k=1}^N\psi (\hat{\mathbf{\theta }})=0.$$

$${\hat{\mathbf{\theta }}}_{\text{FRB}}^{(b)}=\hat{\mathbf{\theta }}-{\left[\nabla {\psi }_N(\hat{\mathbf{\theta }})\right]}^{-1}\cdot {\psi }_N^{(b)}(\hat{\mathbf{\theta }}).$$

$${\psi }_N(\hat{\mathbf{\theta }})=\sum _{i=1}^n\sum _{j=1}^{n-i+1}\left[\left({\psi }_c(r_{ij})-E\left[{\psi }_c(r_{ij})\right]\right)w(x_{ij})\sqrt{{\mu }_{ij}}{x}_{ij}\right].$$

$$\frac{\partial }{\partial \mathbf{\theta }}\left(\sqrt{{\mu }_{ij}}\right)=\frac{1}{2}{\mu }_{ij}^{1/2}{x}_{ij}\text{ and }\frac{\partial }{\partial \mathbf{\theta }}\left({\psi }_c(r_{ij})\right)=-\frac{1}{2}{I}_{|r_{ij}|⩽c}\left(\frac{y_{ij}+{\mu }_{ij}}{{\mu }_{ij}^{1/2}}\right)x_{ij},$$

$$\nabla {\psi }_N(\hat{\mathbf{\theta }})=\sum _{i=1}^n\sum _{j=1}^{n-i+1}w(x_{ij})x_{ij}$$

$$+{\left({\psi }_c(r_{ij})-E\left[{\psi }_c(r_{ij})\right]\right)\frac{1}{2}{\mu }_{ij}^{1/2}{x}_{ij}]}^t.$$

$$\{\begin{array}{l}{Y}_{ij}\sim F_{i,j}\forall (i,j)\ne (l,m)\\ Y_{l,m}\sim (1-\epsilon )F_{l,m}+\epsilon {\Delta }_z\end{array}$$

$$\text{mbox}IF([z,l,m],T,F)=\underset{\epsilon \downarrow 0}{\lim }\frac{T(F_{l,m,\epsilon ,z})-T(F)}{\epsilon }.$$

$$\text{mbox}RESIF([z,l,m],T_{\mathbf{\theta }}^{nr},F_{{\hat{\mathbf{\theta }}}^r})={\left[\text{mbox}IF([z,l,m],T_{\mathbf{\theta }}^{nr},F_{{\hat{\mathbf{\theta }}}^r}{)}^t\cdot \text{mbox}{V}^{-1}(T_{\mathbf{\theta }}^{nr},F_{{\mathbf{\theta }}^r})\cdot \text{mbox}IF([z,l,m],T_{\mathbf{\theta }}^{nr},F_{{\hat{\mathbf{\theta }}}^r})\right]}^{\frac{1}{2}}$$

$$\text{mbox}V(T_{\mathbf{\theta }}^{nr},F_{{\hat{\mathbf{\theta }}}^r})=\text{mbox}{E}_{F_{{\hat{\mathbf{\theta }}}^r}}\left[\text{mbox}IF([z,l,m],T_{\mathbf{\theta }}^{nr},F_{{\hat{\mathbf{\theta }}}^r})\cdot \text{mbox}IF([z,l,m],T_{\mathbf{\theta }}^{nr},F_{{\hat{\mathbf{\theta }}}^r}{)}^t\right].$$

$$\text{mbox}{RESIF}_{ij}=\text{mbox}RESIF([y_{ij},i,j],T_{\mathbf{\theta }}^{nr},F_{{\hat{\mathbf{\theta }}}^r})$$

$$w_{ij}$$

Full text

Keywords: Claims reserving, Generalized linear models, Outliers, Inference, Solvency II

1 Introduction

One of the main challenges for insurers is the task of estimating the future reserves needed to handle the liabilities related to current insurance contracts. The recent Solvency II directive provides new guidelines for the calculation of these provisions. A comprehensive discussion on the Solvency II directive and its implications may be found in Dreksler et al., (2015). They report that technical provisions under Solvency II consist of the following three components: claims provisions, premium provisions and a risk margin. Whereas claims provisions relate to claim events that have already occurred, the premium provisions relate to all future claim events covered by current insurance and reinsurance obligations. The risk margin is described as the amount needed to make the value of the technical provisions equivalent to the amount a (re)insurance undertaking would be expected to require in order to take over and meet (re)insurance obligations.

In this article we focus on the claims provision component. Although micro-level claims reserving models (see for example Norberg, (1993); Antonio and Plat, (2014)) seem to be very promising, Dreksler et al., (2015) state that the current norm is to use deterministic actuarial techniques to asses the undiscounted element of the earned claim reserves and that this is likely to continue in the immediate future. The most popular and widely used technique is the chain-ladder method. A substantial amount of research has already been devoted to exploring the relationship between the chain-ladder method and various stochastic models (see for example Mack, (1993, 1994), Renshaw and Verrall, (1998), Verrall, 1991a ; Verrall, (2000) and Schiegl, (2015)).

To approach claims reserving in a stochastic manner, the chain-ladder method is typically reformulated within the context of generalized linear models (GLM). However, assessing the variability of claims reserving predictions and the construction of upper limits with a certain confidence level still remains a difficult task. Analytical expressions for standard errors in GLM depend on strong model assumptions (Mack et al., (2006); Venter, (2006) and Wüthrich et al., (2008)). An alternative approach to statistical inference that requires less assumptions is the bootstrap technique (see for example Efron, (1981), Hall et al., (1995) and Chatterjee and Lahiri, (2011)). Also in the claims reserving framework, the bootstrap technique is very popular and has been studied by various authors such as Ashe, (1986), Lowe, (1994), Pinheiro et al., (2003), England and Verrall, (2006) and Barnett and Zehnwirth, (2008).

When applying a claims reserving method in practice it may occur that one or a couple of observations deviate from the main pattern. One might naively expect that if the model assumptions hold for almost all available data, then the results of the classical method also hold approximately. This is unfortunately not the case and it is well known in the robust statistics literature that outlier(s) may heavily influence classical statistical techniques. Verdonck and Debruyne, (2011) have shown theoretically that the classical chain-ladder method is very sensitive to outliers and several robust methods have already been developed in the claims reserving framework (see e.g. Brazauskas et al., (2009), Brazauskas, (2009), Verdonck et al., (2009), Verdonck and Debruyne, (2011), Verdonck and Van Wouwe, (2011) and Pitselis et al., (2015)). Robust methods provide estimates for the claim provisions which resemble the classical estimates that would have been obtained if there were no outliers in the data, while they do not model the outlier generating process. As a consequence of fitting the majority of the data well, robust methods also provide a reliable method to detect outliers. Observations which are flagged as outliers can then be examined in detail by experts to understand their origin. If these flagged observations are errors or one-time events, then the robust method indeed yields reliable estimates of the claim provisions as desired. If the flagged observations are not errors and such deviations may reoccur in the future, then the robust reserve may need to be adjusted to cover such potentially large deviations in the future. We believe that expertise from practitioners in the field of the data is needed to achieve this, possibly combined with techniques from extreme value theory. The idea of combining robust statistics (which typically downweights atypical points) and extreme value statistics (which models the extremes) has also been advocated in Dell’Aquila and Embrechts, (2006) and Hubert et al., (2013). It is however outside the scope of this paper to investigate estimators for the extra amount that needs to be added to the robust reserve estimate in such case.

To our knowledge, robust methods in claims reserving have only focused on reliable point estimates and/or outlier detection. However, confidence intervals may also become very unreliable (incorrect coverage) and/or uninformative (very large interval length) when the data are contaminated. In this article we investigate the influence of outliers on the variability of point estimators and we present some state-of-the-art robust bootstrap procedures in the claims reserving framework. These bootstrap techniques provide a complete approximate predictive distribution, from which characteristics such as the mean, variance and quantiles are easily derived. Note that we do not advocate replacing the classical methodology by the robust counterpart, but we advise to always apply both techniques. If the classical and robust method give approximately the same results (and no outliers are detected by the robust method), then it seems safe to continue the analysis with the traditional reserving techniques that will give the most efficient estimates. However, if both methods yield different results, then the robust procedure is helpful to gain insight into the data and will yield more reliable reserve estimates if the outliers are indeed isolated events.

This paper is structured as follows. Section 2 reviews the chain-ladder method and its formulation as a Poisson GLM (which is the default approach to obtain prediction under Solvency II). In this section, we also present the bootstrap procedure in the context of claims reserving. The robust chain-ladder method is briefly described in Section 3 and it is shown that the classical bootstrap method may produce unreliable results when the data contains outlier(s), even when the classical chain-ladder method is replaced by its robust alternative. Therefore, we present several state-of-the-art robust bootstrap methods for GLM within the context of claims reserving in section 4. Their performance is compared by means of a simulation study in section 5 and on a real dataset from a non-life business line in section 6. We also illustrate how the methodology can be applied to investigate one year risk measures and to detect outliers automatically. Some concluding remarks and potential directions for further research are given in section 7.

2 Classical chain-ladder method and bootstrapping

In this section we introduce the classical chain-ladder method and briefly explain how this method can be derived from a basic GLM. Then, we describe the classical bootstrap procedure to obtain inference corresponding to the estimates in this GLM. Finally, we briefly describe the robust GLM estimator.

2.1 Classical chain-ladder method

Insurers have to build up reserves enabling them to pay outstanding claims and to meet claims arising in the future on the written contracts. The chain-ladder method is a very popular algorithm to compute these reserve estimates. Let $Y_{ij}$ and $C_{ij}$ denote respectively the incremental and cumulative claim amount corresponding to accident year $i$ and development year $j$ for $1\leqslant i\leqslant I$ and $1\leqslant j\leqslant J$. We suppose that $I=J=n$ for ease of notation, but all methodology remains valid when $I\neq J$. Variables for which $i+j\leqslant n+1$ correspond to past claims data and are used to predict the future claims which correspond to the variables with $i+j>n+1$. We define the total or overall reserve $R$ as the sum over all future (incremental) claims111Note that we define the total reserve $R$ as the total amount of the incremental claims in the lower right triangle (as in key references England and Verrall, (2002) and Wüthrich and Merz, (2008)), but in practice the final reserve is an amount set by the insurer based on this value.:

[TABLE]

The chain-ladder method uses cumulative data and assumes the existence of development factors $f_{j}$ for $j=2,\ldots,n$ that are estimated by

[TABLE]

This leads to the following estimates for the future claim amounts:

[TABLE]

The chain-ladder method is typically specified in the GLM framework using the following multiplicative model for the incremental claim amounts:

[TABLE]

with $\zeta_{i}$ the exposure (i.e. expected claims up to the latest observed development year) of accident year $i$ and $\xi_{j}$ the expected claims pattern over the different development periods $j$. Since $\zeta_{i}$ and $\xi_{j}$ can only be determined up to a constant factor, the condition $\sum_{j=1}^{J}\xi_{J}=1$ is often imposed. Formulation (2) has a multiplicative structure for the mean, i.e. $\mu_{ij}=E[Y_{ij}]=\zeta_{i}\xi_{j}$ and hence it is straightforward to use a log link function such that

[TABLE]

This structure now has a parameter for each row $i$ and each column $j$ and typically the constraint $\alpha_{1}=\beta_{1}=0$ is used (to obtain a more convenient interpretation of the parameters). Formulations (2) and (3) are equivalent and are reparameterisations of the same structure. The parameters of the first have physical interpretations, whereas the statistical analysis of the latter is more straightforward. The relationship between the parameters and the proof of the equivalence between both models is given in Kremer, (1982) and Verrall, 1991b . The full parameter vector is thus $\boldsymbol{\theta}=\left(\tau,\alpha_{2},\ldots,\alpha_{I},\beta_{2},\ldots,\beta_{I}\right)$ and the number of parameters to be estimated equals $2n-1$. The predictor variable $x_{ij}$ consists of the constant $1$ (to include an intercept) and indicators corresponding to each of the rows and columns of the claims triangle except for the first. By inserting parameter estimates $\hat{\boldsymbol{\theta}}$ into equation (3) and exponentiating, estimates for the future claim amounts and the total reserve are then obtained. When estimating the parameters $\boldsymbol{\theta}$ by maximum likelihood, the results are identical to the estimates obtained by the deterministic chain-ladder method using development factors as shown in Kremer, (1982). For more details we refer to England and Verrall, (2002), Hoedemakers et al., (2005) and Wüthrich and Merz, (2008).

2.2 Classical bootstrap procedure

Point estimates are very useful of course, but the natural next step is to also estimate their precision. We consider the bootstrap methodology to derive inference corresponding to the point estimates. To create bootstrap samples we follow the procedure described in Kaas et al., (2009), which is based on England and Verrall, (1999, 2002).

First, the Pearson residuals, defined as

[TABLE]

are calculated by replacing the $\mu_{ij}$ by their estimates $\hat{\mu}_{ij}=e^{x_{ij}^{t}\hat{\boldsymbol{\theta}}}$. When the number of observations is not very large, the residuals suffer from a small sample bias. To remedy this, England and Verrall, (1999) and England, (2002) propose to adjust the residuals $\boldsymbol{r}=\{r_{ij}|1\leqslant i,j\leqslant n;i+j\leqslant n+1\}$ by multiplying them by a correction factor:

[TABLE]

where $N=n(n+1)/2$ is the sample size and $p=2n-1$ is the number of fitted parameters. Pinheiro et al., (2003) argue that it is better to use an individual adjustment factor for each residual $\{r_{k}|k=1,\ldots,N\}=\{r_{ij}|1\leqslant i,j\leqslant n;i+j\leqslant n+1\}$ as opposed to a global correction factor. In correspondence with the classical linear regression model they propose to use the hat matrix $H$ of the model to standardize the Pearson residuals as follows

[TABLE]

where $h_{kk}$ is the corresponding element on the diagonal of the hat matrix. For a Poisson GLM, the hat matrix is given by $H=W^{1/2}X(X^{T}WX)^{-1}X^{T}W^{1/2}$, where $W$ is a diagonal matrix with elements $\{\mu_{ij}|1\leqslant i,j\leqslant n;i+j\leqslant n+1\}$ on the diagonal (see e.g. McCullagh and Nelder, (1983) for more details). Cordeiro, (2004) recently provided general expressions for the expectation and variance of Pearson residuals in the GLM framework up to the first order $O(N^{-1})$. Using equations (13) and (14), calculated in Appendix 8.2, we obtain that

[TABLE]

where ‘diag’ indicates that we extract the diagonal elements of a matrix. These correction terms lead to the following adjusted residuals

[TABLE]

Replacing the parameters $\mu_{ij}$ in the expressions of $H$ and $W$ by their estimates $\hat{\mu}_{ij}$ yields the expectation and variance correction term of the sample Pearson residuals.

The bootstrap methodology can be applied by using the original Pearson residuals or any of the above explained adjustments. From now on, we use the general notation $\boldsymbol{r}^{*}$ for the residuals of choice, i.e. $\boldsymbol{r}^{*}$ is any of $\boldsymbol{r},\boldsymbol{r}^{E},\boldsymbol{r}^{P}$ or $\boldsymbol{r}^{C}$. In the next sections, all results shown are obtained for $\boldsymbol{r}^{*}=\boldsymbol{r}^{C}$, but the conclusions also hold for the other adjustment options.

After obtaining the residuals, the following bootstrap steps are performed many times, e.g. for $b=1,\,\ldots,$ 10,000:

Resample with replacement from the residuals $\boldsymbol{r}^{*}$ leading to a new vector of residuals $\boldsymbol{r}^{(b)}$. Since there is only one observation for the last accident year and the last development year, the corresponding residuals (on the cornerpoints of the run-off triangle) will be exactly zero when fitting a GLM. As suggested by England, (2002) and Pinheiro et al., (2003), these two residuals are not resampled (whereas all other residuals have equal resampling probabilities). 2. 2.

Backtransform the new residuals $\boldsymbol{r}^{(b)}$ in order to obtain a new pseudo-history:

[TABLE] 3. 3.

Refit the GLM on the new pseudo-history $\boldsymbol{y}^{(b)}$ in order to obtain a new total future reserve estimate $R^{(b)}$.

3 Robust chain-ladder method and classical bootstrapping

The robust chain-ladder method of Verdonck and Debruyne, (2011) is briefly described and the effect of outliers on the total reserve estimates obtained by the classical and robust chain-ladder method is illustrated on the well studied claims data set of Taylor and Ashe, (1983) as well as on a simulated data set. For both estimators we also evaluate the effect of outliers on an estimated upper limit for the total reserve obtained by using the classical bootstrap procedure.

3.1 Robust chain-ladder method

Verdonck and Debruyne, (2011) already showed the sensitivity of the classical chain-ladder method to outliers and therefore proposed a robust alternative, which is based on a general class of $M$-estimators of Mallow’s type proposed by Cantoni and Ronchetti, (2001), with the following set of estimating equations:

[TABLE]

The vector $\boldsymbol{\theta}=\left(c,\alpha_{2},\ldots,\alpha_{I},\beta_{2},\ldots,\beta_{I}\right)$ as in equation (3) and $\mu_{ij}^{\prime}$ denotes the derivative of $\mu_{ij}$ with respect to $\boldsymbol{\theta}$. The constant $a(\boldsymbol{\boldsymbol{\theta}})=\frac{2}{n(n+1)}\sum_{i=1}^{n}\sum_{j=1}^{n-i+1}E\left[\psi_{c}(r_{ij})\right]\omega(x_{ij})\frac{1}{V^{1/2}(\mu_{ij})}\mu_{ij}^{\prime}$ ensures Fisher consistency and

[TABLE]

is the Huber function. The weight function $w(x_{ij})$ gives less weight to observations with outlying predictors, whereas the Huber function limits the effect of observations with an outlying response by truncating the residual. Therefore, the influence function of deviations on the response and on the predictors is bounded separately. Note that for a log-link Poisson GLM one obtains $V^{1/2}(\mu_{ij})=\sqrt{\mu_{ij}}$. Cantoni and Ronchetti, (2001) showed that the estimating equations may equivalently be written as

[TABLE]

with

[TABLE]

and $j_{1}=\lfloor\mu_{ij}-c\sqrt{\mu_{ij}}\rfloor$ and $j_{2}=\lfloor\mu_{ij}+c\sqrt{\mu_{ij}}\rfloor$.

3.2 Illustration of the effect of outlier(s)

In this section we illustrate the classical and robust techniques on the widely studied data set of Taylor and Ashe, (1983), shown in table 2.

The classical and robust total reserve estimates equal respectively $18,680,856$ and $18,562,327$. To study the effect of an outlier on the classical and robust chain-ladder method, one observation $y_{ij}$ is adjusted to $\kappa y_{ij}$ for $\kappa\in[0,10]$. We calculate the total reserve estimates before and after this adjustment, respectively denoted by $\hat{R}$ and $\hat{R}^{\kappa}$. For observations $y_{41}$, $y_{43}$ and $y_{64}$, we plot the ratio $\frac{\hat{R}^{\kappa}}{\hat{R}}$ in Figure 1.

It is immediately clear that the influence of an outlier on the robust chain-ladder method (orange dashed line) is very limited in all cases, even if $\kappa$ becomes very large. On the other hand, the outlier has a large, unbounded effect on the classical chain-ladder method (blue solid line). Depending on the location of the outlier in the run-off triangle, the estimated reserve $\hat{R}^{\kappa}$ can largely increase or decrease. For example, an unusually large claim amount in the past ($\kappa>1$) cannot only result in a much bigger reserve estimate, but may also result in a reserve estimate that is much smaller even when no anomalies occur in the future anymore. This illustrates the sensitivity of the classical methodology and hence, there is no guarantee that the estimated classical reserve suffices to cover future claim amounts.

We now investigate the effect of an outlier on upper limits for the total reserve estimated by the classical bootstrap procedure. For the classical chain-ladder the distribution of the obtained classical bootstrap total reserve estimates is shown in orange (light color) in Figure 2. The 99.5$\%$ quantile of this bootstrap distribution corresponds to a total reserve estimate of approximately $27.8$ million (indicated by a vertical line). To illustrate the effect of an outlier, we multiply observation $y_{27}$ in Table 2 by ten. From the corresponding bootstrap distribution, which is shown in blue (dark color) in Figure 2, it is seen that now a $99.5\%$ quantile of roughly 55 million is obtained. This clearly indicates the sensitivity of the classical chain-ladder method in combination with the classical bootstrap procedure.

To limit the effect of outlying observations in the obtained bootstrap distribution, we now combine the robust chain-ladder method with the traditional bootstrap procedure and apply it to the contaminated data set of Taylor and Ashe, (1983) (i.e. $y_{27}$ is multiplied by ten). Note that for the adjusted residuals $\boldsymbol{r}^{P}$ and $\boldsymbol{r}^{C}$, we need the hat matrix corresponding to the robust Poisson GLM estimator of Cantoni and Ronchetti, (2001). This expression is derived in Appendix 8.1. To compare the distribution of the obtained total reserves for the contaminated data with the previously obtained original reserve estimates (i.e. the reserves obtained by bootstrapping the classical chain-ladder method on the non-contaminated data), the quantile-quantile plot is given in Figure 3. As expected, a large range of the obtained quantiles are now relatively close to the corresponding original quantiles (since we plugged in the robust chain-ladder method). However, the upper quantiles of the robust chain-ladder method on contaminated data deviate significantly from their original counterparts. This is caused by the resampling with replacement in the bootstrap process. During the resampling process, it might occur that the outliers are selected so many times in a bootstrap sample that the robust method can no longer withstand this huge amount of contamination. This will happen several times if we generate 10,000 resamples, yielding unreliable bootstrap estimates that affect the bootstrap distribution, especially if one is interested in an upper quantile that lies far in the tail of the distribution.

This sensitivity of the bootstrap technique, even in combination with robust point estimators, is well-known in the robust statistics literature and was first discussed by Stromberg, (1997). This seminal paper gave rise to the development of several robust alternatives for the classical bootstrap. In the next section we discuss the robust bootstrap methods that have been proposed (up to our knowledge) and we adapt them to the claims reserving framework.

To end this section, we compare the outlier detection capacity of the classical and robust chain ladder on the simulated data triangle in Table 3 which was generated from the benchmark model (see section 5.1.1). In this triangle we replaced observations $y_{16},y_{36},y_{61},y_{65}$ and $y_{24}$ respectively by 33,000, 35,000, 120,000 and 65,000 and 7000. The classical and robust total reserve estimates are $314,240$ and $155,086$, respectively. Hence, inspection of the data is needed to understand this difference. In practice, one often tries to detect outliers using diagnostics starting from a fit obtained by the classical estimation method. The Pearson residuals obtained from the classical chain-ladder method are given in Table 4. Besides both corner points and $r_{51}$, all residuals have an absolute value much larger than 3 (which is typically used as a benchmark in a classical outlier detection rule). Apparently, the classical method is so heavily affected by the outliers that the fitted model does not allow to detect the deviating observations anymore. For example, all observations in row two and six and column six have become inflated, making it impossible to flag observation $y_{24}$ as an outlier. This is called the masking effect. On the other hand, observation $y_{62}$ which is a regular observation yields the largest residual and hence appears to be the most outying observation. This effect is known as swamping. To avoid these adverse effects, the goal of robust statistics is to obtain a fit which is similar to the fit we would have found without outliers in the data. An important benefit is that it then becomes possible again to detect the outliers by their large deviation from the robust fit. The robust chain-ladder method automatically provides a robustness weight for each observation, which is given by

[TABLE]

If a residual $r_{ij}$ is large compared to its Huberized value $\psi_{c}(r_{ij})$, then the observation receives a low weight and is identified as a potential outlier. The robustness weights obtained by fitting the robust chain-ladder method to the data are shown in Table 5. The outlying cells are now correctly flagged by the robustness weights and can be studied in more detail.

4 Robust bootstrap methods

To obtain robust inference, Stromberg, (1997) proposed to use robust measures of location and scale of the classical bootstrap distribution corresponding to the robust estimators. Typically, the median and the median absolute deviation (see Hampel, (1974)) are applied instead of the classical mean and standard deviation. Obviously, we cannot use this method to calculate quantiles of the bootstrap distribution. Stromberg, (1997) also proposed the deleted-$d$ jackknife procedure, in which a bootstrap sample consists of $N-d$ data points instead of $N$ data points to increase the probability of sampling a clean subset. However, in our framework this poses problems as the resampled triangle would then be incomplete.

On the other hand, Singh, (1998) proposed to robustify the bootstrap procedure itself by resampling from a winsorized version of the original sample. In a regression context, one may resample either the residuals or the original observations themselves. Due to the specific structure of a claims triangle, we resample the residuals. To winsorize the residuals the quantiles $Q_{c}$ and $Q_{1-c}$ of the residuals are determined for $c\in\left]0,0.5\right[$. The $c$ percent smallest residuals are then replaced by $Q_{c}$ whereas the $c$ percent largest ones are replaced by $Q_{1-c}$.

More recently, Salibian-Barrera and Zamar, (2002) and Salibian-Barrera et al., (2008) proposed a fast and robust bootstrap method that is asymptotically consistent under reasonable regularity conditions. The corresponding R-package (Van Aelst and Willems, (2013)) considers three multivariate settings: principal component analysis, Hotelling tests and multivariate regression. Another recent alternative is proposed by Amado and Pires, (2004) and Amado et al., (2014). They suggest to resample observations with different probabilities such that potentially harmful observations have low sampling probability. They propose to determine these probabilities based on the influence function of the estimator. We now discuss both state-of-the-art robust bootstrap methods in more detail and describe their implementation in the claims reserving framework.

4.1 Fast and Robust Bootstrap

For each pseudo-history, most robust bootstrap procedures need to recalculate the robust estimator, which in itself is already computationally intensive compared to the classical estimator. An important aim of the fast and robust bootstrap (FRB) is to reduce this computational load (Salibian-Barrera and Zamar,, 2002; Salibian-Barrera et al.,, 2008).

Within the context of claims reserving the FRB procedure can be described as follows. Let $\boldsymbol{\theta}\in\mathbb{R}^{2n-1}$ be the parameter of interest. Note that $\boldsymbol{\theta}$ consists of the intercept $c$ and the parameters $\alpha_{i}$ and $\beta_{j}$ for $2\leqslant i,j\leqslant n$, as before. Now, suppose that $\hat{\boldsymbol{\theta}}$ is a solution of a fixed-point equation $\hat{\boldsymbol{\theta}}=g_{N}(\hat{\boldsymbol{\theta}})$ where $g_{N}$ depends on the sample of size $N=n(n+1)/2$ composed of the historical claims data. Similarly, for a resampled pseudo-history the corresponding estimator $\hat{\boldsymbol{\theta}}^{(b)}$ solves the equation $\hat{\boldsymbol{\theta}}^{(b)}=g^{(b)}_{N}(\hat{\boldsymbol{\theta}}^{(b)})$ where the function $g^{(b)}_{N}$ now depends on the bootstrap sample. As a first step one may consider $\tilde{\boldsymbol{\theta}}^{(b)}=g^{(b)}_{N}(\hat{\boldsymbol{\theta}})$ as an approximation for the bootstrap replicate $\hat{\boldsymbol{\theta}}^{(b)}$ of the estimator $\hat{\boldsymbol{\theta}}$. This simple procedure avoids having to recompute the robust estimator for each bootstrap sample, which may be computationally hard. However, these approximate bootstrap replicates typically underestimate the true variability since all bootstrap approximations $\tilde{\boldsymbol{\theta}}^{(b)}$ are based on the same initial estimate $\hat{\boldsymbol{\theta}}$. Therefore, Salibian-Barrera et al., (2008) proposed a linear correction resulting in the FRB bootstrap replicates

[TABLE]

where $\nabla g_{N}(\hat{\boldsymbol{\theta}})$ corresponds to the gradient of $g_{N}$ evaluated at $\hat{\boldsymbol{\theta}}$ and $I_{p}$ denotes the identity matrix of dimension $p=2n-1$.

This procedure was later extended by Camponovo et al., (2012) to accommodate $M$-estimators defined by the solution of equations of the following type

[TABLE]

Fast and robust bootstrap estimates are then obtained as follows:

[TABLE]

The FRB paradigm can be used for the Poisson GLM estimator of Cantoni and Ronchetti, (2001). Based on the estimating equations (5) we define

[TABLE]

To calculate $\nabla\psi_{N}(\hat{\boldsymbol{\theta}})$ we readily obtain

[TABLE]

with $I$ the indicator function. It then follows that

[TABLE]

The term $E\left[\psi_{c}(r_{ij})\right]$ however is not differentiable when using its exact expression. Therefore, we use the $\chi^{2}$ distribution as an approximation for the Poisson distribution and afterwards apply the Wilson-Hilferty approximation. The detailed derivation can be found in Appendix 8.3.

4.2 Influence Function Bootstrap

Amado and Pires, (2004) propose to determine resampling probabilities for the observations by using the influence function corresponding to the classical estimator. The definition of the influence function is based on the concept of a statistical functional, which we recall first. Let $T_{N}$ denote a $p$-dimensional estimator defined for any sample size $N$, with $N=(n+1)n/2$ and $p=2n-1$ in our setting. A statistical functional $T$ corresponding to the estimator $T_{N}$ is a map defined on $p$-variate distributions $G$ on $\mathbb{R}^{p}$ such that $T(G_{N}(\{y_{ij}|1\leqslant i,j\leq n;i+j\leqslant n+1\})=T_{N}(\{y_{ij}|1\leqslant i,j\leq n;i+j\leqslant n+1\})$ where $G_{N}$ is the empirical distribution of the sample $\{y_{ij}|1\leqslant i,j\leq n;i+j\leqslant n+1\}$.

Let $\{Y_{ij}\,|\,1\leqslant i,j\leqslant n;i+j\leq n+1\}\sim F$ and assume that the marginal distribution of $Y_{ij}$, denoted by $F_{\mu_{ij}}=F_{i,j}$, is Poisson distributed with mean $\mu_{ij}$, see equation (3). For any $1\leqslant l,m\leqslant n$ with $l+m\leqslant n+1$ and $z\geqslant 0$, we define $F_{l,m,\varepsilon,z}$ such that

[TABLE]

where $\Delta_{z}$ is a Dirac measure putting all its mass at value $z$.

The influence function of the functional $T$ at a distribution $F$ is then defined as

[TABLE]

Intuitively, the influence function measures the change of the estimator when the model distribution is perturbed by an infinitesimal small amount of contamination at location $z$.

Amado and Pires, (2004) propose to measure the effect of each data point on the parameter estimates by estimating the value of the population influence function at each observation using estimates of the model parameters based on the sample. However, if a non-robust estimator is employed then an outlier may already have such a high effect on the sample estimates themselves that it results in a low value for the influence function evaluated at the outlier. This phenomenon is known as masking in the robust statistics literature. Conversely, the influence function of a robust estimator will also have small values at the outliers, since outliers should have a low effect on robust estimators (Pison and Van Aelst,, 2004). Therefore, Amado and Pires, (2004) consider the influence function of the classical non-robust estimator (denoted by $T^{nr}_{\boldsymbol{\theta}}$) using parameter estimates $\hat{\boldsymbol{\theta}}^{r}$ obtained by a robust estimator, i.e. the influence function is evaluated at $F_{\hat{\boldsymbol{\theta}}^{r}}$. By combining robust parameter estimates with the influence function of the non-robust estimator, its values for outliers should be very high compared to regular observations. This idea is formalized by the Robust Empirical Standardized Influence Function (RESIF) that is defined as:

[TABLE]

where

[TABLE]

Based on the RESIF function the bootstrap procedure is adjusted. We now discuss this procedure in more detail in the context of claims reserving.

Calculate the values:

[TABLE]

for each of the observations in the sample $\{y_{ij}|1\leqslant i,j\leq n;i+j\leqslant n+1\}$. 2. 2.

Choose constants $d,\gamma$ and $c$ and calculate for each observation $y_{ij}$ a weight $w_{ij}$ as follows:

[TABLE]

where

[TABLE]

and $I_{[0,c]}$ indicates the indicator function on the interval $[0,c]$. 3. 3.

Compute the resampling probability

[TABLE]

for each observation.

The bootstrap procedure is then initiated by resampling the residuals $\boldsymbol{r}^{*}$ using these resampling probabilities. Note that the residuals of the two corner points are still not resampled.

The RESIF for the Poisson GLM estimator of Cantoni and Ronchetti, (2001) is derived in Appendix 8.4 and the result is given in Theorem 1.

Theorem 1**.**

The Robust Empirical Standardized Influence Function of the M-estimator defined by equation (8) for the Poisson GLM is given by

[TABLE]

Amado and Pires, (2004) propose to apply the classical GLM estimator to calculate the parameter estimates corresponding to the bootstrap samples. This drastically reduces the computation time that is needed for the bootstrap. The robust initial estimate combined with the weights based on the influence function protect against outliers in the bootstrap samples. The resulting procedure is both fast and robust. Note that the consistency of this methodology has not been formally proven.

5 Simulation Study

To compare the different bootstrap procedures, their performance is evaluated in an extensive simulation study. Three different models to generate complete run-off triangles (lower and upper part) are considered. We start with a simple benchmark model already allowing us to study several properties of the various bootstrap methods. Then, we turn our attention to more realistic simulation models from the claims reserving literature. To simulate run-off triangles in the claims reserving framework we focus on the models proposed by Schiegl, (2004) and Cowell, (2009). Before discussing the simulation setup and results, we briefly describe the three data generating models and their corresponding parameters. Values for these parameters were chosen to obtain realistically looking data sets (based on the proposed values in Schiegl, (2004) and Cowell, (2009)). Note that other constants than below were also studied, but they led to similar results and conclusions and hence were omitted.

5.1 Data generating models

5.1.1 Benchmark model

The benchmark model is a simplification of the model proposed by Schiegl, (2004). The entries in cell $i,j$ ($i=1,\ldots,10$ and $j=1,\ldots,10$) are simulated from a Poisson distribution with parameter $\lambda_{ij}$ defined by

[TABLE]

where $\lambda_{0}=$10,000 gives the expected value of claim amounts at the first accident and development year, $\eta_{1}=0.3$ indicates the reduction of claim amounts at development period $\frac{n}{2}$ compared to the claim amount in the first development year and $\eta_{2}=0.05$ is the rate of increase for the claim amounts over consecutive accident years.

5.1.2 Schiegl method

This method is proposed in Schiegl, (2004) and the simulation of the run-off triangles (containing 10 accident and development years) is based on the collective model, where a different distribution is used for the number of claims and for a single claim amount. The number of claims $N_{ij}$ (for $i=1,\ldots,10$ and $j=1,\ldots,10$) follows a Poisson distribution with mean $\lambda_{ij}$ as defined by equation (11) of the benchmark model. Individual claim amounts are then generated using the Gamma distribution with mean $r=$ 2,000 such that the final data are obtained as follows

[TABLE]

5.1.3 Kaishev method

Cowell, (2009) noted that, given the same set of parameter values, the distribution of reserves in the Schiegl method is the same regardless of the sample drawn for the upper triangle (since each $C_{ij}$ value is sampled independently from the others). Therefore, Cowell, (2009) presents a slightly more complex model (proposed by Kaishev through personal communication) so that, even if all parameters are kept fixed, different simulations will lead to different distributions of reserves.

As before, the method is based on the collective model and we again generate individual claim amounts using a Gamma distribution with mean $r=$2,000. The total number of claims $N_{i}$ for accident year $i=1,\ldots,10$ is generated from a Poisson distribution with mean $\lambda_{i}$. Each claim (for $k=1,\ldots,N_{i})$ follows a Gamma distribution (yielding a certain claim amount $Z_{k}$) and this amount is then assigned randomly to a certain development year $j=1,\ldots,10$. This random allocation is done using a $\mbox{Beta}(\alpha=1,\beta=1.6)$ distribution (as suggested in Cowell, (2009)). Similar as in the Schiegl method we defined $\lambda_{i}$ as

[TABLE]

where $\lambda_{0}$ and $\eta_{2}$ were also chosen as before.

Note that plugging in a Pareto or a Lognormal distribution for the individual claim size (instead of a Gamma distribution) in the Schiegl or Kaishev model yielded similar conclusions.

5.2 Simulation setup

During the simulation study we consider the following bootstrap procedures:

• •

TCL: standard approach using the classical chain-ladder method and standard bootstrap

• •

CB: robust chain-ladder method combined with the classical bootstrap

• •

WB: winsorized bootstrap

• •

FRB: fast and robust bootstrap procedure

• •

IFB: influence function bootstrap

Note that for IFB, the weight function to determine the resampling probability of each residual contains three tuning parameters, see equation (10). These values need to be chosen to guarantee a sufficiently high efficiency at the uncontaminated case whilst offering robustness when contamination is present. Taking into account this tradeoff between efficiency and robustness, the parameter $c$ is calculated as the $90\%$ quantile of the $\mbox{RESIF}_{ij}$ values, whereas the parameters $d$ and $\gamma$ were set to $30$ and $10$ respectively. For the winsorized bootstrap the $10\%$ highest and $10\%$ smallest residuals are winsorized.

The bootstrap procedures are tested for both clean and contaminated data. Two contamination settings are shown. Firstly, contamination was generated by multiplying observations $Y_{23}$ and $Y_{32}$ in the original clean data by a factor of five. Secondly, observations $Y_{23}$ and $Y_{41}$ in the original clean data are multiplied by two and a half. We also studied other locations and multiplication factors to create outliers, but these results led to similar conclusions. The corresponding tables are not included in the paper, since they do not add new insights in the proposed methodology.

5.3 Simulation results

5.3.1 Benchmark model

For each bootstrap procedure, we investigated the quantiles obtained from 10,000 bootstrap runs. Since we simulate complete data using the benchmark model, the true total reserve can be calculated from the lower triangle according to equation (1). Hence, for 10,000 generated data sets we recorded whether its true total reserve is smaller than the calculated bootstrap quantile, for different probability levels. These coverage percentages are shown in Table 7 for different bootstrap procedures and confidence levels. If the bootstrap procedure is able to approximate the $(1-\alpha)$-quantiles with $\alpha\in(0,0.5)$ well, one may expect the true total reserve to be smaller than the bootstrap quantile in $(1-\alpha)\%$ of the cases.

For the uncontaminated data, it is clear that the WB procedure leads to significantly lower coverage percentages than the other methods. The TCL, CB and IFB methodologies give comparable results whereas the FRB procedure yields slightly lower coverage probabilities (but still much better than the WB method).

To study the effect of outliers, we repeated the simulation study using the first (i.e. multiplying $Y_{23}$ and $Y_{32}$ by five) and second (i.e. multiplying $Y_{23}$ and $Y_{41}$ by two and a half) contamination setting. In both situations, the FRB and IFB clearly outperform the other bootstrap procedures. We immediately see that the TCL and CB procedures are highly influenced by the outliers as almost all the coverage percentages equal $100\%$ (this is because the obtained bootstrap quantiles became very large due to the outliers), whereas the WB bootstrap quantiles are too small leading to the smallest coverage percentages. It can be seen that in contamination scenario two the estimates of the $95\%$ and $99.5\%$ quantiles obtained by the IFB bootstrap procedure are still affected by the outliers. In this setting the outliers are closer to the real data and therefore it is more difficult to detect them.

To study the effect of the residual adjustment on the coverage percentages, Table 7 shows the percentages obtained by using the original Pearson residuals. These coverage percentages are lower and further away from the target than those obtained with the adjusted Pearson residuals. The coverage percentages based on the adjusted residuals using the England and Verrall, (1999) adjustment factor ($\boldsymbol{r}^{E}$) and the correction proposed by Pinheiro et al., (2003) ($\boldsymbol{r}^{P}$) were comparable to the results found in Table 7 (using $\boldsymbol{r}^{C}$) and are therefore omitted.

To compare the speed of the different bootstrap procedures, the average computation time for a bootstrap run consisting of 10,000 bootstrap samples for a run-off triangle, simulated using the benchmark model, is shown in Table 8. These computations were performed on an Intel Ivy Bridge Xeon E5-2680V2 2,8 GHz CPU and averaged across 100 runs.

The FRB method clearly outperforms the other alternatives in terms of computational efficiency. Given the high time complexity and bad performance of the CB and WB methods in the (contaminated) benchmark model, we discard these procedures in further studies. The other procedures are now compared using the more advanced simulation models.

5.3.2 Schiegl and Kaishev model

We have applied the remaining bootstrap techniques on clean and contaminated data sets simulated by the Schiegl and Kaishev models. Boxplots of the obtained $99.5\%$ quantiles based on 10,000 bootstrapped estimates for the total reserve are shown in Figures 5 and 5. The boxes in the boxplots are white for the clean data, orange for the first contamination setting (i.e. multiplying $y_{23}$ and $y_{32}$ by 5) and red for the second contamination setting (i.e. multiplying $y_{23}$ and $y_{41}$ by 2.5).

It is immediately clear from Figures 5 and 5 that TCL breaks down as soon as the data is contaminated, whereas the other methods give results that are almost the same as those obtained by the classical method on the initial clean data without outliers. The presence of the outliers thus only has a minor effect on these methods. In order to better compare the robustness properties of IFB and FRB, we have removed the results for TCL on contaminated data in Figures 7 and 7. When no outliers are included, the boxplots for FRB are nearly identical to the ones obtained with TCL on the clean data, whereas the results for IFB are clearly lower. These differences are even more obvious when looking at the contaminated settings: the boxplots for FRB hardly change when outliers are added, whereas the boxplots for IFB are getting further away from the target (i.e. the results obtained with TCL on the clean data).

The same conclusions are obtained from the results in Tables 10 and 10, where the 10,000 bootstrapped estimates for the total reserve are compared with the true reserve (obtained by summing up the lower triangle of the simulated data set). The percentages in the tables again represent in how many of the 10,000 runs the true reserve was covered by the corresponding confidence interval for different confidence levels. For the uncontaminated data the coverage percentages for TCL are fairly close to their expected values with FRB percentages being only slightly smaller. The IFB method achieves smaller coverage percentages. When outliers are added, TCL clearly breaks down: the bootstrap quantiles become very large, which explains the high coverage rates in almost all situations. However, for data simulated from the Schiegl model and the second contamination setting, the $75\%$ bootstrap quantiles never succeed in covering the real future reserve. Similar observations can be made for the Kaishev model. On the other hand, the coverage percentages for the robust methods always remain close to the target percentages. The results for FRB are again more stable than for IFB.

5.3.3 Conclusion

From this simulation study we can conclude that TCL and CB are very sensitive to outliers and that WB consistently underestimates the total reserve quantile. The IFB and FRB procedures yield efficient and robust results and FRB is generally less influenced by outliers than IFB. Taking computational speed into account as well, we advice to combine the robust chain-ladder method with the FRB approach for prediction purposes in the claims reserving framework.

6 Real data application

We now apply our optimal robust methodology (i.e. robust chain-ladder method combined with the FRB approach) to a run-off triangle from practice, which is shown in Table 11.

This real data set is extracted from the schedule P

• Analysis of Losses and Loss Expenses in the National Association of Insurance Commissioners (NAIC) database and is available on the website of the Casualty Actuarial Society 222http://www.casact.org/research/index.cfm?fa=loss_reserves_data. This database contains data on several lines of business for a variety of businesses in the United States and here we focus on the other liability line of business for the Rockford Mutual Insurance Company. The true total reserve is 1812 (the sum of the observations in the lower triangle of Table 11).

We recommend to always apply both the classical and robust procedures to a given dataset and compare their results. For this example, the classical chain-ladder method leads to an total reserve estimate of 2837 whereas the robust alternative leads to a total claim reserve estimate of 2304, which is a difference of roughly 20$\%$. The obtained bootstrap $99.5\%$ quantile for the total reserve is 4781 for the classical methodology and 3285 for the robust procedure based on FRB. This clear discrepancy may indicate the presence of possible outlier(s) in the data. The bootstrap distributions are shown in figure 8. The true total reserve is marked with a red dashed line.

An advantage of the robust chain-ladder method is that the most influential points can easily be detected using the weights of the robust chain-ladder method defined in equation (7). From Table 12 it is clear that the sixth and seventh development year of accident year 1991 have a very low weight and should be examined in more detail. If these observations are actually errors, then the robust estimate is a reliable total reserve estimate, otherwise the total reserve estimate should be adjusted accordingly (taking into account the possibility of observing such atypical claims in future).

For comparison, the Pearson residuals obtained using the classical chain-ladder method are shown in Table 13. Development year six and seven of accident year 1991 have the largest residuals, but other observations also have large residuals (e.g. above 5 or even above 10 in absolute value). It is thus not obvious to determine how many observations have an unusually large residual and it is also not clear what adjustments should be made to these observations or how they should be treated in the bootstrap process On the other hand, our robust techniques yield an automatic, data driven and objective methodology to treat and detect possible outliers.

Under the new Solvency II regulations insurers are also compelled to investigate one year risk measures (Dreksler et al.,, 2015). Wüthrich et al., (2009) proposed to use the Claims Development Result (CDR) for this purpose. Let $D_{k}=\left\{Y_{ij}\,|\,i+j\leqslant k+1\right\}$ denote the upper triangle of stochastic variables up to calender year k with $n\leqslant k\leqslant 2n-1$. Note that for $k=n$ this corresponds to the past claims data and for $k=2n-1$ the full square of data is available. The true CDR for claims in accident year $i\in\left\{1,\ldots,n\right\}$ during the next calender year $(k,k+1]$ is defined as

[TABLE]

where $\mbox{TR}_{i}=\sum_{j=k-i+1}^{n}{Y_{ij}}$ denotes the total claims amount for accident year $i$. The aggregated CDR is then defined as

[TABLE]

To estimate the CDR the stochastic quantities are replaced by their observations and estimates respectively. Several authors already considered this problem including Diers, (2008) and Ohlsson and Lauzeningks, (2009). Lacoume, (2008), Boisseau, (2006) and Boumezoued et al., (2012) complement these papers and propose to add the following steps to the bootstrap procedure (explained in Section 2.2) in every run $b=1,\ldots,B$:

Denote the fitted values obtained from the GLM in step 3 of the bootstrap loop as $\hat{Y}^{(b)}_{ij}$, see Table 14(b). The bootstrap reserve $R^{(b)}$ then corresponds to $\sum_{i=2}^{n}{\sum_{j=n-i+2}^{n}{\hat{Y}^{(b)}_{ij}}}$. An extended bootstrap history is constructed by adding the first fitted anti-diagonal $\hat{Y}^{(b)}_{ij}$ with $i+j=n+2$ to the bootstrapped historic observations. Based on this extended bootstrap history the GLM is refitted leading to the fitted values $\widetilde{Y}^{(b)}_{ij}$, see Table 14(c). This leads to a new total future reserve estimate $R^{(b)}_{ext}=\sum_{i=3}^{n}{\sum_{j=n-i+3}^{n}{\widetilde{Y}^{(b)}_{ij}}}$. 2. 5.

A bootstrapped value for the CDR may then be obtained by computing :

[TABLE]

The bootstrapped CDR is shown in Figure 9.

The $99.5\%$ quantile corresponds to 1,677 for the classical procedure and 1,160 for the robust procedure. As the Value at Risk under the Solvency II regulations is directly related to this quantile, this example again shows the impact of outliers on standard claims reserving results and the good performance of the proposed robust alternative.

7 Conclusions

The bootstrap procedure is a powerful and popular technique to obtain information from a single sample of data, which would usually be obtained using analytic techniques. It complements sample estimates with measures of accuracy (e.g. standard error, confidence interval, prediction error). In claims reserving, the bootstrap procedure is typically used to produce a standard error of prediction and an approximate predictive distribution for the random future losses. In this paper we studied the effect of outliers on these procedures and it is shown that it is not sufficient to plug a robust estimator into the classical bootstrap procedure. Therefore, we implemented and compared all existing (up to our knowledge) robust bootstrap methods that could be adapted to the claims reserving framework. The Fast and Robust Bootstrap (FRB) approach appears to be a good choice in the sense that it yields very reliable results and is computationally fast.

The presence of outliers may reveal that the data are more heterogeneous than assumed and cannot be handled by the original statistical model. Outliers can be isolated or may come in clusters, indicating that there are subgroups in the population that behave differently. We believe that a robust analysis together with some exploratory data analysis techniques can reveal structures in the data that would remain hidden in a traditional analysis. Robust methods try to fit the model followed by the majority of the data: if the data contain no outliers, they give approximately the same results as the classical method, while if the data is contaminated they give approximately the same results as the classical method applied to the outlier-free data. The best estimates obtained by the robust methodology is hence the result that one would obtain when applying the traditional methodology on the data without outliers. Note that we do not know the proportion of outliers in advance.

Using our robust methodology, insurers obtain estimates that are similar to the results that would have been obtained from outlier-free data. Moreover, the procedure also detects outlying observations (if any). Of course, it is important to examine the detected outliers and to understand the reasons for their atypical behaviour. In this way the proposed robust methodology is helpful to gain insight in the data. If the outliers are errors or due to causes that are not expected to re-occur in the future, then the robust estimates allow to build up more realistic reserve estimates and confidence intervals. However, if such atypical observations are likely to happen again in future, it is necessary to model also their process (outside the scope of this paper) and to predict how much extra reserve (besides the robust total reserve estimate) is needed to cope with such outlying observations in future years. In such a case, the final estimate may for instance be equal to the robust total reserve estimate plus a safe margin for the outliers.

Acknowledgments

This work was supported by the Internal Funds KU Leuven under Grants $C16/15/068$ and $C24/15/001$; the Flemish Science Foundation (FWO) under Grant $1523915N$. The computational resources and services used in this work were provided by the VSC (Flemish Supercomputer Center), funded by the Research Foundation - Flanders (FWO) and the Flemish Government – department EWI.

8 Appendix

8.1 Hat matrix of the robust GLM estimator

To calculate the adjusted Pearson residuals $\boldsymbol{r}^{C}$ an expression for the hat matrix corresponding to the robust estimator of Cantoni and Ronchetti, (2001) is needed. To this end we rewrite the estimator as an iteratively weighted least squares problem. Using a first order Taylor expansion of equation (4), we obtain

[TABLE]

where $\boldsymbol{\overline{\psi}}(\boldsymbol{y},\boldsymbol{\mu})=\sum_{i=1}^{n}\sum_{j=1}^{n-i+1}{\boldsymbol{\psi}(y_{ij},\mu_{ij})}=X^{t}D1_{N}$ with $1_{N}$ a unit vector of length $N$ and D the diagonal matrix with elements

[TABLE]

We then replace $-\frac{\partial}{\partial\boldsymbol{\theta}}\boldsymbol{\overline{\psi}}(\boldsymbol{y},\boldsymbol{\mu})$ by $-\mbox{E}\left[-\frac{\partial}{\partial\boldsymbol{\theta}}\boldsymbol{\overline{\psi}}(\boldsymbol{y},\boldsymbol{\mu})\right]$. Cantoni and Ronchetti, (2001) computed this to be equal to

[TABLE]

where $B$ is the diagonal matrix with elements $\mbox{E}\left[\psi_{c}(r_{ij})\frac{Y_{ij}-\mu_{ij}}{\mu_{ij}}\right]\frac{w(\boldsymbol{x}_{ij})}{\mu^{3/2}_{ij}}$ and

[TABLE]

Equation (12) then reduces to

[TABLE]

which can be rewritten as

[TABLE]

assuming $B$ and $(X^{t}BX)$ are invertible. This expression may be recognized as that of a weighted least squares regression where $\left[B^{-1}D1_{N}+X\boldsymbol{\theta}\right]$ is regressed onto $X$ with weight matrix $B$. The hat matrix for the robust estimator of Cantoni and Ronchetti, (2001) is thus

[TABLE]

8.2 Calculation of the correction terms proposed by Cordeiro (2004)

Based on the general expression given by Cordeiro, (2004), the following vector of expected Pearson residuals is obtained for the Poisson GLM model

[TABLE]

whereas the vector of variances of the Pearson residuals is given by

[TABLE]

For these expressions it holds that

[TABLE]

where $z$ is the diagonal of $Z$ and 1 is a unit vector. Note that $(QHJ-T)z$ corresponds to the diagonal of $(QHJ-T)Z$ and since

[TABLE]

we obtain

[TABLE]

Therefore we obtain

[TABLE]

Since the diagonal of $W\left(X(X^{T}WX)^{-1}X^{T}\right)$ equals the diagonal of the hat matrix $H=W^{1/2}\left(X(X^{T}WX)^{-1}X^{T}\right)W^{1/2}$, we finally obtain that

[TABLE]

Rewriting the expression for $E[\boldsymbol{r}]$ we obtain that

[TABLE]

In practice the matrices $H$ and $W$ may be replaced by their sample estimates.

8.3 Calculations for Fast and Robust bootstrap for GLM

Since the term $E\left[\psi_{c}(r_{ij})\right]$ is not differentiable using its exact expression, we calculate this term by approximating the Poisson distribution with a $\chi^{2}$ distribution and then using the Wilson-Hilferty approximation, see (Johnson et al.,, 1993, chapter 4, p. 162). This leads to the following lemma:

Lemma 1**.**

If $Y_{ij}\sim\mbox{Poisson}(\mu_{ij})$ it holds that

[TABLE]

where

[TABLE]

Using Lemma (1) we may rewrite equation (6) as follows:

[TABLE]

with

[TABLE]

[TABLE]

and $j_{1}=\mu_{ij}-c\mu_{ij}^{1/2}$ and $j_{2}=\mu_{ij}+c\mu_{ij}^{1/2}$. Note that in turn $j_{1}$ and $j_{2}$ are approximated by considering their non-truncated versions. However $j_{1}$ must be bounded from below by zero to ensure all square roots may be calculated without problems. Clearly we need to calculate the terms

[TABLE]

To calculate the first term note that

[TABLE]

where we used the generalized Leibniz integration rule. Moreover,

[TABLE]

This leads to

[TABLE]

Differentiating $\widetilde{z}_{l}$ and $\widehat{z}_{l}$ yields

[TABLE]

Moreover,

[TABLE]

with $\frac{\partial}{\partial\boldsymbol{\theta}}j_{1}=x_{ij}(\mu_{ij}-\frac{c}{2}\mu_{ij}^{1/2})$ and $\frac{\partial}{\partial\boldsymbol{\theta}}j_{2}=x_{ij}(\mu_{ij}+\frac{c}{2}\mu_{ij}^{1/2})$.

8.4 Proof of theorem 1

The estimator by Cantoni and Ronchetti, (2001) may be generally written as

[TABLE]

where

[TABLE]

can be seen as the score function and $a\left(\boldsymbol{\boldsymbol{\theta}}\right)=\frac{1}{N}\sum_{i=1}^{n}\sum_{j=1}^{n-i+1}E\left[\nu\left(y_{ij},\mu_{ij}\right)\right]w\left(x_{ij}\right)\mu_{ij}^{\prime}$. This is the structure of an $M$-estimator which leads to the following influence function:

[TABLE]

with $M\left(\psi,F\right)=-E\left[\frac{\partial}{\partial\boldsymbol{\theta}}\psi\left(z,\mu_{lm}\right)\right]$. More details may be found in e.g. Hampel et al., (1986). The asymptotic variance can then be found as

[TABLE]

where $Q\left(\psi,F\right)=E\left[\psi\left(z,\mu_{lm}\right)\psi\left(z,\mu_{lm}\right)^{t}\right]$. Using the above expressions, the standardised influence function of the $M$-estimator becomes

[TABLE]

where the arguments of the function have been dropped for lightning the notation.

Choosing $\nu\left(z,\mu_{lm}\right)=\psi_{c}(r_{lm})\frac{1}{V^{1/2}(\mu_{lm})}$ leads to the robust estimator of Cantoni and Ronchetti, (2001) as detailed in section 3.1. However, for the influence function bootstrap procedure we need the expression for the influence function of the classical estimator. Therefore, we choose $\nu\left(z,\mu_{lm}\right)=\frac{r_{lm}}{V^{1/2}(\mu_{lm})}$ and set $w\left(\cdot\right)=1$ which are the functions corresponding to the Poisson GLM yielding the same estimates as the classical chain-ladder. Thus (9) reduces to

[TABLE]

It is easily checked that $E[r_{lm}]=0$, which leads to $a\left(\boldsymbol{\theta}\right)=0$. Hence, $\psi\left(z,\hat{\mu}_{lm}\right)=r_{lm}\hat{\mu}^{1/2}_{lm}$ and thus (15) becomes

[TABLE]

This immediately simplifies to

[TABLE]

which proves the theorem.