Stochastic loss reserving with mixture density neural networks

Probabilistic forecasts of outstanding claims are essential for capital allocation and optimising the risk margins set when reporting liabilities and pricing products. In this paper, we explain how such forecasts can be obtained with the Mixture Density Network (“MDN”) in a flexible way, working with aggregate loss triangles. The MDN is a design developed by Bishop (1994) (see also Bishop, 2006), with applications found in many fields, such as financial modelling (Ormoneit and Neuneier, 1996), acoustics (Zen and Senior, 2014) and electrical engineering (Vossen, Feron and Monti, 2018). The essential idea is that the parameters of the distribution will be estimated by the NN architecture, allowing for heterogeneity. Related works in the (non mixture) Gaussian setting can be found in Nix and Weigend (1994) and Lakshminarayanan et al. . The MDN assumes that the target output follows a mixture (typically Gaussian) distribution, allowing effectively a flexible distribution to be fit to the data. In this paper, we will assume that incremental claims follow a mixture Gaussian or mixture Log-Gaussian distribution.

A special case of an MDN was used by Kuo (2020) who mixed a shifted lognormal with a degenerate distribution to represent (positive) and zero cashflows for individual claims, respectively. In contrast, in our application we utilise a mixture of Gaussian distributions, whereby the number of components as well as all parameters are kept flexible.

In loss reserving, previous literature has fitted members of the exponential family using NNs. For instance, Gabrielli, Richman and Wüthrich (2020) considered aggregate claims, and Delong et al. (2021); Gabrielli (2021) considered individual claims. However, they typically focused on the central estimate rather than distributional properties.

A typical approach is to utilise the NN to replace the linear component in a GLM setup (see also Denuit and Trufin, 2020, for additional discussions), including the GLM-based assumption of a homogeneous mean-variance relationship. Such an assumption is not as flexible for the modelling of volatility, as compared to the mixture Gaussian fit by the MDN, which allows for heterogeneity for each random variable. In addition, the mixture Gaussian, given sufficient components, can approximate any distribution within a desired accuracy (see Nguyen and McLachlan, 2019, for details). As the MDN considers a wide range of distributions of varying volatility and shape, the training and fitting process of the network is akin to distribution selection. Another benefit of MDNs is that the central estimates can be derived directly from the fitted mixture Gaussian parameters, which means that location and shape are fitted simultaneously under the single model.

In addition to analysing central estimates, we use qualitative and quantitative measures to assess the distributional and quantile accuracy of the fitted MDN forecasts; see Section 4.4. Overall, we show that the MDN’s flexibility yields more accurate probabilistic forecasts when compared to the cross-classified over-dispersed Poisson (“ccODP”), in both simulated and real data; see also Section 1.2.4.

Neural networks are highly flexible in their design. Since different designs produce varying results, it is vital to have a clear methodology for testing between these different model designs.

Commonly, the data is partitioned into training, validation and testing sets. The network is trained on the training set until the validation loss is minimised, then projected onto the testing set. The model producing the lowest test error is preferred.

To test between different models, a loss triangle must be split into training, validation and testing sets. Tashman (2000) and Bergmeir and Benítez (2012) explore two popular methods; fixed origin and rolling origin, used to partition sequential data. Balona and Richman (2020) recently apply the rolling origin method to loss triangles, in order to compare between different traditional reserving models.

This paper contributes by performing the rolling origin data partition exclusively within the aggregate loss triangle, using it to select neural network designs. Taking the latest calendar periods for testing has been done by Kuo (2019); Ramos-Pérez, Alonso-González and Núñez-Velázquez (2019), which was facilitated by simultaneously training the NN on multiple triangles. This paper performs model testing, selection and fitting using only one triangle at a time. Furthermore, the neural network model testing framework implemented in this paper outlines a validation set, which is essential in training NNs.

Furthermore, this paper contributes by implementing a model searching and selection algorithm to loss triangle reserving, which methodically searches and selects different network features such as the number of layers, nodes, and components. This algorithm extends the methodology implemented by Gabrielli, Richman and Wüthrich (2020) by also searching for the best-performing regularisation coefficients. This algorithm also makes NN design selection more methodical, requiring less expert input and assisting the popularisation of NNs in practice.

While the current neural network applications to loss reserving have shown their accuracy and flexibility in modelling, there are practical and technical obstacles which have hindered the acceptance and use of neural networks by the actuarial community (Rossouw and Richman, 2019; Wüthrich and Merz, 2019).

