Combined modelling of micro-level outstanding claim counts and individual claim frequencies in non-life insurance

Actuarial departments in insurance companies are responsible for many different tasks related to risk modelling of insurance operations. Two of these tasks are closely interlinked: reserving actuaries need to compute accurate predictions of incurred liabilities for insurance portfolios. Pricing actuaries use these results to develop risk models, which allow them to estimate expected losses for a range of policy parameters and thus to determine the prices at which these policies are sold. While developing models on per-policy data (subsequently referred to as the micro-level) is common practice for risk modelling, micro-level methods for reserving are still not widely adopted in the industry. Adoption is slow in part due to the requirements necessary for micro-level reserving (such as fine-grained claims development information, the necessary data warehousing and compute infrastructure as well as advanced modelling skills such as machine learning in combination with actuarial expertise) and the strong governing bodies that will require tried and tested reserving methods for many business processes in the foreseeable future.

In the research literature on reserving in non-life insurance, there are two main approaches to micro-level reserving. On the one hand, there are methods that operate in discrete time (such as development years), which have similarities with macro-level reserving methods like Bornhuetter–Ferguson. Micro-level methods of this kind typically deal with claim-level information which is only available for reported claims, and therefore work on the reserve of claims that have been reported but not settled (RBNS). Examples of this kind are [6, 9, 11, 12, 24], among others. On the other hand, there are continuous time claim development models which are based on the assumption that claim development is regarded as a point process (e.g., a position-dependent marked Poisson process); respective model parameters are then estimated from observed data. This strain of research goes back to [17, 18]. More recent papers studying a continuous time claim development model are [4, 5, 19, 21], among others.

The current paper aims to combine the task of predicting IBNR claim counts on a policy level with the development of a matching micro-level claim intensity model, thus addressing the two tasks of reserving and risk modelling simultaneously in the hope of improving accuracy in both disciplines. Two main approaches are proposed to achieve this goal. Our first approach is based on an explicit micro-level claim occurrence process model inspired by [17]. Suitable methods are proposed for estimating all model parameters, which results in a fully fitted continuous-time process model that includes a claim intensity model. We thereby extend a related approach in [7], where a submodel for reporting delays was studied. Secondly, we propose a new method that uses classical triangle-level reserving methods as input to the estimation of a micro-level claim intensity model. This approach allows for micro-level allocation of IBNR claim counts that is consistent with the triangle. Due to the nature of the underlying triangle-based reserving methods, the approach only allows discrete time steps.

For each of the two approaches, we discuss the design and estimation of suitable (parametric) sub-models involving classical multilayer perceptron (MLP) neural networks, see [14]. For the estimation, it is taken into account that the observed data is subject to random truncation (indeed, a claim is only observed at a given calendar time if the sum of its associated (random) reporting delay and its (random) accident time does not exceed the given calendar time). Predictors for the number of IBNR claims on a per-policy level are derived from the estimated models. They are compared with each other as well as with classical chain ladder methods in an extensive simulation study as well as on a real dataset. It is found that the new predictors provide accurate predictions as well as appropriate fitted models for claim intensities even in simulation scenarios involving non-homogeneous portfolios and in the real-life example. Moreover, in contrast to classical factor-based reserving methods, the predictor based on the first approach may yield a non-zero number of claims even for policies without already reported claims. This is natural as it allows for interpreting the prediction as the expected number of unreported claims for that particular policy.

For learning the neural networks, we utilize the TensorFlow framework [1] with custom loss functions that take random truncation into account. The implementation is written using the R language and its interface binding packages keras and tensorflow [2, 10] to utilize TensorFlow. Core functionality used for estimation and prediction is freely available as an R package reservr [20].

The remaining parts of the paper are organized as follows. In Sect. 2, we introduce the micro-level claim process used throughout as well as the observation setting. Modelling and fitting of the individual components of the claim process are discussed in Sect. 3. Based on a completely estimated micro-level model, we derive a corresponding micro-level predictor for the IBNR claim count in Sect. 4. Section 5 introduces an alternative micro-level predictor based on similar ideas, but using a triangle-based global reserving method and discrete time steps. After defining performance metrics in Sect. 6, we present results on a large-scale simulation study in Sect. 7. An application to a real dataset from a motor legal insurance portfolio is presented in Sect. 8. Finally, Sect. 9 concludes.

