2.1. The probabilistic setup

Let $J=\left(J_z\right)_{z \geqslant 0}$ be a non-explosive pure jump process on a finite state space $\mathcal{S}$ and take $(\Omega, \mathcal{F}, \mathbb{P})$ to be the underlying probability space. Here, we let $\mathcal{S}=\{1, \ldots, k\}$ , with $k\in \mathbb{N}$ . Denote by Y the possibly infinite absorption time of J. Let X be a random variable with values in $\mathbb{R}^d$ equipped with the Borel $\sigma-$ algebra. We assume that the distribution of X admits a density g with respect to the Lebesgue measure $\lambda$ . Absolute continuity in case of known atoms may be relaxed as discussed in Bladt and Furrer (Reference Bladt and Furrer2023b).

We introduce a multivariate counting process N with components $N_{j h}=\left(N_{j h}(z)\right)_{z \geqslant 0}$ given by

\begin{align*}N_{j h}(z)=\#\left\{s \in(0, z]: J_{s-}=j, J_s=h\right\}, \quad z \geqslant 0,\end{align*}

for $j, h \in \mathcal{S}, j \neq h$ . Let us introduce the following assumption.

Assumption 1. For every $z\geqslant 0$ and $j,h\in\mathcal{S},j\neq h$ there exists some $\delta > 0$ such that

\begin{align*}\mathbb{E}\left[N_{jh}(z)^{1+\delta}\right]< \infty.\end{align*}

Taking the largest z and the minimum across the whole states provides a delta which is independent of j,h,z.

For any x in the support of X define the conditional occupation probabilities according to

\begin{align*}p_j(z \mid x)=\mathbb{E}\left[\unicode{x1D7D9}_{\left\{J_z=j\right\}} \mid X=x\right] .\end{align*}

We denote by $p(z \mid x)$ the row vector with elements $p_j(z \mid x)$ . We now introduce the cumulative conditional transition rates and cast the conditional occupation probabilities as a product integral of these cumulative rates. To this end, let

\begin{align*}p_{j h}(z \mid x)=\mathbb{E}\left[N_{j h}(z) \mid X=x\right],\end{align*}

and define the cumulative conditional transition rates via

\begin{align*}\Lambda_{j h}(z \mid x)=\int_0^z \frac{1}{p_j(s-\mid x)} p_{j h}(\mathrm{d} s \mid x).\end{align*}

A cumulative transition hazard matrix is defined by having the j,j entry $\Lambda_{j,j}(z \mid x)=\sum_{j\neq}\Lambda_{j,h}(z \mid x)$ for $j, h \in \mathcal{S}, j \neq h$ .

Let $q(\cdot \mid x)$ be the interval function associated with the conditional transition probabilities:

\begin{align*}q_{j h}(s, z \mid x)=\unicode{x1D7D9}_{\left\{p_j(s|x)> 0\right\}} \frac{\mathbb{E}\left[\unicode{x1D7D9}_{\left\{J_z=h\right\}} \unicode{x1D7D9}_{\left\{J_s=j\right\}} \mid X=x\right]}{p_j(s|x)}+\unicode{x1D7D9}_{\left\{p_j(s|x)=0\right\}} \unicode{x1D7D9}_{\{j=h\}}.\end{align*}

In Overgaard (Reference Overgaard2019), it is shown that this assumption ensures that the product integral

is well defined. Therefore,

(2.1)

We are interested in the situation where the observation of J is right-censored. To this end, we introduce a strictly positive random variable W describing right censoring, so that we actually observe the triplet

\begin{align*}\left(X,\left(J_z\right)_{0 \leqslant z \leqslant W}, Y \wedge W \right) .\end{align*}

Whenever the Markov property fails for the jump process J, the following assumption is key for our estimation procedure.

Assumption 2. Right censoring is conditionally entirely random:

\begin{align*} W {\perp \!\!\! \perp}\ J \mid X.\end{align*}

Other censoring schemes are possible if J is Markovian, though we do not make that assumption in the sequel.

We introduce

\begin{align*}p_j^c(z \mid x) =\mathbb{E}\left[\unicode{x1D7D9}_{\left\{J_z=j\right\}} \unicode{x1D7D9}_{\{z< W\}} \mid X=x\right],\end{align*}

and

\begin{align*}p_{j h}^c(z \mid x) =\mathbb{E}\left[N_{j h}(z \wedge W) \mid X=x\right] .\end{align*}

from which we may conclude that

(2.2) \begin{equation}\Lambda_{j h}(z \mid x)=\int_0^z \frac{1}{p_j^c(s-\mid x)} p_{j h}^c(\mathrm{d} s \mid x),\end{equation}

within the interior of the support of W.

2.1.1. Estimators

Consider the i.i.d. replicates $\left(X^{i},\left(J_z^{i}\right)_{0 \leqslant z \leqslant W^{i}}, Y^{i} \wedge W^{i}\right)_{i=1}^n$ . Define $\delta^i:=\unicode{x1D7D9}_{\left\{Y^i\leqslant W^i\right\}}$ , which equals 1 when the observation is absorbed (closed claim), and zero otherwise. In a reserving context, this distinction is usually reflected by specifying two types of claims: closed claims and RBNS (open) claims, and we write, respectively, $n=n^{\texttt{Closed}}+n^{\texttt{RBNS}}$ .

We now construct an estimator for the conditional occupation probabilities $p_j(z\mid x)$ , but from Equation (2.1), it essentially suffices to estimate the cumulative conditional transition rates $\Lambda_{j h}(z \mid x)$ , which we do using Equation (2.2), Let $K_b$ be kernel functions of bounded variation and bounded support, with $b=1,\ldots,d$ and let $(a_n)$ be a bandwidth sequence. We impose the following standard assumption, confering with Stute (Reference Stute1986).

Assumption 3. The band sequence satisfies:

1. 1. $a_n \rightarrow 0$ as $n \rightarrow \infty$ .

2. 2. $n a_n^d \rightarrow \infty$ as $n \rightarrow \infty$ .

3. 3. $\sum_{n=1}^{\infty} c_n^r< \infty$ for some $r \in(1,1+\delta)$ where $c_n:=\left(n a_n^d\right)^{-1} \log n$ .

For instance, one may take, for $(1+\delta)^{-1}< \eta< 1$ , any sequence satisfying

\begin{align*}a_n^d \sim \frac{\log n}{n^{1-\eta}}.\end{align*}

We introduce the usual kernel estimator for the density g:

\begin{align*}{\unicode{x1D558}}^{(n)}(x)=\frac{1}{n} \sum_{i=1}^n {\unicode{x1D558}}^{(n, i)}(x)=\frac{1}{n} \sum_{i=1}^n \prod_{b=1}^d \frac{1}{a_n} K_b\left(\frac{x_b-X_b^{i}}{a_n}\right) .\end{align*}

It is well-known (cf. Ghosh, Reference Ghosh2018, Ch. 1.4) that

\begin{align*}{\unicode{x1D558}}^{(n)}(x) \stackrel{\text { a.s. }}{\rightarrow} g(x), \quad n \rightarrow \infty,\end{align*}

and, for any random variable V with i.i.d. replicates $(V^i)$ satisfying $\mathbb{E}[|V|^{1+\delta}]< \infty$ that

(2.3) \begin{align}\frac{1}{n} \frac{\sum_{i=1}^n V^{i} {\unicode{x1D558}}^{(n, i)}(x)}{{\unicode{x1D558}}^{(n)}(x)} \stackrel{\operatorname{a.s.}}{\rightarrow} \mathbb{E}[V \mid X=x], \quad n \rightarrow \infty .\end{align}

Based on the kernel estimator for g, we form the following Nadaraya–Watson type kernel estimators:

(2.4) \begin{align}& \mathbb{N}_{j h}^{(n)}(z \mid x)=\frac{1}{n} \sum_{i=1}^n \mathbb{N}_{j h}^{(n, i)}(z \mid x)=\frac{1}{n} \sum_{i=1}^n \frac{N_{j h}^{i}\left(z \wedge W^i\right) {\unicode{x1D558}}^{(n, i)}(x)}{{\unicode{x1D558}}^{(n)}(x)} \end{align}

(2.5) \begin{align}& \mathbb{I}_j^{(n)}(z \mid x)=\frac{1}{n} \sum_{i=1}^n \mathbb{I}_j^{(n, i)}(z \mid x)=\frac{1}{n} \sum_{i=1}^n \frac{\unicode{x1D7D9}_{\left\{z< W^{i}\right\}} \unicode{x1D7D9}_{\left\{J_z^{i}=j\right\} }{\unicode{x1D558}}^{(n, i)}(x)}{{\unicode{x1D558}}^{(n)}(x)} \end{align}

Remark 2. Due to Assumption 1, it follows from Equation (2.3) that

\begin{align*}\mathbb{N}_{j h}^{(n)}(z \mid x) \stackrel{\text { a.s. }}{\rightarrow} p_{j h}^c(z \mid x), \quad n \rightarrow \infty.\end{align*}

Furthermore, it holds that

\begin{align*}\mathbb{I}_j^{(n)}(z \mid x) \stackrel{\text { a.s. }}{\rightarrow} p_j^c(z \mid x), \quad n \rightarrow \infty .\end{align*}

Definition 1. The conditional Nelson–Aalen estimator for the cumulative conditional transition rates is defined as

\begin{align*}\mathbb{\Lambda}_{j h}^{(n)}(z \mid x)=\int_0^z \frac{1}{\mathbb{I}_j^{(n)}(s-\mid x)} \mathbb{N}_{j h}^{(n)}(\mathrm{d} s \mid x) .\end{align*}

Definition 2. The conditional Aalen–Johansen estimator for the conditional occupation probabilities is defined according Equation (2.1), but with the cumulative conditional transition rates replaced by the conditional Nelson–Aalen estimator, that is,

where ${\unicode{x1D561}}_j^{(n)}(0 \mid x)=\mathbb{I}_j^{(n)}(0 \mid x).$ The estimator ${\unicode{x1D561}}_j^{(n)}(0 \mid x)$ for the occupation probabilities $p_j(z|x)$ denotes the jth entry of the row vector ${\unicode{x1D561}}^{(n)}(z \mid x)$ , an estimator for $p(z|x)$ .

Remark 3. We refer to Bladt and Furrer (Reference Bladt and Furrer2023b) for the proof that the estimators are pathwise strongly consistent and asymptotically Gaussian distributed.

When all the subjects start in state 1 at inception, we can use the estimator of the occupation probabilities of the absorbing state k, ${\unicode{x1D561}}_k^{(n)}$ to model the cumulative density function of the individual claim size.

Remark 4. The Volterra integral equation defines the product integral of the matrix function $\Lambda(\mathrm{d} s \mid x)$ as the unique solution $p(\cdot \mid x)$ to the equation

\begin{align*} p(t \mid x)=\operatorname{Id}+\int_{(0,t]}p(s-\mid x) \Lambda(\mathrm{d} s \mid x). \end{align*}

However, when a piecewise constant estimator of $\Lambda$ is given by $\mathbb{\Lambda}^{(n)}$ , as is described in detail in the previous section, the computation of the product integral simplifies. More precisely, denote $0=t_0< t_1< t_2< \cdots$ the sequence of time points describing all observed jump times in the sample. Then the recursion

\begin{align*}{\unicode{x1D561}}^{(n)}_j(t_{\ell+1} |x)-{\unicode{x1D561}}^{(n)}_j(t_{\ell} |x)=\sum_{\substack{k \in \mathcal{Z} \\ k \neq j}} {\unicode{x1D561}}^{(n)}_k(t_{\ell} | x)\Delta\mathbb{\Lambda}_{kj}^{(n)}(t_{\ell +1} | x) -{\unicode{x1D561}}^{(n)}_j(t_{\ell} | x) \sum_{\substack{k \in \mathcal{Z} \\ k \neq j}} \Delta\mathbb{\Lambda}_{jk}^{(n)}(t_{\ell +1} | x)\end{align*}

and starting value ${\unicode{x1D561}}^{(n)}(0 | x) = \mathbb{I}^{(n)}(0 | x)$ completely characterize the product integral ${\unicode{x1D561}}^{(n)}_j(\cdot |x)$ at jump times (and constant otherwise).

2.2. Practical considerations for claims development

In this section, we offer practical insights into our modeling approach, with the aim of assisting the reader in understanding our multi-state setup for claims reserving. Refer to Figure 1, which conceptually illustrates our model featuring a state space $\mathcal{S}$ with $k=5$ states.

Figure 1.

Example of multi-state model for claims reserving, with $k=5$ . We interpret time spent in each state as increasing claims size clock, instead of calendar time. The states represent the development periods (DP’s) of the individual claims.

In our framework, we represent the active state as IBNR claims, while we treat RBNS claims as transitioning through intermediate states. Any claim that reaches state k is considered closed in our system due to the presence of an absorbing state. Consequently, we do not allow claims to reopen in our framework. Notably, unlike a traditional survival study, we do not explicitly model the time component, and the number of IBNR claims, along with their characteristics, remains unknown at the time of evaluation. In Section 3, we delve into how to model RBNS and IBNR claims, showcasing how our model can accommodate such scenarios. Furthermore, we have designed our state space so that state DP3+ at step $k-1$ in Figure 1 is the final relevant DP that we explicitly model.

Finally, it is worth mentioning that we have not introduced development triangles yet, which are the standard aggregate datasets commonly used in claims reserving applications. The choice of parameter k aligns with the data that, in a conventional modelling approach, would be aggregated into development triangles of varying sizes. Since our proposed framework does not explicitly model the time component, we assume that payments close at the development year following the last payment. This assumption results in a development triangle with $k-1$ development times. This aspect will be particularly relevant for our data application, where we compare our model to the benchmark chain ladder method outlined in Mack (Reference Mack1993), and we need to construct a development triangle for this purpose.

3.1. A comparative literature review

Research papers on loss reserving can be broadly categorized into two main streams of research: models based on development triangles (aggregate loss reserving models) and models based on individual data (individual loss reserving models). The relationship between individual data and aggregate modelling has been extensively studied in Miranda et al. (Reference Miranda, Nielsen and Verrall2012), Hiabu (Reference Hiabu2017), and Bischofberger et al. (Reference Bischofberger, Hiabu and Isakson2020), utilizing different survival analysis tools.

The existing body of literature on aggregate reserving has already provided the building blocks for predicting the distribution of loss reserves, with a specific emphasis on models that seek to emulate the empirical chain-ladder algorithm. In this section, we briefly mention some of these models and defer to the comprehensive discussion in Hess and Schmidt (Reference Hess and Schmidt2002) for their detailed categorization. For a distribution-independent approach to computing the mean squared prediction error of the loss reserve, one can refer to Mack (Reference Mack1993), where the author decomposes the prediction error of the loss reserve into its elemental components, namely process variance and estimation variance, as expressed by $\mathbb{E}[(\hat R- R)^2]=\mathrm{Var}(R)+\mathbb{E}[(\hat R- \mathbb{E} [\hat R])^2]$ . The model has received several statistical critiques through the years but remains one of the best-performing reserving models in an actuary’s toolkit.

In England and Verrall (Reference England and Verrall1999), the authors replicate the point estimates of the chain-ladder method and ingeniously employ bootstrap and generalized linear models (GLM) to compute $\mathbb{E}[(\hat R- R)^2]$ . For an extension of the model, allowing for negative payments as well as an exploration of the model assumptions underpinning it, one may consult Verrall (Reference Verrall2000). An interesting comparison of the fundamental assumptions underpinning the aforementioned models is presented in Mack and Venter (Reference Mack and Venter2000), where the authors conclude that only the approach outlined in Mack (Reference Mack1993) is the true underlying model of the chain-ladder algorithm. In a more recent development, in Al-Mudafer et al. (Reference Al-Mudafer, Avanzi, Taylor and Wong2022), the approach in England and Verrall (Reference England and Verrall1999) is integrated into the framework of mixture density networks (Bishop, Reference Bishop1994).

Shifting our focus to non-parametric models for individual reserving, the current discussion primarily centers on providing precise point forecasts for the loss reserve. In contrast, Delong and Wüthrich (2023) discuss the role of the variance function in regression problems for pricing.

In another recent development, Wüthrich (Reference Wüthrich2018) incorporates neural networks into the chain-ladder model, introducing features to the existing methodology. Further related work can be found in Lopez et al. (Reference Lopez, Milhaud and Thérond2019) and Lopez and Milhaud (Reference Lopez and Milhaud2021), where the authors employ a tree-based algorithm for modeling RBNS claims using censored data. In particular, in Lopez and Milhaud (Reference Lopez and Milhaud2021), the use of bootstrapping for estimating the uncertainty in the claims reserve is discussed.

Our manuscript introduces a novel approach to the individual reserving problem, allowing us to explicitly model curves of individual claim sizes denoted as $p_k(z|x)$ within our framework.

3.2. Prediction of RBNS and IBNR claim costs

We can calculate the final cost of RBNS claims numerically using the conditional Aalen–Johansen estimate ${\unicode{x1D561}}_k^{(n)}(z \mid x)$ for $p_k(z|x)$ . The predictor of the final total cost of RBNS claims is defined by

\begin{equation*}\hat Y^{\texttt{RBNS}}=\sum^n_i\unicode{x1D7D9}_{\left\{\delta^i = 0\right\}} \left(W^{i} + \widehat{\mathbb{E}[Y| Y> W^{i}, X=x]}\right)\end{equation*}

with

(3.1) \begin{equation}\widehat{\mathbb{E}\left[Y| Y> W^{i}, X=x\right]} = \frac{1}{1-{\unicode{x1D561}}_k^{(n)}(W^i \mid x)}\int^{+\infty}_{W^{i}} (1-{\unicode{x1D561}}_k^{(n)}(y \mid x)) \mathrm{d} y.\end{equation}

We remark that the general formula for the mth moment is explicitly given by

\begin{align*}\widehat{\mathbb{E}\left[Y^m| Y> W^{i}, X=x\right]} = \frac{1}{1-{\unicode{x1D561}}_k^{(n)}(W^{i} \mid x)} \int_{W^{i}}^{+\infty} m y^{m-1}(1-{\unicode{x1D561}}_k^{(n)}(y \mid x)) \mathrm{d} y.\end{align*}

Strictly speaking, from a reserving standpoint, we also need to propose a strategy for reserving for the (unknown) IBNR claims. We denote the total number of IBNR claims by $n^{\texttt{IBNR}}$ , which is a random variable with values in $\mathbb{N}_0$ . The individual cost of the IBNR claims is described by a sequence of iid copies of the random variable $\tilde Y$ , namely $\left\{ \tilde Y^{i^\prime} \right\}_{i^\prime \in \mathbb{N}}$ . Under the framework of Bladt and Furrer (Reference Bladt and Furrer2023b), $\tilde Y$ has cumulative distribution function $p_k(z)$ . As we will describe shortly, we specify a model that does not consider the features, since in this case we do not have access to future exposures by feature. We propose to correct the model for the random presence of IBNR using a compound distribution, under the assumptions in Klugman et al. (Reference Klugman, Panjer and Willmot2012, p. 148). We describe the total cost of IBNR claims as follows:

\begin{align*}Y^{\texttt{IBNR}}=\sum^{n^{\texttt{IBNR}}}_{i^\prime=1} \tilde Y^{i^\prime}.\end{align*}

The expression for $Y^{\texttt{IBNR}}$ is well known in non-life mathematics and is often referred to as the collective risk model, Parodi (Reference Parodi2014). The moments of $Y^{\texttt{IBNR}}$ can be calculated in closed form. In particular, the expected cost of IBNR claims is predicted by

