Neural networks for insurance pricing with frequency and severity data:

a benchmark study from data preprocessing to technical tariff
a,c]Freek Holvoet
a,b,c]Katrien Antonio
a]Roel Henckaerts
[a]Faculty of Economics and Business, KU Leuven, Belgium.
[b]Faculty of Economics and Business, University of Amsterdam, The Netherlands.
[c]LRisk, Leuven Research Center on Insurance and Financial Risk Analysis, KU Leuven, Belgium.

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

§ 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 Wüthrich 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. Wüthrich, 2019 compares the performance of neural networks and GLMs on a frequency case study. Both Wüthrich, 2019 and Schelldorfer and Wüthrich, 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.

code-before = [HTML]FFFFFF1,3,5,7,9,11,13
Contribution     Categorical treatment     Model architecture     # Data sets     Case study     Interpretation tools

height 1.0pt
Dugas et al., 2003     $-$     LR, GLM, DT, NN, SVM     $1$     Tech. tariff     $-$

height 0.3pt
Yang et al., 2018     $-$     TDBoost     $1$     Tweedie compound     PDP, VIP

height 0.3pt
Henckaerts et al., 2018     $-$     GLM     $1$     Freq, sev     $-$

height 0.3pt
Wüthrich, 2019      Dummy encoding, embedding layers     GLM, NN, CANN     $1$     Freq     Avg. neuron activation

height 0.3pt
Schelldorfer and Wüthrich, 2019      Embedding layers     CANN     $1$     Freq     $-$

height 0.3pt
Ferrario et al., 2020     One-hot encoding     Boosted trees, NN     $1$     Freq     $-$

height 0.3pt
Noll et al., 2020     Dummy encoding, empirical means, one-hot encoding     GLM, DT, Boosted trees, NN     $1$     Freq     Loss-per-label

height 0.3pt
Henckaerts et al., 2021     $-$     DT, RF, GBM     $1$     Freq, sev, tech. tariff     PDP, VIP, ICE

height 0.3pt
Kuo and Richman, 2021     One-hot encoding, embedding layers, attention layers     GLM, NN, Transformer, TabNET     $1$     Sev     $-$

height 0.3pt
Delong and Kozak, 2023      Autoencoder     NN     $1$     Freq     $-$

height 0.3pt
Meng et al., 2022     Convolutional autoencoder     GLM     $1$     Freq     $-$

height 0.3pt
Henckaerts et al., 2022     $-$     GBM     $6$     Freq     PDP, SHAP, Surrogates

height 0.3pt
This paper     Autoencoder     GLM, GBM, NN, CANN     $4$     Freq, sev, tech. tariff     PDP, VIP, Surrogates, Shapley