The general model for the claim arrivals in a given insurance portfolio is the same as in [7], and builds on the notion of (position-dependent) marked Poisson processes [17]. More precisely, we consider an insurance portfolio containing \({{\mathcal {I}}}\in \mathbb {N}\) independent risks. Each risk is described by a coverage period \(C = [t_{\text {start}}, t_{\text {end}}]\), and by risk features \({\bar{x}} \in {\bar{\mathfrak {X}}}\), where \({\bar{\mathfrak {X}}}\) is a feature space containing both discrete and continuous features; for example, information on the insured product and chosen options such as deductibles. We write \(x = (C, {\bar{x}}) \in \mathfrak {X}= \{\text {intervals on } [0,\infty )\} \times {\bar{\mathfrak {X}}}\), and assume that x is constant over the course of the contract. In practice, risk features do change over time, but not very often, whence such a contract could be modelled as two separate risks.

Each risk can potentially incur claims during its coverage period, which will formally be modelled by a claim arrival process. Note that (marked) claim arrival processes provide a natural mathematical model for random event analysis, which have been proposed for actuarial risk modelling in [17]; see also the textbook [16]. Each claim in the claim arrival process occurs at some (calendar) accident time \(t_{\text {acc}} \in [t_\text {start}, t_{\text {end}}]\), and will be associated with several claim features \(y \in \mathfrak {Y}\) like the type of the claim or the value of the cars involved in some accident; here, \(\mathfrak {Y}\) denotes some suitable feature space. Moreover, the claim will not be immediately known to the insurer, but it will rather be reported at (calendar) reporting time \(t_{\text {report}} \in [t_{\text {acc}}, \infty )\). Formally, both the claim features y and the reporting delay \(d_{\text {report}} {:}{=}t_{\text {report}} - t_{\text {acc}}\) will be assumed to be a mark on the claim arrival process.

Let \(\delta _z\) denote the dirac-measure at z. Following the notation in [15], we arrive at the following definition of a claim arrival process. The process is illustrated in Fig. 1 for the case that the claim feature is binary.

Note that the claim intensity \(\lambda (x^{(i)}, t)\) controls the expected number of claims per exposure (the claim frequency), i.e., the expected amount of claims of the ith risk in the period \(A\subset C(x^{(i)})\) is given by \(\int _{A} \lambda (x^{(i)}, t) \mathrm {\,d}t\).

In practice, the three building blocks of Definition 2.1, i.e., \(\lambda (x,t), P_{Y|x,t}\) and \(P_{D|x,t,y}\), are unknown and must be estimated based on observational data that is available to the insurer at some given calendar time \(\tau > 0\) of observing the portfolio. More precisely, the insurer observes data from Observation Scheme 2.2, which is again illustrated in Fig. 1.

Clearly, estimating the building blocks of Definition 2.1 can only be feasible if additional model assumptions are made. Those assumptions, as well as respective fitting procedures, will be described in the next section.

Modelling and estimating the claim intensity \(\lambda (x,t)\), the claim feature distribution \(P_{Y|x,t}\) and the reporting delay distribution \(P_{D|x,t,y}\) from Definition 2.1 will be done iteratively, starting with the latter. We discuss each aspect in a separate subsection.

The models and estimators from the previous section can be used in various ways to define predictors for IBNR claim numbers; see [7] for an example that only involves the reporting delay model. Throughout this section, we describe a predictor that is based on the full (estimated) claim arrival model. Alternative intermediate predictors will be defined in the simulation study.