\begin{align*}\hat Y^{\texttt{IBNR}} = \widehat{\mathbb E[Y^{\texttt{IBNR}}]}= \hat n^{\texttt{IBNR}}\widehat{\mathbb{E} [\tilde Y]}, \end{align*}

with the m-th moment of $\tilde Y$ predicted by

\begin{align*}\widehat{ \mathbb{E} [\tilde Y^m]} = \int^{\infty}_0 m y^{m-1}(1-{\unicode{x1D561}}_k^{(n)}(y)) \mathrm{d} y,\end{align*}

fitting the unconditional Aalen–Johansen to the sample n to obtain ${\unicode{x1D561}}_k^{(n)}(z)$ . To project the number of IBNR claims, we simply use the chain ladder model on the triangle of reporting delays $\mathcal{D}^{(k-1)}$ and obtain $\hat n^{\texttt{IBNR}}$ .

Similarly, the variance of the IBNR claims is predicted by

\begin{align*}\widehat {\mbox{Var}( Y^{\texttt{IBNR}})}= \hat n^{\texttt{IBNR}} \left[\widehat{\mathbb{E} [\tilde Y^2]}- \widehat{\mathbb{E} [\tilde Y]}^2\right]+\widehat{\mathbb{E} [\tilde Y]}^2 \widehat{\mbox{Var}(n^{\texttt{IBNR}})}.\end{align*}

We denote by $\widehat{\mbox{Var}(n^{\texttt{IBNR}})}$ a predictor for the variance of $n^{\texttt{IBNR}}$ . In this manuscript, we consider the predictor for the process variance discussed in Mack (Reference Mack1993).

The total size of claims is $Y^{\texttt{TOT}}= Y^{\texttt{Closed}}+Y^{\texttt{RBNS}}+Y^{\texttt{IBNR}}$ and we estimate it with $\hat Y^{\texttt{TOT}}= Y^{\texttt{Closed}} + \hat Y^{\texttt{RBNS}}+\hat Y^{\texttt{IBNR}}$ . We remark that in several lines of business the contribution of the true cost of IBNR claims ( $ Y^{\texttt{IBNR}}$ ) to $Y^{\texttt{TOT}}$ is minor compared to the cost of RBNS claims ( $ Y^{\texttt{RBNS}}$ ). This is the case for those lines of business that have data with mostly single payments and a short time lag between reporting and settlement. See for instance the results presented in Friedland (Reference Friedland2010, p. 80) on the US market data, comparing for different lines of insurance the size of the estimated cost of IBNR claims to paid and outstanding claims. As we will show in Section 5, our real dataset exhibits exactly this behaviour and the IBNR correction will have a minor impact on the estimated final cost.

4.2. Evaluating the models performance

Using the approach introduced in Section 4.1, we present two simulated scenarios: Alpha and Beta. In scenario Alpha, we simulate directly with unconditional transition intensities that we have fitted to the real data. In the Beta scenario, we introduce the effect of the accident period characteristics and the occupation intensities. In each scenario we simulate, for the different choices of $k=4,5,6,7$ , a total number of 20 datasets to smooth possible fluctuations due to the random generator and report the average results that we obtain. Each dataset in scenario Alpha contains $1200(k-1)$ observations. The dataset in scenario Beta contains $1200(k-1)-100(k-2)$ observations. As motivated in Section 4.1, we let the volumes decrease in recent accident periods. We chose arbitrarily to decrease the volume by 100 in every accident period subsequent to the first one.

The choice of measure of model performance can be an ardouos task when models are very different, as are AJ and CL. The first performance measure that we propose is the error incidence,

\begin{align*}\texttt{EI}= \frac{\hat Y^{\texttt{TOT}}}{Y^{\texttt{TOT}}}-1,\end{align*}

