Neural networks for insurance pricing with frequency and severity data:
  a benchmark study from data preprocessing to technical tariff
Insurers usually turn to generalized linear models for modeling claim
frequency and severity data. Due to their success in other fields, machine
learning techniques are gaining popularity within the actuarial toolbox. Our
paper contributes to the literature on frequency-severity insurance pricing
with machine learning via deep learning structures. We present a benchmark
study on four insurance data sets with frequency and severity targets in the
presence of multiple types of input features. We compare in detail the
performance of: a generalized linear model on binned input data, a
gradient-boosted tree model, a feed-forward neural network (FFNN), and the
combined actuarial neural network (CANN). The CANNs combine a baseline
prediction established with a GLM and GBM, respectively, with a neural network
correction. We explain the data preprocessing steps with specific focus on the
multiple types of input features typically present in tabular insurance data
sets, such as postal codes, numeric and categorical covariates. Autoencoders
are used to embed the categorical variables into the neural network, and we
explore their potential advantages in a frequency-severity setting. Model
performance is evaluated not only on out-of-sample deviance but also using
statistical and calibration performance criteria and managerial tools to get
more nuanced insights. Finally, we construct global surrogate models for the
neural nets' frequency and severity models. These surrogates enable the
translation of the essential insights captured by the FFNNs or CANNs to GLMs.
As such, a technical tariff table results that can easily be deployed in
practice.
Freek Holvoet, Katrien Antonio, Roel Henckaerts
# Neural networks for insurance pricing with frequency and severity data:

###### Abstract

Insurers usually turn to generalized linear models for modelling claim frequency and severity data. Due to their success in other fields, machine learning techniques are gaining popularity within the actuarial toolbox. Our paper contributes to the literature on frequency-severity insurance pricing with machine learning via deep learning structures. We present a benchmark study on four insurance data sets with frequency and severity targets in the presence of multiple types of input features. We compare in detail the performance of: a generalized linear model on binned input data, a gradient-boosted tree model, a feed-forward neural network (FFNN), and the combined actuarial neural network (CANN). Our CANNs combine a baseline prediction established with a GLM and GBM, respectively, with a neural network correction. We explain the data preprocessing steps with specific focus on the multiple types of input features typically present in tabular insurance data sets, such as postal codes, numeric and categorical covariates. Autoencoders are used to embed the categorical variables into the neural network and we explore their potential advantages in a frequency-severity setting. Finally, we construct global surrogate models for the neural nets' frequency and severity models. These surrogates enable the translation of the essential insights captured by the FFNNs or CANNs to GLMs. As such, a technical tariff table results that can easily be deployed in practice.

**Practical applications summary:** This paper explores how insights captured with deep learning models can enhance the insurance pricing practice. Hereto we discuss the required data preprocessing and calibration steps, and we present a work flow to construct GLMs for frequency and severity data by leveraging the insights obtained with a carefully designed neural network.

**JEL classification:** G22

**Key words:** property and casualty insurance, pricing, neural networks, embeddings, interpretable machine learning

## 1 Introduction

One of the central problems in actuarial science is the technical pricing of insurance contracts. Premiums are determined at the time of underwriting, while the actual cost of the contract
is only known when claims are processed. The technical premium is defined as the expected loss on a contract. In property and casualty insurance (P&C), expected losses are often estimated by independently modelling the frequency and severity of claims in function of policy and policyholder information. Hence, the modelling of historical data sets, with policyholder characteristics and the observed claim frequency and severity, is key in the design of predictive models. These historical data sets are of tabular structure, containing both numerical, categorical and spatial variables.

Industry-standard is the use of generalized linear models (GLM), introduced by Nelder and Wedderburn (1972), as a predictive modelling tool for claim frequency and severity. Haberman and Renshaw (1996), De Jong and Heller (2008), Ohlsson and Johansson (2010) and Denuit et al. (2019) apply GLMs for non-life insurance pricing. Frees and Valdez (2008) and Antonio et al. (2010) convert the numerical inputs to categorical format for use in a frequency GLM. Henckaerts et al. (2018) present a data-driven method for constructing both a frequency and severity GLM on categorized input data, by combining evolutionary trees and generalized additive models to convert the numerical inputs to categorical variables.

In recent years, machine learning techniques for actuarial purposes have been rising in popularity because of their strong predictive powers. Both Wuthrich and Buser (2021) and Denuit et al. (2020) detail the use of tree-based models in an actuarial context. Liu et al. (2014) use Adaboost for claim frequency modelling. Henckaerts et al. (2021) compare the performance of decision trees, random forests, and gradient boosted trees for modelling claim frequency and severity. Moreover, their paper studies a range of interpretational tools to look under the hood of these predictive models and compares the resulting technical tariffs with managerial tools. Instead of modelling the claim frequency and severity independently, the total loss random variable can be modelled directly via a gradient boosting model with Tweedie distributional assumption, see Yang et al. (2018) and Hainaut et al. (2022). Henckaerts and Antonio (2022) combine tabular contract and policyholder specific information with telematics data in a gradient boosting model for usage-based pricing. Henckaerts et al. (2022) construct a surrogate model on top of a gradient boosting model (GBM) to translate the insights captured by a GBM into a tariff table. A benchmark study on six data sets then examines the robustness of the proposed strategy.

Deep learning methods have been popular in the field of machine learning for many years. An early study of deep learning in an actuarial context is Dugas et al. (2003), comparing the performance of a GLM, decision tree, neural network and a support vector machine for the construction of a technical insurance tariff. Ferrario et al. (2020) use neural networks for frequency modelling and discuss various preprocessing steps. Wuthrich (2019) compares the performance of neural networks and GLMs on a frequency case study. Both Wuthrich (2019) and Schelldorfer and Wuthrich (2019) propose a combined actuarial neural network (CANN) for claim frequency modelling. The CANN starts with a GLM and builds a neural network adjustment on top of the GLM predictions, via a skip connection between input and output layer.

Categorical or factor data must be transformed into numerical representations in order to be utilized by neural networks (Guo and Berkhahn, 2016). This transformation is known in the literature as embedding, which maps categorical variables into numerical vectors. The choice of embedding technique can significantly impact the neural network's performance; see, for example, the claim severity study by Kuo and Richman (2021) where embedding layers are used in both a feed-forward neural network and a transformer network. Embedding layers allow a neural network to learn meaningful representations from the categorical inputs during the training of the neural network. Delong and Kozak (2023) suggest using autoencoders as an alternative method for categorical embedding. An autoencoder is a type of neural network that
learns to compress and to reconstruct data in an unsupervised manner. Using an autoencoder, a compact, numerical representation of the factor input data results that can then be used in both frequency as well as severity modelling. Delong and Kozak (2023) compare different setups of the autoencoder for claim frequency modelling and highlight the importance of normalization of the resulting numerical representation before using it in a feed-forward neural network. Meng et al. (2022) use the same technique in a claim frequency case study with telematic input data and extend the autoencoder with convolutional layers to process input data in image format.

Table 1 gives an overview of the discussed literature on deep learning for insurance pricing. We list the treatment techniques applied to categorical input data, the model architectures used and the extent of the case studies covered by these papers. Lastly, we summarize the interpretation tools used by the authors to extract insights from the model architectures.

Historical claim data sets are often of tabular structure, meaning they can be represented in matrix notation, with each column representing an input variable and each row representing a

\begin{table}
\begin{tabular}{c|c c c c c} \hline \hline
vector of policyholder information. Several papers recently questioned the performance of neural networks on tabular data. Borisov et al. (2022) compare 23 deep learning models on five tabular data sets and show how different tree-based ensemble methods outperform them. They highlight the predictive powers of techniques that combine gradient boosting models with neural network models, such as DeepGBM (Ke et al., 2019), which combines a GBM and a neural network for, respectively, numerical and categorical input features and TabNN (Ke et al., 2018), which bins the input features based on a GBM and uses the resulting bins in a neural network. Shwartz-Ziv and Armon (2022) analyze eight tabular data sets and compare five ensemble methods with four deep learning methods, concluding that the best performer combines gradient boosted trees and a neural network. Grinsztajn et al. (2022) compare the performance of gradient boosted trees, a random forest and different neural network structures on 45 different tabular data sets, highlighting the importance of data normalization and categorical treatment for deep learning models.