More precisely, for each given period \(I_p(k)=[k p, (k + 1) p)\) of length \(p>0\) and each claim feature set \(\mathfrak {Y}' \subset \mathfrak {Y}\) and each reporting interval \((\tau _0, \tau _1] \subset [0,\infty ]\), we derive a predictor for the number of claims policy i has incurred within period \(I_p(k)\) with claim features in \(\mathfrak {Y}'\) and with a reporting time in \((\tau _0, \tau _1]\). For that purpose, let \(S'(k) {:}{=}I_p(k) \times \mathfrak {Y}' \times [0, \infty )=\{(t,y,d): t \in I_p(k)\}\) denote the set of claims that occurred in the kth period with claim features in \(\mathfrak {Y}'\) and letdenote the set of claims reported between times \(\tau _0\) and \(\tau _1\), see Fig. 2 for an illustration. For completeness, let \(R_{\tau _0:\tau _1}=\varnothing \) if \(\tau _0 \ge \tau _1\).

Note that the target number of claims for the ith policy can then be written asand that we observe, under Observation Scheme 2.2, the respective number of reported claims \(\xi ^{(i)}_r(S'(k) \cap R_{\tau _0:\tau _1}) = \xi ^{(i)}\bigl (S'(k) \cap R_{\tau _0:\min (\tau _1, \tau )}\bigr )\), which is zero if \(\tau _0 > \tau \).

Now, if Assumption 3.4 is met for the given period length \(p>0\), we obtain that, by the restriction theorem (Theorem 5.2 in [15]),Here, \(\xi _{nr}^{(i)} = \xi ^{(i)}-\xi _{r}^{(i)}\) denotes the unknown number of unreported claims, \(I_p(x,k)\) is the coverage time of the policy associated with x within the kth period, see (8), and \(I_{\tau _0:\tau _1}(\tau , s) {:}{=}( \max (\tau , \tau _0) - s, \tau _1 - s ]\), with the convention that the interval is the empty set if \(\max (\tau , \tau _0)> \tau _1\). As is well-known, if \(\lambda (x,t), P_{Y|x,t}\) and \(P_{D|x,t,y}\) were known, this would be the best \(L^2\)-predictor for \(\xi ^{(i)}(S'(k) \cap R_{\tau _0:\tau _1})\) (the total number of claims in \(S'(k)\) that are reported between \(\tau _0\) and \(\tau _1\)) in terms of the observed counterpart of reported claims \(\xi _r^{(i)}(S'(k) \cap R_{\tau _0:\tau _1})\). Replacing the unknown objects on the right-hand side by suitable estimators as in the previous sections, we arrive at the predictorIn contrast to classical factor-based reserving methods, this predictor may yield a non-zero expected number of claims even for policies without already reported claims. This allows for the individual-level count predictions to have a meaningful interpretation as the expected number of unreported claims for that particular policy.

We propose an alternative estimator for the claim intensity \(\lambda \), which is similar to the estimator from Sect. 3.3. However, instead of being based on preliminary estimators of the claim feature and reporting delay distributions, the new estimator is based on classical chain ladder factors. In a second step, the estimator is used to define a new predictor for IBNR claims, similarly as in Sect. 4.

Given a partition \(\mathfrak {Y}= \mathfrak {Y}_1 \cup \mathfrak {Y}_2 \cup \cdots \cup \mathfrak {Y}_M\) into groups of claim features and a development period length p where \(\tau = \bar{\tau } p\) for some \({\bar{\tau }}\in \mathbb {N}_{\ge 2}\), we make the following classical chain ladder assumption: for any group index \(m \in \{1, \dots , M\}\) and any development period \(j\in \{1, \dots , {\bar{\tau }}-1\}\), there exists a factor \(f_j^{\text {CL}, \mathfrak {Y}_m}\) called chain ladder factor such that, for any policy i and any accident period \(k\in \{0, \dots , {\bar{\tau }}-1\}\),where \(S_m(k) {:}{=}I_p(k)\times \mathfrak {Y}_m \times [0,\infty )\) denotes the set of claims with claim features from \(\mathfrak {Y}_m\) occurring in the kth period and \(R_{0:kp}\) from (11) denotes the set of claims reported until time kp. Iterating the equation for fixed k and with \(j={\bar{\tau }}-1, \dots , {\bar{\tau }}-k\), we obtain thatwhereis the chain ladder factor-to-ultimate. Equation (13) has the following interpretation: the expected number of claims from \(\mathfrak {Y}_m\) with accident period k that are reported within the next k periods from today is equal to the \(\mathfrak {Y}_m\)-specific factor-to-ultimate multiplied with the expected number of claims from \(\mathfrak {Y}_m\) with accident period k that have been reported until today (which is observable). Under the additional assumption that every claim is developed within at most \({\bar{\tau }} \) periods, we have that \(\xi ^{(i)} (S_m(k) \cap R_{0:(k+{\bar{\tau }} ) p}) = \xi ^{(i)} (S_m(k))\). Hence, if Assumption 3.4 is met for \(p>0\) specified above, the left-hand side of (13) can be written asOn the other hand, for the expression on the right-hand side of (13), we observe that \(\xi ^{(i)} (S_m(k) \cap R_{0:{\bar{\tau }} p}) = \xi _r^{(i)}(S_m(k))\) is the reported number of claims from \(\mathfrak {Y}_m\) with accident period k. Hence, combining the previous equations with (13), we obtain thatIn view of the basic Poisson assumption on \(\xi ^{(i)}\) from Definition 2.1 (and hence on the reported and unreported counterparts \(\xi ^{(i)}_r\) and \(\xi ^{(i)}_{nr}=\xi ^{(i)}-\xi _{r}^{(i)}\) from Observation Scheme 2.2), we obtain thatThis can be used to estimate the unknown claim intensities on \(\mathfrak {Y}_m\), i.e.,Indeed, let \({\mathcal {G}}\) denote a set of MLPs \(g:\mathfrak {X}\times [0,\infty ) \rightarrow [0, \infty )\) as in Model 3.5. We assume that the claim intensity on \(\mathfrak {Y}_m\) satisfies, for some \(g^{\mathfrak {Y}_m} \in {\mathcal {G}}\),where \(k_t=\lfloor \tfrac{t}{p}\rfloor \) indicates that the \(k_t\)-th period of length p contains t and where \(t_k= kp+\tfrac{p}{2}\) denotes the mid-point of the kth period.

As in (9), we arrive at the per-triangle losswhereAs in (10), \( \textrm{ex}^{\text {CL}, \mathfrak {Y}_m}(x^{(i)}, k) \) may now be interpreted as the expected number of claims that have occurred in the accident period \(I_p(k)\) and are reported by calendar time \(\tau \) for a policy with constant claim intensity \( \lambda ^{\mathfrak {Y}_m}(x^{(i)}, t_k) = 100\%\). In practice, the chain ladder factors within the loss must be estimated, for which we apply the well-known estimatorsNote that in contrast to the micro-level approach from the previous sections, there is no explicit model for the distribution of claim features on \(\mathfrak {Y}\). If \(g^{\mathfrak {Y}_m}=\hat{\lambda }^{\text {CL}, \mathfrak {Y}_m}(x, t)\) are maxima of the per-triangle losses \(\mathcal {L}^{\text {CL}, \mathfrak {Y}_m}\), the triangle-level claim intensity estimates can be aggregated to a common intensity estimatorFinally, exploiting \(\mathbb {E}[\xi ^{(i)}(S) \mid \xi _r^{(i)}(S)] = \xi _r^{(i)}(S) + \mathbb {E}[\xi _{nr}^{(i)}(S)]\) similar as in Sect. 4, we may define an IBNR-predictor as follows: recalling \(S_m(k) = I_p(k)\times \mathfrak {Y}_m \times [0,\infty )\), the claims from \(\mathfrak {Y}_m\) occurring in period k, letRecalling the notation \(R_{\tau _0:\tau _1} {:}{=}\{(t, y, d): \tau _0 < t + d \le \tau _1 \}\), similar derivations show that this predictor can also be extended to a predictor for reporting times \(\tau _0 = {\bar{\tau }}_0 p\) to \(\tau _1 = {\bar{\tau }}_1 p\) with \({\bar{\tau }}_0<{\bar{\tau }}_1 \in \mathbb {N}_0 \cup \{+\infty \}\):Here, we define the empty product as 1, i.e., \(\widehat{\textrm{FtU}}^{\text {CL}, \mathfrak {Y}_m}_j {:}{=}1\) for all \(j \le 0\). Finally, it is worthwhile to mention that in contrast to the predictor described in Sect. 4, it is not possible to use the chain ladder networks to obtain predictions for arbitrary \(\tau _0, \tau _1\) that are not whole multiples of p.

The quality of competing predictors may be assessed by suitable error measures. In this section, we define two such measures: an individual mean squared prediction error, and an aggregated global mean squared prediction error.

We start by considering the individual error measure. For \(S = [0, \tau ) \times \mathfrak {Y}' \times [0,\infty ) \subset [0,\infty ) \times \mathfrak {Y}\times [0,\infty )\) and \(0 \le \tau _0 < \tau _1 \le \infty \), let \({\hat{N}}^{(i)}_{\tau _0:\tau _1}(S)\) denote individual claim count predictions for \(N_{\tau _0:\tau _1}^{(i)}(S) = \xi ^{(i)}(S \cap R_{\tau _0:\tau _1})\), the number of claims in S incurred by policy \(x^{(i)}\) that are reported between \(\tau _0\) and \(\tau _1\); recall \(R_{\tau _0:\tau _1}=\{(t, y, d): \tau _0 < t + d \le \tau _1 \}\). Let \(q>0\) denote an evaluation period length (for instance, \(q=365\) corresponding to a year; note that there should be no confusion with the running index q used in Sect. 3.2), which is assumed to be a divisor of the total observation length \(\tau \) from Observation Scheme 2.2, and let \(\mathfrak {Y}' \subset \mathfrak {Y}\) denote an evaluation set of claim features. We then definewhere, for \(\ell \in \mathbb {N}_0\), recalling the notation \(I_q(x, \ell ) = C(x) \cap [\ell q, (\ell +1) q )\) for the covering time of policy x within the \(\ell \)th period of length q,Note that in practice the measure can only be calculated for \(\tau _1 \le \tau \) (with \(\tau \) the most recent date for which data is available) and on selected tests sets (for instance, in a back-testing approach). In controlled simulation experiments, see Sect. 7, we may and will use \(\tau _1=\infty \), thereby aiming at predicting the total number of unreported claims for each policy. Moreover, for using the error measures with the predictors from Sect. 4 and 5, the evaluation period q must be a multiple of the homogeneity period length p (unless one is willing to use the extension discussed in Remark 4.1).

The quality of individual claim count predictors may alternatively be assessed by first aggregating the individual predictions and then using standard global error measures; the predictors may then even be compared with classical methods for aggregated data like the standard chain ladder approach. Aggregated predictions are obtained from individual predictions straightforwardly: for \(A = \mathfrak {X}' \times S\) with policies from \(\mathfrak {X}' \subset \mathfrak {X}\) and claims from S as above, letwhich is to be considered a predictor for the aggregated claim numberFor q and \(\mathfrak {Y}'\) as in (15), we then definewhere \(A_{q, \ell , \mathfrak {Y}'} {:}{=}\mathfrak {X}\times [\ell q, (\ell + 1) q) \times \mathfrak {Y}' \times [0, \infty )\) comprises all claims in the portfolio from \(\mathfrak {Y}'\) that have occurred in the \(\ell \)th period of length q. Note that this measure has also been used in [7], Formula (20). Its application is limited to the constraints mentioned above for \(\tau _1\) and q.

In this section, we will study the performance of the new estimators and predictors within nine different simulation scenarios taken from [7]. We start by restating a brief, partially verbatim summary of the simulation models taken from the last-named paper:

The underlying portfolios build upon the car insurance data set described in Appendix A in [23]. The latter data set provides claim counts for 500,000 insurance policies, where each policy is associated with the risk features

which correspond to age of driver, age of car, power of car, fuel type of car, brand of car, and area code, respectively; see also (A.1) in [23] for further details. Next to that, the data set also provides the variable \(\texttt{truefreq}\), which corresponds to the claim intensity \(\lambda (x)\) in our model.

Each portfolio is considered over ten periods of 365 days, that is, the portfolio coverage period is the interval [0, 3650]. The different scenarios are as follows:

The baseline scenario. The baseline scenario/portfolio is characterized by a homogeneous exposure as well as position-independent claim intensity, occurrence process, and reporting process. It may be considered the vanilla portfolio that practitioners often aim at by careful selection of considered risks and suitable transformations, e.g., adjustment for inflation. More precisely:Eight non-homogeneous scenarios.

Eight non-homogeneous scenarios are obtained by altering a single element of the baseline scenario: Figure 3 illustrates the effect of the different scenarios on exposure, claim counts and reporting delays. The precise functional relationships are documented in Appendix D in [8].

The simulation study was conducted with 50 data seeds for each of the 9 scenarios, i.e., with 450 simulated portfolio datasets in total.

In this section, we will apply the different methods to a large real dataset containing motor legal insurance claims provided by a German insurance company. The dataset is described in Sect. 8.1. More detail on the prediction methods and estimation procedure can be found in Sect. 8.2. Due to the nature of real world data, observations are only available for a limited time frame. Therefore, model performance metrics cannot use \(\infty \) as the time of evaluation, but must instead use a finite cutoff date. We examined two artificial truncation points, \(\tau = \textrm{31st December 2017} and \tau = \textrm{31st December 2018}\) and evaluate predictions for one year into the future, i.e. \(\text {RMSE}_{\tau :\tau + 365}(\mathfrak {Y}, 365)\) and \(\text {RMSE}^{\text {expo}}_{\tau :\tau + 365}(\mathfrak {Y}, 365)\). Results of this examination are presented and discussed in Sect. 8.3.

Two new methods for joint prediction of micro-level IBNR claim counts and claim frequencies have been developed and applied to real and simulated data. Results show promising accuracy on an exposure level compared to the theoretical optimum under laboratory conditions. The new methods also permit assigning IBNR claim count predictions on a policy level such that policies without claims can receive a non-zero IBNR prediction, which is an advantage for analysis of small sub-portfolios where applying chain ladder estimation factors—even with parameters estimated on a larger dataset - produce highly volatile estimates. The presented case studies uncover several opportunities for further research: