Continuous chain-ladder with paid data

TL;DR

This paper develops a continuous-time, non-parametric framework for predicting outstanding liabilities using hazard functions and kernel smoothing, demonstrating consistency and improved estimation methods.

Contribution

It introduces a continuous-time chain-ladder model with histogram and kernel estimators, extending traditional methods with new smoothing techniques and theoretical consistency results.

Findings

The proposed methods are consistent under increasing claim data and decreasing aggregation levels.

Kernel-based estimators outperform traditional development factors in simulations.

Real-data application confirms the effectiveness of the new estimators.

Abstract

We introduce a continuous-time framework for the prediction of outstanding liabilities, in which chain-ladder development factors arise as a histogram estimator of a cost-weighted hazard function running in reversed development time. We use this formulation to show that under our assumptions on the individual data chain-ladder is consistent. Consistency is understood in the sense that both the number of observed claims grows to infinity and the level of aggregation tends to zero. We propose alternatives to chain-ladder development factors by replacing the histogram estimator with kernel smoothers and by estimating a cost-weighted density instead of a cost-weighted hazard. Finally, we provide a real-data example and a simulation study confirming the strengths of the proposed alternatives.

Tables

Table 1. Table 1: Forecasted claim amount for future years and by accident year estimated by chain-ladder (CL), manually corrected local linear (LL) and local constant estimators (LC).

future year	2014	2015	2016	2017	2018	2019	2020	2021	2022	2023	total
CL	3972072	2997371	2241199	1613228	1157272	699522	401116	190240	49779	0	13321800
LL	4606662	2676251	1944796	1396172	940474	562684	296735	120538	19231	10	12563553
LC	4706341	2734479	2014402	1473544	1020410	639236	357675	172538	49401	3420	13171446
accident year	2004	2005	2006	2007	2008	2009	2010	2011	2012	2013	total
CL	0	31977	216901	496947	830549	1503785	1778076	2221387	2387308	3854869	13321800
LL	18	28133	179632	440042	823458	1323252	1734902	2055000	2364939	3614176	12563553
LC	5144	62617	248154	499977	906861	1403358	1792929	2100946	2372382	3779076	13171446

Table 2. Table 2: Scenarios in the simulation study.

Scenario	$T$	$U$	$Z$
1	decreasing beta	truncated mixed normal	moderately decreasing
2	decreasing beta	truncated mixed normal	heavily decreasing
3	decreasing beta	boundary challenge	moderately decreasing
4	decreasing beta	boundary challenge	heavily decreasing
5	mixture of betas	truncated mixed normal	moderately decreasing
6	mixture of betas	truncated mixed normal	heavily decreasing
7	mixture of betas	boundary challenge	moderately decreasing
8	mixture of betas	boundary challenge	heavily decreasing

Table 3. Table 3: Median, mean and standard deviation of the squared relative errors err 2 superscript err 2 \operatorname{err}^{2} in the reserve estimate for 100, 1000, 10,000, and 100,000 observations. The statistics are taken over 1000 simulation runs.

		LL		LC		CL
		Median	Mean (s.d.)	Median	Mean (s.d.)	Median	Mean (s.d.)
1	100	0.2937	0.3897 (0.4169)	0.2821	0.4008 (0.4046)	0.3525	0.5465 (0.4169)
	1000	0.0999	0.1198 (0.0911)	0.0945	0.1170 (0.0952)	0.1230	0.1459 (0.0911)
	10000	0.0376	0.0439 (0.0326)	0.0353	0.0411 (0.0303)	0.0425	0.0501 (0.0326)
	1e+05	0.0253	0.0259 (0.0155)	0.0242	0.0251 (0.0146)	0.0162	0.0186 (0.0155)
2	100	0.3140	0.3892 (0.3650)	0.2473	0.3743 (0.4307)	0.3626	0.5370 (0.3650)
	1000	0.1094	0.1317 (0.1032)	0.1164	0.1328 (0.0956)	0.1266	0.1514 (0.1032)
	10000	0.0464	0.0541 (0.0380)	0.0656	0.0694 (0.0434)	0.0537	0.0583 (0.0380)
	1e+05	0.0302	0.0308 (0.0168)	0.0429	0.0431 (0.0171)	0.0256	0.0267 (0.0168)
3	100	0.3350	0.3871 (0.3078)	0.2773	0.3209 (0.2608)	0.4443	0.6169 (0.3078)
	1000	0.1129	0.1309 (0.0979)	0.1147	0.1317 (0.0943)	0.1397	0.1581 (0.0979)
	10000	0.0510	0.0589 (0.0407)	0.1056	0.1055 (0.0510)	0.0962	0.0969 (0.0407)
	1e+05	0.0392	0.0393 (0.0192)	0.1000	0.1001 (0.0189)	0.0824	0.0821 (0.0192)
4	100	0.3241	0.3552 (0.3022)	0.4216	0.4199 (0.2195)	0.4687	0.6987 (0.3022)
	1000	0.1198	0.1393 (0.1025)	0.2368	0.2424 (0.1193)	0.1979	0.2078 (0.1025)
	10000	0.0721	0.0776 (0.0497)	0.1941	0.1939 (0.0537)	0.1669	0.1664 (0.0497)
	1e+05	0.0475	0.0478 (0.0203)	0.1650	0.1644 (0.0196)	0.1445	0.1443 (0.0203)
5	100	0.2871	4.189e-01 (0.8021)	0.1877	2.831e-01 (0.3990)	0.3198	1.663e+11 (0.8021)
	1000	0.1255	0.1399 (0.0952)	0.1161	0.1286 (0.0869)	0.1512	0.2437 (0.0952)
	10000	0.0490	0.0575 (0.0437)	0.0660	0.0724 (0.0465)	0.0763	0.0841 (0.0437)
	1e+05	0.0270	0.0293 (0.0195)	0.0479	0.0479 (0.0210)	0.0288	0.0319 (0.0195)
6	100	0.3169	3.943e-01 (0.4260)	0.2303	2.800e-01 (0.2365)	0.3452	1.028e+10 (0.4260)
	1000	0.1214	0.1385 (0.1017)	0.1228	0.1384 (0.0949)	0.1502	0.2113 (0.1017)
	10000	0.0461	0.0539 (0.0408)	0.0636	0.0694 (0.0454)	0.0680	0.0758 (0.0408)
	1e+05	0.0295	0.0320 (0.0194)	0.0522	0.0519 (0.0199)	0.0301	0.0334 (0.0194)
7	100	0.5853	7.398e-01 (1.0259)	0.5016	5.832e-01 (0.6756)	0.6462	6.066e+11 (1.0259)
	1000	0.2089	0.2954 (0.5423)	0.2858	0.3258 (0.3315)	0.3892	0.4983 (0.5423)
	10000	0.0981	0.1134 (0.1183)	0.2109	0.2157 (0.1070)	0.1802	0.1846 (0.1183)
	1e+05	0.0648	0.0675 (0.0414)	0.1699	0.1702 (0.0568)	0.1489	0.1478 (0.0414)
8	100	0.4721	5.814e-01 (1.0504)	0.4595	4.445e-01 (0.4085)	0.6020	1.307e+11 (1.0504)
	1000	0.1661	0.2012 (0.1496)	0.2717	0.2816 (0.1633)	0.3565	0.4672 (0.1496)
	10000	0.1119	0.1217 (0.0822)	0.2317	0.2303 (0.0908)	0.2104	0.2059 (0.0822)
	1e+05	0.0934	0.0923 (0.0423)	0.1962	0.1948 (0.0487)	0.1819	0.1796 (0.0423)

Table 4. Table 4: Distribution, basis price, and price factor of phone brands and models.

Brand	Probability	Basis price	Model	Probability	Price factor
Brand 1	0.45	$600	0	0.05	1
Brand 2	0.30	$550	1	0.10	1.15
Brand 3	0.15	$300	2	0.35	${1.15}^2$
Brand 4	0.10	$150	3	0.50	${1.15}^3$

Table 5. Table 5: Incidence hazard rates and parameters of the beta distribution, which determines the proportion of the phone price that is paid for a claim.

Incident	Yearly hazard rate	$\alpha$	$\beta$
Breakage	0.15	2	5
Oxidation	0.05	5	3
Theft	$0.05\times \text{model}$	5	0.5

Equations

$$N_i(t)=I(t\ge T_i),$$

$$T_i\le today-U_i,$$