In light of these recent papers questioning the performance of deep learning architectures on tabular data, this paper aims to explore the added value of deep learning for non-life insurance pricing using tabular frequency and severity data. For this, we extend the analyses performed in Henckaerts et al. (2021) to deep learning models. Our study is an extension of the existing literature in five directions. First, we extend the CANN model architecture from Schellorfer and Wuthrich (2019) by combining a GBM baseline with neural network adjustments. Moreover, we study both trainable and non-trainable adjustments. Second, we compare a neural network, the proposed CANN structures and two benchmark models, a GLM and a GBM, by considering predictive accuracy and using interpretation tools. The GLM is constructed on categorized input data, following the approach outlined in Henckaerts et al. (2018); the GBM follows the setup from Henckaerts et al. (2021). Third, we study the autoencoder embedding technique from Delong and Kozak (2023) and highlight its importance in frequency-severity modelling. Because the autoencoder is trained in an unsupervised setting, the embedding can be learned on the frequency data and transferred to the severity setting, where we typically have fewer data points. Fourth, our case study is not limited to frequency modelling only but studies both frequency and severity modelling. We use four different insurance data sets to study the impact of sample size and the composition of the input data. Lastly, we use a set of interpretation techniques to capture insights from the constructed frequency and severity models and extend the surrogate GLM technique from Henckaerts et al. (2022) to deep learning architectures. We compare the resulting technical tariffs based on their Lorenz curves and look at the balance achieved by each model at portfolio level. This allows us to get a robust look at the possibilities of neural networks for frequency-severity pricing, from preprocessing steps to technical tariff.

## 2 Technical insurance pricing: notation and set-up

This paper assumes access to an insurance data set with tabular structure, meaning the data can be written in matrix notation, with each column representing a variable and each row representing a data point. We denote a data set as \(\mathcal{D}=\left(\mathbf{x}_{i},y_{i}\right)_{i=1}^{n}\), where each \(\mathbf{x}_{i}\) is a \(p\)-dimensional data point with response \(y_{i}\). Each data point \(\mathbf{x}_{i}\) can be written as a vector \(\left(x_{i,1},\ldots,x_{i,p}\right)\), where each entry \(x_{i,j}\) represents the value of input variable \(j\) for data point \(i\). When not referencing a specific observation \(i\), we often omit the subscript \(i\) and write \(\left(\mathbf{x},y\right)\), with \(\mathbf{x}=\left(x_{1},\ldots,x_{p}\right)\), each \(x_{j}\) representing a variable in our data set \(\mathcal{D}\).

The variables in our data sets can be either numerical or categorical. Assuming \(c\) categorical
variables, we order the variables in \(\mathcal{D}\) as follows:

\[\mathcal{D}=\big{(}\underbrace{x_{1},\ldots,x_{p-c}}_{\text{numerical variables}},\underbrace{x_{p-c+1},\ldots,x_{p}}_{\text{categorical variables}},\underbrace{y}_{\text{response variable}}\big{)}.\]

Insurance data sets can also contain spatial information. A spatial variable is either numerical, i.e., latitude and longitude coordinates, or categorical, i.e., postal code of residence. We do not denote spatial variables separately, but count them as a numerical or categorical variable. When introducing a data set, we specify how the spatial information is encoded.

For frequency-severity modelling, we work with a frequency data set \(\mathcal{D}^{\text{freq}}\) and a severity data set \(\mathcal{D}^{\text{sev}}\), where a data point \(\mathbf{x}_{i}\) represents information about policyholder \(i\). In \(\mathcal{D}^{\text{freq}}\), the response \(y_{i}\) is the number of claims reported by policyholder \(i\). The severity data set \(\mathcal{D}^{\text{sev}}\) consists of the policyholders from \(\mathcal{D}^{\text{freq}}\) who had at least one claim. In \(\mathcal{D}^{\text{sev}}\) we use the average claim size over all reported claims as the response \(y\). Because \(\mathcal{D}^{\text{sev}}\subseteq\mathcal{D}^{\text{freq}}\), but with a different response, we often omit the superscript freq and use \(\mathcal{D}\) and \(\mathcal{D}^{\text{sev}}\). We denote the number of observations as \(n_{f}\) for \(\mathcal{D}\) and \(n_{s}\) for \(\mathcal{D}^{\text{sev}}\), with \(n_{s}\leq n_{f}\). Note that for both \(\mathcal{D}\) and \(\mathcal{D}^{\text{sev}}\), we have the same variables \(x_{1},\ldots,x_{p}\), except that we add the extra variable _exposure-to-risk_\(e\) to the frequency data set. Exposure is the fraction of the year the insurance covered the policyholder. This is only relevant for frequency modelling, hence we do not add this variable to the severity data set. In \(\mathcal{D}^{\text{sev}}\) we do take into account the observed number of claims for each data point, to be used as a weight in the loss function.

For a regression model \(f(\cdot)\) with input covariates \((x_{1},\ldots,x_{p})\) and target response \(y\), we write the model prediction for data point \(i\) as \(f\left(x_{i,1},\ldots,x_{i,p}\right)=\hat{y}_{i}\). We train \(f\) on a training set, denoted as \(\mathcal{D}^{\text{train}}\subset\mathcal{D}\), by choosing the model-specific parameters that minimize a chosen loss function \(\sum_{\mathbf{x}_{i}\in\mathcal{D}^{\text{train}}}\mathscr{L}\left(\hat{y}_{i},y_ {i}\right)\). The out-of-sample performance of a trained model is calculated on the test set \(\mathcal{D}^{\text{test}}=\mathcal{D}\backslash\mathcal{D}^{\text{train}}\) as \(\sum_{\mathbf{x}_{i}\in\mathcal{D}^{\text{test}}}\mathscr{L}\left(\hat{y}_{i},y_ {i}\right)\). We follow the loss functions proposed by Wuthrich and Buser (2021) and Henckaerts et al. (2021) for modelling claim frequency and severity. For claim frequency modelling, where the claim count is typically assumed to be Poisson distributed, we use the Poisson deviance:

\[D_{\text{Poisson}}(f(\mathbf{x}),\mathbf{y})=\frac{2}{n_{f}}\sum_{i=1}^{n_{f}}\left(y_ {i}\ln\frac{y_{i}}{f(\mathbf{x}_{i})}-(y_{i}-f(\mathbf{x}_{i}))\right). \tag{1}\]

Note that when using the exposure-to-risk \(e\) in the frequency model, we replace each prediction \(f(\mathbf{x}_{i})\) with \(e_{i}\cdot f(\mathbf{x}_{i})\) in the Poisson loss function above. Claim severity data are often assumed to be long-tailed and right-skewed, so we use the gamma deviance given by

\[D_{\text{gamma}}(f(\mathbf{x}),\mathbf{y})=\frac{2}{n_{s}}\sum_{i=1}^{n_{s}}\alpha_{i }\left(\frac{y_{i}-f(\mathbf{x}_{i})}{f(\mathbf{x}_{i})}-\ln\frac{y_{i}}{f(\mathbf{x}_{i} )}\right), \tag{2}\]

where the weight \(\alpha_{i}\) is the observed number of claims for data point \(i\).

## 3 Deep learning architectures and preprocessing steps

### Neural network architectures

Feed-forward neural networkA feed-forward neural network (FFNN) is a type of machine learning model that utilizes interconnected layers, represented by \(\mathbf{z}^{(m)}\) with \(m=0,\ldots,M+1\)
The input layer, represented by \(\mathbf{z}^{(0)}\), provides the network with input data, while the output layer, represented by \(\mathbf{z}^{(M+1)}\), gives the network's prediction. Between the input and output layers, there can be one or more hidden layers, represented by \(\mathbf{z}^{(1)},\dots,\mathbf{z}^{(M)}\). When there are two or more hidden layers, we call the neural network a deep learning model. Each layer \(\mathbf{z}^{(m)}\) consists of \(q_{m}\) nodes, so it can be expressed as a vector \(\mathbf{z}^{(m)}=\left(z_{1}^{(m)},\dots,z_{q_{m}}^{(m)}\right)\).

Each node in a layer, excluding the input layer, is connected to all nodes in the previous layer through weights, represented by \(W_{m}\in\mathbb{R}^{q_{m}\times q_{m-1}}\), and a bias term, represented by \(\mathbf{b}_{m}\in\mathbb{R}^{q_{m}}\). An activation function \(\sigma^{(m)}\), \(m=1,\dots,M+1\), adds non-linearity to the network and allows it to learn complex relationships between inputs and outputs. The activation function is applied to the weighted sum of inputs to a node, along with its bias. Each layer \(\mathbf{z}^{(m)}\) can be written in function of the previous layer as follows:

\[\mathbf{z}^{(m)}=\sigma^{(m)}\left(W_{m}\cdot\mathbf{z}^{(m-1)}+\mathbf{b}_{m}\right). \tag{3}\]

Calculating the output of the FFNN in function of the input consists of performing a matrix multiplication for each layer and applying the activation functions. The value of a layer \(\mathbf{z}^{(m)}\) for input \(\mathbf{x}_{i}\) is denoted as \(\mathbf{z}^{(m)}_{i}\) and the value of a specific node \(j\) as \(z^{(m)}_{ij}\). When referencing a node without a specific input, we omit the subscript \(i\) and write \(z^{(m)}_{j}\).