Overview of the literature on deep learning for insurance pricing. Categorical treatment describes the preprocessing steps taken for categorical input variables. Model architecture lists the different models compared in the paper. We give the number of benchmark data sets
in # Data sets and the focus of the case study in Case study. Lastly, Interpretational tools list the tools used to look under the hood of the fitted models. The used abbreviations are: linear regression (LR), generalized linear model (GLM), decision tree (DT), neural network (NN), support vector machine (SVM), Tweedie boosted tree model (TDBoost), combined actuarial neural network (CANN), random forest (RF), gradient boosting model (GBM), frequency case study (freq), severity case study (sev), technical tariff structure (tech. tariff), partial dependency plot (PDP), variable importance plot (VIP), individual conditional expectation (ICE) and Shapley additive explanations (SHAP).

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 <ref> 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 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 [Schelldorfer and Wüthrich, 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.

§ 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 $\D = \left(\x_i,y_i\right)^{n}_{i=1}$, where each $\x_i$ is a $p$-dimensional data point with response $y_i$. Each data point $\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 $(\x,y)$, with $\x=(x_1,\ldots,x_p)$, each $x_j$ representing a variable in our data set $\D$.

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

\[\D = \big(\underbrace{x_1,\ldots,x_{p-c}\strut}_\text{numerical variables},\underbrace{x_{p-c+1},\ldots,x_p\strut}_\text{categorical variables},\underbrace{y\strut}_\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 $\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 $\text{freq}$ and use $\mathcal{D}$ and $\mathcal{D}^{\text{sev}}$. We denote the number of observations as $n_f$ for $\D$ and $n_s$ for $\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 $\left(x_1,\ldots,x_p\right)$ 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 $\D^{\text{train}}\subset\D$, by choosing the model-specific parameters that minimize a chosen loss function $\sum_{\x_i\in\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 $\D^{\text{test}} = \D\backslash\D^{\text{train}}$ as $\sum_{\x_i\in\D^{\text{test}}} \mathscr{L}\left(\hat{y}_i,y_i\right)$. We follow the loss functions proposed by Wüthrich 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:

\begin{equation}
D_{\text{Poisson}}(f(\boldsymbol{x}),\boldsymbol{y}) = \frac{2}{n_f} \sum_{i=1}^{n_f}\l y_i\ln\frac{y_i}{f(\boldsymbol{x}_i)}-(y_i-f(\boldsymbol{x}_i)) \r.
\label{eq_devpoisson}
\end{equation}

Note that when using the exposure-to-risk $e$ in the frequency model, we replace each prediction $f(\boldsymbol{x}_i)$ with $e_i\cdot f(\boldsymbol{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

\begin{equation}
D_{\text{gamma}}(f(\boldsymbol{x}),\boldsymbol{y}) = \frac{2}{n_s}\sum_{i=1}^{n_s} \alpha_i\l\frac{y_i-f(\boldsymbol{x}_i)}{f(\boldsymbol{x}_i)}-\ln\frac{y_i}{f(\boldsymbol{x}_i)} \r,
\label{eq_devgamma}
\end{equation}

where the weight $\alpha_i$ is the observed number of claims for data point $i$.

§ DEEP LEARNING ARCHITECTURES AND PREPROCESSING STEPS

§.§ Neural network architectures

Feed-forward neural network
A feed-forward neural network (FFNN) is a type of machine learning model that utilizes interconnected layers, represented by $\boldsymbol{z}^{(m)}$ with $m=0,\ldots,M+1$. The input layer, represented by $\boldsymbol{z}^{(0)}$, provides the network with input data, while the output layer, represented by $\boldsymbol{z}^{(M+1)}$, gives the network's prediction. Between the input and output layers, there can be one or more hidden layers, represented by $\boldsymbol{z}^{(1)}, \ldots,\boldsymbol{z}^{(M)}$. When there are two or more hidden layers, we call the neural network a deep learning model. Each layer $\boldsymbol{z}^{(m)}$ consists of $q_m$ nodes, so it can be expressed as a vector $\boldsymbol{z}^{(m)} = \left(z^{(m)}_1,\ldots,z^{(m)}_{q_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 $\boldsymbol{b}_m\in\mathbb{R}^{q_{m}}$. An activation function $\sigma^{(m)}$, $m=1,\ldots,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 $\boldsymbol{z}^{(m)}$ can be written in function of the previous layer as follows:

\begin{equation}
\boldsymbol{z}^{(m)} = \sigma^{(m)}\left(W_m\cdot\boldsymbol{z}^{(m-1)} + \boldsymbol{b}_m\right).
\label{eq_nnlayer}
\end{equation}

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 $\boldsymbol{z}^{(m)}$ for input $\x_i$ is denoted as $\boldsymbol{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 $\D$, the dimension $q_0$ of the input layer is equal to the number of variables $p$ in the data set[The dimension of the input layer can be larger than $p$ when using an encoding technique, such as one-hot encoding.]. We write the input layer as $\left(x_1,\ldots,x_p\right)$ 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 <ref> gives a schematic overview of a feed-forward neural network.

Structure of a feed-forward neural network with $p$-dimensional input layer, hidden layers $\boldsymbol{z}^{(1)}, \ldots,\boldsymbol{z}^{(M)}$, with $q_1,\ldots,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 networks
Wüthrich, 2019 and [Schelldorfer and Wüthrich, 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 <ref> shows the structure of the CANN model.

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 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:

\begin{equation}
\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).
\label{eq_cannoutput}
\end{equation}

The case study in [Schelldorfer and Wüthrich, 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 <ref> 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 variables

We 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:

\begin{equation}
\label{eq_normalization}
\end{equation}

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 $\D^{\text{train}}\subset\D$ to train the model, we calculate the $\mu_{x_j}$ and $\sigma_{x_j}$ only on the data set $\D^{\text{train}}$ to avoid data leakage.

Categorical variables
The 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 <ref> shows the general structure of such an autoencoder. It consists of an input layer of dimension $c$, one hidden layer $\boldsymbol{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 $\boldsymbol{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 $\boldsymbol{b}_{\text{dec}}\in\mathbb{R}^c$. For input $\boldsymbol{x}_i\in\mathbb{R}^c$, the encoded and decoded representations are calculated as

\begin{equation}
\begin{aligned}
& \boldsymbol{z}_i^{\text{enc}} = \sigma^{(\text{enc})}\left(W_{\text{enc}}\cdot\boldsymbol{x}_i+\boldsymbol{b}_{\text{enc}}\right), \\
& \boldsymbol{x}_i^{\text{dec}} = \sigma^{(\text{dec})}\left(W_{\text{dec}}\cdot\boldsymbol{z}_i^{\text{enc}}+\boldsymbol{b}_{\text{dec}}\right). \\
\end{aligned}
\label{eq_ae_calculation}
\end{equation}

General structure of an autoencoder neural network. The connections between the input layer and the encoded layer are called the encoder, highlighted in red, the connections between encoded layer and decoded output layer are called the decoder. After training, the encoded vector is a $d$-dimensional representation of the $c$-dimensional input.

The CANN model where the last $c$ input nodes represent the one-hot encoding of the categorical variables in the data set. Using the encoder part of the pretrained autoencoder, shown with the red box, we embed the categorical variables into a $d$-dimensional numerical representation, and connect the embeddings to the CANN model structure.

Our proposed network structure combines the autoencoder embedding technique from [Delong and Kozak, 2023] and the CANN structure from
[Schelldorfer and Wüthrich, 2019].

The autoencoder is trained on all data points $\boldsymbol{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_{\boldsymbol{x}_i\in\mathcal{D}}^{}\mathscr{L} \left( \boldsymbol{x}^{\text{dec}}_i, \boldsymbol{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^{(j)}_1,\ldots,x^{(j)}_{L_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 (<ref>).

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^{(j)}_1,\ldots,x^{(j)}_{L_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^{(j)}_1,\ldots,x^{(j)}_{L_j}\right)$ and for each $h\in\{1,\ldots,L_j\}$, the softmax transformation of the output node $x^{(j,\text{dec})}_h$ is defined as:

\begin{equation}
x^{(j,\text{dec})}_h\mapsto\tilde{x}^{(j,\text{dec})}_h = \frac{\exp\left(x^{(j,\text{dec})}_h\right)}{\exp\left(x^{(j,\text{dec})}_1+\ldots+x^{(j,\text{dec})}_{L_j}\right)},\qquad h=1,\ldots,L_j.
\label{eq_softmaxtranform}
\end{equation}

The use of the softmax activation function ensures that the values of the decoded vectors $\left(x^{(j,\text{dec})}_1,\ldots,x^{(j,\text{dec})}_{L_j}\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 $\x_i^{\text{OH}}$ the one-hot encoding of all categorical variables for policyholder $i$ and $\tilde{\boldsymbol{x}}^{\text{dec}}_i$ the values of the autoencoder's output layer for policyholder $i$, the cross-entropy loss function is defined as:

\begin{equation} \mathscr{L}^{\text{CE}} \left( \tilde{\boldsymbol{x}}^{\text{dec}}_i, \boldsymbol{x}_i^{\text{OH}} \right) = - \sum_{j=1}^{c} \sum_{h=1}^{L_j} x_{ih}^{(j)}\cdot\log\left(\tilde{x}^{(j,\text{dec})}_{ih}\right) .
\label{eq_bincrossentropy}
\end{equation}

After training the autoencoder and applying the trained autoencoder on each policyholder $\x_i\in\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^{\text{enc}}_{i1},\ldots,z^{\text{enc}}_{id})$ as calculated by Equation (<ref>). We call the vector $\z^{\text{enc}}_i$ the embedding of the categorical inputs of $\x_i$. To use the embedding together with the numerical features of $\x_i$, we normalize the values in the nodes $z^{\text{enc}}_1,\ldots,z^{\text{enc}}_d$ by scaling the weight matrix $W_{\text{enc}}$ and bias vector $\boldsymbol{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^{\text{enc}}_{i1},\ldots,z^{\text{enc}}_{id}$ for all $\x_i\in\D^{\text{train}}$, we scale the weight matrix $W_{\text{enc}}$ and bias vector $\boldsymbol{b}_{\text{enc}}$ of the pre-trained encoder as follows:

\begin{equation}
W_{\text{enc}}\mapsto \tilde{W}_{\text{enc}} = \begin{pmatrix}
\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_{dc}}{\sigma_d}\\
\end{pmatrix} , \boldsymbol{b}_{\text{enc}}\mapsto \tilde{\boldsymbol{b}}_{\text{enc}} = \begin{pmatrix}
\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 - \mu_d}{\sigma_{d}}\\
\end{pmatrix}.
\label{eq_normmatrix}
\end{equation}

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 <ref> and <ref> with the encoding part of the trained autoencoder, as shown in Figure <ref> for the CANN. For clarity, we omit the one-hot encoding notation of each variable in Figure <ref>. 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(<ref>) or Equation (<ref>).

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 $\left[t_{k,\text{min}},t_{k,\text{max}}\right]$. The search space $\mathcal{S}$ is the space consisting of all possible values for all tuning parameters:

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

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 <ref>, 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.

Example of random grid search with two tuning parameters $t_1$ and $t_2$. The search space $\mathcal{S} = \left[t_{1,\text{min}},t_{1,\text{max}}\right]\times\left[t_{2,\text{min}},t_{2,\text{max}}\right]$ 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.

We use the extensive cross-validation scheme proposed by Henckaerts et al., 2021, as sketched in Figure <ref>. 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.

Representation of the $6$ times $5$-fold cross-validation scheme, figure from Henckaerts et al., 2021.

§ PERFORMATIVE COMPARISON BETWEEN BENCHMARK MODELS AND DEEP LEARNING ARCHITECTURES

Section <ref> 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 models introduced in Section <ref>. Section <ref> 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 <ref>. Lastly, Section <ref> compares the autoencoder embedding against using the one-hot encoding when used in the deep learning models under consideration.

§.§ Data sets

The used data sets are an Australian, Belgian, French[The French data set in the  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.] and Norwegian MTPL data set, available through the R packages  [Dutang and Charpentier, 2019] and  [Henckaerts and Antonio, 2022, Henckaerts, 2021]. Table <ref> gives an overview of the number of records for each data set and the number of continuous, categorical and spatial variables.

code-before = [HTML]FFFFFF1,2,4,5,7,9
Australian MTPL     Belgian MTPL     French MTPL     Norwegian MTPL

height 0.3pt

6lNumber of observations

height 0.3pt
Frequency     $67\,856$      $163\,212$     $668\,897$      $183\,999$

Severity     $4\,624$      $18\,276$     $24\,944$     $8\,444$

height 0.3pt

6lCovariates: number and type

height 0.3pt
Continuous     1     4     2     0

Categorical     4     5     5     3

Spatial     0     1     2     1

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.

The spatial variables are listed separately in Table <ref>. 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 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 (<ref>) and Equation (<ref>) are used for, respectively, frequency and severity modelling. We follow the repeated $5$-fold cross-validation scheme as described in Section <ref>. 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 (<ref>). 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 <ref>.

code-before = [HTML]FFFFFF1,2,5
Tuning parameter     Range

height 0.3pt
Activation function for hidden layers     ReLU, sigmoid, softmax

Batch size     $[10\,000,50\,000]$     Frequency

$[200,10\,000]$     Severity

Number of hidden layers     $[\num{1},\num{4}]$

Nodes per hidden layer     $[\num{10},\num{50}]$

Dropout rate     $[\num{0},\num{0.1}]$

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.

The softmax activation function is applied on the collection of all nodes in a hidden layer, rather than individually at each node within the layer.

The cross-validation scheme is shown in Algorithm <ref>, 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 (<ref>) and (<ref>). 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.

model class () and corresponding tuning grid $\mathcal{R}$
data $\mathcal{D}$ with $6$ disjoint stratified subsets $\mathcal{D}_1,\ldots,\mathcal{D}_6$
$\ell=1,\ldots, 6$
leave out $\mathcal{D}_{\ell}$ as test set
pre-processing steps
continuous variable $x_j\in\mathcal{D}$
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 (<ref>)

categorical variable $x_j\in\mathcal{D}$
construct one-hot encoding $\left(x^{(j)}_1,\ldots,x^{(j)}_{L_j}\right)$ of variable $x_j$

train an autoencoder $f^{\text{AE}}_{d}$ 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}^{\text{BCE}}(\cdot,\cdot)$

select $f^{\text{AE}}_{d}$ with the lowest $d$ while $\mathscr{L}^{\text{BCE}}(\cdot,\cdot) < 0.001$
calculate the scaled $\tilde{W}_{\text{enc}}$ and bias vector $\tilde{\boldsymbol{b}}_{\text{enc}}$ for the encoder part of $f^{\text{AE}}_{d}$
add the encoder part of $f^{\text{AE}}_{d}$, with $\tilde{W}_{\text{enc}}$ and $\tilde{\boldsymbol{b}}_{\text{enc}}$, to the model class

cross-validation steps
tuning parameter point $s\in\mathcal{R}$
$k \in \left\lbrace 1,\ldots,6\right\rbrace  \setminus \ell$
train a model $f_{\ell k}$ of  on $\mathcal{D} \setminus \left\lbrace \mathcal{D}_{\ell},\mathcal{D}_k \right\rbrace $
evaluate the model performance on $\mathcal{D}_k$ using loss function $\mathscr{L}(\cdot,\cdot)$
$\text{valid\_error}_{\ell k} \leftarrow \frac{1}{|\mathcal{D}_k|}\underset{\boldsymbol{x}_i \in \mathcal{D}_k}{\sum} \mathscr{L}\{y_i,f_{\ell k}(\boldsymbol{x}_i)\}$

$\text{valid\_error}_\ell \leftarrow \frac{1}{5}\sum_{k \in \left\lbrace 1,\ldots,6\right\rbrace  \setminus \ell}\text{valid\_error}_{\ell k}$

optimal parameter point $s^{*}_{\ell}\in\mathcal{R}$ minimizes $\text{valid\_error}_\ell$

out-of-sample performance
rep $\in\{1,2,3\}$
select a validation set $\mathcal{D}_{\text{val}}$ containing $20\%$ of the records in $\mathcal{D} \setminus \mathcal{D}_\ell$
train a model $f_{\ell,\text{rep}}$ of  on $\mathcal{D} \setminus \{ \mathcal{D}_\ell, \mathcal{D}_{\text{val}} \}$ using the optimal parameter point $s^{*}_{\ell}$ and using $\mathcal{D}_{\text{val}}$ for early stopping
evaluate the model performance on $\mathcal{D}_\ell$ using loss function $\mathscr{L}(\cdot,\cdot)$
$\text{test\_error}_{\ell,\text{rep}} \leftarrow \frac{1}{|\mathcal{D}_\ell|}\underset{\boldsymbol{x}_i \in \mathcal{D}_\ell}{\sum} \mathscr{L}\{y_i,f_{\ell,\text{rep}}(\boldsymbol{x}_i)\}$

$\text{test\_error}_{\ell} \leftarrow \frac{1}{3}\sum_{\text{rep}} \text{test\_error}_{\ell,\text{rep}}$

optimal tuning parameters + performance measure for each of the six folds.

Pseudocode to sketch the pipeline for calculating out-of-sample performances with the neural network structures from Section <ref>, 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 <ref>. We have out-of-sample deviances for each withheld test set, measured in Poisson deviance (<ref>) or gamma deviance (<ref>). We show frequency on the left-hand side and severity on the right-hand side.

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].

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$.

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.

Figure <ref> 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.

§ 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

\begin{equation}
\text{VIP}_{x_j} = \sum_{\x_i\in\D} \text{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),
\end{equation}

where $\tilde{x}_{i,j}$ is a random permutation of the values observed for $x_j$ in the data set $\D$. A large value for $\text{VIP}_{x_j}$ indicates that the variable significantly influences the model output and is therefore considered important. Figure <ref> 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

\begin{equation}
\overline{\text{VIP}}_{x_j} = \frac{\text{VIP}_{x_j}}{\sum_{x_j\in\mathcal{D}} \text{VIP}_{x_j}},\qquad\text{where the sum runs over all over all variables}\,\,x_j.
\label{eq_vip}
\end{equation}

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.

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

\begin{equation}
\text{PD}_{x_j} = \left\{\frac{1}{|\mathcal{D}|}\sum_{\x_i\in\mathcal{D}}f\left(x_{i,1},\ldots,X_{o,j},\ldots,x_{i,p}\right)\,;\, \forall X_{o,j} \in X_j\right\},
\label{eq_pdp}
\end{equation}

where $\left(x_{i,1},\ldots,X_{o,j},\ldots,x_{i,p}\right)$ is the data point $\x_i\in\mathcal{D}$ with element $x_{i,j}$ replaced by the value $X_{o,j}\in X_j$. The vector $\text{PD}_{x_j}$ can be seen as the average prediction on $\mathcal{D}$, while letting the 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 (<ref>) can be extended to a two-way interaction partial dependence effect by letting two variables range over their respective variable spaces.

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.

Figure <ref> 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.

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 <ref> displays the partial dependence effect of the policyholder's age when calibrated on the claim severity data. Similar to the effects portrayed in Figure <ref>, 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 for the French data, the FFNN model reveals an almost flat effect across all age groups.

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.

Figure <ref> 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.

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.

We consider the partial dependence effect of the postal code in the Belgian frequency data set in Figure <ref>. 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 construction
[Henckaerts 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 <ref> and calibrated in Section <ref>. 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 <ref>. 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 GLM on the segmented frequency and severity data leads to a tabular premium structure incorporating the insights captured by the CANN architectures.

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.

Figure <ref> shows the partial dependence effects of the CANN GBM flexible for the bonus-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.

code-before = [HTML]FFFFFF1,3,5,6
Binned GLM     Surrogate GLM

height 0.3pt
VehValue     VehValue

DrivAge     DrivAge

VehValue $\times$ DrivAge

height 0.3pt
$\mathbf{0.3817}$     $\mathbf{0.3805}$

Australian data set

code-before = [HTML]FFFFFF1,3,5,7,9,11,12
Binned GLM     Surrogate GLM

height 0.3pt
ageph     ageph

bm     bm

power     power

coverage      coverage

fuel     fuel

postcode     postcode

ageph $\times$ power

height 0.3pt
$\mathbf{0.5314}$     $\mathbf{0.5308}$

Belgian data set

code-before = [HTML]FFFFFF1,3,5,7,9,11,13,15,17,19
Binned GLM     Surrogate GLM

height 0.3pt
VehValue     VehPower

VehAge     VehAge

DrivAge     DrivAge

BonusMalus     BonusMalus

VehGas     VehGas

Density     Density

VehBrand $\times$ Region

Density $\times$ Region

VehPower $\times$ Region

VehPower $\times$ VehBrand

DrivAge $\times$ Region

BonusMalus $\times$ Region

VehBrand $\times$ Density

BonusMalus $\times$ VehBrand

DrivAge $\times$ BonusMalus

height 0.3pt
$\mathbf{0.2761}$     $\mathbf{0.2738}$

French data set

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$.

Table <ref> 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 profiles

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

code-before = [HTML]FFFFFF1,3,5,7,9,10,12
Variables     Low risk     Medium risk     High risk

height 0.3pt
Vehicle power     $4$     $6$     $9$

Vehicle age     $3$     $2$      $1$

Policyholder age     $[21,26[$     4 $[30,40[$     $\ge70$

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

height 0.3pt

4lPredicted number of claims

height 0.3pt
Surrogate GLM     $0.020$     $0.106$     $0.361$

CANN GBM flexible     $0.021$     $0.101$     $0.519$

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.

We construct a low, medium, and high-risk profile based on the frequency surrogate GLM for the French data set. Table <ref> 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 <ref>. For the GLM, we 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.

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).

§ 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 <ref> 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, and the tariff structures shown here are based on out-of-sample predictions. The GBM and the CANN model deviate slightly from perfect balance.

code-before = [HTML]FFFFFF1,2,4,5,7,9,11

Observed     GLM     GBM     CANN GBM flex     Surrogate GLM

height 0.3pt

7lObserved and predicted losses

height 0.3pt
Australia (AU$)     $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$     -

height 0.3pt

7lRatio of predicted losses over observed losses

height 0.3pt
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$     -

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.

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^f_i$ is the evaluation of $F_n$ in $f(\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{equation}
r^f_i = F_n\left\{f^{\text{freq}}(\x_i) \times f^{\text{sev}}(\x_i)\right\}.
\label{eq_riskscore}
\end{equation}

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

\begin{equation}
LC^f(s) = \frac{\sum_{i=1}^{n} L_i\,\mathbbm{1}\{r_i^f\le s\}}{\sum_{i=1}^{n} L_i},
\label{eq_lorenz}
\end{equation}

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^f_i$ 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.

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 equility.

Figure <ref> 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.

§ 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 W$001021$N. 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>.

An R Markdown demonstration is available on that github page, called NNforFreqSevPricing.nb.html.

[Antonio et al., 2010]
Katrien Antonio, Edward W Frees, and Emiliano A Valdez.
A multilevel analysis of intercompany claim counts.
ASTIN Bulletin: The Journal of the International Actuarial
Association, 400 (1):0 151–177, 2010.

[Bergstra and Bengio, 2012]
James Bergstra and Yoshua Bengio.
Random search for hyper-parameter optimization.
Journal of machine learning research, 130
(2):0 281–305, 2012.

[Borisov et al., 2022]
Vadim Borisov, Tobias Leemann, Kathrin Sessler, Johannes Haug, Martin
Pawelczyk, and Gjergji Kasneci.
Deep neural networks and tabular data: A survey.
IEEE transaction on neural networks and learning systems,
PP:0 1–21, 2022.

[De Jong and Heller, 2008]
Piet De Jong and Gillian Z. Heller.
Generalized linear models for insurance data.
Cambridge University Press, 2008.

[Delong and Kozak, 2023]
Łukasz Delong and Anna Kozak.
The use of autoencoders for training neural networks with mixed
categorical and numerical features.
ASTIN Bulletin: The Journal of the ternational Actuarial
Association, page 1–20, 2023.

[Denuit et al., 2019]
Michel. Denuit, Donatien Hainaut, and Julien Trufin.
Effective Statistical Learning Methods for Actuaries I: GLMs
and Extensions.
Springer Actuarial Lecture Notes. Springer International Publishing,

[Denuit et al., 2020]
Michel. Denuit, Donatien Hainaut, and Julien Trufin.
Effective Statistical Learning Methods for Actuaries II:
Tree-Based Methods and Extensions.
Springer Actuarial Lecture Notes. Springer International Publishing,

[Dugas et al., 2003]
Charles Dugas, Yoshua Bengio, Nicolas Chapados, Pascal Vincent, Germain
Denoncourt, and Christian Fournier.
Statistical learning algorithms applied to automobile insurance
In CAS Forum, volume 1, pages 179–214. Citeseer, 2003.

[Dutang and Charpentier, 2019]
Christophe Dutang and Arthur Charpentier.
CASdatasets: insurance datasets.
<http://cas.uqam.ca>, 2019.
R package version 1.0.10.

[Ferrario et al., 2020]
Andrea Ferrario, Alexander Noll, and Mario V Wüthrich.
Insights from inside neural networks.
Available at SSRN 3226852, 2020.

[Frees and Valdez, 2008]
Edward W Frees and Emiliano A Valdez.
Hierarchical insurance claims modeling.
Journal of the American Statistical Association, 1030
(484):0 1457–1469, 2008.

[Friedman, 2001]
Jerome H Friedman.
Greedy function approximation: a gradient boosting machine.
The Annals of statistics, 29:0 1189–1232, 2001.

[Gao et al., 2022]
Guangyuan Gao, Shengwang Meng, and Mario V. Wüthrich.
What can we learn from telematics car driving data: A survey.
Insurance: Mathematics and Economics, 104:0 185–199,

[Gielis, 2020]
Simon Gielis.
Neural networks for non-life insurance pricing, 2020.
[Master's thesis, KU Leuven].

[Goodfellow et al., 2016]
Ian Goodfellow, Yoshua Bengio, and Aaron Courville.
Deep Learning.
MIT Press, 2016.

[Grinsztajn et al., 2022]
Leo Grinsztajn, Edouard Oyallon, and Gael Varoquaux.
Why do tree-based models still outperform deep learning on typical
tabular data?
In Thirty-sixth Conference on Neural Information Processing
Systems Datasets and Benchmarks Track, 2022.

[Guo and Berkhahn, 2016]
Cheng Guo and Felix Berkhahn.
Entity embeddings of categorical variables.
arXiv preprint, 2016.
URL <https://arxiv.org/abs/1604.06737>.

[Haberman and Renshaw, 1996]
Steven Haberman and Arthur E. Renshaw.
Generalized linear models and actuarial science.
Journal of the Royal Statistical Society. Series D (The
Statistician), 450 (4):0 407–436, 1996.

[Hainaut et al., 2022]
Donatien Hainaut, Julien Trufin, and Michel Denuit.
Response versus gradient boosting trees, GLMs and neural networks
under Tweedie loss and log-link.
Scandinavian Actuarial Journal, 20220 (10):0
841–866, 2022.

[Hastie et al., 2009]
Trevor Hastie, Robert Tibshirani, and Jerome H Friedman.
The elements of statistical learning: data mining, inference,
and prediction, volume 2.
Springer, 2009.

[Henckaerts, 2021]
Roel Henckaerts.
maidrr: Model-Agnostic Interpretable Data-driven
<http://github.com/henckr/maidrr>, 2021.
R package version 1.0.0.

[Henckaerts and Antonio, 2022]
Roel Henckaerts and Katrien Antonio.
The added value of dynamically updating motor insurance prices with
telematics collected driving behavior data.
Insurance: Mathematics and Economics, 105:0 79–95,

[Henckaerts et al., 2018]
Roel Henckaerts, Katrien Antonio, Maxime Clijsters, and Roel Verbelen.
A data driven binning strategy for the construction of insurance
tariff classes.
Scandinavian Actuarial Journal, 20180 (8):0
681–705, 2018.

[Henckaerts et al., 2021]
Roel Henckaerts, Marie-Pier Côté, Katrien Antonio, and Roel Verbelen.
Boosting insights in insurance tariff plans with tree-based machine
learning methods.
North American Actuarial Journal, 250 (2):0
255–285, 2021.

[Henckaerts et al., 2022]
Roel Henckaerts, Katrien Antonio, and Marie-Pier Côté.
When stakes are high: balancing accuracy and transparency with
Model-Agnostic Interpretable Data-driven suRRogates.
Expert Systems with Applications, 202:0 117230–,

[Ke et al., 2018]
Guolin Ke, Jia Zhang, Zhenhui Xu, Jiang Bian, and Tie-Yan Liu.
TabNN: a universal neural network solution for tabular data.
In URL https://openreview. net/forum, 2018.

[Ke et al., 2019]
Guolin Ke, Zhenhui Xu, Jia Zhang, Jiang Bian, and Tie-Yan Liu.
DeepGBM: A deep learning framework distilled by GBDT for online
prediction tasks.
In Proceedings of the 25th ACM SIGKDD International Conference
on Knowledge Discovery & Data Mining, pages 384–394, 2019.

[Keskar et al., 2016]
Nitish Shirish Keskar, Dheevatsa Mudigere, Jorge Nocedal, Mikhail Smelyanskiy,
and Ping Tak Peter Tang.
On large-batch training for deep learning: Generalization gap and
sharp minima.
arXiv preprint arXiv:1609.04836, 2016.

[Kingma and Ba, 2014]
Diederik P Kingma and Jimmy Ba.
Adam: A method for stochastic optimization.
arXiv preprint arXiv:1412.6980, 2014.
URL <https://arxiv.org/abs/1412.6980>.

[Kuo and Richman, 2021]
Kevin Kuo and Ronald Richman.
Embeddings and attention in predictive modeling.
arXiv preprint, 2021.
URL <https://arxiv.org/abs/2104.03545>.

[Liu et al., 2014]
Yue Liu, Bing-Jie Wang, and Shao-Gao Lv.
Using multi-class Adaboost tree for prediction frequency of auto
Journal of Applied Finance and Banking, 40
(5):0 45, 2014.

[Lorenz, 1905]
Max O Lorenz.
Methods of measuring the concentration of wealth.
Publications of the American statistical association,
90 (70):0 209–219, 1905.

[Meng et al., 2022]
Shengwang Meng, He Wang, Yanlin Shi, and Guangyuan Gao.
Improving automobile insurance claims frequency prediction with
telematics car driving data.
ASTIN Bulletin: The Journal of the International Actuarial
Association, 520 (2):0 363–391, 2022.

[Nelder and Wedderburn, 1972]
John Ashworth Nelder and Robert WM Wedderburn.
Generalized linear models.
Journal of the Royal Statistical Society: Series A (General),
1350 (3):0 370–384, 1972.

[Noll et al., 2020]
Alexander Noll, Robert Salzmann, and Mario V Wüthrich.
Case study: French motor third-party liability claims.
Available at SSRN 3164764, 2020.

[Ohlsson and Johansson, 2010]
Esbjörn Ohlsson and Björn Johansson.
Non-life insurance pricing with generalized linear models,
volume 174.
Springer, 2010.

[Olden et al., 2004]
Julian D. Olden, Michael K. Joy, and Russell G. Death.
An accurate comparison of methods for quantifying variable importance
in artificial neural networks using simulated data.
Ecological modelling, 1780 (3-4):0 389–397,

[Schelldorfer and Wüthrich, 2019]
Jürg Schelldorfer and Mario V Wüthrich.
Nesting classical actuarial models into neural networks.
Available at SSRN 3320525, 2019.

[Shapley et al., 1953]
Lloyd S Shapley et al.
A value for n-person games.

[Shwartz-Ziv and Armon, 2022]
Ravid Shwartz-Ziv and Amitai Armon.
Tabular data: Deep learning is not all you need.
Information Fusion, 81:0 84–90, 2022.

[Sola and Sevilla, 1997]
J. Sola and J. Sevilla.
Importance of input data normalization for the application of neural
networks to complex industrial problems.
IEEE Transactions on Nuclear Science, 440
(3):0 1464–1468, 1997.

[Srivastava et al., 2014]
Nitish Srivastava, Geoffrey Hinton, Alex Krizhevsky, Ilya Sutskever, and Ruslan
Dropout: a simple way to prevent neural networks from overfitting.
The Journal of Machine Learning Research, 150
(1):0 1929–1958, 2014.

[Wang and Song, 2011]
Haizhou Wang and Mingzhou Song.
Ckmeans. 1d. dp: optimal k-means clustering in one dimension by
dynamic programming.
The R journal, 30 (2):0 29, 2011.

[Wüthrich, 2019]
Mario V Wüthrich.
From generalized linear models to neural networks, and back.
Available at SSRN 3491790, 2019.

[Wüthrich and Buser, 2021]
Mario V Wüthrich and Christoph Buser.
Data analytics for non-life insurance pricing.
Swiss Finance Institute Research Paper, 0 (16-68),

[Yang et al., 2018]
Yi Yang, Wei Qian, and Hui Zou.
Insurance premium prediction via gradient tree-boosted Tweedie
compound Poisson models.
Journal of Business & Economic Statistics, 360
(3):0 456–470, 2018.