$$N_i^R(t)=I(t\ge T_i^R),T_i^R=\mathcal{T}-T_i,$$

$$\mathcal{F}^{R}_{it}=\sigma\mathopen{}\mathclose{{}\left(\bigg{\{}T^{R}_{i}\leq s:\ s\leq t\bigg{\}}\cup\bigg{\{}U_{i}\leq s:\ s\leq t\bigg{\}}\cup\mathcal{N}}\right),$$

$$\lambda_{i}^{R}(t)=\lim_{h\downarrow 0}h^{-1}E\mathopen{}\mathclose{{}\left[N^{R}_{i}\mathopen{}\mathclose{{}\left\{(t+h)-}\right\}-N^{R}_{i}(t-)|\ \mathcal{F}^{R}_{it-}}\right]=\alpha^{R}(t)Y_{i}^{R}(t),$$

$${\alpha }^R(t)$$

$$Y_i^R(t)$$

$${\tilde{\lambda }}_i^R(t)$$

$${\tilde{\alpha }}_{\ast }^R(t)$$

$${\tilde{Y}}_i^R(t)$$

$${\mathcal{R}}_n(t)$$

$$\widetilde{f}_{T}(t)=\frac{E\mathopen{}\mathclose{{}\left[Z|\ T=t}\right]}{E[Z]}f_{T}(t).$$

$$\widetilde{\alpha}^{R}(t)=\frac{\widetilde{f}^{R}_{T}(t)}{\widetilde{S}^{R}_{T}(t)}=\frac{E\mathopen{}\mathclose{{}\left[Z|\ T^{R}=t}\right]}{E[Z|T^{R}\geq t]}\frac{f_{T}^{R}(t)}{\int_{t}^{\mathcal{T}}f_{T}^{R}(s)\mathrm{d}s}.$$

$$\underset{t\in [0,\mathcal{T}]}{\sup }{\mathcal{R}}_n(t)\to 1.$$

$${\tilde{\alpha }}_{\ast }^R(t)={\tilde{\alpha }}^R(t).$$

$$r_n=\frac{{\int }_0^{\mathcal{U}}{\int }_{1-u}^{\mathcal{T}}{\tilde{f}}_{T,U}(t,u)\mathrm{d}t\mathrm{d}u}{{\int }_0^{\mathcal{U}}{\int }_0^{\mathcal{T}-u}{\tilde{f}}_{T,U}(t,u)\mathrm{d}t\mathrm{d}u}\sum _{i=1}^n{Z}_i,$$

$$\lim_{\varepsilon\downarrow 0}\sum_{i=1}^{n}\int\mathopen{}\mathclose{{}\left[\mathopen{}\mathclose{{}\left\{\frac{1}{\varepsilon}\int_{s}^{s+\varepsilon}\mathrm{d}\widetilde{N}^{R}_{i}(w)-\widehat{\widetilde{\alpha}^{R}}(t)\widetilde{Y}^{R}_{i}(s)}\right\}^{2}-\xi(\varepsilon)}\right]W_{h}(s,t)\frac{1}{\widetilde{Y}_{i}^{R}(s)}\,\mathrm{d}s,$$

$$W_h^H(s,t)=I\{s\in [t_l,t_{l+1})\}.$$

$${\hat{\tilde{\alpha }}}_{H,h}^R(l)=\frac{{\sum }_{i=1}^n{\int }_0^{\mathcal{T}}I\{s\in [t_l,t_{l+1})\}\mathrm{d}{\tilde{N}}_i^R(s)}{{\sum }_{i=1}^n{\int }_0^{\mathcal{T}}I\{s\in [t_l,t_{l+1})\}{\tilde{Y}}_i^R(s)\mathrm{d}s},l=1,\dots ,m.$$

$$\gamma (t)=S_T^R(t)F_U(\min (\mathcal{U},t)).$$

$$(nh)^{1/2}\mathopen{}\mathclose{{}\left\{\widehat{\widetilde{\alpha}}^{R}_{H,h}(l)-\widetilde{\alpha}^{R}(t)-B_{H}(t)}\right\}\rightarrow N\mathopen{}\mathclose{{}\left\{0,\sigma_{H}^{2}(t)}\right\},$$

$$B_H(t)$$

$${\sigma }_H^2(t)$$

$$\widehat{\widetilde{S}}{}_{T}^{R}(t)=\prod_{s\leq t}\mathopen{}\mathclose{{}\left\{1-\Delta\widehat{\widetilde{A}}{}^{R}(s)}\right\},$$

$$\left(\begin{array}{c}{\hat{\theta }}_0\\ ⋮\\ {\hat{\theta }}_p\end{array}\right)=\arg \underset{\theta \in {\mathbb{R}}^{p+1}}{\min }\underset{\epsilon \downarrow 0}{\lim }\sum _{i=1}^n\int$$

$$\times K_h(t-s)\frac{1}{{\tilde{Y}}_i^R(s)}\mathrm{d}s.$$

$$2\theta \sum _{i=1}^n{\int }_0^{\mathcal{T}}{K}_h(t-s){\tilde{Y}}_i^R(s)\mathrm{d}s=2\sum _{i=1}^n{\int }_0^{\mathcal{T}}{K}_h(t-s)\hat{\tilde{S}}{}_T^R(s)\mathrm{d}{\tilde{N}}_i^R(s),$$

$${\hat{\tilde{f}}}_T^{R,0,h}(t)=\frac{{\sum }_{i=1}^n{\int }_0^{\mathcal{T}}{K}_h(t-s)\hat{\tilde{S}}{}_T^R(s)\mathrm{d}{\tilde{N}}_i^R(s)}{{\sum }_{i=1}^n{\int }_0^{\mathcal{T}}{K}_h(t-s){\tilde{Y}}_i^R(s)\mathrm{d}s}.$$

$${\hat{\tilde{f}}}_T^{0,h}(t)={\hat{\tilde{f}}}_T^{R,0,h}(\mathcal{T}-t).$$

$${\mu }_j(K)={\int }_0^{\mathcal{T}}{s}^{j}K(s)\mathrm{d}s,R(K)={\int }_0^{\mathcal{T}}{K}^2(s)\mathrm{d}s.$$

$$(nh)^{1/2}\mathopen{}\mathclose{{}\left\{\widehat{\widetilde{f}}_{T}^{0,h}(t)-\widetilde{f}_{T}(t)-B_{0}(t)}\right\}\rightarrow N\mathopen{}\mathclose{{}\left\{0,\sigma_{0}^{2}(t)}\right\},$$

$$B_0(t)$$

$${\sigma }_0^2(t)$$

$$G_j(t)$$

Full text

Keywords: chain-ladder method; general insurance; granular reserving; nonparametric estimation; survival analysis.

1 Introduction

The classical run-off triangle used for the prediction of outstanding liabilities can be explained as a two-way ANOVA arrangement, where data is organized on a two-dimensional plane of (cohort, age) with cohort being the accident date or underwriting date of a claim, and age the time from that date to a payment. Developed at least in the beginning of the last century, the chain-ladder method is still the industry standard for estimating the future cost of outstanding liabilities from these run-off triangles. However, as a deterministic algorithm, chain-ladder does not specify the assumptions that it is based on, nor the uncertainty of the estimation.

Stochastic models around the chain-ladder method are the Mack Model (Mack, 1993) and multiplicative models in Kremer (1982); Verrall (1991); Renshaw & Verrall (1998) and Kuang et al. (2009) among many others. A comprehensive review is given in England & Verrall (2002). The drawback of these papers is that they do not discuss how the data arises as aggregation from individual data. This is needed when one wants to understand the underlying assumptions of the model. Taylor (1986) coined the term macro-models to describe these previous models and defined models that begin on an individual level as micro-models. Macro models have assumptions which are hard to justify once a data generating process with individual payments is considered. The assumptions of the most widely used Mack model can hardly be justified if one considers that the cells in the classical run-off triangle are aggregations of individual payments. Under Mack’s assumptions, not a single future payment can be independent of the past. This is because the conditional expectation of the next cell within a row of the run-off triangle is a multiple of all previous observations in the same row. The other big class of models is those of Kremer (1982); Verrall (1991); Renshaw & Verrall (1998) and Kuang et al. (2009). They assume that the expected claim amount in one cell is the product of a row factor and a column factor — representing underwriting/accident date and payment delay, respectively. This multiplicative structure implies that there is no interaction effect between rows and columns working on the expected claim amount. In Hiabu (2017) it has been shown that this non-interaction assumption generally does not hold because the cells in the run-off triangle are aggregated as parallelograms as illustrated in Figure 1. These parallelograms will generally introduce interdependencies, which violate the multiplicative structure assumption leading to an interaction effect. Hence, assuming a multiplicative structure produces a bias that grows with the level of aggregation. Therefore, as done in this paper, consistency of payment predictions can only hold in a continuous framework where the level of aggregation is understood to converge to zero with increasing number of observations.

Recent literature connects the chain-ladder method and its data to counting process theory in survival analysis. Hiabu et al. (2016) introduced a statistical model including the data generating process which is built on the continuous model of Martínez-Miranda et al. (2013). The sampling technique of the chain-ladder method is different from other sampling techniques used in classical (bio-)statistical literature. Individuals or policies are only followed if a failure, i.e., a claim occurs. This has the advantage that less data is required than in classical survival data. Truncation occurs when cohort plus age is greater than the date of data collection. However, Martínez-Miranda et al. (2013) and Hiabu et al. (2016) only considered claim counts, and ignored its associated payments.

In this paper, we introduce a micro-model in continuous time in which chain-ladder development factors, applied on a paid triangle, are a histogram estimator of a cost-weighted hazard function running in reversed development time. We establish new assumptions under which consistency of the development factors is achieved. Consistency is understood in the sense that both the number of observed claims grows to infinity and the level of aggregation tends to zero. Finally, we improve on chain-ladder estimation by replacing the histogram estimator with kernel smoothers and by estimating a cost-weighted density instead of a cost-weighted hazard.

There is also a growing literature on micro-models for estimating outstanding liabilities in non-life insurance that is not based on the chain-ladder idea. Arjas (1989) and Norberg (1993) formulated models in a classical bio-statistical setup with a non-homogeneous marked Poisson process. A strong case study in this setting has been developed in Antonio & Plat (2014) and the models have been further studied in Huang et al. (2015) and Huang et al. (2016). These models are more complex than chain-ladder models. They model each delay component in the claims process separately and require full inference on the marked point process, for instance the distribution of the mark/cost. They also require additional information about the exposure, i.e., information about the number of policies underwritten. The assumptions of this paper can be used to decide whether this additional complexity is beneficial noting that additional complexity introduces bias and is only advisable if a significantly better fit can be obtained. Additional complexity might also be necessary if claims with different accident dates, e.g., due to calendar time effects, are not independent, as investigated in Shi et al. (2012); Merz et al. (2013); Lee et al. (2015); Badescu et al. (2016); Avanzi et al. (2016); Lee et al. (2017) and Crevecoeur et al. (2019). If more complexity is justified, the estimator presented in this paper can be used as a building block in those and other more complex models. However, this is beyond the scope of this present paper.

This paper is structured as follows. Section 2 describes the mathematical model and Section 3 links chain-ladder development factors to that framework by identifying them as histogram estimator of a hazard function. Section 4 proposes improvements on chain-ladder development factors by replacing the histogram with kernel smoothers and by estimating a density function instead of a hazard function. We provide a data application and simulation studies in Sections 5 and 6. All proofs can be found in the appendix.

2 Mathematical framework

We start by putting the unique sampling scheme of chain-ladder into a micro-structure framework. We observe counting processes $(N_{i}(t))_{t\in[0,\mathcal{T}]}$, $\mathcal{T}>0$, for claims $i=1,\dots n$ and call $t$ development time. Each counting process starts with value zero at the underwriting date underlying its claim. It jumps, with jump-size one, whenever a payment is made. Additionally to every jump, we observe a mark indicating the size of the payment made. The number of counting processes, $n$, varies over calendar-time: We follow retrospectively only those claims for which at least one payment has been observed, i.e., we do not follow every claim in the policy book. In this paper, we make the following assumptions.

1. [M1]

All claims are independent. 2. [M2]

Every claim consists of only one payment.

Assumptions [M1] and [M2] are rather strong but are made to simplify the mathematical derivations yielding a first and clean step towards a better understanding of chain-ladder on a micro-structure level. Possible ways to relax these assumptions are weak dependency (instead of [M1]) and a Markov process structure where every jump triggers a new state (instead of [M2]). This, however, is beyond the scope of the present paper.

The jump-time in development direction corresponding to the payment for claim $i$ is denoted by $T_{i}$. Thus, we get

[TABLE]

where $I(\cdot)$ denotes the indicator function. As pointed out in Hiabu et al. (2016), statistical inference on the counting process $N_{i}$ is not directly feasible. We only follow a claim once we have observed at least one payment. Therefore, by design it holds

[TABLE]

where $U_{i}$ is the underwriting date or accident date of claim $i$. Hence, by not following every policy, we are exposed to a right-truncation problem instead of a right-censoring problem. In the sequel, for notational convenience, we parameterize the dates such that $today=\mathcal{T}$ which yields $T_{i}\leq\mathcal{T}-U_{i}$.

A solution to the right-truncation problem is to reverse the time of the counting process leading to a tractable left-truncation problem (Ware & DeMets, 1976). To this end we consider the counting processes

[TABLE]

each with respect to the filtration

[TABLE]

satisfying the usual conditions (Andersen et al., 1993, p. 60), and where $\mathcal{N}$ is the set of all zero probability events. It is well known (Andersen et al., 1993, Theorem II.4.2), that the intensity process of $N_{i}^{R}$ is

[TABLE]

where

[TABLE]

which is a product of a deterministic function and a predictable function. This structure is called Aalen’s multiplicative intensity model (Aalen, 1978), and enables nonparametric estimation and inference on the deterministic factor $\alpha^{R}(t)$, which is done in Hiabu (2017).

Let $Z_{i}$ denote the payment size of claim $i$ and consider the process $\widetilde{N}^{R}_{i}(t)=Z_{i}N_{i}^{R}(t)$. Ignoring for now the necessary regularity conditions, it is straight forward to see that

[TABLE]

which asymptotically satisfies Aalen’s multiplicative intensity structure, if $\mathcal{R}_{n}(t)$ converges to 1. This convergence will be verified below and it is sufficient to apply the well developed techniques for counting processes, which we do in this paper. In the next section we will show that chain-ladder development factors (which are defined for instance in Taylor (1986)) are a nonparametric histogram estimator of $1+\widetilde{\alpha}_{\ast}^{R}(t)$.

When the goal is to predict outstanding liabilities, one is interested in the untruncated versions of the truncated observations. We will indicate these variables by suppressing the subscript, i.e., three-dimensional random variable $(T,U,Z)$ has the same distribution as $(T_{i},U_{i},Z_{i})$, for every $i=1,\dots,n$, if conditioned on the event $\{T\leq\mathcal{T}-U\}$. We make the following assumptions on the untruncated objects.

1. [M3]

The random variables $T$ and $U$ have, respectively, strictly positive continuous density functions $f_{T}$ with support $[0,\mathcal{T}]$ and $f_{U}$ with support $[0,\mathcal{U}]$, $\mathcal{U}\leq\mathcal{T}$, each with respect to the Lebesgue measure. Moreover, the continuous joint density $g$ of $(T,U,Z)$ with respect to the Lebesgue measure exists and $E[|Z|]<\infty$ 2. [CLM1]

The random variables $T$ and $U$ are independent. 3. [CLM2]

There exist functions $m_{1},m_{2}$ such that $E[Z|T,U]=m_{1}(T)m_{2}(U)$.

Assumption [M3] ensures that the intensity $\widetilde{\lambda}_{i}^{R}$ is well defined. Note that [CLM1] is a statement about the untruncated objects and does not imply that $U_{i}$ and $T^{R}_{i}$ are independent, noting that $U_{i}\leq T^{R}_{i}$. The second part of Assumption [CLM2] means that there is no interaction effect between development time and underwriting date on the claim amount.

To align the density of $T$ with our setting, we define the cost-weighted density of $T$ as

[TABLE]

Conditional expectations to point events with probability zero here and below are well defined through the continuous density $f_{T}$. Analogously, we define the cost-weighted density in reversed time as $\widetilde{f}_{T}^{R}(t)=\widetilde{f}_{T}(\mathcal{T}-t)$ and $f_{T}^{R}(t)=f_{T}(\mathcal{T}-t)$. Moreover, the underlying hazard rate in reversed time can be derived from the above definition as

[TABLE]

Proposition 2.1

Given [M1]–[M3], for $n\to\infty$, it holds

[TABLE]

If additionally Assumptions [CLM1] and [CLM2] hold, then