The inputs of the neural network are the data points in a data set \(\mathcal{D}\), the dimension \(q_{0}\) of the input layer is equal to the number of variables \(p\) in the data set1. We write the input layer as \((x_{1},\dots,x_{p})\) to indicate that each node in the input layer represents an input variable from the data set. The target variable \(y\) in our insurance data sets is one-dimensional, so the output layer \(z^{(M+1)}\) has only one node and \(q_{M+1}=1\). We write the output node as \(\hat{y}\). Figure 1 gives a schematic overview of a feed-forward neural network.

Footnote 1: The dimension of the input layer can be larger than \(p\) when using an encoding technique, such as one-hot encoding.

Figure 1: Structure of a feed-forward neural network with \(p\)-dimensional input layer, hidden layers \(\mathbf{z}^{(1)},\dots,\mathbf{z}^{(M)}\), with \(q_{1},\dots,q_{M}\) nodes, respectively. The network has a single output node \(\hat{y}\).
When modelling claim frequency and severity data with GLMs, an actuary typically relies on a Poisson GLM with a log-link function for frequency and a gamma GLM with a log-link function for severity modelling. To mimic this log-link relationship between covariates and output in our FFNN, we use an exponential activation function for the output layer in both the frequency and the severity model. As such, we obtain strictly positive predictions for claim counts and claim amounts.

Combined actuarial neural networksWuthrich (2019) and Schelldorfer and Wuthrich (2019) propose a combination of a GLM with a FFNN, called the Combined Actuarial Neural Network (CANN). A CANN model calibrates a neural network adjustment on top of the GLM prediction. We refer to the GLM prediction as the _initial model prediction_, denoted as \(\hat{y}^{\text{IN}}\). We use \(\hat{y}^{\text{IN}}\) as an input node in a FFNN but do not connect this node to the hidden layers. Instead, \(\hat{y}^{\text{IN}}\) directly connects to the output node of the FFNN via a so-called skip connection. The adjustment made by the neural network on the initial model prediction is called the _adjustment model prediction_ and denoted as \(\hat{y}^{\text{NN}}\). The combination of the initial model prediction and the adjustment calibrated by the neural net is the resulting CANN model prediction, denoted as \(\hat{y}\). Figure 2 shows the structure of the CANN model.

The output node of the CANN model, \(\hat{y}\), is only connected to the initial model input \(\hat{y}^{\text{IN}}\) and the neural network adjustment \(\hat{y}^{\text{NN}}\). We use the exponential activation function in the output layer to ensure the log-link relationship between inputs and the predicted output. Because \(\hat{y}^{\text{IN}}\) is a prediction at the level of the response, we apply a log transform on the initial model predictions. The output of the CANN model is then calculated as:

\[\hat{y}=\exp\left(w_{\text{NN}}\cdot\hat{y}^{\text{NN}}+w_{\text{IN}}\cdot \ln\left(\hat{y}^{\text{IN}}\right)+b\right). \tag{4}\]

Figure 2: Structure of a Combined Actuarial Neural Network (CANN). The initial model prediction \(\hat{y}^{\text{IN}}\) is connected via a skip-connection to output node of the FFNN.
The case study in Schelldorfer and Wuthrich (2019) fixes the weights and bias in the output of the CANN as follows

\[w_{\text{NN}}=1,\,w_{\text{IN}}=1\,\text{and}\,b=0.\]

Following Gielis (2020), we call this the _fixed_ CANN, as the output weights are fixed and not trainable. In our case study, we also run experiments with trainable weights in the output layer and refer to this model as the _flexible_ CANN. This flexibility allows the training of the neural network to put more, or less, weight on the initial model prediction. This can potentially improve the predictive accuracy of the flexible CANN compared to the fixed CANN. Moreover, the initial model input is not restricted to GLM predictions and we will also run experiments in Section 4.4 with an input prediction established with a carefully trained GBM. According to Henckaerts et al. (2021) the GBMs are capable of achieving a higher predictive accuracy compared to a GLM. Using the GBM predictions as initial model input can therefore potentially increase the performance of the CANN model, compared to a CANN using the GLM predictions.

### Preprocessing steps

Continuous variablesWe normalize the continuous input variables to ensure that each variable in the input data has a similar scale. This is important because most neural network training algorithms use gradient-based optimization, which can be sensitive to the scale of the input data (Sola and Sevilla, 1997). For a continuous variable \(x_{j}\) in the input data \(\mathcal{D}\), we use normalization around zero as a scaling technique. Hereto, we replace each value \(x_{i,j}\) as follows:

\[x_{i,j}\mapsto\tilde{x}_{i,j}=\frac{x_{i,j}-\mu_{x_{j}}}{\sigma_{x_{j}}}, \tag{5}\]

where \(\mu_{x_{j}}\) and \(\sigma_{x_{j}}\) are the mean and standard deviation of the variable \(x_{j}\) in the data set \(\mathcal{D}\). When using a subset \(\mathcal{D}^{\text{train}}\subset\mathcal{D}\) to train the model, we calculate the \(\mu_{x_{j}}\) and \(\sigma_{x_{j}}\) only on the data set \(\mathcal{D}^{\text{train}}\) to avoid data leakage.

Categorical variablesThe FFNN and CANN models generate output by performing matrix multiplications and applying activation functions. Therefore, all inputs must be in numerical format. So-called embedding techniques convert categorical input variables to a numerical format. In this study, we utilize the _autoencoder embedding_ proposed by Delong and Kozak (2023). Autoencoders are neural networks commonly used for dimensionality reduction (Goodfellow et al., 2016). They consist of two components: an encoder and a decoder. The encoder maps a numerical input vector to a lower-dimensional representation, while the decoder reconstructs the original input from this representation. During training, the autoencoder minimizes the difference between the original and reconstructed inputs, resulting in an encoder that captures the most important characteristics of the data.

Figure 2(a) shows the general structure of such an autoencoder. It consists of an input layer of dimension \(c\), one hidden layer \(\mathbf{z}^{\text{enc}}\) of dimension \(d\) and an output layer of the same dimension as the input layer. The encoding layer is defined by the activation function \(\sigma^{\text{(enc)}}\), weights matrix \(W_{\text{enc}}\in\mathbb{R}^{d\times c}\) and bias vector \(\mathbf{b}_{\text{enc}}\in\mathbb{R}^{d}\). Similarly, the output layer is defined by activation function \(\sigma^{\text{(dec)}}\), weight matrix \(W_{\text{dec}}\in\mathbb{R}^{c\times d}\) and bias vector \(\mathbf{b}_{\text{dec}}\in\mathbb{R}^{c}\). For input \(\mathbf{x}_{i}\in\mathbb{R}^{c}\), the encoded and decoded representations are calculated as

\[\mathbf{z}_{i}^{\text{enc}} =\sigma^{\text{(enc)}}\left(W_{\text{enc}}\cdot\mathbf{x}_{i}+\mathbf{b}_{ \text{enc}}\right), \tag{6}\] \[\mathbf{x}_{i}^{\text{dec}} =\sigma^{\text{(dec)}}\left(W_{\text{dec}}\cdot\mathbf{z}_{i}^{\text {enc}}+\mathbf{b}_{\text{dec}}\right).\]
The autoencoder is trained on all data points \(\mathbf{x}_{i}\) in a data set \(\mathcal{D}\) by adjusting the weight matrices and bias vectors in order to minimize a chosen loss function \(\sum_{\mathbf{x}_{i}\in\mathcal{D}}\mathscr{L}\left(\mathbf{x}_{i}^{\text{dec}},\mathbf{x}_ {i}\right)\).

Our study employs an autoencoder to construct an embedding for multiple categorical input variables. First, we construct the one-hot encoded representation of each categorical variable (Ferrario et al., 2020). One-hot encoding maps a categorical variable \(x_{j}\) with \(L_{j}\) levels to a binary vector \(x_{j}^{\text{OH}}=\left(x_{1}^{(j)},\ldots,x_{L_{j}}^{(j)}\right)\) in the space \(\{0,1\}^{L_{j}}\). If we have \(c\) categorical variables, the dimension of all one-hot representations together equals \(\sum_{j=1}^{c}L_{j}\).

Second, we train an autoencoder using the combined one-hot representations of the categorical variables as input nodes. As such, the input layer has a dimension of \(\sum_{j=1}^{c}L_{j}\). The input layer is connected to an encoded layer of dimension \(d\), which is then connected back to the output layer of dimension \(\sum_{j=1}^{c}L_{j}\). We use the identity function as activation function for both \(\sigma^{\text{(enc)}}\) and \(\sigma^{\text{(dec)}}\) in Equation (6).

Following the construction in Delong and Kozak (2023), we apply a softmax transformation on the output layer of the autoencoder after the activation function \(\sigma^{\text{(dec)}}\). For each categorical variable \(x_{j}\), exactly one value in the input nodes \(\left(x_{1}^{(j)},\ldots,x_{L_{j}}^{(j)}\right)\) is one and the rest of the input nodes takes the value zero. Therefore, we apply the softmax activation function to the output layer of the autoencoder for each group of nodes corresponding to the one-hot encoding of a categorical variable. For each categorical variable \(x_{j}^{\text{OH}}=\left(x_{1}^{(j)},\ldots,x_{L_{j}}^{(j)}\right)\) and for each