Producing interpretable forecasts is important for justifying business decisions to stakeholders. The neural network’s lack of interpretability has contributed to its lack of acceptance by actuaries. Wüthrich and Merz (2019) tackle this by adopting the resnet design from He et al. (2016) to form the GLM-NN hybrid CANN (Combined Actuarial Neural Network) architecture.

In this paper, we adapt the MDN to the above approach, to create the ResMDN. While Gabrielli, Richman and Wüthrich (2020), Poon (2019) and Richman and Wuthrich (2021) have adopted a similar methodology in loss reserving, the ResMDN provides the additional distributional flexibility of a (heterogeneous) mixture Gaussian distribution assumption. It boosts all parameters of the mixture Gaussian simultaneously and directly within one network, while maintaining a fixed GLM initialisation for interpretability.

We show that the ResMDN can successfully boost the embedded GLM, the ccODP, improving structural deficiencies in both mean and volatility estimates where visible, while maintaining the interpretable GLM backbone in its forecasts.

As a further consideration, the NN’s black box modelling can cause long range forecasts to become unstable (Rossouw and Richman, 2019). In other words, it may not be adequate to project the behaviour of the highly flexible function fit to the upper triangle to the lower triangle. In this paper, central estimate projections are able to be explicitly constrained by the actuary, forcing the neural network to only fit functions which produce reasonable mean forecasts. This framework thus also allows actuarial judgement to be explicitly incorporated into the modelling process.

Recent applications of machine learning to loss reserving have been mainly focused on more granular claims datasets, such as individual claims data. This article, however, works exclusively with aggregate loss triangles. Developing a neural network model that is applied exclusively to loss triangles is highly relevant, as it enhances comparability and interpretability of results (as opposed to traditional reserving), and some insurers may still have limited data, preventing them from applying individual loss reserving methodologies. Furthermore, it is not obvious that machine learning methods can replicate their documented accuracy when presented with limited data, which is an interesting research question in itself.

In this paper, we demonstrate how neural networks can improve on traditional models, even with limited triangular data:

• —

  Data scarcity: The MDN was applied by Kuo (2020), but it was used in Individual Claims modelling. The current paper shows that the MDN also finds success when applied to loss triangles.

• —

  Model selection: Rossouw & Richman (2019) outline a fixed origin data partition (see Bergmeir & Benitez (2012) for more detail), but use more granular data, hence this methodology hasn’t been tested on a sparse loss triangle. While Gabrielli et al (2020) works with loss triangles, individual claims data are used to partition into training and validation sets. Balona & Richman (2020) recently apply the rolling origin model validation methodology (see Bergmeir & Benitez (2012) for more details) exclusively to loss triangles. However, as only non-machine learning models were tested, no validation set was constructed in the paper.

  With the rolling origin methodology and model searching algorithm, the neural network designs chosen in this paper produce robust, smooth, and accurate central and distributional forecasts, showcasing their combined effectiveness. Given that neural networks can perform unreliably with small datasets, these results show that the rolling origin partition is feasible for loss triangles of a certain size (40 $\times$ 40 triangles were used in this paper), which should encourage further neural network modelling with loss triangles.

• —

  Data variety: Additionally, for NNs to prove their practicality and reliability, they must provide accurate results in a variety of different triangles (of instance, in a range of claim development situations which would invalidate the chain ladder assumptions). The current paper contributes to the literature by testing the MDN on a variety of environments of varying complexity and specifications, using both simulated and real data. Some of the environments tested are derived in concept from triangles used by Harej et al. (2017); Gabrielli et al. (2020). The MDN achieved superior probabilistic forecasts relative to the ccODP in all environments tested (see Section 3.3 for details).

Incremental claims $X_{i,j}$ are assumed to follow a mixture Gaussian distribution

$$f_{\hat{X}_{i,j}}(x)=\sum_{k=1}^{K}\alpha_{i,j,k}\phi({x|\mu_{i,j,k},{\sigma_{i,j,k}}})$$

(2.3)

With that distributional assumption, the output layer of the MDN estimates the parameters of the mixture distribution, $(\boldsymbol{\alpha},\boldsymbol{\mu},\boldsymbol{\sigma}$), which are used to form a mixture Gaussian density. A Negative Log Likelihood (NLL) loss function

$$NLLLoss(\textbf{X},\hat{\textbf{X}}|\textbf{w})=-\frac{1}{|\textbf{X}|}\sum_{i,j:X_{i,j}\in Train}ln(f_{\hat{X}_{i,j}}(X_{i,j}|\textbf{w}))$$

(2.4)

is used to train the MDN, where X is the set of cells $X_{i,j}$ in the training set, $|\textbf{X}|$ is the cardinal of X, $\hat{\textbf{X}}$ is the set of predicted distributions of $X_{i,j}$ and w is the set of weights in the MDN.

The MDN isn’t structurally restricted to fitting mixture Gaussians to the response; it can estimate the parameters of any desired distribution so long as the loss function is specified accordingly. In this paper, only the mixture Gaussian framework was considered, with the output layer estimating the $(\boldsymbol{\alpha},\boldsymbol{\mu},\boldsymbol{\sigma})$ parameters.

As Bishop (1994) notes, the mixture Gaussian distribution, given a sufficient number of components and hidden layers, is capable of approximating any desired distribution within a desired accuracy. Mixture densities with more components will certainly be more flexible, however, practical obstacles such as over-parametrisation and data insufficiency will limit the range of distributions fit by the MDN.

Alternatively, while maintaining a mixture Gaussian output layer, a mixture Log-Gaussian can also be fit to $X_{i,j}$ by fitting a mixture Gaussian to $ln(X_{i,j})$—This holds by the definition of a mixture random variable; see Appendix A for details). This distribution helps to address the practical limitations to the flexibility of the mixture Gaussian by providing a positive, heavier-tailed option. Furthermore, taking the log of incremental claims linearises the data, which can make training simpler and more efficient. Both the mixture Gaussian and mixture Log-Gaussian distributions achieved impressive results, analysed in Section 5.

Figure 1 provides a basic visualisation of the MDN’s design. Mixture Density Networks differ from other neural networks due to their output layer, which estimates a mixture distribution, commonly mixture Gaussian, to the response variable. For a detailed overview of neural networks, their design, mechanism and terminology, see, for instance, Richman (2020). In the case of aggregate loss triangles that we consider, the input variables of the MDN are $i$ and $j$, the accident and development periods, respectively. These variables are passed through a fully connected hidden layer, which consists of hidden layers, each layer consisting of neurons. Each neuron in a hidden layer takes a weighted sum of the previous layer’s output, before passing it through an activation function. The final hidden layer’s output is then passed to the output layer, which produces the desired distribution parameters.

Let $w_{a,c}^{l}$ be the weight parameter connecting the $a^{th}$ neuron in the $l^{th}$ layer to the $c^{th}$ neuron in the $(l+1)^{th}$ layer. Define $b_{a}^{l-1}$ as the bias term added to the $a^{th}$ neuron in the $l^{th}$ layer. Weighted sums of the inputs, $(i,j)$, are passed into the first hidden layer, before an activation function $g$ yields such output for the $p^{th}$ node in that layer:

$$\textbf{z}_{i,j,p}^{1}=g(i\times w_{1,p}^{0}+j\times w_{2,p}^{0}+b_{p}^{0}).$$

(2.5)

We assume $L$ hidden layers in the MDN, with $D$ nodes in each layer. Each node in successive hidden layers take a weighted sum of the output from nodes in the previous layer, such that $\textbf{z}_{i,j,p}^{L}$, the $p^{th}$ node in the final hidden layer, is calculated as

$$\textbf{z}_{i,j,p}^{L}=g(\sum_{d=1}^{D}w_{d,p}^{L-1}\textbf{z}_{i,j,d}^{L-1}+b_{p}^{L-1}).$$

(2.6)

The output layer is split into three sections, each with $K$ nodes, $K$ being the number of components in the Mixture Density. Let’s call these sections the alpha, mu and sigma sections, respectively. Similarly to the hidden layers, each node in the output layer takes a weighted sum of the output of all nodes in the last hidden layer, Layer $L$. The weighted sums are then passed through different activation functions for each section to yield the final output of the MDN, $(\boldsymbol{\alpha},\boldsymbol{\mu},\boldsymbol{\sigma})$. Specifically:

Alpha:

The output of the $k^{th}$ nodes of the alpha section

$$\textbf{z}_{i,j,k}^{\alpha}=\sum_{d=1}^{D}w_{d,k}^{L,\alpha}\textbf{z}_{i,j,k}^{L}+b_{k}^{L,\alpha}\text{, for }k=1,2,..K$$

(2.7)

leads to

$$\alpha_{i,j,k}=\frac{e^{\textbf{z}_{i,j,k}^{\alpha}}}{\sum_{k=1}^{K}e^{\textbf{z}_{i,j,k}^{\alpha}}}.$$

(2.8)

Note that the output $\textbf{z}_{i,j,k}^{\alpha}$ was passed through a Softmax activation function, which ensures that $\sum_{k=1}^{K}\alpha_{i,j,k}=1$.

Mu:

Similarly,

$$\textbf{z}_{i,j,k}^{\mu}=\sum_{d=1}^{D}w_{d,k}^{L,\mu}\textbf{z}_{i,j,k}^{L}+b_{k}^{L,\mu}\text{, for }k=1,2,..K$$

(2.9)

leads to

$$\mu_{i,j,k}=\textbf{z}_{i,j,k}^{\mu}$$

(2.10)

as there are no constraints on the mu layer which would require an activation function. Bishop (1994) notes that such a design represents an ‘un-informative prior’ on $\mu$, which befits the lack of constraints on the mean.

Sigma:

Finally, the sigma output

$$\textbf{z}_{i,j,k}^{\sigma}=\sum_{d=1}^{D}w_{d,k}^{L,\sigma}\textbf{z}_{i,j,k}^{L}+b_{k}^{L,\sigma}\text{, for }k=1,2,..K$$

(2.11)

is passed through an exponential function,

$${\sigma_{i,j,k}}=e^{\textbf{z}_{i,j,k}^{\sigma}},$$

(2.12)

which ensures the standard deviation is always positive (Hjorth and Nabney, 2000).

Here, $w_{d,k}^{L,\alpha},w_{d,k}^{L,\mu},w_{d,k}^{L,\sigma}$ are the weights connecting the output of node $d$ in layer $L$ to node $k$ in the alpha, mu and sigma layers, respectively. Thus, for each input cell $(i,j)$, a unique combination of parameters,

$$(\alpha_{i,j,1},\alpha_{i,j,2}.....\alpha_{i,j,K},\mu_{i,j,1},\mu_{i,j,2}.....\mu_{i,j,K},\sigma_{i,j,1},\sigma_{i,j,2}.....\sigma_{i,j,K}),$$

is produced in the output layer of the MDN, which then generates the probability density for $\hat{X}_{i,j}$

$$f_{\hat{X}_{i,j}}(x)=\sum_{k=1}^{K}\alpha_{i,j,k}\phi(x|\mu_{i,j,k},\sigma_{i,j,k}).$$

(2.13)

The ResMDN embeds an approximation of the Cross-Classified Over-Dispersed Poisson (ccODP) model (see Section 2.2 for more detail), which follows the distribution outlined in (2.1). In neural network terminology, an embedding is the mapping of a discrete or categorical input variable into a numerical vector, which is then fed into the network (Richman, 2020). Since the output of the ResMDN takes the form of parameters for a mixture Gaussian distribution, the GLM initialisation’s density is approximated by a mixture Gaussian, spread evenly over $K$ components. The parameters of this approximation are fed into the ResMDN as embeddings of the GLM. The distribution $f_{\hat{X}_{i,j}^{ccODP}}(x)$ of $X_{i,j}$ as estimated by the ccODP is approximated by

$$f_{\hat{X}_{i,j}^{ccODP}}(x)\approx\sum_{k=1}^{K}\alpha^{GLM}_{i,j,k}\phi_{i,j,k}(x|{\mu^{GLM}_{i,j,k},\sigma^{GLM}_{i,j,k}}),$$

(2.14)

where

$$\alpha_{i,j,k}^{GLM}=\frac{1}{K},\quad\mu_{i,j,k}^{GLM}=E[\hat{X}_{i,j}^{ccODP}]=A_{i}B_{j},\quad\sigma_{i,j,k}^{GLM}=\sqrt{Var[\hat{X}_{i,j}^{ccODP}]}=\sqrt{DA_{i}B_{j}}.$$

(2.15)

The ResMDN’s structure resembles the MDN very closely. The Input Layer of the ResMDN consists of the accident and development periods, $(i,j)$, as well as a unique categorical integer, $c_{i,j}=40*(i-1)+j$, which allows the ResMDN’s embedding layer to identify the specific cell $(i,j)$ and produce the corresponding GLM loss estimate for that cell as output. Hence, each cell $(i,j)$ is assigned a number from $(0,1599)$. The variables $(i,j)$ are passed through an MDN which excludes the activations of the final output layer. An embedding layer takes the categorical input $c_{i,j}$ and produces as output:

$$\bigg{(}ln(\boldsymbol{\alpha^{GLM}_{i,j}}),\boldsymbol{\mu^{GLM}_{i,j}},ln(\boldsymbol{\sigma^{GLM}_{i,j}})\bigg{)}$$

(2.16)

The outputs from both the fully connected component and embedding layer are added together, before the Softmax and exponential activations are applied to the alpha and sigma additions, respectively. The final output of the ResMDN consists of the mixture Gaussian parameter estimates, $(\boldsymbol{\alpha^{ResMDN}},\boldsymbol{\mu^{ResMDN}},\boldsymbol{\sigma^{ResMDN}})$. Figure 2 provides a visualisation of the ResMDN model. The key design features of the ResMDN are outlined below:

• 1.

  Fully Connected Module: Let $\textbf{z}_{i,j,k}^{\alpha},\textbf{z}_{i,j,k}^{\mu},\textbf{z}_{i,j,k}^{\sigma}$ be as described in (2.7), (2.9) and (2.11), that is, the MDN’s output before the output layer activations are applied. The fully connected module of the ResMDN performs the function

  $$(i,j)\mapsto\{(\textbf{z}_{i,j,k}^{\alpha},\textbf{z}_{i,j,k}^{\mu},\textbf{z}_{i,j,k}^{\sigma}),k=1,2,..K\},$$

  (2.17)

• 2.

  Embedding Layer: The embedding layer weights are pre-set to provide the mapping

  $$c_{i,j}\mapsto\bigg{(}ln(\boldsymbol{\alpha^{GLM}_{i,j}}),\boldsymbol{\mu^{GLM}_{i,j}},ln(\boldsymbol{\sigma^{GLM}_{i,j}}\bigg{)}.$$

  (2.18)

  The log of the $\boldsymbol{\alpha}$ and $\boldsymbol{\sigma}$ parameters are produced in the embedding layer, since the Softmax and exponential Activation functions will take the exponent in the output layer nodes.

• 3.

  Addition and Final Activation: The output from the embedding and fully connected modules are added together element-wise:

  	Addition:	$\displaystyle\bigg{(}ln(\alpha^{GLM}_{i,j,k}),\mu^{GLM}_{i,j,k},ln(\sigma^{GLM}_{i,j,k}),\textbf{z}_{i,j,k}^{\alpha},\textbf{z}_{i,j,k}^{\mu},\textbf{z}_{i,j,k}^{\sigma}\bigg{)}$
  		$\displaystyle\mapsto\bigg{(}ln(\alpha^{GLM}_{i,j,k})+\textbf{z}_{i,j,k}^{\alpha},\mu^{GLM}_{i,j,k}+\textbf{z}_{i,j,k}^{\mu},ln(\sigma^{GLM}_{i,j,k})+\textbf{z}_{i,j,k}^{\sigma}\bigg{)}$
  		$\displaystyle\quad=\bigg{(}ln(\alpha^{GLM}_{i,j,k})+\textbf{z}_{i,j,k}^{\alpha},\mu^{ResMDN}_{i,j,k},ln(\sigma^{ResMDN}_{i,j,k})\bigg{)},$		(2.19)

  before the Softmax and exponential activations are applied to the alpha and sigma layers, respectively:

  	Final Activations:	$\displaystyle\bigg{(}ln(\alpha^{GLM}_{i,j,k})+\textbf{z}_{i,j,k}^{\alpha},\mu^{ResMDN}_{i,j,k},ln(\sigma^{ResMDN}_{i,j,k})\bigg{)}$
  		$\displaystyle\mapsto\left(\frac{{\alpha^{GLM}_{i,j,k}e^{\textbf{z}_{i,j,k}^{\alpha}}}}{\sum_{k=1}^{K}\alpha^{GLM}_{i,j,k}e^{\textbf{z}_{i,j,k}^{\alpha}}},\mu^{ResMDN}_{i,j,k},e^{ln(\sigma^{ResMDN}_{i,j,k})}\right)$
  		$\displaystyle\quad=\left(\alpha^{ResMDN}_{i,j,k},\mu^{ResMDN}_{i,j,k},\sigma^{ResMDN}_{i,j,k}\right).$		(2.20)

  Hence, the boosted mixture Gaussian parameters are produced in the Output Layer.

• 4.

  Initialisation: The activations, $\textbf{z}_{i,j,k}^{\alpha},\textbf{z}_{i,j,k}^{\mu},\textbf{z}_{i,j,k}^{\sigma}$, are generated using the parameters, $(\textbf{w}_{L},\textbf{b}_{L})$, defined in (2.7)–(2.12). These parameters, representing the weights in the final hidden layer, are initialised at 0, such that

  $$\textbf{z}_{i,j,k}^{\alpha},\textbf{z}_{i,j,k}^{\mu},\textbf{z}_{i,j,k}^{\sigma}=0,\alpha_{i,j,k}^{ResMDN}=\alpha_{i,j,k}^{GLM},\mu_{i,j,k}^{ResMDN}=\mu_{i,j,k}^{GLM},\sigma_{i,j,k}^{ResMDN}=\sigma_{i,j,k}^{GLM},$$

  hence producing the GLM approximation in the Output Layer at the initialisation of the ResMDN. This initialisation follows the methodology of Gabrielli, Richman and Wüthrich (2020) closely.

During training, the embedding layer maintains constant output, while the fully connected module adjusts its weights to capture non-linearities which the GLM has missed. The ResMDN’s overall function at the termination of training is

	$\displaystyle(i,j,c_{i,j})$	$\displaystyle\mapsto\left(\frac{{\alpha^{GLM}_{i,j,k}e^{\textbf{z}_{i,j,k}^{\alpha}}}}{\sum_{k=1}^{K}\alpha^{GLM}_{i,j,k}e^{\textbf{z}_{i,j,k}^{\alpha}}},\mu^{GLM}_{i,j,k}+\textbf{z}_{i,j,k}^{\mu},{\sigma^{GLM}_{i,j,k}e^{\textbf{z}_{i,j,k}^{\sigma}}}\right)$
		$\displaystyle\quad=\left(\alpha^{ResMDN}_{i,j,k},\mu^{ResMDN}_{i,j,k},\sigma^{ResMDN}_{i,j,k}\right)\text{, for }k=1,2,...,K.$		(2.21)

The NN boosting terms are relatively easy to analyse in relation to the GLM fit, especially the mean and volatility terms. Furthermore, the black box neural network modelling is only applied to the residuals, meaning the lack of interpretability is restricted to that domain only. Hence the ResMDN improves the interpretability of the model compared to the MDN.

Thanks to the flexibility offered by SynthETIC (see Avanzi, Taylor, Wang and Wong, 2021, for a description of the R simulation package), four different claim environments were simulated, of various features and complexities. Simulating different environments was done to test the MDN’s versatility and ability to capture complex trends and produce accurate forecasts in a variety of controlled, challenging environments. As a loss triangle is a collection of random variables, it is important to run the MDN on a large sample of triangles to gain a better understanding of its accuracy and also test its ability to provide consistent results. Hence, for each simulated environment, 50 independent triangles were simulated (of size $40\times 40$), leading to 200 triangles. Namely, the scenarios, or environments, were:

1. 1.

   Environment 1 - Simple, short tail claims: This environment simulates short tail claims which are homogeneous in composition for all accident quarters. The reporting and settlement delays have been approximately calibrated to show similar characteristics to the simulator developed by Gabrielli and Wüthrich (2018). Figure 4 plots the incremental claims for Environment 1. A spike in claim payment in DQ2 complicates this dataset, but such a feature is not uncommon in practice. This environment was a preliminary test to the feasibility of MDNs in modelling 40x40 triangles and producing reasonable results.

   

2. 2.

   Environment 2 - Increase in claim processing speed: A gradual shift from long tail to short tail claims along accident quarters is simulated, that is, an increase in claims processing speed. Initially, there are more long tail claims, however, the proportion of these claims decreases, while the proportion of short claims increases. From the incremental claim plot shown below, later AQs see higher losses early on, due to the increasing proportion of short tail claims. Figure 5 plot the incremental claims of environment 2. The main question to be answered in testing this dataset is, given the systematic volatility in the claims data, can the MDN accurately distinguish between systematic and unsystematic volatility and capture the distribution of data points accurately? That is, will the MDN learn that claims are getting shorter, or will it attribute the trend to noise?

   

3. 3.

   Environment 3 - Inflation shock: Superimposed inflation is changed instantly from 0% to 8% per annum, starting at AQ30. The 8% inflation remains constant in the lower triangle. This environment tests the ability of the MDN to recognise changes in calendar effects and adapt projections accordingly.

   Only the last 10 calendar quarters in the upper triangle contain information regarding the inflation shock, which increases the difficulty for the MDN. A further complication is that the rolling origin partition contains little training points featuring the change in inflation, hence this environment assesses its ability to capture recent trends. Figure 6 plots the incremental claims.

   

4. 4.

   Environment 4 - High Systematic Complexity: This environment is the default triangle generated by the SynthETIC simulator, which was designed to mimic features seen in real data (see Avanzi, Taylor, Wang and Wong, 2021, for details). Complex dependencies exist between claim size, reporting and settlement delay, and superimposed inflation. Settlement delay, which depends on claim size, declines over the first 20 AQs. Superimposed inflation is up to 30%, but declines for larger claims. A legislative change at AQ20 causes small claims to face a reduction in size and settlement speed until AQ30. The general trend can be summarised by slow development, high volatility and high superimposed inflation. The volatility is primarily caused by the low claim frequency and highly volatile severity. Hence, it is normal for claims in one AQ to follow a completely different pattern (reporting, settlement, volume, development pattern) than claims in the adjacent AQ. Figure 7 plots the incremental claims. This environment assesses the MDN’s ability to produce accurate forecasts in a volatile environment.

   

We apply the MDN to a dataset obtained through a collaborative project between Allianz, University of New South Wales (UNSW), Suncorp and Insurance Australia Group (IAG). This forms the AUSI acronym, which we will use to refer to this dataset. The Auto Bodily Injury line of business is used, which features slow claim development and high volatility. The AUSI dataset consists of transactional data for individual claims, which we aggregated into quarterly triangles. We use quarterly data from January 2005 to December 2014, which provides a 36$\times$36 upper triangle and a 4 quarter forecasting period, which will be used to compare the MDN’s results to the ccODP. Using the individual claims data, each claim was randomly allocated to one of ten triangles, meaning that ten aggregate triangles of roughly equal size were created using this dataset. As mentioned earlier, running the MDN on multiple triangles better assessed the model’s consistency. Figure 8 plots the incremental claims for this dataset. This environment aims to assess the MDN’s ability to provide accurate forecasts for real data.

The MDN produced excellent central estimate projections in all environments. Where the environment had more structural heterogeneity, i.e. where the ccODP assumptions are not satisfied, the MDN decisively outperformed the ccODP in all metrics. Several key observations can be made:

• 1.

  In environment 1, the data (by design) satisfies the ccODP assumptions well, hence the ccODP was very competitive. Nevertheless, the MDN slightly outperformed in all quantitative metrics when measuring the accuracy of incremental claims, $X_{i,j}$. This can be attributed to the smooth function fit by the MDN.

• 2.

  In environment 2, the MDN successfully learned that claims processing speed is increasing, predicting a sharper spike in claim payments in the later AQs. Figure 13 plots the results for this environment. The ccODP, assuming homogeneity in claim development, approximated claims as medium-tailed, leading to a clear over-estimation of claims in later AQs.

• 3.

  For environment 3, the MDN accurately captured the inflation shock at calendar quarter 30 (CQ30) onwards. Figure 14 plots the results. The ccODP did not keep pace with the increased inflation due to its limited ability to handle heterogeneity, leading to its under-estimation of claims from CQ30 onwards.

• 4.

  In both environment 4 and the AUSI environment, the MDN handled volatile data well and provided accurate central estimates, outperforming the ccODP. This accuracy is shown in Figure 15, which plots the MDN’s central estimates against the empirical mean based on 10 simulations.

Figure 16 provides boxplots of the MDN’s % reduction of the RMSE relative to the ccODP for 50 triangles in each of the environments tested (10 for the AUSI environment). The boxplots show that the MDN achieved a positive reduction in the RMSE (lower RMSE) relative to the ccODP for the majority of triangles in each environment, further showing the MDN’s higher forecasting accuracy.

Despite the MDN’s success, it showed weaknesses in several areas, some of which were mitigated:

• 1.

  Using the NLL loss function alone can encourage the MDN to over-estimate the volatility when its central estimate is inaccurate. This was seen in environment 2, where a MDN with a NLL loss function under-estimates claims in the (40,2) cell, leading to an excessively high volatility estimate for that region. In addition to sigma activity regularisation (see Section 3.2, adding an MSE term to the loss function helped to resolve this issue, as it encouraged the MDN to achieve more accurate central estimates. Hence, an MSE term was added to the loss function for environments 1 and 2, and is recommended for loss triangles with sharp shifts in claims development.

• 2.

  Taking the log of aggregate claims linearised the data, which often led to faster and more accurate modelling. In this paper, environments 1 and 3 were fit with a mixture Log-Gaussian, as they showed significantly more accurate results, especially capturing the claims decay in later DQs.

The Mixture Density Network produced very smooth and accurate volatility estimates. Where noise in the data was low, the MDN projected low volatility, and vice versa. Overall, it outperformed the ccODP qualitatively and quantitatively in estimating the volatility of individual cells.

In relation to individual cells, the MDN’s risk margin estimates at the 25th and 75th percentile were more accurate overall than the ccODP’s margins in almost all environments tested. The ccODP’s variance is a function of its mean, and hence it failed where the central estimates failed. For example, in environment 2 (Figure 17), the ccODP over-estimates claims in later AQs, which led to it over-estimating margins in that same period. In environment 3, the ccODP under-estimated claims in later calendar Quarters (CQs), as it didn’t effectively capture the inflation shock. This led to volatility estimates also being too low in those periods, as Figure 18 illustrates. The MDN dealt with these structural issues more effectively, leading to more accurate dispersion estimates of claims. In the volatile AUSI environment, the MDN produced smooth and accurate margin estimates, as Figure 19 illustrates.

While some correlation was found in the results between central and volatility forecasts, the MDN allows a large degree of independence between the $\mu$ and $\sigma$ parameters estimated, allowing it to fit a much wider range of distributions than the ccODP.

Figure 16 provides boxplots of the MDN’s increase in the log score relative to the ccODP for 50 triangles in each of environments 1,2,3 and 4. The AUSI boxplot was based on 10 triangles. The boxplots show that the MDN achieved a higher log score relative to the ccODP for the majority of triangles in each environment, indicating a more accurate probabilistic forecast.

With the MDN’s success, there are some weaknesses which must be addressed:

• 1.

  The MDN still showed signs of attributing noise to systematic trends. In environment 2, even though the DQ2 spike was fixed, the volatility was still too high for AQ30 and AQ40.

• 1.

  Similar to central estimates, the MDN often over-estimates volatility in later DQs, also due to a lack of data in that region.

The MDN provided more accurate 75% and 95% quantiles for all environments in the majority of triangles run for each environment. These results follow from the MDN’s ability to provide more accurate central estimates and volatility estimates. The quantile analysis was mainly quantitative, using the quantile scores. Table 1 confirms that in all environments, the MDN reduces the 75% and 95% quantile scores for the majority of triangles, indicating more accurate quantile estimates at the 75% and 95% levels.

The MDN and ccODP models were run on fifty triangles of each of environments 1,2,3,4, and the ten triangles partitioned from the AUSI data. The quantitative metrics are calculated for the 50 triangles and averaged, with results between the MDN and ccODP models compared in Table 1. As the table shows, the MDN, on average, had a lower RMSE and Quantiles Scores and had a higher log score for each environment, which is a significant out-performance by the MDN. Table 2 further reinforces these results by showing the percentage of triangles in which the MDN outperformed the ccODP for each quantitative metric. In each environment, the MDN outperforms the ccODP in each metric for the majority of triangles.

The MDN, in all environments except environment 1, showed more accurate central and dispersion estimates of total reserves compared to the ccODP estimate (the dispersion accuracy is qualitatively measured through visualising reserve density plots, Figure 20 plots these results). An empirical distribution of total reserves, based on many simulations, is used as an estimate of the actual distribution of $R$. This out-performance follows from the MDN modelling the mean and volatility of individual cells more accurately than the ccODP. Because the claims in environment 1 are homogeneous in development, the ccODP provides highly competitive results and the MDN didn’t outperform. However, for the more complicated environments 2,3 and 4, the MDN had more accurate 75% and 95% quantiles of total reserves compared to the ccODP. Table 3 calculated the quantitative metrics for both models for total reserve estimates, $\hat{R}$. The qualitative analysis (Figure 20) also supports these results.

While the ResMDN was able to detect the ccODP’s systematic errors, it generally failed to outperform the non-embedded MDN in the above examples as it only partially corrected the errors (see Figure 22 for an illustration). Nevertheless, the ResMDN showed great function and potential in correcting an embedded GLM’s mean and volatility estimates, while maintaining a fair level of interpretability. While the current implementation is not as accurate as the MDN, if this is criteria is essential then one possibility is to consider a more sophisticated GLM backbone.

Finally, it should also be noted that an additional benefit of the ResMDN is that the ResMDN’s training time was noticeably faster than the MDN. For environment 2, when a few triangles were analysed, the ResMDN finished training in 2218 epochs on average, compared to 3868 for the MDN. This 43% reduction in training time is due to the ResMDN’s GLM initialisation being more accurate than the MDN’s random starting fit, meaning that less parameter adjustment is required. This observation is in parallel to the findings of Gabrielli, Richman and Wüthrich (2020), who also noted faster training times under the CANN structured model.