where $\hat Y^{\texttt{TOT}} $ is the predicted final claim size and $Y^{\texttt{TOT}}$ is the actual final claim size, which is available from the simulation. The EI is an interesting benchmark that provides insight into the relative accuracy of the AJ compared to the CL, though it does not highlight the fact that AJ provides a much more nuanced description of the claim size than point estimates. In Figure 3, we provide a box plot of the EI for the Alpha scenario (left side) and the Beta scenario (right side) for the different values of k. Compared to the CL, the AJ model performs favourably for each choice of k and in both scenarios.

Figure 3.

Box plots of EI for AJ and CL over the 20 simulations, for each value of k, in the Alpha scenario (left column) and the Beta scenario (right column).

While the EI is an interesting benchmark, to better understand the performance of the models compared to the data at hand, we would ideally like to use a proper scoring measure (Gneiting and Raftery, 2007). In principle, a scoring metric with the property of being proper is able to evaluate the better model by assigning better scores to better models. In this manuscript we consider the continuously ranked probability score (CRPS) as a proper scoring metric:

\begin{align*}\operatorname{CRPS}\left({\unicode{x1D561}}_k^{(n)}(z \mid x), y\right) &= \int_{0}^{+\infty} \left({\unicode{x1D561}}_k^{(n)}(z \mid x) - \unicode{x1D7D9}_{\{y \leqslant z\}}\right)^2 \mathrm{d} z \\&= \int^{y}_0 \left({\unicode{x1D561}}_k^{(n)}(z \mid x)\right)^2 \mathrm{d}z + \int^{+\infty}_y \left(1-{\unicode{x1D561}}_k^{(n)}(z \mid x)\right)^2 \mathrm{d}z,\end{align*}

where ${\unicode{x1D561}}_k^{(n)}(z \mid x)$ is the predicted curve of the individual claims sizes and $y_i$ is the (true) value of the final claim size that is available from the simulation. The interested reader can refer to the established statistical literature on proper scoring for a detailed discussion on this topic Gneiting and Raftery, Reference Gneiting and Raftery2007; Gneiting and Ranjan, Reference Gneiting and Ranjan2011).

4.3. Empirical analysis

The average results of our 20 simulations for the Alpha scenario are reported in Table 1. Each column of the table represents one of the AJ models fitted for the different choices of k (column one). The target $Y^{\texttt{TOT}}$ (average simulated total cost) assuming $Y^{\texttt{IBNR}}=0$ is in column four. The predicted final cost is $\hat Y^{\texttt{TOT}} = Y^{\texttt{Closed}}+ \hat{Y^{\texttt{RBNS}}}$ . The EI results show that, on average, our models outperform the CL for every choice of k except $k=5$ (columns five and six). For all choices of k, we find that the relative variability of our model is less than that provided by the CL (columns seven and eight). We compute $\mbox{sd}(Y^{\texttt{TOT}})=\sqrt{\sum^{n^{\texttt{RBNS}}}_i \mbox{Var}(Y^i)}$ . In this manuscript, we will not discuss the estimation error, and the results we present with respect to the variability of the reserve refer only to the process variance. In column nine, we report the (average) CRPS.

Table 1.

Results for scenario Alpha. For each value of k (column one) we present the average results over the 20 simulations. Each row of the table corresponds to a different AJ model specification (column three). The table includes the (average) actual Y $^{\texttt{TOT}}$ simulated total cost (column four) and the error incidence for the AJ and the CL (columns five and six). In columns seven and eight, we reported the coefficients of variation of Y $^{\texttt{TOT}}$ . The results for the CRPS are reported in column nine.

Table 2.

Results for scenario Beta. For each value of k (column one), we present the average results over the 20 simulations. Each row of the table corresponds to a different AJ model specification (column three). The table includes the (average) actual simulated total cost (column four) and the error incidence for the AJ and the CL (columns five and six). In columns seven and eight, we reported the coefficients of variation of Y $^{\texttt{TOT}}$ . The results for the CRPS, relative to the AJ model with features, are reported in column nine.

Table 2 reports the average results for the 20 simulations of the Beta scenario. For each simulation, we fit two models, with and without the use of U as a feature ( $\checkmark$ and ✗ in the table). The model with U is correctly specified assuming the Beta scenario. The aim of this experiment is to show that the CRPS is able to judge which is the better model in a context where we know the data generation process. As expected, we obtain a lower average CRPS for the models that include features. The CRPS results in Table 2 are relative to the AJ model with feature U. For example, for $k=6$ we see that the AJ model without the feature U shows (on average) a CRPS $3.9 \%$ higher than the CRPS of the AJ model with feature U. For the purpose of an easier interpretation, the CRPS results will be shown relative to the AJ with feature model also in the following tables. The (average) target $Y^{\texttt{TOT}}$ is reported in the fourth column. The results for EI indicate that the AJ model generally outperforms the CL model. The relative variability results show comparable results for AJ and CL.

5.1. Model comparison on different datasets

Consider the first application, where we censor our data at different depths and fit an AJ model using the largest possible state space. Using our model, it is possible to explore and compare the individual claim size curves for the different values of claim_type, as shown in Figure 5.

The results are reported in Table 4, and they show that in general the AJ models outperform the CL in terms of EI, except for the scenario with $k=7$ . For each scenario, we are able to select a model with or without features using the CRPS measure. The rows with a $\checkmark$ (column two) refer to models that include the claim_type feature in the fit. The target $Y^{\texttt{TOT}}$ is shown in column three. While in other scenarios, the CRPS for the AJ model that includes the claim_type feature is close to or better than the CRPS of the model without the feature, in $k=5$ the model without the feature has a much lower CRPS. The relative variability results show that AJ and CL have comparable results. In fact, only for $k=4$ does CL have a lower relative variability.

Table 4.

For different choices of k (column one), we fit a model with and without $\texttt{claim_type}$ (column two). The target Y $^{\texttt{TOT}}$ is shown in column three. The EI is shown in columns four and five. The CV of Y $^{\texttt{TOT}}$ are displayed in columns seven and eight for the AJ and CL respectively. The CRPS is shown in column nine.

Figure 5.

For each different dimension $k=4,5,6,7,8$ , we provide the individual claim size curves for our observations by covariate value $\texttt{claim_type}$ . The x-axis represents Z and we scale it by $\texttt{log2}$ to ease the plot visualization.

For each choice of k, we can simply select the best-performing model in terms of CRPS and compute the claims reserve, the results are displayed in Table 5.

5.2. Model comparison on a single dataset

In the second application, we censor our data after 5 accident period and construct the AJ model for $k=4,5,6$ . For each choice of k, we fit a model both with and without the feature claim_type. Table 6 shows that in terms of CRPS, the model with $k=4$ using the claim_type feature is the best model (minimum CRPS). This model is also the best in terms of EI.

Using the CRPS, we are able to select the model with $k=4$ and that uses claim_type as a feature to compute the claims reserve. The model provides a reserve of $11.485$ millions and a standard deviation of $ 2.0385$ millions. The strategy presented in this last section is of particular interest to reserving actuaries. While the empirical study in Section 5.1 is interesting for understanding the behaviour of our model on different data, in practice the maximum depth of the data is selected using indicators such as claims settlement speed or using the expert judgement of an experienced actuary Wüthrich and Merz (Reference Wüthrich and Merz2008, p. 12). The depth of the data is then known and the main interest is simply to select the best model, as shown here. Selecting the most appropriate model can be very difficult, see again Wüthrich and Merz (Reference Wüthrich and Merz2008, p. 13), and the strategy presented in this section provides a mathematically sound approach to doing so.

Table 5.

We selected for each data k (column one), the best-performing model in terms of CRPS and present the reserve (column one) and the standard deviation of the reserve (column two).

Table 6.

For the dataset with depth 5 accident periods, we fit the AJ model both including and excluding the $\texttt{claim_type}$ feature (column two). The target Y $^{\texttt{TOT}}$ is displayed in column three. The EI for the AJ and the CL can be found in columns four and five. The coefficient of variation of Y $^{\texttt{TOT}}$ is displayed in columns six and seven. The CRPS is displayed in column eight.