Figure 3: Our proposed network structure combines the autoencoder embedding technique from Delong and Kozak (2023) and the CANN structure from Schelldorfer and Wüthrich (2019).
\(h\in\{1,\ldots,L_{j}\}\), the softmax transformation of the output node \(x_{h}^{(j,\text{dec})}\) is defined as:

\[x_{h}^{(j,\text{dec})}\mapsto\tilde{x}_{h}^{(j,\text{dec})}=\frac{\exp\left(x_{h }^{(j,\text{dec})}\right)}{\exp\left(x_{1}^{(j,\text{dec})}+\ldots+x_{L_{j}}^{ (j,\text{dec})}\right)},\qquad h=1,\ldots,L_{j}. \tag{7}\]

The use of the softmax activation function ensures that the values of the decoded vectors \(\left(x_{1}^{(j,\text{dec})},\ldots,x_{L_{j}}^{(j,\text{dec})}\right)\) sum up to one for each variable \(x_{j}\).

To train the autoencoder, we use the cross-entropy loss function, which is suitable because of the 0/1 values in the input data. With \(\mathbf{x}_{i}^{\text{OH}}\) the one-hot encoding of all categorical variables for policyholder \(i\) and \(\tilde{\mathbf{x}}_{i}^{\text{dec}}\) the values of the autoencoder's output layer for policyholder \(i\), the cross-entropy loss function is defined as:

\[\mathscr{L}^{\text{CE}}\left(\tilde{\mathbf{x}}_{i}^{\text{dec}},\mathbf{x}_{i}^{\text {OH}}\right)=-\sum_{j=1}^{c}\sum_{h=1}^{L_{j}}x_{ih}^{(j)}\cdot\log\left(\tilde {x}_{ih}^{(j,\text{dec})}\right). \tag{8}\]

After training the autoencoder and applying the trained autoencoder on each policyholder \(\mathbf{x}_{i}\in\mathcal{D}\), the vector of categorical inputs \((x_{i,p-c+1},\ldots,x_{i,p})\) is accurate, compact and numerically represented in the vector \((z_{i1}^{\text{enc}},\ldots,z_{id}^{\text{enc}})\) as calculated by Equation (6). We call the vector \(\mathbf{z}_{i}^{\text{enc}}\) the embedding of the categorical inputs of \(\mathbf{x}_{i}\). To use the embedding together with the numerical features of \(\mathbf{x}_{i}\), we normalize the values in the nodes \(z_{1}^{\text{enc}},\ldots,z_{d}^{\text{enc}}\) by scaling the weight matrix \(W_{\text{enc}}\) and bias vector \(\mathbf{b}_{\text{enc}}\) of the trained encoder. With \(\mu_{1},\ldots,\mu_{d}\) the means, and \(\sigma_{1},\ldots,\sigma_{d}\) the standard deviations, of the values \(z_{i1}^{\text{enc}},\ldots,z_{id}^{\text{enc}}\) for all \(\mathbf{x}_{i}\in\mathcal{D}^{\text{train}}\), we scale the weight matrix \(W_{\text{enc}}\) and bias vector \(\mathbf{b}_{\text{enc}}\) of the pre-trained encoder as follows:

\[W_{\text{enc}}\mapsto\tilde{W}_{\text{enc}}=\left(\begin{array}{cccc}\frac{ w_{11}}{\sigma_{1}}&\frac{w_{12}}{\sigma_{1}}&\ldots&\frac{w_{1c}}{\sigma_{1}} \\ \frac{w_{21}}{\sigma_{2}}&\frac{w_{12}}{\sigma_{2}}&\ldots&\frac{w_{1c}}{\sigma_ {2}}\\ \vdots&\vdots&\ddots&\vdots\\ \frac{w_{d-1,1}}{\sigma_{d-1}}&\frac{w_{d-1,2}}{\sigma_{d-1}}&\ldots&\frac{w_ {d-1,c}}{\sigma_{d-1}}\\ \frac{w_{d1}}{\sigma_{d}}&\frac{w_{d2}}{\sigma_{d}}&\ldots&\frac{w_{d}}{\sigma _{d}}\\ \end{array}\right),\mathbf{b}_{\text{enc}}\mapsto\tilde{\mathbf{b}}_{\text{enc}}= \left(\begin{array}{c}\frac{b_{1}-\mu_{1}}{\sigma_{1}}\\ \frac{b_{2}-\mu_{2}}{\sigma_{2}}\\ \vdots\\ \frac{b_{d-1}-\mu_{d-1}}{\sigma_{d-1}}\\ \frac{b_{d-1}-\mu_{d-1}}{\sigma_{d-1}}\\ \frac{b_{d-1}-\mu_{d}}{\sigma_{d}}\\ \end{array}\right). \tag{9}\]

Having access to the trained and scaled autoencoder, we now add the encoder part to the FFNN and the CANN structures by replacing the input nodes of the categorical variables in Figure 1 and 2 with the encoding part of the trained autoencoder, as shown in Figure 2(b) for the CANN. For clarity, we omit the one-hot encoding notation of each variable in Figure 3. We say the autoencoder is _pre-trained_ because we perform a first training and scaling of the autoencoder before training the neural network architectures with the added encoder. Adding the encoder to the network allows the network to finetune the weights and biases of the pre-trained encoder with respect to the considered regression task and its applicable loss function as in Equation(1) or Equation (2).

Autoencoders used to embed categorical variables provide several advantages over one-hot encoding (Delong and Kozak, 2023). Firstly, they allow for a significantly smaller dimension of the encoding compared to the dimension resulting from one-hot encoding. Secondly, autoencoders enable the encoding of all categorical variables together, capturing interactions between variables more effectively than variable specific encoding does. Lastly, autoencoders prove advantageous in multi-task scenarios such as frequency-severity modeling. Learning to encode
categorical variables solely on the severity dataset can be problematic due to its smaller size. Since autoencoders are unsupervised learning methods, we can train the autoencoder using all data available, and add the resulting pre-trained encoder to both frequency and severity models.

### Training and tuning neural networks

We train the FFNN and CANN models using the Adam optimization algorithm. Adam, introduced by Kingma and Ba (2014), is a stochastic gradient descent algorithm with an adaptive learning rate. Iteratively, the Adam algorithm changes the weights and biases in the network to minimize the loss between predictions \(\hat{y}\) and the observed responses \(y\). We use batches of training data for each training iteration to speed up optimization; see Keskar et al. (2016). The size of the batches is a parameter that needs to be tuned. The network size is also tuned; the number of hidden layers \(M\), and the number of nodes in each layer \(q_{1},\ldots,q_{M}\) are tuning parameters. We use a drop-out rate (Srivastava et al., 2014) to avoid overfitting, and consider this rate to be a tuning parameter as well. The drop-out rate is the percentage of nodes in each layer that are disconnected from the next and previous layer during each iteration of the Adam algorithm. The last tuning parameter is the choice of activation functions \(\sigma^{(1)},\ldots,\sigma^{(M)}\). To simplify the tuning process, we use layers of equal sizes, \(q_{1}=\ldots=q_{M}=q\), and apply the same activation function for all hidden layers, \(\sigma^{(1)}=\ldots=\sigma^{(M)}=\sigma\). Hence, only the value for \(q\) and the activation function \(\sigma\) are tuned and applied to each hidden layer.

We deploy a random grid search, introduced by Bergstra and Bengio (2012), to determine the optimal value for each tuning parameter. For each tuning parameter \(t_{k}\), with \(k=1,\ldots,K\), we define a range of possible values \([t_{k,\min},t_{k,\max}]\). The search space \(\mathcal{S}\) is the space consisting of all possible values for all tuning parameters:

\[\mathcal{S}=[t_{1,\min},t_{1,\max}]\times\ldots\times[t_{K,\min},t_{K,\max}]\,.\]

The _random grid_\(\mathcal{R}\subset\mathcal{S}\) consists of randomly drawn points in the search space \(\mathcal{S}\). Each point \(s\in\mathcal{R}\) represents a set of candidate tuning parameter values. Out of the random grid \(\mathcal{R}\), we select the optimal point \(s^{*}\) with a cross-validation scheme. In Figure 4, we give an example of a search space defined by two tuning parameters and a random grid of size nine sampled in the search space.

We use the extensive cross-validation scheme proposed by Henckaerts et al. (2021), as sketched in Figure 5. We divide the data set \(\mathcal{D}\) in six disjoint and stratified subsets \(\mathcal{D}_{1},\ldots,\mathcal{D}_{6}\). We define six data folds; in data fold \(\ell\), for \(\ell=1,\ldots,6\), we select a hold-out test set \(\mathcal{D}_{\ell}\) and use five-fold cross-validation (Hastie et al., 2009) on the data set \(\mathcal{D}\backslash\mathcal{D}_{\ell}\). Each cross-validation loop uses four out of the five data subsets in \(\mathcal{D}\backslash\mathcal{D}_{\ell}\) to train the neural network. The fifth subset is used both for early stopping and to calculate the validation error. The cross-validation error is the average validation error over the five validation sets. We then determine the optimal point \(s^{*}_{\ell}\in\mathcal{R}\) which minimizes the cross-validation error for data fold \(\ell\). We use the six optimal tuning parameter sets \(s^{*}_{1},\ldots,s^{*}_{6}\) to determine out-of-sample performance on the test set \(\mathcal{D}_{1},\ldots,\mathcal{D}_{6}\) of each data fold \(\ell=1,\ldots,6\). As such, we obtain an out-of-sample prediction for every point in the data set.

## 4 Performance comparison between benchmark models and deep learning architectures

Section 4.1 introduces four data sets that are used in our benchmark study. In this study we compare the performance of the deep learning architectures against two benchmark mod
els introduced in Section 4.2. Section 4.3 covers the used tuning parameter grid for both the autoencoder and the deep learning architectures. We compare the statistical out-of-sample performance of the models under study in Section 4.4. Lastly, Section 4.5 compares the autoencoder embedding against using the one-hot encoding when used in the deep learning models under consideration.

Figure 4: Example of random grid search with two tuning parameters \(t_{1}\) and \(t_{2}\). The search space \(\mathcal{S}=[t_{1,\min},t_{1,\max}]\times[t_{2,\min},t_{2,\max}]\) is shown in the figure by the dotted square. The random grid \(\mathcal{R}\) consists of nine randomly drawn points \(s_{1},\ldots,s_{9}\) from \(\mathcal{S}\). The optimal point \(s^{*}\in\mathcal{R}\) is then selected via a cross-validation scheme.

Figure 5: Representation of the 6 times 5-fold cross-validation scheme, figure from Henckaerts et al. (2021).
### Data sets

The used data sets are an Australian, Belgian, French2 and Norwegian MTPL data set, available through the R packages CASdatasets(Dutang and Charpentier, 2019) and maidrr(Henckaerts and Antonio, 2022; Henckaerts, 2021). Table 2 gives an overview of the number of records for each data set and the number of continuous, categorical and spatial variables.

Footnote 2: The French data set in the CASdatasets package contains 35 560 claims, but only 24 000 claims have a claim amount. We exclude the policies with claims but without claim amount from our study.

The spatial variables are listed separately in Table 2. The Belgian spatial variable is the postal code, which is converted to two continuous variables, the latitude and longitude coordinates of the center of that postal code. The French data includes two spatial variables: the French district, which is categorical, and the logarithm of the population density of the place of residence, which is a continuous variable. The Norwegian data has one spatial variable denoting the population density of the region of residence as a categorical variable.

### Benchmark models

To enable an assessment of the predictive performance of the neural network and CANN structures we construct two benchmark models; a generalized linear model (GLM) and a gradient boosting model (GBM), for both frequency as well as severity. Predictions from these benchmark models are then also used as the initial model inputs in the CANN models. For the Belgian data set we use the GLM constructed in Henckaerts et al. (2018) and the GBM from Henckaerts et al. (2021). For the other data sets, we follow the construction methods outlined in the mentioned papers.

For the construction of the GLM, we follow the strategy proposed in Henckaerts et al. (2018) and start from a generalized additive model (GAM), including interaction effects between continuous variables. Based on the insights from the GAM, we bin the continuous variables using a regression tree. On the binned input data, we construct a GLM. We repeat the construction of the GLM six times, each time withholding a subset \(\mathcal{D}_{\ell}\), \(\ell=1,\ldots,6\). This way we obtain GLM based out-of-sample predictions for all observations in the data set \(\mathcal{D}\).

GBM is an ensemble method combining multiple decision trees (Friedman, 2001). A GBM has two tuning parameters; the number of trees and the depth of each tree. We use a tuning grid

\begin{table}
\begin{tabular}{c c c c} \hline \hline Australian MTPL & Belgian MTPL & French MTPL & Norwegian MTPL \\ \hline \hline \multicolumn{3}{c}{**Number of observations**} \\ \hline Frequency & 67 856 & 163 212 & 668 897 & 183 999 \\ Severity & 4 624 & 18 276 & 24 944 & 8 444 \\ \hline \multicolumn{3}{c}{**Covariates: number and type**} \\ \hline Continuous & 1 & 4 & 2 & 0 \\ Categorical & 4 & 5 & 5 & 3 \\ Spatial & 0 & 1 & 2 & 1 \\ \hline \hline \end{tabular}
\end{table}
Table 2: Overview of the structure of the data sets used in the benchmark study. The number of records for both frequency and severity modelling is given, as well as the number of different input variables per type.
with the following values:

\[\text{Number of trees:}\,\{100,300,500,\ldots,5\,000\},\] \[\text{Depth of each tree:}\,\{1,2,3,4,5,6,7,8,9,10\}.\]

We use three hyper parameters of which the values are not tuned: shrinkage = 0.01, bagging fraction = 0.75 and minimum observations per node is 0.75% of the number of records in the training data. The loss functions in Equation (1) and Equation (2) are used for, respectively, frequency and severity modelling. We follow the repeated 5-fold cross-validation scheme as described in Section 3.3. With optimal tuning parameters, we fit a GBM for each data fold and look at the prediction for the observations in the corresponding test set. As such, we obtain an out-of-sample prediction for every data point in the portfolio.

### Neural network models

The pre-training of the autoencoder uses the Nadam optimizer algorithm, a batch size of 1 000 and a randomly selected validation set of 20% of the frequency data set \(\mathcal{D}\) for early stopping. The number of nodes \(d\) in the encoding layer is tuned, testing across the values \(\{5,10,15\}\), selecting the lowest value of \(d\) while the loss \(\mathscr{L}^{\text{CE}}(\cdot,\cdot)<0.001\), as calculated with Equation (8). After the autoencoder is trained and scaled, the encoder is used in each FFNN and CANN structure, for both frequency and severity modelling.

For both the FFNN and the CANN models, a random grid \(\mathcal{R}\) of size 40 is sampled from the search space \(\mathcal{S}\) defined by the tuning parameters and their respective ranges as shown in Table 3.

The cross-validation scheme is shown in Algorithm 1, starting with the pre-processing steps, and resulting in out-of-sample performances for each holdout test set \(\mathcal{D}_{1},\ldots,\mathcal{D}_{6}\). For \(\ell=1,\ldots,6\), we train a network on the data \(\mathcal{D}\backslash\mathcal{D}_{\ell}\), choosing a random validation set consisting of 20% of the training data for early stopping. With this model, we construct out-of-sample predictions on the test set \(\mathcal{D}_{\ell}\) and calculate the out-of-sample loss using the loss functions in Equation (1) and (2). Because optimization in a neural network is dependent on the random initialization of the weights, we train the model three times, and use the average out-of-sample loss over the three trainings. This ensures an objective out-of-sample loss evaluation, without the risk of accidentally getting well or badly initialised weights.

\begin{table}
\begin{tabular}{l l c} \hline \hline
**Tuning parameter** & **Range** & \\ \hline Activation function for hidden layers & ReLU, sigmoid, softmax\({}^{3}\) & \\ Batch size & \([10\,000,50\,000]\) & Frequency \\  & \([200,10\,000]\) & Severity \\ Number of hidden layers & \([1,4]\) & \\ Nodes per hidden layer & \([10,50]\) & \\ Dropout rate & \([0,0.1]\) & \\ \hline \hline \end{tabular}
\end{table}
Table 3: Collection of tuning parameters and their respective ranges for the random grid search tuning strategy. This range is used for both the FFNN and CANN structures.
```
Input: model class (mclass) and corresponding tuning grid \(\mathcal{R}\) data \(\mathcal{D}\) with 6 disjoint stratified subsets \(\mathcal{D}_{1},\ldots,\mathcal{D}_{6}\); for\(\ell=1,\ldots,6\)do  leave out \(\mathcal{D}_{\ell}\) as test set; foreach continuous variable \(x_{j}\in\mathcal{D}\)do  calculate mean \(\mu_{x_{j}}\) and standard deviation \(\sigma_{x_{j}}\) on the data \(\mathcal{D}\setminus\mathcal{D}_{\ell}\);  normalize the variable \(x_{j}\) in the data set \(\mathcal{D}\) with Equation (5); foreach categorical variable \(x_{j}\in\mathcal{D}\)do  construct one-hot encoding \(\left(x_{1}^{(j)},\ldots,x_{L_{j}}^{(j)}\right)\) of variable \(x_{j}\); for\(d\in\{5,10,15\}\)do  train an autoencoder \(f_{d}^{\mathrm{AE}}\) using the one-hot encoding of all categorical variables in \(\mathcal{D}\setminus\mathcal{D}_{\ell}\) as input and \(d\) nodes in the encoding layer;  evaluate the model performance on \(\mathcal{D}\setminus\mathcal{D}_{\ell}\) using loss function \(\mathscr{L}^{\mathrm{BCE}}(\cdot,\cdot)\);  select \(f_{d}^{\mathrm{AE}}\) with the lowest \(d\) while \(\mathscr{L}^{\mathrm{BCE}}(\cdot,\cdot)<0.001\);  calculate the scaled \(\tilde{W}_{\mathrm{enc}}\) and bias vector \(\tilde{\mathbf{b}}_{\mathrm{enc}}\) for the encoder part of \(f_{d}^{\mathrm{AE}}\);  add the encoder part of \(f_{d}^{\mathrm{AE}}\), with \(\tilde{W}_{\mathrm{enc}}\) and \(\tilde{\mathbf{b}}_{\mathrm{enc}}\), to the model class mclass; foreach tuning parameter point \(s\in\mathcal{R}\)do for\(k\in\{1,\ldots,6\}\setminus\ell\)do  train a model \(f_{\ell k}\) of mclass on \(\mathcal{D}\setminus\{\mathcal{D}_{\ell},\mathcal{D}_{k}\}\);  evaluate the model performance on \(\mathcal{D}_{k}\) using loss function \(\mathscr{L}(\cdot,\cdot)\);  valid_error\({}_{\ell k}\leftarrow\frac{1}{|\mathcal{D}_{k}|}\underset{\mathbf{x}_{i}\in\mathcal{D}_{k}} {\sum}\mathscr{L}\{y_{i},f_{\ell k}(\mathbf{x}_{i})\}\);  valid_error\({}_{\ell}\leftarrow\frac{1}{5}\sum_{k\in\{1,\ldots,6\}\setminus\ell}\) valid_error\({}_{\ell k}\);  optimal parameter point \(s_{\ell}^{*}\in\mathcal{R}\) minimizes valid_error\({}_{\ell}\); for\(\mathrm{rep}\in\{1,2,3\}\)do  select a validation set \(\mathcal{D}_{\mathrm{val}}\) containing 20% of the records in \(\mathcal{D}\setminus\mathcal{D}_{\ell}\);  train a model \(f_{\ell,\mathrm{rep}}\) of mclass on \(\mathcal{D}\setminus\{\mathcal{D}_{\ell},\mathcal{D}_{\mathrm{val}}\}\) using the optimal parameter point \(s_{\ell}^{*}\) and using \(\mathcal{D}_{\mathrm{val}}\) for early stopping;  evaluate the model performance on \(\mathcal{D}_{\ell}\) using loss function \(\mathscr{L}(\cdot,\cdot)\);  test_error\({}_{\ell,\mathrm{rep}}\leftarrow\frac{1}{|\mathcal{D}_{\ell}|}\underset{\mathbf{x}_{i} \in\mathcal{D}_{\ell}}{\sum}\mathscr{L}\{y_{i},f_{\ell,\mathrm{rep}}(\mathbf{x}_{i })\}\);  test_error\({}_{\ell}\leftarrow\frac{1}{3}\sum_{\mathrm{rep}}\) test_error\({}_{\ell,\mathrm{rep}}\); Output: optimal tuning parameters + performance measure for each of the six folds.
```

**Algorithm 1**Pseudocode to sketch the pipeline for calculating out-of-sample performances with the neural network structures from Section 3.1, including the data pre-processing steps, the cross-validation scheme as outlined in Henckaerts et al. (2021) with the random grid search methodology and the repeated out-of-sample loss calculation to avoid local minima solutions.
### Out-of-sample performances

The two benchmark models enable a comparison between their performance and those obtained with the proposed neural network architectures. Moreover, they serve as initial model input for the CANN models. We investigate the predictive performance of seven models for each data set: GLM, GBM, neural network, the CANN with GLM input, with fixed and flexible output layer, and the CANN with GBM input, with both fixed and flexible output layer. We compare the out-of-sample performances of the benchmark models and the neural network structures in Figure 6. We have out-of-sample deviances for each withheld test set, measured in Poisson deviance (1) or gamma deviance (2). We show frequency on the left-hand side and severity on the right-hand side.

Among the four data sets analyzed, the combination of a neural network and a gradient boosting model (CANN GBM flex) consistently yields the lowest deviance when modelling claim frequency. This aligns with recent research highlighting the predictive performance of combining a gradient boosting model with a neural network (Borisov et al., 2022; Ke et al., 2019; Shwartz-Ziv and Armon, 2022). However, for the Norwegian data set, which has few input variables, the impact was less pronounced, with similar performance observed across all models except for the feed-forward neural network, which exhibits slightly higher deviance. Regarding claim severity modelling, no single model consistently achieves the lowest deviance across all data sets and test sets. For the Australian and Norwegian data sets, all models perform comparably in terms of deviance. The CANN models with GBM input demonstrate the lowest deviance for the Belgian data set, while for the French data set, the CANN model with GLM input achieves the best results. Notably, the CANN models with a flexible output layer structure outperform those with a fixed output layer in most cases, for both frequency and severity modelling. This suggests that the more adaptable combination of the initial model input and the neural network adjustment leads to reduced deviance.

### Comparison of categorical embedding methods

We investigate the impact of the autoencoder embedding compared to directly utilizing one-hot encoded categorical variables. For each data set, we train a FFNN and a CANN model using the one-hot encoding of each categorical variable and a FFNN and a CANN model using the autoencoder embedding. We do not tune the models but choose a set of tuning parameters that we apply to all shown models. This means the deviance of each model is not relevant here, only the difference in deviance between the model with one-hot encoding and the model with autoencoder embedding. This approach allows us to isolate the embedding technique's effect on a model's predictive performance. Each model is trained on \(\mathcal{D}\setminus\mathcal{D}_{1}\), and the out-of-sample performance is calculated on the out-of-sample test set \(\mathcal{D}_{1}\).

Figure 7 displays the predictive accuracy of each model under consideration, with the frequency models in the top and the severity models in the bottom row. For frequency modeling, the autoencoder embedding has the most pronounced effect on the performance of the FFNNs, leading to a lower deviance compared to models utilizing one-hot encoding. However, the impact on CANN models appears to be negligible. In the case of severity modeling, both FFNNs and CANNs demonstrate an improved predictive performance when using the autoencoder embedding. Only the FFNN on Australian data set and the CANN model on Belgian data set perform similarly when comparing the one-hot encoding and the autoencoder embedding for claim severity modelling. The reduced deviance in most severity models highlights the benefits of unsupervised learning through the autoencoder approach.
Figure 6: Out-of-sample performance comparison between the different models for each data set. The left-hand side shows the performance of the frequency models and the right-hand side for the severity models. From top to bottom, we show the results on the Australian, Belgian, French and Norwegian data sets. The deviances for the GLM and GBM on the Belgian data correspond to the results reported in, respectively, Henckaerts et al. (2018) and Henckaerts et al. (2021).
## 5 Looking under the hood: interpretation tools and surrogate models

First, we use two model interpretation tools to look under the hood of the constructed models. Second, we translate the model insights into a tariff structure by constructing GLM surrogates along the work flow presented in Henckaerts et al. (2022). All results shown in this section are calculated using data fold one, meaning the models are trained using data subsets \(\mathcal{D}_{2},\ldots,\mathcal{D}_{6}\) and the shown results are calculated on the test set \(\mathcal{D}_{1}\).

### Variable importance

We measure variable importance using the permutation method from Olden et al. (2004). Hereby, we consider the average change in predictions when a variable is randomly permuted. For a trained model \(f\), we measure the importance of a variable \(x_{j}\) by calculating

\[\mathrm{VIP}_{x_{j}}=\sum_{\mathbf{x}_{i}\in\mathcal{D}}\mathrm{abs}\big{(}f\left( x_{i,1},\ldots,x_{i,j},\ldots,x_{i,p}\right)-f\left(x_{i,1},\ldots,\tilde{x}_{i,j}, \ldots,x_{i,p}\right)\big{)}, \tag{10}\]

where \(\tilde{x}_{i,j}\) is a random permutation of the values observed for \(x_{j}\) in the data set \(\mathcal{D}\). A large value for \(\mathrm{VIP}_{x_{j}}\) indicates that the variable significantly influences the model output and is therefore considered important. Figure 8 shows the variable importance of each variable in the four data sets for both frequency and severity modelling. For clarity, we show the relative VIP of each variable, calculated as

\[\overline{\mathrm{VIP}}_{x_{j}}=\frac{\mathrm{VIP}_{x_{j}}}{\sum_{x_{j}\in \mathcal{D}}\mathrm{VIP}_{x_{j}}},\qquad\text{where the sum runs over all over all variables $x_{j}$.} \tag{11}\]

Figure 7: Comparison of one-hot encoding and autoencoder embedding on the out-of-sample performance of both the FFNN and the CANN model. Top row shows the effect on frequency modelling and bottom row on severity modelling.
By comparing the variable importance of the GBM with the CANN model, we can evaluate the impact of the neural network adjustment component within the CANN. In general, most variables show similar importance in both the GBM and the CANN GBM flexible models, indicating that the adjustment calibrated by the neural network does not substantially alter the importance of the relationships between input variables and the response variable. However, notable changes are observed for certain variables, such as the postal code in the frequency model for the Belgian data set and the vehicle age and brand in the frequency model for the French data set. When we compare the variable importance of the GBM and CANN GBM with the FFNN, we observe more substantial changes, particularly in claim severity modelling. This shows that the FFNN models a significantly different relationship between the input variables and the output variable, compared to the GBM and CANN GBM flexible model.

### Partial dependence effects

We consider partial dependence effects (Hastie et al., 2009; Henckaerts et al., 2021) to explore the relationship between an input variable and the model output. Let the variable space \(X_{j}\) be the vector containing all possible values for variable \(x_{j}\). For a trained model \(f\), the partial dependency effect of a variable \(x_{j}\) is the vector calculated as

\[\mathrm{PD}_{x_{j}}=\left\{\frac{1}{|\mathcal{D}|}\sum_{\mathbf{x}_{i}\in\mathcal{ D}}f\left(x_{i,1},\ldots,X_{o,j},\ldots,x_{i,p}\right)\text{ ; }\forall X_{o,j}\in X_{j}\right\}, \tag{12}\]

where \((x_{i,1},\ldots,X_{o,j},\ldots,x_{i,p})\) is the data point \(\mathbf{x}_{i}\in\mathcal{D}\) with element \(x_{i,j}\) replaced by the value \(X_{o,j}\in X_{j}\). The vector \(\mathrm{PD}_{x_{j}}\) can be seen as the average prediction on \(\mathcal{D}\), while letting the

Figure 8: Relative variable importance in the GBM, the FFNN and the CANN GBM flexible. Top row shows the effects for the frequency models, the bottom row for the severity models.
variable \(x_{j}\) range over all possible values in the variable space \(X_{j}\). A partial dependence plot is the plotted effect between \(X_{j}\) and \(\text{PD}_{x_{j}}\). Equation (12) can be extended to a two-way interaction partial dependence effect by letting two variables range over their respective variable spaces.

Figure 9 shows the partial dependence effect between the policyholder's age and the predicted claim frequency across the four data sets in the benchmark study. We compare the effects of the benchmark GBM, the FFNN and the CANN GBM flexible. The effect in all three models is similar for the Australian, French and Norwegian data. However, for the Belgian data set, the GBM and CANN GBM flexible show a similar partial dependence effect, while the FFNN shows a very different pattern. The partial dependence effect of this FFNN shows a less complex, less nuanced relationship between age of the policyholder and claim frequency. Across the four data sets, the average predicted claim frequency decreases with age, which is an expected relationship between age and claim frequency. For the Belgian and French data sets, we observe an increasing effect for the older ages.

Figure 10 displays the partial dependence effect of the policyholder's age when calibrated on the claim severity data. Similar to the effects portrayed in Figure 9, the three models applied to the Australian and Norwegian data sets exhibit a comparable effect. For the Belgian and French data sets, the FFNN showcases a notably distinct partial dependence effect. Specifically

Figure 10: Partial dependence effect of the policyholder’s age across the four data sets, claim severity models. We compare the benchmark GBM, the FFNN and the CANN GBM flexible.

Figure 9: Partial dependence effect of the policyholder’s age across the four data sets, claim frequency models. We compare the benchmark GBM, the FFNN and the CANN GBM flexible.
for the French data, the FFNN model reveals an almost flat effect across all age groups.

Figure 11 shows the partial dependence effect of the bonus-malus score for the Belgian and French frequency data sets. For both data sets, the three models show an increasing relation between the level occupied in the bonus-malus scale and the expected claim frequency. According to the FFNN, the partial dependence is a distinctly smoother effect compared to the effect calibrated by the GBM and CANN GBM flexible, showing again the less complex, less nuanced relationships captured by the FFNN.

We consider the partial dependence effect of the postal code in the Belgian frequency data set in Figure 12. We compare the partial dependence effect with the empirical claim frequency in the Belgian data, calculated as the number of claims per postal code divided by the sum of the exposure for that postal code. The effect in the GBM and CANN GBM flexible is very similar, with a higher expected claim frequency around the capital of Belgium. The effect in the FFNN also shows a higher expected number of claims in the capital but the calibrated spatial effect is much smoother. This aligns with the smoother partial dependence effects for the policyholder age and bonus-malus in the Belgian frequency FFNN model. Empirically, we see a higher concentration of claims per unit of exposure in and around the capital and for some postal codes in the west and east of Belgium. This effect is visible for the GBM and CANN model but not for the FFNN.

### Surrogate models for practical applications

Surrogate model constructionHenckaerts et al. (2022) present a workflow for constructing a surrogate GLM by leveraging insights obtained with a black box model. In our study, we apply this technique to the CANN model with GBM input, as discussed in Section 3 and calibrated in Section 4. To create the surrogates, we first calculate the partial dependence effect for each individual variable and for interactions between any two variables, as discussed in Section 5.2. Next, we use the dynamic programming algorithm introduced by Wang and Song (2011) to segment the input data into homogeneous groups based on these partial dependence effects. On the resulting binned data set, we fit a generalized linear model. Constructing the surrogate

Figure 11: Partial dependence effect of the bonus-malus score for the Belgian data set, claim frequency model (left) and claim severity model (right). We compare the benchmark GBM, the FFNN and the CANN GBM flexible.
GLM on the segmented frequency and severity data leads to a tabular premium structure incorporating the insights captured by the CANN architectures.

Figure 13 shows the partial dependence effects of the CANN GBM flexible for the bonus

Figure 12: Partial dependence relationship of the spatial variable and the expected number of claims in the GBM, FFNN and CANN GBM flexible for the Belgian data. We compare the modelled effects with the empirical claim frequency in the Belgian data set.

Figure 13: Partial dependence plots for three variables in the French data set; left to right: bonus-malus scale, policyholder age and region. In color, we show the binning of the input data with the frequency surrogate GLM. Each color represents one bin of the input variable.
malus score, policyholder age and the region variable from the French data set, with respect to frequency modelling. In color, we show the obtained data segmentation. The so-called surrogate GLM is fitted on the segmented input data. The benchmark GLM, constructed via the approach in Henckaerts et al. (2018) (hereafter referred to as binned GLM), is also fitted on binned data, therefore, it is insightful to compare both the predictive accuracy and the selected variables as obtained with both techniques. To avoid data leakage in the comparison between two models, we compare the predictive accuracy on a withheld test set.

Table 4 shows the variables included in the binned GLMs and the surrogate GLMs for the Australian, Belgian and French data set. The out-of-sample performance of these models are evaluated on the withheld data set \(\mathcal{D}_{1}\). We excluded the Norwegian data set from the surrogate fitting, as this data set only consists of categorical variables. The surrogate technique selects more variables and performs better on the out-of-sample test set than the binned GLM. This finding is consistent across all three data sets. Hence, the surrogate GLM benefits from the insights learned from the neural network adjustments in the CANN compared to the direct construction of the binned GLM.

Identification of risk profilesWe estimate the number of claims and the claim severities using the surrogate models constructed for the frequency and severity CANN GBM flexible, respectively.

We construct a low, medium, and high-risk profile based on the frequency surrogate GLM for the French data set. Table 5 compares these profiles via their expected claim frequency according to the surrogate GLM and the CANN GBM flexible model. We compare the influence of each variable on the assessed risk using two local interpretation tools in Figure 14. For the GLM, we

\begin{table}

\end{table}
Table 4: Comparison between the benchmark GLM and the surrogate GLM for frequency modelling on the Australian, Belgian and French data sets. The surrogate GLM is constructed from the CANN with GBM input and flexible output layer. The last row shows the Poisson deviance of both GLMs on the out-of-sample data set \(\mathcal{D}_{1}\).
show the fitted coefficients on the response scale. A value lower (higher) than one means the feature's value leads to a lower (higher) prediction than the baseline prediction obtained with the intercept of the GLM. The uncertainty of each contribution is shown with the 95% confidence interval. Shapley values (Shapley et al., 1953) are used to compare the feature contributions of the GLM to the influences in the CANN model. A positive (negative) Shapley value indicates that this feature's value leads to a higher (lower) than average prediction. The effect in the GLM and CANN model mostly align. We see a strong impact of the variables region, driver age and bonus-malus score on the predicted number of claims. The variable area was not selected in the surrogate GLM construction, and its Shapely value is negligible in all tree risk profiles.

## 6 Managerial insights: a comparison between technical tariff stuctures

We now combine the predictions for claim frequency and severity to a technical tariff. For each data set and each fold, we make predictions for all observations in the test set \(\mathcal{D}_{\ell}\), using the model trained on the data subsets \(\mathcal{D}\setminus\mathcal{D}_{\ell}\). As such, we obtain out-of-sample predictions for both the expected number of claims and the expected claim severities for each policyholder in the data set. The predicted loss, or technical tariff, for each policyholder is the expected number of claims times the expected claim severity.

Table 6 shows the total predicted loss next to the total loss observed in each data set. We compare the results from the benchmark GLM and GBM, the CANN GBM flexible and the surrogate GLM. We also show the ratio of predicted losses over the observed losses. A ratio of one means the model has perfect balance at portfolio level. For the Norwegian data set, the predicted losses are very close to the observed losses for all models. For the Australian and Belgian data, both GLM models are close to balance, meaning the predicted losses are close to the observed losses. Although a canonical link GLM satisfies the balance property (Nelder and Wedderburn, 1972), our severity models use a gamma distribution with non-canonical log-link,

\begin{table}
\begin{tabular}{l c c c} \hline \hline
**Variables** & **Low risk** & **Medium risk** & **High risk** \\ \hline Vehicle power & 4 & 6 & 9 \\ Vehicle age & 3 & 2 & 1 \\ Policyholder age & \([21,26[\) & \(4\,[30,40[\) & \(\geq 70\) \\ Bonus-malus scale & 50 & 70 & 190 \\ Vehicle brand & B12 & B5 & B11 \\ Fuel type & Regular & Regular & Diesel \\ Population density of area & 2.71 & 665.14 & \(22\,026.47\) \\ District of residence & Midi-Pyrenees & Basse-Normandie & Corse \\ \hline
**Predicted number of claims** & & & \\ \hline Surrogate GLM & 0.020 & 0.106 & 0.361 \\ CANN GBM flexible & 0.021 & 0.101 & 0.519 \\ \hline \hline \end{tabular}
\end{table}
Table 5: Example of a low, medium and high risk profile for the French data set, using the surrogate based on the CANN model with GBM input and flexible output layer, withholding test set one. We compare the predicted number of claims for each profile.
and the tariff structures shown here are based on out-of-sample predictions. The GBM and the CANN model deviate slightly from perfect balance.

To compare tariff structures, we follow the methodology from Henckaerts and Antonio (2022) using risk scores. For a model \(f\), let \(F_{n}\) be the empirical cumulative distribution function of the predictions made by the model \(f\). For each policyholder \(i\), the risk score \(r_{i}^{f}\) is the evaluation of \(F_{n}\) in \(f(\mathbf{x}_{i})\). For frequency-severity modelling, with a frequency model \(f^{\text{freq}}\) and a severity model \(f^{\text{sev}}\), the risk score is calculated as

\begin{table}
\begin{tabular}{c c c c c} \hline \hline \multicolumn{1}{c}{**Observed and predicted losses**} & \multicolumn{1}{c}{**GDM**} & \multicolumn{1}{c}{**GBM**} & \multicolumn{1}{c}{**CANN**} & \multicolumn{1}{c}{**Surrogate GLM**} \\ \hline Australia (AU8) & 9 314 604 & 9 345 113 & 9 136 324 & 9 154 467 & 9 355 718 \\ Belgium (€ & 26 464 970 & 26 399 027 & 26 079 709 & 25 720 143 & 26 345 969 \\ France (€ & 58 872 147 & 56 053 341 & 56 207 993 & 58 629 584 & 57 048 375 \\ Norway (NOK) & 206 649 080 & 206 634 401 & 206 475 980 & 206 494 683 & - \\ \hline \multicolumn{1}{c}{**Ratio of predicted losses over observed losses**} & & & & \\ \hline Australia & - & 1.00 & 0.98 & 0.98 & 1.00 \\ Belgium & - & 1.00 & 0.99 & 0.97 & 1.00 \\ France & - & 0.95 & 0.95 & 1.00 & 0.97 \\ Norway & - & 1.00 & 1.00 & 1.00 & - \\ \hline \hline \end{tabular}
\end{table}
Table 6: Comparison between the total observed losses for each data set and the total predicted losses for the GLM, GBM, CANN GBM flex and the surrogate GLM. We also show the ratio of total predicted losses over observed losses.

Figure 14: Comparison between the low, medium and high risk profiles in the frequency models on the French data set according to the Shapley values of the CANN GBM flexible model (top row) and the fitted coefficients in the surrogate GLM (bottom row).
\[r_{i}^{f}=F_{n}\left\{f^{\text{freq}}(\mathbf{x}_{i})\times f^{\text{sev}}(\mathbf{x}_{i}) \right\}. \tag{13}\]

We compare the risk scores of multiple models using Lorenz curves (Lorenz, 1905). For a model \(f\), the Lorenz curve evaluated in \(s\in[0,1]\) is

\[LC^{f}(s)=\frac{\sum_{i=1}^{n}L_{i}\,\mathbbm{1}\{r_{i}^{f}\leq s\}}{\sum_{i=1} ^{n}L_{i}}, \tag{14}\]

with \(L_{i}\) the observed loss for policyholder \(i\). We visualize the Lorenz curve by plotting the pairs \(\left(s,LC^{f}(s)\right)\), for \(s\in[0,1]\). The Lorenz curves shows the accumulation of losses, ordered by risk score \(r_{i}^{f}\) obtained with model \(f\). A model with a better risk classification accumulates losses slower for low risk scores and faster for high risk scores. A Lorenz curve further away from the equality line at 45%, represents a better risk classification.

Figure 15 shows the Lorenz curves for the four data sets in the benchmark study where the top row compares the benchmark GBM with the CANN GBM flexible model and the bottom row the benchmark GLM with the surrogate GLM. For the Australian and French data sets, the tariff structure from the CANN model is (slightly) preferred to that of the benchmark GBM according to the Lorenz curve. For the Belgian and Norwegian data sets, the CANN model is also preferred, but the two curves are very similar. For the Australian, Belgian and French data sets, the Lorenz curves of the benchmark GLM and the surrogate GLM show a very similar pattern, with the surrogate model having a preferable curve over the benchmark GLM, showing that the higher predictive accuracy of the surrogate model also results in a slightly better risk classification in the tariff structure.

Figure 15: Lorenz curve comparison between the GBM benchmark model and the CANN GBM flexible model in the top row. Bottom row compares the GLM benchmark with the surrogate GLM. A dashed line is added to show line of equality.
## 7 Conclusion

This paper explores the potential of deep learning models for the analysis of tabular frequency and severity data in non-life insurance pricing. We detail a benchmark study using extensive cross-validation on multiple model architectures. Categorical input features are embedded by use of an autoencoder of which the encoder part is integrated into the FFNN and CANN structures. Our results demonstrate the performance gains achieved using the autoencoder embedding technique for categorical input variables, especially in modelling claim severities. Interpretation techniques are applied to both frequency and severity modelling.

The literature often questions the value created when analyzing tabular data with deep learning models. Indeed, our feed-forward neural network does not improve upon a carefully designed GBM. Combining gradient-boosted trees and neural networks leads to higher accuracy for frequency modelling. This aligns with what we see in other fields, where GBM and neural network combinations outperform the corresponding stand-alone models on tabular data. In modelling the severity data, out-of-sample deviances are relatively similar across benchmark models and deep learning architectures. This suggests that the added value created by using the deep learning approach is limited when applied to these datasets. Data sets with a high dimensional set of input variables and/or complex input features might benefit more from using deep learning models to model claim frequency and severity data. Gao et al. (2022), for instance, use deep learning models to analyze high-frequency time series of telematics data.

The end result of our study is a technical tariff structure on categorized input data, using a GLM that is carefully designed as a global surrogate for the deep learning model. Hence, this surrogate GLM leverages the insights from the deep learning models. The surrogate GLMs lead to a lower out-of-sample deviance than the benchmark GLMs and have a better risk differentiation. The workflow to construct a GLM as a global surrogate for a deep learning model can potentially be of interest to insurance companies aiming to harvest refined insights from their available data while aiming for an interpretable, explainable technical tariff. The latter consideration is of much interest in light of the GDPR algorithmic accountability clause. Further research could look into data sets with high dimensional feature sets, including features with images, text or times series, and explore the value of deep learning architectures via a carefully designed benchmark study, as outlined in this paper.

## Acknowledgements

Katrien Antonio gratefully acknowledges funding from the FWO and Fonds De La Recherche Scientifique - FNRS (F.R.S.-FNRS) under the Excellence of Science (EOS) program, project ASTeRISK Research Foundation Flanders [grant number 40007517]. The authors gratefully acknowledge support from the Ageas research chair on insurance analytics at KU Leuven, from the Chaire DIALog sponored by CNP Assurances and the FWO network W001021N. The authors thank Simon Gielis for his contributions (as MSc student) in the early stage of the research.

## Declaration of interest statement

The authors declare no potential conflict of interests.
## Supplemental material

The results in this paper were obtained using R. All code is available through github: [https://github.com/freekholvoet/NNforFreqSevPricing](https://github.com/freekholvoet/NNforFreqSevPricing). An R Markdown demonstration is available on that github page, called NNforFreqSevPricing.nb.html.