[TABLE]

Proof 2.2**.**

See Appendix A.1.

Proposition 2.1 explains why we gave the last two assumptions a ‘CLM’-prefix. The convergence in (1) ensures that Aalen’s multiplicative intensity model is approximately satisfied and chain-ladder development factors are histogram estimates of $1+\widetilde{\alpha}_{\ast}^{R}$. But it is only via equation (2) that chain-ladder development factors also approximate $1+\widetilde{\alpha}^{R}$. For the latter to be true, we assume [CLM1] and [CLM2]. Under these two additional conditions, the chain-ladder algorithm predicts the right object and leads to a sensible quantification of the outstanding liabilities.

We close this section with some further remarks.

Remark 2.3** (Remark 1).**

(Exposure). As in traditional chain-ladder and in contrast to classical survival data, because all failures (claims) are observed, there is no additional information needed about the number of individuals under risk (i.e. exposure in form of the number of underwritten policies) in order to estimate future claim amounts. The unique sampling leads to a right-truncation which is solved by reversing the time. This is different to the approaches described in Arjas (1989) and Norberg (1993).

Remark 2.4** (Remark 2).**

(Cost-weighted density). The density $\widetilde{f}_{T}$ is the continuous analogue to the column parameter of chain-ladder, which is often called $\beta$ and which is considered in Kremer (1982); Verrall (1991); Renshaw & Verrall (1998) and Kuang et al. (2009). Moreover, integrating to one, $\widetilde{f}_{T}$ is indeed a density function.

Remark 2.5** (Remark 3).**

(Predicting outstanding liabilities). An estimator of the cost-weighted hazard, $\widetilde{\alpha}^{R}(t)$, in conjunction with a chain-ladder algorithm can be used to predict outstanding liabilities. Alternatively, as proposed in this paper, one can employ estimators of the cost-weighted densities $\widetilde{f}_{T},\widetilde{f}_{U}$. If the maximum development time of a claim is $\mathcal{T}$, then the expected outstanding liabilities for claims underwritten in $[0,\mathcal{U}]$, $\mathcal{U}\leq\mathcal{T}$, is given as

[TABLE]

where \widetilde{f}_{T,U}(t,u)=E[Z]^{-1}E\mathopen{}\mathclose{{}\left[Z|\ T=t,U=u}\right]f_{T,U}(t,u) is the cost-weighted density of $(T,U)$. The total amount of payments until today is given by $\sum_{i=1}^{n}Z_{i}$ and the fraction in $r_{n}$ gives the expected ratio between outstanding payments and past payments. Note that under Assumptions [CLM1] and [CLM2], the cost-weighted joint density factorizes into $\widetilde{f}_{T,U}(t,u)=\widetilde{f}_{T}(t)\widetilde{f}_{U}(u)$. In Section 4 we propose estimators for $\widetilde{f}_{T}$. Due to symmetry, the component $\widetilde{f}_{U}$ can be estimated by swapping the roles of $T$ and $U$. Outstanding liabilities are estimated by replacing $\widetilde{f}_{T}$ and $\widetilde{f}_{U}$ in $r_{n}$ with their estimates. Developing estimation theory for $r_{n}$ is rather involved because of the non-trivial integrals and the ratio-structure in (3). We only consider a simulation study for the prediction performance of $r_{n}$ in Section 6.

Remark 2.6** (Remark 4).**

(Assumptions). While the model is built around the observation of independent claims with one single payment each (Assumptions [M1] and [M2]), allowing for claim clusters and multiple payments per claim is feasible and would only require some further assumptions.

Assumption [CLM1] is analogue to the usual multiplicity assumption found for example in Kremer (1982); Verrall (1991); Renshaw & Verrall (1998) and Kuang et al. (2009). The difference is that [CLM1] refers to claim counts and not to claim amounts. However, [CLM1] and [CLM2] together imply the multiplicity of aggregated expected claim amounts as assumed in the literature — if ignoring the potential bias arising from aggregation. Chain-ladder development factors can be biased if the cost-weighed development delay with density $\widetilde{f}_{T}$ is neither exponentially distributed nor uniformly distributed within each development period (Hiabu, 2017).

3 Chain-ladder development factors

We now discuss how hazard rates can be estimated in the framework of Section 2. In the setting of Proposition 2.1, the intensity of $\widetilde{N}_{i}^{R}$ at $t$ is asymptotically equal to $\widetilde{\alpha}^{R}(t)\widetilde{Y}_{i}^{R}(t)$. We use this fact to construct a least squares criterion to estimate $\widetilde{\alpha}^{R}(t)$. Given a smoothing parameter, $h>0$, and a weight function $W_{h}(\cdot,\cdot)$, we look for estimators $\widehat{\widetilde{\alpha}}^{R}(t)$ that minimize

[TABLE]

where the expression $1/{\widetilde{Y}_{i}^{R}(s)}$ is understood as being zero whenever $\widetilde{Y}_{i}^{R}(s)$ is zero. The term $\xi(\varepsilon)=\{\varepsilon^{-1}\int_{s}^{s+\varepsilon}\mathrm{d}\widetilde{N}^{R}_{i}(w)\}^{-2}$ is a vertical shift subtracted to make the integral well-defined. Since $\xi(\varepsilon)$ does not depend on $\widehat{\widetilde{\alpha}}^{R}(t)$, the estimator $\widehat{\widetilde{\alpha}}^{R}(t)$ is defined by a local weighted least squares criterion. We understand the integral with respect to $\mathrm{d}\widetilde{N}_{i}$ as a Stieltjes integral.

Let $0=t_{1}<...<t_{m}=T$ be an equidistant partition of the interval $[0,\mathcal{T}]$ with bin-width $h$ and some integer $m$. For $t\in[t_{l},t_{l+1})$, we set

[TABLE]

The first order condition minimizing (4) under the weighting $W^{H}_{h}(s,t)$ leads to the histogram estimator

[TABLE]

Analogue to Hiabu (2017) it can be shown that, up to lower order terms, chain-ladder development factors equal $1+\widehat{\widetilde{\alpha}}^{R}_{H,h}(l)$. Therefore, the following proposition can also be interpreted as a central limit theorem for development factors. We make the following assumptions.

• [S1]

The bandwidth $h=h(n)$ satisfies $h\rightarrow 0$ and $n^{1/4}h\rightarrow\infty$ for $n\rightarrow\infty$.

• [S2]

The density $f_{T}$ is two times continuously differentiable.

• [S3]

The function $l(t)=E[Z_{1}|\ Y_{1}^{R}(t)=1]$ is continuously differentiable.

We define the following quantity.

[TABLE]

Proposition 3.1**.**

Under Assumptions [M1]–[M3], [CLM1], [CLM2], and [S1]–[S3], for $t\in(t_{l},t_{l+1})$, $n\rightarrow\infty$, it holds

[TABLE]

in distribution, where

[TABLE]

Proof 3.2**.**

See Appendix A.5.

Hence, apart from the usual regularity condition, under [CLM1] and [CLM2], consistency of the development factors is achieved if both the number of observations, $n$, goes to infinity and the level of aggregation, $h$, tends to zero.

In the next section, we propose two improvements to chain-ladder development factors. Firstly, we replace the histogram weighting, $W_{h}(\cdot,\cdot)$, with local polynomial kernel smoothers leading to a reduced bias. Secondly, we will work with the density function instead of the hazard function because we expect estimation of the density to be more robust. This is because the hazard function, due to the bounded support, usually increases heavily at the right boundary whereas the shape of the density is less explosive. A simulation study in Bischofberger et al. (2019) confirms this heuristic.

4 Local polynomial density estimation

In this section we introduce two nonparametric estimators of the one-dimensional cost-weighted density $\widetilde{f}_{T}$: the local constant estimator and the local linear estimator. The idea of local polynomial fitting is quite old and might originate from early time series analysis (Macaulay, 1931). It has been adapted to the regression case in Stone (1977) and Cleveland (1979). A general overview of local polynomial fitting can be found in Fan & Gijbels (1996).

Note that $\widetilde{f}_{U}$ can be estimated analogously by inverting the roles of $T$ and $U$ and adapting the definitions of $N_{i}$, $Y_{i}$ etc. The joint cost-weighted density $\widetilde{f}_{T,U}$ is then estimated by $\widehat{\widetilde{f}}_{T}\times\widehat{\widetilde{f}}_{U}$ in line with Remark 3.

We first define the cost-weighted Kaplan-Meier product-limit estimator of the survival function \widetilde{S}_{T}^{R}(t)=\int_{t}^{\infty}\widetilde{f}^{R}_{T}(s)\mathrm{d}s=\{E\mathopen{}\mathclose{{}\left[Z|\ T^{R}\geq t}\right]/E[Z]\}\int_{t}^{\infty}f^{R}_{T}(s)\mathrm{d}s as

[TABLE]

where \widehat{\widetilde{A}}{}^{R}(t)=\sum_{i=1}^{n}\int_{0}^{t}Z_{i}\mathopen{}\mathclose{{}\left\{\sum_{j=1}^{n}Z_{j}Y_{j}^{R}(s)}\right\}^{-1}\mathrm{d}N^{R}_{i}(s) is motivated by the Aalen estimator, estimating \widetilde{A}^{R}(t)=\int_{0}^{t}E\mathopen{}\mathclose{{}\left[Z|\ T^{R}=s}\right]\{E\mathopen{}\mathclose{{}\left[Z|\ T^{R}\geq s}\right]\}^{-1}\alpha^{R}(s)\mathrm{d}s. Here the product can be understood as simple finite product because of the finite number of jump points of $\widehat{\widetilde{A}}{}^{R}(t)$ as explained in (Andersen et al., 1993, p. 89). Let $q_{p}(z)=\sum_{i=0}^{p}\theta_{i}z^{i}$ denote a polynomial of degree $p$. For $t\in[0,\mathcal{T}]$, we define the local polynomial estimator of degree $p$, $\widehat{\widetilde{f}}_{T}^{R,p,h}(t)$ of $\widetilde{f}^{R}_{T}(t)$ as the minimizer $\widehat{\theta}_{0}$ in the equation

[TABLE]

For a kernel $K$ and bandwidth $h>0$, we set $K_{h}(t)=h^{-1}K(t/h)$ as usual. The expression $\xi_{f}(\varepsilon)=\{\varepsilon^{-1}\int_{s}^{s+\varepsilon}\widehat{\widetilde{S}}{}_{T}^{R}(w)\mathrm{d}\widetilde{N}^{R}_{i}(w)\}^{-2}$ is needed to make (6) well defined. Since $\xi_{f}(\varepsilon)$ does not depend on $q_{p}$, $\widehat{\theta}_{0}$ is defined by a local weighted least squares criterion.

In the sequel we will only consider the cases $p=0,1$, i.e., the local constant and local linear case. While a higher degree in conjunction with higher order kernels improves the asymptotic properties, finite sample studies show that improvements are only visible with unrealistically big sample sizes. In the local constant case of (6) we derive the first order condition

[TABLE]

and conclude the local constant estimator

[TABLE]

The final estimator in non-reversed time is then simply defined as

[TABLE]

We add the following assumption.

• [S4]

*The kernel $K$ is symmetric, has bounded support and has finite second moment, and it holds $\int K(u)\mathrm{d}u=1$. *

Other kernels can also be used but they will require a more complex estimator. Moreover, we introduce the following notation. For every kernel $K$ and $j\geq 0$, let

[TABLE]

Proposition 4.1**.**

Under Assumptions [M1]–[M3], [CLM1], [CLM2], and [S1]–[S4], for $t\in(0,\mathcal{T})$, $n\rightarrow\infty$, it holds

[TABLE]

in distribution, where

[TABLE]

Proof 4.2**.**

See Appendix A.3.

For the local linear case, we introduce the following quantities. For $t\in[0,\mathcal{T}]$, set

[TABLE]

The first order condition for $p=1$ then reads

[TABLE]

Hence, the solution $\widehat{\theta}_{0}$ is given by

[TABLE]

where

[TABLE]

If $K$ is a second-order kernel, then

[TABLE]

so that $\overline{K}_{t,h}$ can be interpreted as a second-order kernel with respect to the measure $\mu$, which is defined via $\mathrm{d}\mu(s)=n^{-1}\sum_{i=1}^{n}\widetilde{Y}_{i}^{R}(s)\mathrm{d}s$.

The local linear estimator in non-reversed time is defined as

[TABLE]

Proposition 4.3**.**

Under Assumptions [M1]–[M3], [CLM1], [CLM2], and [S1]–[S4], for $t\in(0,\mathcal{T})$, $n\rightarrow\infty$, it holds

[TABLE]

in distribution, where

[TABLE]

for $R(K)=\int K^{2}(s)\mathrm{d}s$.

Proof 4.4**.**

See Appendix A.4.

One alternative to estimate the cost-weighted density is to use a semiparametric asymmetric kernel density estimator which better accounts for the tail (Gustafsson et al., 2009). We chose not to do so in this paper, since a nonparametric estimation technique is more in the spirit of the chain-ladder technique as explained in the previous section.

5 Data Application: Estimating outstanding liabilities

We apply our estimator on a data set from a motor insurance in Cyprus which was collected between 2004 and 2013. The data contains $n=\,$ 51,216 closed claims $(T_{i},U_{i},Z_{i})$, $i=1,\dots,n$ consisting of their payment delay until the final payment $T_{i}$, their accident dates $U_{i}$, and the total claim amount $Z_{i}\geq 0$. First, we estimate the marginal cost-weighted densities $\widetilde{f}^{T}$ and $\widetilde{f}^{U}$ of $T$ and $U$, respectively, and in particular we forecast the outstanding claim amount $r_{n}$ consisting of all claims for accidents that have already incurred but have not been paid yet (see Remark 3 in Section 2). Afterwards, we illustrate our model assumptions [CLM1] and [CLM2] on the data set.

5.1 Estimation and forecasting

For the estimation of outstanding liabilities, we calculate the components $\widehat{\widetilde{f}}^{0,h}=\widehat{\widetilde{f}}_{T}^{0,h}\widehat{\widetilde{f}}_{U}^{0,h}$ and $\widehat{\widetilde{f}}^{1,h}=\widehat{\widetilde{f}}_{T}^{1,h}\widehat{\widetilde{f}}_{U}^{1,h}$ using the Epanechnikov kernel $K(s)=0.75(1-s^{2})I(|s|\leq 1)$. For data-driven bandwidth selection, we use cross-validation (Rudemo, 1982; Hall, 1983; Bowman, 1984). The score function $Q_{T}^{j}$ is motivated by the minimization problem which lead to the local polynomial estimators introduced in Chapter 4. For the estimation of $\widetilde{f}_{T}^{R}$ we want to minimize $\sum_{i=1}^{n}\int_{0}^{1}\{\widehat{\widetilde{f}}_{T}^{R,j,h}(t)-\widetilde{f}_{T}^{R}(t)\}^{2}\widetilde{Y}^{R}_{i}(t)\mathrm{d}t$ in $h$ for $j=0,1$. Since $\widetilde{f}_{T}$ is unknown, we select the bandwidth $h_{j}^{T}$ as the minimizer of

[TABLE]

in $h$ instead (Nielsen et al., 2009). The “leave-one-out” terms are given as

[TABLE]

While this bandwidth selection works well for the estimators of the weighted density of the accident date $U$, we get unrealistic estimates for $\widetilde{f}_{T}$. We decided to adjust the bandwidth manually and calculated $\widehat{\widetilde{f}}_{T}^{j,h}$ for a small bandwidth $h_{T}^{1}=10$ days for delays shorter than $1.5$ years (= 548 days) and we used a large bandwidth $h_{T}^{1}=500$ days to estimate $\widetilde{f}_{T}$ for $t>548$ days. The optimal bandwidths for $\widetilde{f}_{U}$ by cross-validation are $h_{U}^{1}=471$ and $h_{U}^{0}=280$ days, respectively. We remark that a full investigation of local bandwidth selection is beyond the scope of this paper.

The results are given in Figure 2. Since most claims were paid off after 1.5 years, our density estimators for $\widetilde{f}_{T}$ are almost zero for $t>1.5$ years. Big outliers in that area are oversmoothed, which reflects the possibility of large payments with high delays better than a small number of sharp local maxima of the density at the positions of the outliers and a density of 0 elsewhere.

For cost-weighted density estimators $\widehat{\widetilde{f}}_{T}$ and $\widehat{\widetilde{f}}_{U}$, we estimate the reserve by

[TABLE]

The reserve estimate $\hat{r}_{n}$ is motivated by the representation of the reserve $r_{n}$ in equation (3) in Remark 3. Estimates for outstanding claim payments per future year, per accident year, and in total are given in Table 1. We compare the estimators with local bandwidth correction with the results obtained through the classical chain-ladder method with quarterly aggregated data. Whereas all three total reserve forecasts are very similar, one can see differences for very short and very large delays. Furthermore, it is striking that both smoothed estimators forecast a non-zero claim amount for 2023 but chain-ladder estimates it to be 0. The difference between chain-ladder and smoothed density estimators for very short and mainly for large delays has already been observed for non-cost-weighted estimators in Hiabu et al. (2016). Moreover, as explained in Hiabu (2017), chain-ladder tends to overestimate the total reserve whereas the estimate from the local linear estimator is asymptotically unbiased. The local constant estimator is known to suffer from bias at boundaries, i.e., weaker performance than the local linear one for very short and very large delays (Fan & Gijbels, 1996; Wand & Jones, 1994). The undersmoothed densities estimators with bandwidths obtained from cross-validation yield similar estimates for the reserve although the shape of the density estimates is unrealistically rough for larger delays (13,030,459 in the local linear case and 13,268,768 in the local constant case).

We are aware that these forecasts are just point estimates for the reserve. We investigate variation in the forecast under a controlled setting in the simulation study in the next chapter.

5.2 Illustration of assumptions

For Assumption [CLM1], the independence between $T$ and $U$ could not be assured by an independence test based on Conditional Kendall’s tau for truncated data (Austin & Betensky, 2014; Martin & Betensky, 2005). To get more insight we aim to visualize the underlying dependency. We aggregate the data into three-month bins $q_{1},\dots,q_{40}$. Then we introduce a triangle with aggregated observations $\mathcal{N}^{Q}_{r,s}=\sum_{i=1}^{n}I(U_{i}\in q_{r},T_{i}\in q_{s})$, $r,s,=1,\dots,40$, $r+s\leq 40$, and calculate the development factors $\alpha(r,s)=\sum_{l=1}^{s+1}\mathcal{N}^{Q}_{r,l}/\sum_{l=1}^{s}\mathcal{N}^{Q}_{r,l}$, for development quarter $s$ and accident quarter date $r$. The values of $\alpha$ for the first six development quarters are given in Figure 3a. Under the assumption of independence between $T$ and $U$ (Assumption [CLM1]), the function $\alpha$ is independent of the accident date $r$ and hence each plot should show points scattered around a horizontal line yielding a flat regression line.

The $p$-values for the linear regression slope parameters were only significant at 5%-level in the first two quarters. For comparison, Figure 3b shows the development factors on independently simulated variables for which linear trends were insignificant at 10%-level in every quarter. Except for the first plot in Figure 3a, indeed none of the plots and linear fits from Figure 3a are visually distinguishable from any plot in Figure 3b. Clearly, this is not a sound method to prove independence but it illustrates that the dependence in our data might come from the first quarter only. Advisable extensions of our model, that handle possible dependence in our data, are e.g. seasonal effects as considered in Lee et al. (2015) or operational time (Lee et al., 2017). We do not consider these approaches in this paper.

For Assumption [CLM2], we have to verify that the expected cost conditioned on $T$ and $U$ is multiplicatively separable, i.e., that there exist functions $m_{1}$, $m_{2}$ such that $E[Z|T,U]=m_{1}(T)m_{2}(U)$. Similarly to the above, we use a visual approach on aggregated data to illustrate this setting in the first six quarters $q_{1},\dots,q_{6}$ of the accident years. Assumption [CLM2] is satisfied if for all observations $i$ with $U_{i}$ in quarter $q_{k}$ it holds $Z_{i}=c_{k}m_{1}(T_{i})+\varepsilon_{i}$ for a quarter dependent constant $c_{k}>0$, and a mean-zero error $\varepsilon_{i}$. Figure 4 shows the claim cost given the delay for claims in the first six quarters. Under Assumption [CLM2], the points in each plot should be generated by the same regression function after normalizing with the accident date quarter dependent factor $c_{k}$. We use a linear interpolation to compare the structure of the observations. All but the third plot show very similar development of claim costs. In the third plot, the claim costs increase much faster due to some outliers that are not visible in the plot. For comparison, we generate 50,000 observations from Scenario 5 in Section 6 where $E[Z|T,U]=m_{1}(T)m_{2}(U)$ with $m_{1}(t)=(t+0.75)$, $m_{2}(u)=(u-0.25)^{2}+1$.

We conclude that while the data does not fully follow our assumptions, it is suitable enough for the illustration purpose of this paper.

6 Simulation study

This chapter shows the performance of our new estimators on simulated data. Our first finding is that the local linear estimator outperforms the local constant one at boundaries (Section 6.1). Secondly, when estimating the reserve, our local constant estimator is best for small sample sizes and the local linear one for large sample sizes whereas the performance of chain-ladder varies (Section 6.2). Last, in a micro model in Section 6.3, we see that reliable monthly forecasts can only be obtained with our density estimators and not with chain-ladder.

6.1 Weighted density estimation

We perform a simulation study to show the performance of our estimators for a selection of distributions if the optimal bandwidth is chosen. We simulate truncated observations $(T_{i},U_{i},Z_{i})$, $i=1,\dots,n$ on $\mathcal{I}=\{(t,u):0\leq u,t\leq 200,u+t\leq 200\}$. With the true weighted densities $\widetilde{f}_{T}$ and $\widetilde{f}_{U}$ being known, we calculate the local constant and local linear estimators $\widehat{\widetilde{f}}^{j,h}_{A}$, $j=0,1$, with the best bandwidth $h^{j}_{A}$ with respect to the integrated squared error

[TABLE]

for $A\in\{T,U\}$. We choose eight different settings for the distributions of $T$, $U$, and $Z$. The choice of the distributions is motivated by empirical distributions on the one hand and challenging estimation settings for the distribution of $T$ and $U$ are added on the other hand. The observations of $T$ and $U$ are simulated independently and truncated on $[0,1]$. The probability density functions for $T$ and $U$ are shown in Figure 5a–d and the values of $Z$ given one choice of $(T,U)$ are illustrated as histograms in Figure 5e and f. For simulated conditional claim costs $Z$ given $T$ and $U$, we take gamma distributions with shape parameter $k=1$ and different scale parameters $\theta=(T+0.75)((U-0.25)^{2}+1)$ and $\theta=T(U^{2}-U+1)$. Note that Assumption [CLM2] holds because of the identity $E[Z|T,U]=k\theta$ of the gamma distribution.

We take all combinations of these distributions and label the eight scenarios as given in Table 2. For each scenario 1000 random samples of sizes 100, 1000, 10,000 and 100,000 are generated.

We investigate 32 cases arising from eight different scenarios and four different sample sizes. The exact results are omitted here. We only give our main conclusion and focus on the reserve estimate in more detail in the next section. The local constant and the local linear estimators perform similarly in terms of empirical mean integrated squared error (eMISE) with the local linear estimators being more stable. In 26 out of 32 cases, the eMISE of $\hat{f}_{1,h_{1}^{X},K}^{X}$ is more than 25% lower than the of $\hat{f}_{0,h_{0}^{X},K}^{X}$. In 4 cases it even improves the eMISE by more than 75%. On the other hand, there are only two cases where the local linear estimator leads to an increase in the eMISE by 0.7% and 6.7%, respectively. In the other covariate, both density estimators perform equally well. However, the local linear estimator performs better in scenarios with the boundary challenge distribution for $U$. This reflects aforementioned weakness of the local linear kernel density estimator close to boundary regions. The difference is biggest in Scenarios 7 and 8 where the local linear estimator is able to make up for the lack of observations in the corner.

6.2 Estimates for stimulated outstanding liabilities

Next, we compare the reserve estimates

[TABLE]

with the true outstanding claim amount $r_{n}$ defined in equation (3).

Table 3 contains the mean, standard deviation and the median of the errors in the estimation of the squared relative errors

[TABLE]

The results are compared to the estimation through chain-ladder applied on the triangle arising from the aggregation of the simulated observations of $T$ and $U$ into 20 bins each. This aggregation is comparable to quarterly aggregation on real data.

First, we want to note that there was a complete breakdown of the chain-ladder algorithm for too small numbers of observations which resulted in an invalid estimate in our implementation. Moreover, in most cases of the simulation study, our local polynomial density estimators outperform chain-ladder. The reserve estimates from the local linear estimators were strikingly better in the boundary challenge Scenarios 3, 4, 7 and 8 for numbers of observations larger than $n=1000$. For $n=100$ the local constant reserve estimate was best in six out of eight scenarios. In boundary challenge scenarios, chain-ladder not only lead to invalid results for small sample sizes but it also resulted in extreme outliers. An illustration of the results is given in Figure 6. It shows a scenario in which the local linear density estimator is the only one that estimates the altitude of the maximum in the joint density almost correctly.

We conclude that the local linear estimator performs best for $n\geq 1000$ and that the local constant one does for smaller sample sizes. Detailed results can be found in Table 3.

6.3 Simulation of a micro model

In this section, we investigate the performance of our estimators and chain-ladder on simulated data arising from a micro model. We simulate different steps in the underwriting and payment process separately and for different types of policies in one line of business and estimate outstanding payments under different circumstances. We then compare estimates of the reserve with actual future payments.

To create our data set, we follow the “central scenario” simulation in Baudry & Robert (2019). We generate mobile phone insurance policies that were underwritten over two years and estimate the outstanding liabilities at different times throughout the underwriting period and some months after the last policy was sold. We assume that the insurance provider covers damage in the three events breakage, oxidation, and theft. For this purpose the three policy types “breakage”, “breakage and oxidation”, and “breakage, oxidation, and theft” are underwritten with probabilities 0.25, 0.45, and 0.30. Moreover, there are four different mobile phone brands with four different models each specifying the price of the phone. The frequencies of brands and models and their basis prices and model prize factors are given in Table 4.

Following Baudry & Robert (2019), we simulate insurance policies that are underwritten independently between the first day of 2016 and the last day of 2017. The number of underwritten policies per day follows a Poisson point process with constant intensity $\lambda_{0}(t)=700$, independently of policy type, phone brand, or model. Each policy covers exactly the period of the next 360 days after the underwriting day. For claims, we simulate the three incidents through a competing risk model with the constant hazards in Table 5. All events are recorded daily and we identify a year with the grid $\{1/360,2/360,\dots,1\}$. After an incident has happened at time $T_{0}$, the reporting time is generated from the reporting delay hazard

[TABLE]

for $a=0.4$, $b=10$. Hence, we assume a maximum reporting delay of one year. Denoting the reporting day by $T_{1}$, the payment day is then generated from the payment delay hazard

[TABLE]

with $a=7$, $b=7$ and where $m=40/360$ and $d=10/360$. Hence, all claims are settled within 10 to 50 days. Note that both delays are independent of the incident or underwriting day $T_{0}$ and $T_{1}$ and, thus, Assumptions [M3] and [CLM1] are satisfied. Moreover, reporting delay $\alpha_{0}$ and payment delay $\alpha_{1}$ are both independent of phone brand, model, type of policy, and type of incident. Last, we assume that the whole claim is settled in a single payment (Assumption [M2]) which is a random proportion of the phone price following a beta distribution with parameters given in Table 5.

We estimate the outstanding payments at each month from September 2016 to May 2018 with our new estimators and with monthly aggregated chain-ladder and compare it with the simulated future payments. The whole scenario is then repeated 200 times. The average reserves over all 200 simulation runs and their empirical 95% confidence intervals are given in Figure 7a. It shows an increase in payments until February 2017, with new policies being underwritten every day. The payments stabilize afterwards in a balance between new policies and their claims and old policies expiring after 360 days. After December 2017, there is a decrease in payments since no new policies are underwritten anymore after 2017 and remaining policies expire. The medians of the actual outstanding future payments are taken over all 200 simulation runs and labeled as “true” reserves. Moreover, the mean squared error

[TABLE]

for reserve estimate $\hat{r}_{n}^{[k]}$ and true reserve $r_{n}^{[k]}$ and simulation runs $k=1,\dots,200$ is given in Figure 7b. The reserve estimates from our local linear estimators $\widehat{\widetilde{f}}^{1,h}_{T},\widehat{\widetilde{f}}^{1,h}_{U}$ have the lowest bias and variance. Whereas the local constant estimator suffers from bias, chain-ladder suffers heavily from variance (as found in Baudry & Robert (2019)). The latter is due to the monthly aggregation for chain-ladder which, however, is necessary to derive the monthly cash-flow. Our proposed kernel smoothers use larger bandwidths, and thus reduce variance, while still providing a monthly cash-flow.

To compare this setting with the previous section, the marginal cost-weighted distributions of the delay from incident day until payment day and of the incident day are given in Figure 8a and b, respectively. The real distribution of the data is approximated by the average over the empirical distributions from each of the 200 simulation runs. Figure 8 illustrates how the development factors of the chain-ladder method lead to a histogram instead of a smooth kernel estimator as described in Section 3.

Acknowledgments

The authors would like to thank the editors and two anonymous referees for useful comments and suggestions which helped to improve this research article.

Appendix A Proofs

In the proofs below we will use the symbols $o_{p}(1)$ and $O_{p}(1)$ which are the probabilistic counterparts to the Landau symbols $o(1)$ and $O(1)$. A precise definition and explanation can be found in Appendix A of Pollard (2012). Further we will us the short-hand \Delta N^{R}_{i}(t)=\lim_{h\downarrow 0}N^{R}_{i}\mathopen{}\mathclose{{}\left\{(t+h)-}\right\}-N^{R}_{i}(t-).

A.1 Proof of Proposition 2.1

For the proof it suffices to show that

[TABLE]

and

[TABLE]

since $N_{i}$, $Y_{i}$, $Z_{i}$, $i=1,\dots,n$, are $iid$.

For the convergence of $(\ref{proof:prop1:1})$ we note that $\sup_{t\in[0,\mathcal{T}]}|n^{-1/2}\sum^{n}_{i=1}Y_{i}^{R}(t)-E[Y_{i}^{R}(t)]|=o_{p}(\log(n))$ and $\sup_{t\in[0,\mathcal{T}]}|n^{-1/2}\sum^{n}_{i=1}Y_{i}^{R}(t)Z_{i}-E[Y_{i}^{R}(t)Z_{i}]|=o_{p}(\log(n))$. Both statements follow from a strengthened Glivenko-Cantelli Theorem, since we have $Y^{R}_{i}(t)=I(U_{i}<t<T_{i}^{R})=I(U_{i}<t)+I(T_{i}<\mathcal{T}-t)-1$; see (Van der Vaart, 2000, Chapter 19.1) for more details. Next we argue that (9) is equivalent to [CLM2]. We note that

[TABLE]

Now, since $T^{R}$ and $U$ are independent, we get

[TABLE]

Hence equation (9) is equivalent to

[TABLE]

With continuity arguments this holds if and only if $\int_{0}^{\infty}zg(z,s,u)\mathrm{d}z$ is multiplicatively separable in $s$ and $u$, i.e, [CLM2] holds.

A.2 Estimation of the weighted survival function

We begin by investigating the asymptotic behaviour of

[TABLE]

instead of $\widehat{\widetilde{A}}{}^{R}(t)$. Later with Lemma A.5, we show that the difference between the two terms is uniformly of stochastic order $n^{-3/2}$ and hence negligible. One can hence carry this result over to a result for the actual estimators $\widehat{\widetilde{S}}{}^{R}$ and $\widehat{\widetilde{f}}^{0,h}$ or $\widehat{\widetilde{f}}^{1,h}$, respectively.

To start with, we analyze the process $\widehat{\widetilde{A}}{}^{R,*}_{i}(t)=\int_{0}^{t}\{\sum_{j\neq 1}Z_{j}Y^{R}_{j}(s)\}\mathrm{d}\widetilde{N}^{R}_{i}(s)$, where the integral can be understood pathwise in Lebesgue-Stieltjes sense.

We start by deriving the compensator of $\widehat{\widetilde{A}}{}^{R,*}_{i}$:

[TABLE]

Note that the error $O(n^{-2})$ comes from the Taylor expansion of $E[X^{-1}]$ at $E[X]$ for $X=\sum_{j\neq i}Z_{j}Y_{j}^{R}(t)$. Since $Y_{j}^{R}(t)\leq 1$ for all $t$ and $Z_{j}\in L_{1}(\Omega)$, this error is also uniform of order $O(n^{-2})$. Hence,

[TABLE]

is asymptotically a compensator of the uniformly integrable submartingale $\widehat{\widetilde{A}}{}^{R,*}_{i}$. We denote the resulting process by $\widetilde{M}^{R}_{i}=\widehat{\widetilde{A}}{}^{R,*}_{i}-\widetilde{\Lambda}_{i}^{R}$ which is, up to a lower order term of $O(n^{-2})$, a martingale. Since $\widetilde{M}^{R}_{i}$ is cadlag with finite variation, the quadratic variation equals the sum of square differences:

[TABLE]

Note that we used $\Delta\widetilde{M}^{R}_{i}=\Delta\widehat{\widetilde{A}}{}^{R,*}_{i}$, since $\widetilde{\Lambda}_{i}^{R}$ is continuous. As [\widetilde{M}^{R}_{i}(t)]=\mathopen{}\mathclose{{}\left(\widehat{\widetilde{A}}{}^{R,*}_{i}}\right)^{2}, by similar arguments as before we can calculate its compensator to derive the predictable variation process

[TABLE]

Proposition 4.1 is based on the following intermediate result.

Lemma A.1**.**

Under Assumptions [M1]–[M3], [CLM1]–[CLM2] and S1–S4, it holds that

[TABLE]

in distribution in Skorokhod topology sense, where $W$ is a zero mean Gaussian martingale with covariance, $\mathrm{Cov}\{W(s),W(t)\}=\sigma^{2}(s\wedge t)$.

Proof A.2**.**

This follows from a martingale central limit theorem in Rebolledo (1980) as illustrated in Andersen et al. (1993, p. 83). For the assumptions to be satisfied, we verify that

[TABLE]

where we have used that $\langle\widetilde{M}^{R}_{i},\widetilde{M}^{R}_{j}\rangle=0$ for $i\neq j$. For the Lindeberg condition we observe

[TABLE]

where we used that the jumps $t_{i}$ happen at the same time with zero probability. The condition follows from the terms in the sum being $o(n^{-2})$, since $n^{1/2}Z_{i}(\sum_{j\neq i}Z_{j}Y_{j}^{R}(s))^{-1}\rightarrow 0$.

Corollary A.3**.**

Under Assumptions [M1]–[M3], [CLM1]–[CLM2] and S1–S4, it holds that

[TABLE]

Proof A.4**.**

This follows from Lemma A.1 with Lenglart’s inequality and the functional delta method, since $\widehat{\widetilde{S}}{}^{R}$ and $\widetilde{S}^{R}$ are functionals of $\widehat{\widetilde{A}}{}^{R}$ and $\widetilde{A}^{R}$, respectively.

Indeed, it holds $\widetilde{S}^{R}(t)=\prod_{s\leq t}(1-\Delta\widetilde{A}^{R}(s))$ and $\widetilde{A}(t)={(n-1)}{n}^{-1}\sum_{i=1}^{n}\hat{\Lambda}_{i}^{R}(t)+o_{p}(1)$. The conclusion from Lenglart’s inequality (see Andersen et al. (1993)) is

[TABLE]

for every $\delta>0$ and every $\eta>0$. Hence, the fact that $n^{1/2}\langle\sum_{i=1}^{n}\widetilde{M}^{R}_{i}(1)\rangle\to\sigma^{2}(1)$ from the proof of Lemma A.1 implies $n^{1/2}\sup_{s}|\widehat{\widetilde{A}}{}^{R}(s)-\widetilde{A}^{R}(s)|=O_{p}(1)$. Therefore, we get $n^{1/2}\sup_{t}|\widehat{\widetilde{S}}{}^{R}(t)-\widetilde{S}^{R}(t)|=O_{p}(1)$ (Andersen et al., 1993, p. 86) and the conclusion follows.

A.3 Proof Proposition 4.1

We first split the estimation error into a stable part and a martingale part via

[TABLE]

where

[TABLE]

We now discuss the asymptotics of $B_{0}$ and $V_{0}$ separately, starting with $V_{0}$. Defining

[TABLE]

leads to

[TABLE]

where the error is again uniformly bounded (see Proposition 2.1). Hence, since $\overline{M}^{R}_{i}$ is a martingale, $V_{0}$ is, up to lower order terms, a martingale as well. With similar arguments as before we get

[TABLE]

Similarly as in the proof of Lemma A.1, we use the martingale central limit theorem from Rebolledo (1980) to show

[TABLE]

in distribution in Skorohod topology sense, where $W$ is a zero mean Gaussian martingale with covariance $\mathrm{Cov}\{W(s),W(t)\}=\sigma^{2}(s\wedge t)$. The assumptions are satisfied by similar arguments as in Lemma A.1. We only illustrate the derivation of $\sigma^{2}(t)$. It holds

[TABLE]

With the uniform convergences of $\widehat{\widetilde{S}}{}^{R}$, $\sum_{i=1}^{n}Y_{i}^{R}$ and $\sum_{i=1}^{n}Z_{i}Y_{i}^{R}(t-uh)$, we conclude that

[TABLE]

for $n\to\infty$, which coincides with $(\ref{V})$.

We continue with the asymptotics for $B_{0}$. After expanding $\widetilde{f}_{T}^{R,0,h}(t)$ and replacing $\widehat{\widetilde{S}}{}^{R}(s)$ by $\widetilde{S}^{R}(s)+O_{p}(n^{-1/2})$, which we can do with Corollary A.3, we have that

[TABLE]

From the remark in A.1 we can further use that $n^{-1}\sum_{i=1}^{n}Z_{i}Y_{i}^{R}(s)$ converges uniformly to $E[Z_{1}|Y^{R}_{1}=1]\gamma(s)$ and it even holds that $\sup_{s\in[0,\mathcal{T}]}|n^{-1}\sum_{i=1}^{n}\{Z_{i}Y_{i}^{R}(s)-E[Z_{1}|Y^{R}_{1}(s)=1]\gamma(s)\}|=o_{p}(n^{-1/2}\log(n))$. Hence,

[TABLE]

Note that an error term of order $O(n^{-2})$ is likewise of order $o(h^{2})$ by Assumption S1. The proof is concluded by two Taylor expansions in the numerator and one in the denominator and using that $K$ is a second order kernel.

The implication for the estimator built from $\widehat{\widetilde{A}}{}^{R,*}$ instead of $\widehat{\widetilde{A}}{}^{R}$ follows with Lemma A.5.

Lemma A.5**.**

It holds $n^{3/2}\sup_{s\in[0,\mathcal{T}]}|\widehat{A}(s)-\widehat{A}^{*}(s)|=O_{p}(1)$.

Proof A.6**.**

The proof follows from the fact that $Z_{i}\geq 0$ and $Y_{i}^{R}\leq 1$, $i=1,\dots,n$ together with an analogous argumentation as in the last proof.

First, because of the non-negativity of $Z_{i}$ and the boundedness of $Y_{i}^{R}$, it holds

[TABLE]

for $t\in[0,1]$, with the notation $\widehat{\widetilde{A}}{}^{R,*}_{i}(t)$ from Section A.2. From here, a completely analogous argumentation with a martingale \widetilde{M}^{R,*}_{i}(t)=\Big{(}\widehat{\widetilde{A}}{}^{R,*}_{i}\Big{)}^{2}-\Lambda_{i}^{R,*}(t) for the compensator

[TABLE]

leads to the central limit theorem

[TABLE]

The same argument with Lenglart’s inequality as in the proof of Corollary A.3 yields the conclusion.

A.4 Proof of Proposition 4.3

We first introduce the notation

[TABLE]

for every kernel $K$ and remind of the notation

[TABLE]

Since

[TABLE]

one can easily verify that $n^{-1}\sum_{i}\overline{K}_{t,h}(t-s)\widetilde{Y}_{i}^{R}(s)$ converges locally uniform almost surely to $\overline{K}^{*}_{h}(t-s)$, where $\overline{K}^{*}_{h}$ arises from $\overline{K}^{*}$ by replacing $u$ and $K(u)$ with the local versions $h^{-1}u,h^{-1}K(u/h)$ (Nielsen & Tanggaard, 2001). Furthermore, if $K$ is symmetric, then $\overline{K}^{*}(t)=K(t)$.

From equation (11), Assumption S3 and Corollary A.3, we conclude that it is enough to consider the asymptotic behaviour of

[TABLE]

Analogously to the local constant case, we split the estimation error into a stable and a martingale part

[TABLE]

where

[TABLE]

The asymptotic limit of the bias part, $B_{1}$, is now easily derived via a second order Taylor expansion. The martingale part can be concluded with similar arguments as in Appendix A.2.

A.5 Proof of Proposition 3.1

The proof follows along the same lines as the proof of Proposition 4.1 above, just simpler, with the choice

[TABLE]