# A Neural Frequency-Severity Model and Its Application to Insurance Claims

Dong-Young Lim  School of Mathematics, The University of Edinburgh, UK.

###### Abstract

This paper proposes a flexible and analytically tractable class of frequency-
severity models based on neural networks to parsimoniously capture important
empirical observations. In the proposed two-part model, mean functions of
frequency and severity distributions are characterized by neural networks to
incorporate the non-linearity of input variables. Furthermore, it is assumed
that the mean function of the severity distribution is an affine function of
the frequency variable to account for a potential linkage between frequency
and severity. Within our modelling framework, we provide explicit closed-form
formulas for the mean and variance of the aggregate loss. Components of the
proposed model including parameters of neural networks and distribution
parameters can be estimated by minimizing the associated negative log-
likelihood functionals with neural network architectures. Furthermore, we
leverage the Shapely value and recent developments in machine learning to
interpret the outputs of the model. Applications to a synthetic dataset and
insurance claims data illustrate that our method outperforms the existing
methods in terms of interpretability and predictive accuracy. Keywords:
frequency-severity model; dependence modelling;aggregate risk losses;
interpretable machine learning

##  1 Introduction

When we deal with a dataset that involves a joint distribution of frequency
and severity such as insurance claims, health care expenditures and
operational risk losses, frequency-severity models are widely considered. For
example, Chavez-Demoulin et al. [2006], Dahen and Dionne [2010] and Brechmann
et al. [2014] developed frequency-severity models to quantify and analyze
operational losses, and Frees et al. [2011a] used the aggregate loss model for
predicting health care expenditures. In particular, they are a fundamental
tool for risk management and assessment in insurance business. This paper
proposes a novel modelling framework based on neural networks which unifies
neural networks and frequency-severity models to address fundamental empirical
observations. We call this new framework the neural frequency-severity model,
which we abbreviate to NeurFS.

The first observation is that severity is positively(or negatively) related to
frequency. Motivated by the empirical feature, two approaches have been mainly
explored in the generalized linear model(GLM) framework. Gschlö$\beta$l and
Czado [2007], Frees et al. [2011a] and Garrido et al. [2016] used the
frequency variable, which is the number of claims in the context of actuarial
science, as an additional covariate to model the severity distribution,
allowing to incorporate the dependence between frequency and severity. On the
other hand, the traditional GLM has been extended by combining marginal
distributions for frequency and severity via a bivariate copula, for example,
see Czado et al. [2012] and Kramer et al. [2013]. Two ideas are extensively
reviewed and compared to traditional independent models in Shi et al. [2015].

Despite such developments in the literature, the GLM based models have their
obvious limitation in the sense that mean functions of marginal distributions
are restricted to a linear form, which is not enough to accommodate complex
relations in covariates. For example, Kelly and Nielson [2006] found that the
severity against age is ‘U’ shaped in North American auto insurance data,
indicating the strong evidence on the nonlinear relationship between severity
and age. In response to the fact, the second feature this paper considers is
that input variables can be non-linearly related to frequency and/or severity.
Recently, there is a growing body of literature extending frequency-severity
models based on machine learning techniques to capture the nonlinearity. For
example, Yang et al. [2018] introduced a gradient boosting algorithm for the
nonlinear Tweedie model and Zhou et al. [2020] extended it to the zero-
inflated Tweedie model. Huang and Meng [2019] investigated several machine
learning techniques in predicting frequency and severity based on telematics
driving data. However, these models need to assume that frequency and severity
are independent. In addition, Guelman [2012] and Su and Bai [2020] applied the
gradient boosting to modelling frequency and severity distributions that
allows their dependency. But, in their models, it is impossible to compute and
analyze the aggregate loss, which is the most important quantity in analyzing
frequency-severity data. More importantly, these models rely on the feature
importance measure and dependence plot to interpret the outputs of models.
Unfortunately, many researchers have revealed that these tools are
inconsistent and inaccurate, resulting in misleading decision-making when
deploying the model in a real business, see Lundberg and Lee [2017].

The goal of our modelling is to overcome the limitations stated above by
combining the advantages of the existing models and nonlinear regression based
on neural networks. In particular, the use of neural networks allows our model
to automatically extract complex relations in data and efficiently handle
high-dimensional and large scale data. Also, the feasibility of NeurFS is
fully supported by the availability of massive data, increase in computing
power and storage capacity, and advanced machine learning algorithms.

The key contributions of our work are summarized as follows:

  * •

Analytical tractability: We provide the explicit closed-form solutions for
mean and variance of the aggregate loss in the presence of the nonlinearity
and dependence structure in data.

  * •

Flexibility: A wide range of probability distributions can be utilized for
marginal frequency and severity distributions in our modelling framework. This
paper mainly considers the exponential family of distributions due to its
popularity.

  * •

Efficient estimation: We suggest efficient neural network architectures to
estimate parameters of the model. In our training scheme, the parameters of
neural networks and distribution parameters are simultaneously estimated by
solving the associated optimization problems.

  * •

Interpretability: The main components of NeurFS are neural networks, of which
the great expressivity is strength in modelling, but weakness in interpreting
the result. We improve the interpretability of our models by applying the
Shapley value-based explanation method, proposed by Lundberg and Lee [2017].

This paper is organized as follows. In Section 2, we develop NeurFS with some
assumptions. In Section 3, we present some background on neural networks and
notations. Also, we describe the estimation process and the network
architecture. Section 4 demonstrates the accuracy and effectiveness of our
model through an artificial example. In Section 5, we apply NeurFS to real
insurance claim data and explain the results via the SHAP model.

##  2 Neural frequency-severity model

###  2.1 Marginal distributions and assumptions

This subsection introduces two marginal distributions of frequency-severity
models with associated regularities to develop NeurFS. Let $N$ be the number
of claim counts for a fixed time period. We interchangeably use the term,
frequency, to represent the random variable $N$. We assume that there exists
some constant $\delta>0$ such that for $|t|<\delta$,

$M_{N}(t):=E[e^{tN}]<\infty.$

We note that the moment-generating function $M_{N}(t)$ is well-defined within
the interval $\mathcal{I}_{N}=(-\delta,\delta)$ with $0\leq\delta\leq\infty$
and the endpoints may be included. In particular, we have the following useful
identity for the $k$th order derivative of the moment generating function of
$N$:

$M_{N}^{(k)}(\gamma):=\frac{\partial^{k}M_{N}(t)}{\partial
t^{k}}\bigg{|}_{t=\gamma}=E[N^{k}e^{\gamma N}]$ (1)

for $\gamma\in\mathcal{I}_{N}$. A common and popular candidate for $N$ is the
Poisson distribution. In addition, several parametric count distributions have
been discussed to capture the excess zero observations and over-dispersion
observed in frequency data. Without additional effort, our approach can
utilize these distributions as the frequency part. We list some count
distributions and their moment generating functions in Table 6 and 7.

Consider the aggregate loss random variable during the given time period

$S=\sum_{i=1}^{N}Y_{i}$

where $Y_{i}$ is the size of the individual claim(or severity). Define $S=0$
when $N=0$ by convention. We further define the average claim size(or average
severity) as

$\overline{Y}=\left\\{\begin{array}[]{ll}&\frac{1}{N}\sum_{j=1}^{N}Y_{j}\quad\mbox{if
}N>0\\\ &0\quad\mbox{if }N=0,\\\ \end{array}\right.$ (2)

and then we have $S=N\overline{Y}$. Traditionally, many studies assumed that
$N$ is independent of $Y_{i}$’s to exploit the favorable fact that the mean of
the aggregate loss can be computed by the product of $E[N]$ and $E[S]$.
However, the model without taking account of the dependence effect between
frequency and severity is inevitably faced with model risk. Indeed, various
papers have revealed that the average claim size is positively(or negatively)
related to the number of claims in non-life insurances. For example, see Shi
et al. [2015] and Garrido et al. [2016]. NeurFS incorporates the dependence
structure by describing the average severity in terms of covariates as well as
the frequency.

We suppose that $Y_{i}$’s are mutually independent conditional on $N$. Notice
that the assumption does not imply $Y_{i}$ and $\overline{Y}$ are independent
of $N$. While the flexibility of our approach allows us to use a wide spectrum
of distributions for the severity, we mainly consider the exponential family
of distributions throughout this paper due to its generality and popularity in
actuarial science. If the distribution of $Y_{i}$ belongs to the exponential
family of distributions, we write $Y_{i}\sim EDF(\theta,\phi)$ and the
probability density is defined by

$f_{Y_{i}}(y;\theta,\phi)=c(y,\phi)\exp\bigg{\\{}\frac{y\theta-a(\theta)}{\phi}\bigg{\\}}$

where $\theta\in\mathbb{R}$ is the canonical parameter and
$\phi\in\mathbb{R}^{+}$ is the dispersion parameter. The normalizing function
$c(\cdot)$ and the cumulant function $a(\cdot)$ are determined from the actual
probability function. Expectation and variance of the random variable $Y_{i}$
are expressed as

$E[Y_{i}]=a^{\prime}(\theta)\quad$

and

$Var[Y_{i}]=\phi a^{\prime\prime}(\theta)=\phi V(E[Y_{i}])$

where $a^{\prime}(\theta)$ and $a^{\prime\prime}(\theta)$ are the first and
second derivatives of $a(\theta)$ with respect to $\theta$, respectively, and
$V$ is called the variance function.

One useful property of the exponential family of distributions is that they
are closed under the convolution operation. That is, the average claim size
$\overline{Y}$ also belongs to the exponential family of distributions with a
different dispersion parameter. Specifically, if $Y_{i}\sim EDF(\theta,\phi)$,
conditional on $N$, we then have $\overline{Y}|N\sim EDF(\theta,\phi/N)$. This
well-known property can be easily confirmed by using the moment generating
function of the exponential family of distributions. We refer to Frees [2014].
As a result, one can identify the average severity based on the knowledge
about the individual severity distribution that is well documented in the
literature. For example, $V[\overline{Y}|N]$ is given by

$Var[\overline{Y}|N]=\frac{\phi}{N}V(E[\overline{Y}|N])$ (3)

where $V(\cdot)$ is the variance function of $Y_{i}$.

###  2.2 Dependent and nonlinear model

This subsection develops NeurFS and compute the mean and variance of the
aggregate loss random variable. Let
$\textbf{x}_{i}=(x_{i1},\cdots,x_{ip})\in\mathbb{R}^{p}$ be the set of
covariates that explains the $i$-th observation, $N_{i}$ be the number of
claims and $\overline{Y}_{i}$ be the average severity. We model the mean
function of random variable $N_{i}$ in terms of a regression function
$F(\textbf{x}):\mathbb{R}^{p}\rightarrow\mathbb{R}$ with the log link :

$\ln\lambda_{i}:=\ln\lambda(\mathbf{x}_{i},t_{i})=\ln\mathbb{E}[N_{i}|\textbf{x}_{i},t_{i}]=\ln
t_{i}+F(\textbf{x}_{i})$ (4)

or equivalently, $\lambda_{i}=t_{i}\exp{\\{F(\textbf{x}_{i})\\}}$ where
$t_{i}$ is the exposure variable.

Next, the mean function of random variable $\overline{Y}_{i}$ is specified by
a regression function $S(\textbf{x}):\mathbb{R}^{p}\rightarrow\mathbb{R}$ and
$N_{i}$

$\ln\mu_{i}:=\ln\mu(\textbf{x}_{i},N_{i};\gamma)=\ln\mathbb{E}[\overline{Y}|\mathbf{x}_{i},N_{i}]=S(\textbf{x}_{i})+\gamma
N_{i}$ (5)

or equivalently,

$\mu_{i}=\exp{\\{S(\textbf{x}_{i})+\gamma N_{i}\\}}$

where $\gamma\in\mathcal{I}_{N}$ is a parameter that controls the dependency
of the frequency and severity.111 The range of the dependence parameter
$\gamma$ is dependent on the frequency distribution of interest It is worth
noting that we do not impose any linearity or other restrictions on $F(\cdot)$
and $S(\cdot)$ to incorporate the non-linearity of covariates. Instead,
$F(\cdot)$ and $S(\cdot)$ are represented by neural networks that can
approximate any complex functions thanks to the universal approximation
theorem. Henceforth, we call the triplet of
$(F(\textbf{x}_{i}),S(\textbf{x}_{i}),\gamma)$ as the neural frequency-
severity model.

NeurFS is fully analytically tractable in the sense that once we have
estimated $F(\cdot)$ and $S(\cdot)$, we can obtain explicit closed form
solutions for the mean and variance of the aggregate loss of $\textbf{x}_{i}$,
which are important quantities for analyzing frequency-severity data. In
particular, these statistics are essential for managing insurance portfolios.
Given $\textbf{x}_{i}$, the mean of the aggregate loss $S_{i}$ is given by

$\displaystyle\mathbb{E}\big{[}S_{i}|\textbf{x}_{i}\big{]}=\mathbb{E}\big{[}N_{i}\overline{Y}_{i}|\textbf{x}_{i}\big{]}$
$\displaystyle=$
$\displaystyle\mathbb{E}\bigg{[}N_{i}\mathbb{E}[\overline{Y}_{i}|N_{i},\textbf{x}_{i}]\big{|}\textbf{x}_{i}\bigg{]}$
$\displaystyle=$
$\displaystyle\mathbb{E}\bigg{[}N_{i}e^{S(\textbf{x}_{i})+\gamma
N_{i}}\big{|}\textbf{x}_{i}\bigg{]}$ $\displaystyle=$ $\displaystyle
e^{S(\textbf{x}_{i})}\mathbb{E}\bigg{[}N_{i}e^{\gamma
N_{i}}\big{|}\textbf{x}_{i}\bigg{]}$ $\displaystyle=$ $\displaystyle
e^{S(\textbf{x}_{i})}M_{N}^{(1)}(\gamma|\textbf{x}_{i}).$

where $M_{N}^{(1)}(\gamma|\textbf{x})$ is defined in (1). Note that
$M_{N}^{(1)}(\gamma|\textbf{x})$ is always well-defined when
$\gamma\in\mathcal{I}_{N}$.

The variance of the aggregate loss of $\textbf{x}_{i}$ can be also derived
analytically. From the law of total variance and (3), one calculates that

$\displaystyle Var\big{[}S_{i}\big{|}\textbf{x}_{i}\big{]}$ $\displaystyle=$
$\displaystyle\mathbb{E}\bigg{[}Var[S_{i}|N_{i},\textbf{x}_{i}]\bigg{|}\textbf{x}_{i}\bigg{]}+Var\bigg{[}\mathbb{E}[S_{i}|N_{i},\textbf{x}_{i}]\bigg{|}\textbf{x}_{i}\bigg{]}$
$\displaystyle=$
$\displaystyle\mathbb{E}\bigg{[}N_{i}^{2}Var[\overline{Y_{i}}|N_{i},\textbf{x}_{i}]\bigg{|}\textbf{x}_{i}\bigg{]}+Var\bigg{[}N_{i}\mu_{i}\bigg{|}\textbf{x}_{i}\bigg{]}$
$\displaystyle=$ $\displaystyle\mathbb{E}\bigg{[}N_{i}\phi
V(\mu_{i})\big{|}\textbf{x}_{i}\bigg{]}+\mathbb{E}\bigg{[}N_{i}^{2}\mu_{i}^{2}\big{|}\textbf{x}_{i}\bigg{]}-\mathbb{E}\bigg{[}N_{i}\mu_{i}\big{|}\textbf{x}_{i}\bigg{]}^{2}$
$\displaystyle=$
$\displaystyle\phi\mathbb{E}\bigg{[}N_{i}V(\mu_{i})\big{|}\textbf{x}_{i}\bigg{]}+e^{2S(\textbf{x}_{i})}\bigg{[}\mathbb{E}\big{[}N_{i}^{2}e^{2\gamma
N_{i}}\big{|}\textbf{x}_{i}\big{]}-\mathbb{E}\big{[}N_{i}e^{\gamma
N_{i}}|\textbf{x}_{i}\big{]}^{2}\bigg{]}$ $\displaystyle=$
$\displaystyle\phi\mathbb{E}\bigg{[}N_{i}V(\mu_{i})\big{|}\textbf{x}_{i}\bigg{]}+e^{2S(\textbf{x}_{i})}\bigg{[}M_{N}^{(2)}(2\gamma|\textbf{x}_{i})-M_{N}^{(1)}(\gamma|\textbf{x}_{i})^{2}\bigg{]}$

where $M_{N}^{(2)}(2\gamma|\textbf{x})$ is defined in (1), which is also well-
defined for $2\gamma\in\mathcal{I}_{N}$.

The following theorem summarizes the mean and variance of the aggregate loss
for $\textbf{x}_{i}$ within NeurFS.

######  Theorem 2.1.

Given the neural frequency-severity model,
$(F(\textbf{x}_{i}),S(\textbf{x}_{i}),\gamma$), the mean and variance of the
aggregate loss are given by

$E\big{[}S_{i}\big{|}\textbf{x}_{i}\big{]}=e^{S(\textbf{x}_{i})}M_{N}^{(1)}(\gamma|\textbf{x}_{i}),$

and

$Var\big{[}S_{i}\big{|}\textbf{x}_{i}\big{]}=\phi
E\bigg{[}N_{i}V(\mu_{i})\bigg{|}\textbf{x}_{i}\bigg{]}+e^{2S(\textbf{x}_{i})}\bigg{[}M_{N}^{(2)}(2\gamma|\textbf{x}_{i})-M_{N}^{(1)}(\gamma|\textbf{x}_{i})^{2}\bigg{]}$

for $2\gamma\in\mathcal{I}_{N}$ where $M_{N}(\cdot|\textbf{x})$ is the moment-
generating function of the random variable $N$, and
$M_{N}^{(1)}(\cdot|\textbf{x})$ and $M_{N}^{(2)}(\cdot|\textbf{x})$ are the
first and second-order derivative of $M_{N}(\cdot|\textbf{x})$ defined in (1).
$\phi$ and $V$ are the dispersion parameter and variance function of the
individual claim loss $Y_{1}$, respectively.

When $\gamma=0$, the NeurFS model reduces to a nonlinear prediction model
under the independent assumption between frequency and severity distributions.
$M_{N}$, $M_{N}^{(1)}$, $M_{N}^{(2)}$ and $V$ for popular distributions in
actuarial science can be found in Table 6 and 7.

######  Example 2.1.

Suppose that frequency and severity follow the zero-inflated Poisson
distribution and Gamma distribution, respectively. That is, $N_{i}\sim
ZIP(\lambda_{i},\pi)$ and $Y_{i}\sim Gamma(\mu_{i},\nu)$. Then, we obtain

$E\big{[}S_{i}\big{|}\textbf{x}_{i}\big{]}=(1-\pi)\lambda_{i}\exp\\{\lambda_{i}(e^{\gamma}-1)+S(\textbf{x}_{i})+\gamma\\}$

and

$\displaystyle
Var\big{[}S_{i}\big{|}\textbf{x}_{i}\big{]}=(1-\pi)\lambda_{i}e^{2(S(\textbf{x}_{i})+\gamma)}\bigg{\\{}(1+\lambda_{i}e^{2\gamma}+\frac{1}{\nu})e^{\lambda_{i}(e^{2\gamma}-1)}-e^{2\lambda_{i}(e^{\gamma}-1)}\bigg{\\}}.$
(6)

Moreover, we remark that it is straightforward to consider $\gamma$, $\alpha$
and $\phi$ are also dependent on $\textbf{x}_{i}$. However, it is assumed that
these values are constants throughout our analysis to keep our model being
parsimonious.

##  3 Neural network and estimation

This section discusses an estimation procedure for
$(F(\textbf{x}),S(\textbf{x}),\gamma)$ and associated distribution parameters.
Before proceed, we present some background on neural networks. Simply
speaking, the (feed-forward) neural network is constructed by the composition
of affine transformations and non-linear functions, so called the activation
function. One of the key reasons for the success of neural networks lies on
the universal approximation theorem which states that they can approximate any
function provided that the number of layers and neurons is sufficient. In a
statistical model context, the expressive power of neural networks enables the
proposed model to capture complex non-linearity in inputs.

###  3.1 Definitions and notations

Over the two past decades, various types of neural network architecture have
been developed to deal with different data features and tasks arising in
science and engineering. In this paper, we focus on the feed-forward neural
network(also known as the fully connected network), which is one of the
fundamental architectures.

Let $L\in\mathbb{N}$ and
$[N_{0},N_{1},N_{2},\cdots,N_{L},N_{L+1}]\in\mathbb{N}^{L+2}$ be the number of
hidden layers and the number of neurons in each layer, respectively. Here,
$N_{0}$ and $N_{L+1}$ represent the dimension of input and output layers,
respectively and $N_{1},\cdots,N_{L}$ are the dimensions of hidden layers. In
this paper, we will set $N_{0}=p$ and $N_{L+1}=1$ where $p$ is the number of
covariates. Let $\sigma:\mathbb{R}\rightarrow\mathbb{R}$ denote the activation
function and define

$\displaystyle W^{l}:\mathbb{R}^{N_{l-1}}$ $\displaystyle\longrightarrow$
$\displaystyle\mathbb{R}^{N_{l}}$ x $\displaystyle\longmapsto$ $\displaystyle
A^{l}\textbf{x}+b^{l}$

for $l=1,\cdots,L+1$ where
$A^{l}\in\mathbb{R}^{\mathbb{N}_{l}\times\mathbb{N}_{l-1}}$ is the weight
matrix between the $(l-1)$-th layer and $l$-th layer and
$b^{l}\in\mathbb{R}^{N_{l+1}}$ is the bias term. The parameters of the neural
network, the set of trainable parameters of all weight matrices and bias terms
in the network, are denoted by $\textbf{w}\in\Theta\subseteq\mathbb{R}^{q}$
where $\Theta$ is a parameter space of the neural network and
$q=\sum_{i=1}^{L+1}(N_{i-1}N_{i}+N_{i})$. Then, the neural network
$G(\cdot;\cdot):\mathbb{R}^{p}\times\mathbb{R}^{q}\rightarrow\mathbb{R}$ is
defined by

$G(\textbf{x};\textbf{w})=\bigg{(}W_{L+1}\circ F_{L}\circ\cdots\circ
F_{1}\bigg{)}(\textbf{x})$ (7)

with $F_{l}:=\sigma\circ W^{l}$ where the activation function is applied
element-wise. In particular, the case when $L\geq 2$ is called “deep” neural
network.

Theoretical properties of neural networks as function approximators have been
extensively studied in the literature. See Hornik et al. [1989], Hornik [1991]
and Barron [1993]. To make our paper self-contained, we record the following
theorem from Hornik et al. [1989].

######  Theorem 3.1.

(Universal approximation theorem, Hornik et al. [1989]) Suppose the activation
function is bounded and nonconstant then, for any finite measure $\mu$,
multilayer feedforward networks can approximate any function in $L^{p}(\mu)$
(the space of all functions on $\mathbb{R}^{n}$ such that
$\int_{\mathbb{R}^{n}}|f(x)|^{p}d\mu(x)<\infty$) arbitrarily well, provided
that sufficiently many hidden layers are available.

Moreover, if the activation function is continuous, bounded and nonconstant
then, for arbitrary compact subsets $X\subseteq\mathbb{R}^{n}$, multilayer
feedforward networks can approximate any continuous function on $X$ arbitrary
well with respect to uniform norm, provided that sufficiently many hidden
layers are available.

######  Remark 3.1.

The sigmoid and tanh activation functions satisfy the condition in Theorem
3.1. On the other hand, in practice, it is common to use piecewise-linear
functions such as Rectified Linear Units(ReLUs) and its variants as an
activation function. We refer to Pan and Vivek [2016] for the expressiveness
of feed-forward neural networks with ReLu activations.

###  3.2 Estimation of $F(\cdot)$, $S(\cdot)$ and model parameters

Suppose that we construct a dataset from $m$ independent observations:

$\\{\textbf{z}_{i}\\}_{i=1}^{m}:=\\{{\textbf{x}}_{i},n_{i},t_{i},\overline{y}_{i}\\}_{i=1}^{m}$

where $\textbf{x}_{i}$ is the feature set of covariates, $n_{i}$ is the number
of claims, $t_{i}$ is the observed period and $\overline{y}_{i}$ is the
average severity. Let $f_{N}(\cdot;\lambda,\bm{\alpha})$ be the density of the
frequency distribution with the mean function $\lambda$ and a vector of other
parameters $\bm{\alpha}$. For example, $\bm{\alpha}$ is the proportion of
structural zeros for the zero-inflated Poisson distribution and $\bm{\alpha}$
is the empty set for one-parameter distributions like the Poisson
distribution. Also, we denote the conditional density of the average severity
given the number of claims by $f_{\overline{Y}|N}(\cdot;\mu,\phi)$ with the
mean function $\mu$ and the dispersion parameter $\phi$ for the severity
distribution. Given
$\\{{\textbf{x}}_{i},n_{i},t_{i},\overline{y}_{i}\\}_{i=1}^{m}$, then, the
log-likelihood functional of the frequency model is given by

$\displaystyle\ell_{N}(F,\bm{\alpha}|\\{{\textbf{x}}_{i},n_{i},t_{i},\overline{y}_{i}\\}_{i=1}^{m})$
$\displaystyle=$ $\displaystyle\sum_{i=1}^{m}\ln
f_{N}(n_{i}|\lambda_{i},\bm{\alpha})$ $\displaystyle=$
$\displaystyle\sum_{i=1}^{m}\ln
f_{N}(n_{i}|t_{i}\exp\\{F(\textbf{x}_{i};\cdot)\\},\bm{\alpha}).$

where $F(\cdot;\cdot)$ is a neural network for the frequency model, which is
defined in (7). Similarly, the log-likelihood functional of the average
severity model can be written as

$\displaystyle\ell_{\overline{Y}|N}(S,\gamma,\phi|\\{n_{i},\overline{y}_{i},t_{i},{\bf
x}_{i}\\}_{i=1}^{m})$ $\displaystyle=$ $\displaystyle\sum_{i=1}^{m}\ln
f_{\overline{Y}|N}(\overline{y}_{i}|\mu(x_{i},n_{i};\gamma),\phi)$
$\displaystyle=$ $\displaystyle\sum_{i=1}^{m}\ln
f_{\overline{Y}|N}(\overline{y}_{i}|\exp\\{S(\textbf{x}_{i};\cdot)+\gamma
n_{i}\\},\phi).$

where $S(\cdot;\cdot)$ is the second neural network for the severity model.
Then, the mean function for the frequency model is defined as
$F^{*}(\cdot)=F(\cdot;\textbf{w}^{*})$ by minimizing the negative log-
likelihood(NLL) functional:

$\displaystyle(\textbf{w}^{*},\bm{\alpha}^{*})=\operatorname*{arg\,min}_{\textbf{w}\in\Theta_{F},\alpha\in\mathcal{A}}-\sum_{i=1}^{m}\ln
f_{N}(n_{i}|t_{i}\exp{(F(\textbf{x}_{i};\textbf{w}))},\bm{\alpha})$ (8)

where $\Theta_{F}$ is the parameter space of the neural network of the
frequency model and $\mathcal{A}$ is the set of possible values for
$\bm{\alpha}$. Notice that we have used the entire dataset to estimate the
frequency model. On the other hand, the subset of the dataset for which the
number of claims is positive is only used for estimating $S(\cdot)$, $\gamma$
and $\phi$. More specifically, the mean function for the severity model is
defined as $S^{*}(\cdot)=S(\cdot;\textbf{w}^{*})$ where

$\displaystyle(\textbf{w}^{*},\gamma^{*},\phi^{*})=\operatorname*{arg\,min}_{\textbf{w}\in\Theta_{S},2\gamma\in\mathcal{I}_{N},\phi\in\mathbb{R}^{+}}-\sum_{i\in\mathcal{K}}\ln
f_{\overline{Y}|N}(\overline{y}_{i}|\exp\\{S(\textbf{x}_{i};\textbf{w})+\gamma
n_{i}\\},\phi)$ (9)

where $\mathcal{K}=\\{i\in\\{1,\cdots,m\\}\ |\ n_{i}>0\\}$ and $\Theta_{S}$ is
the parameter space of the neural network of the severity model.

We emphasize that our estimating process is more efficient and easier to be
implemented compared to the existing non-linear(but independent) models of
Yang et al. [2018]. In the model of Yang et al. [2018], a prediction function
and other distribution parameters have to be estimated separately. Thus, one
needs to sequentially solve three optimization problems when there are two
distribution parameters. In contrast, the frequency mean estimate
$F^{*}(\cdot)$ and $\alpha^{*}$ are simultaneously obtained by solving the
optimization problem (8) in our learning process. Likewise, the severity mean
estimate $S^{*}(\cdot)$ and two parameters $\gamma^{*}$ and $\phi^{*}$ are
computed at once from the optimization problem (9). This property is extremely
useful when there are many parameters associated with frequency and severity
distributions to be estimated.

In the specification of Example 2.1, the above estimation process is stated as
follows:

$\displaystyle\mbox{Frequency model}:$
$\displaystyle\quad\min_{\textbf{w}\in\Theta_{F},\
\pi\in[0,1)}-\sum_{i=1}^{m}\bigg{\\{}{\bf
1}_{n_{i}=0}\ln\bigg{(}\pi+(1-\pi)\exp{(-t_{i}\exp
F(\textbf{x}_{i};\textbf{w}))}\bigg{)}$ (10) $\displaystyle+$
$\displaystyle{\bf
1}_{n_{i}>0}\bigg{(}\ln(1-\pi)-t_{i}\exp{F(\textbf{x}_{i};\textbf{w})}+n_{i}(\ln
t_{i}+F(\textbf{x}_{i};\textbf{w}))\bigg{)}\bigg{\\}},$
$\displaystyle\mbox{Severity model}:$
$\displaystyle\quad\min_{\textbf{w}\in\Theta_{S},\gamma\in\mathbb{R},\phi\in\mathbb{R}^{+}}-\sum_{i\in\mathcal{K}}\bigg{\\{}-\ln\overline{y}_{i}-\ln\Gamma\bigg{(}\frac{1}{\phi}\bigg{)}$
(11) $\displaystyle+$
$\displaystyle\frac{1}{\phi}\bigg{(}\ln\bigg{(}\frac{\overline{y}_{i}}{\phi}\bigg{)}-(S(\textbf{x}_{i};\textbf{w})+\gamma
n_{i})\bigg{)}-\frac{\overline{y}_{i}}{\phi}\exp{(-S(\textbf{x}_{i};\textbf{w})-\gamma
n_{i})}\bigg{\\}}.$

We give details about training specification and neural network architectures.
For simplicity’s sake, we focus on training the frequency model (10). The
severity model (11) can be similarly treated. Let
$U(\textbf{w},\pi|\\{\textbf{z}\\}_{i=1}^{m})$ denote the objective function
in (10).

To train neural networks, it is common to use stochastic gradient descent(SGD)
algorithm which updates parameters at each iteration as follows: for
$n\in\mathbb{N}$,

$\displaystyle\textbf{w}_{n+1}$ $\displaystyle=$
$\displaystyle\textbf{w}_{n}-\frac{\eta}{b}\sum_{i=1}^{b}\nabla_{\textbf{w}_{n}}U(\textbf{w}_{n},\pi_{n};\textbf{z}_{i})$
$\displaystyle\pi_{n+1}$ $\displaystyle=$
$\displaystyle\pi_{n}-\frac{\eta}{b}\sum_{i=1}^{b}\nabla_{\pi_{n}}U(\textbf{w}_{n},\pi_{n};\textbf{z}_{i})$

where $b$ is mini-batch size and $\textbf{z}_{i}$ is randomly chosen from the
dataset. Along with the great success of nerual networks, several variants of
SGD such as ADAM of Kingma and Ba [2015] and AMSGrad of Reddi et al. [2018]
have been proposed to boost the speed and stability of the convergence. The
network parameter w is initialized according to the rule proposed by He et al.
[2015] and the initial value for $\pi$ is uniformly distributed on $[e^{-2},\
e^{-1}]$. We use the exponential linear unit(ELU) of Clevert et al. [2016]
defined as

$\sigma(x)=\left\\{\begin{array}[]{ll}x&\quad\mbox{if }x>0,\\\
e^{x}-1&\quad\mbox{if }x\geq 0.\end{array}\right.$

for all experiments in this paper.

Distribution parameters are often restricted to some pre-specified domains.
For example, the proportion of structural zeros $\pi$ in the zero-inflated
Poisson distribution has a value between $0$ and $1$ and the dispersion
parameter in the severity distribution has a positive value. To ensure these
parameters to stay within given domains, we apply proper transformations such
as $x\mapsto x^{2}$ for positive parameters and $x\mapsto\frac{1}{1+e^{-x}}$
for $[0,1]$-valued parameters at the last layer. As a result, $\pi\in[0,1]$ is
expressed by $\frac{1}{1+e^{-\overline{\pi}}}$ with some trainable parameter
$\overline{\pi}$. The final networks for the frequency model and severity
model are visualized in Figure 1 and 2.

Figure 1: Neural network architecture for the frequency model in Example 2.1.
The boxes colored in red represent input tensors. Figure 2: Neural network
architecture for the severity model in Example 2.1. The boxes colored in red
represent input tensors.

##  4 Simulation study

This subsection examines the performance of NeurFS by estimating the mean
functions of frequency, severity and aggregate loss, and associated parameters
from a synthetic dataset. To demonstrate the importance of taking into account
non-linearity and dependence in covariates, we compare the proposed algorithm
to two existing models in the literature; 1. linear and independent model(GLM)
2. linear but dependent model(dGLM) of Garrido et al. [2016].

Let $\textbf{x}_{i}=(x_{i1},x_{i2})\in\mathbb{R}^{2}$ be a vector of
covariates for the $i$-th observation. For simplicity, we set $t_{i}=1$ for
all observations. Define $F(\textbf{x}_{i})=(x_{i1}-0.5)^{2}+(x_{i2}-0.5)^{2}$
and $S(\textbf{x}_{i})=x_{i1}^{2}+x_{i2}^{2}$. Moreover, suppose that the
number of claims $N$ follows the zero-inflated Poisson distribution with the
mean function

$\ln\lambda_{i}=F(\textbf{x}_{i})=(x_{i1}-0.5)^{2}+(x_{i2}-0.5)^{2}.$ (12)

and $\pi=0.2$. Also, assume that the average severity $\overline{Y}_{i}|N_{i}$
is generated from the Gamma distribution with the mean function

$\ln\mu_{i}=S(\textbf{x}_{i})+\gamma N_{i}=x_{i1}^{2}+x_{i2}^{2}+\gamma
N_{i},$ (13)

and the dispersion parameter $\phi=1$, and $\gamma=0.5$. Then, from Theorem
2.1 and Example 2.1, we can analytically compute the true mean of the
aggregate loss $S_{i}$

$\displaystyle E\big{[}S_{i}\big{|}\textbf{x}_{i}\big{]}$ $\displaystyle=$
$\displaystyle\exp\\{\lambda_{i}(e^{\gamma}-1)+F(\textbf{x}_{i})+S(\textbf{x}_{i})+\gamma\\}.$

We generate 40,000 samples

$\\{\textbf{x}_{i},n_{i},\overline{y}_{i}\\}_{i=1}^{40,000}$

with $x_{ij}\sim Uni(\\{0,\frac{1}{10},\frac{2}{10},\cdots,1\\})$ for
$i=1,\cdots,40,000$ and $j=1,2.$, i.e., $P(x_{ij}=\frac{k}{10})=1/11$ for
$k=0,1,\cdots,10$.

The dataset exhibits both the nonlinearity of covariates and the dependency
between frequency and severity. Thus, it is able to identify the impact of
model misspecification stemmed from the lack of non-linearity and dependence
by comparing the three different models.

Under GLM, the estimation errors come from both the linear assumption in
covariates and the independent assumption between frequency and severity. It
is expected that dGLM works well when there is a linear relationship among
variables and the average severity is linearly dependent on the frequency.
However, dGLM is still exposed to the presence of non-linearity in covariates
like GLM. In contrast, our model remains valid even when the non-linearity and
the dependency are present in the data.

#### Frequency model

To estimate the mean function $F(\cdot)$, we use a neural network with $2$
hidden layers and $25$ neurons in each layer. The trained mean function
$F^{*}(\cdot)$ is obtained by solving (10) with AMSGrad of Reddi et al. [2018]
for $50$ epochs. For AMSGrad, the initial learning rate $\eta$ is set to
$0.01$ and decays by a constant factor $0.9$ every $5$ epochs. The mini-batch
size $b$ is 512.

Figure 3 shows that the true mean function $\exp(F(\cdot))$ and the estimated
functions $\exp(F^{*}(\cdot))$ obtained from three different models. It is
apparent that the NeurFS model replicates the true mean function accurately,
but GLM and dGLM generate completely different mean functions. Note that GLM
and dGLM produce the equivalent mean function for the frequency model. To
quantify the accuracy of models, we compute mean absolute error(MAE) and root
mean squared error(RMSE) for $\exp(F(\cdot))$, that measures $\ell_{1}$ and
$\ell_{2}$ norm of the prediction error on the grid
$\\{(x_{1},x_{2}):x_{1},x_{2}\in\\{0,\frac{1}{10},\frac{2}{10},\cdots,1\\}\\}$,
respectively. Table 1 displays MAE and RMSE of NeurFS, GLM and dGLM, showing
that the proposed method significantly performs better than other two models
in terms of given measures.

Figure 3: Fitted curves for the mean function $\exp(F)$ of the frequency model defined in (12) | GLM | dGLM | NeurFS
---|---|---|---
Frequency model
MAE | 0.2456 | 0.2456 | $\bm{0.1513}$
RMSE | 0.2920 | 0.2920 | $\bm{0.1894}$
$\pi$ | 0.2591(29.55%) | 0.2591(29.55%) | $\bm{0.2007(0.35\%)}$
Severity model
MAE | 4.1135 | 0.2065 | $\bm{0.0485}$
RMSE | 4.4662 | 0.2867 | $\bm{0.0625}$
$\phi$ | 2.1888(118.88%) | 1.0300(3%) | $\bm{0.9915(0.85\%)}$
$\gamma$ | N/A | 0.5190(3.8%) | $\bm{0.5022(0.44\%)}$
Aggregate loss
MAE | 3.3289 | 2.7832 | $\bm{0.6797}$
RMSE | 5.1093 | 4.5908 | $\bm{1.0706}$
Table 1: Results of GLM, dGLM and the NeurFS for the frequency model, severity
model and aggregate loss. The true value of $\pi$ is $0.2$. The true value of
$\phi$ is $1$ and the true value of $\gamma$ is 0.5. NA means “not available”.
The percentage of relative error is shown in parenthesis.

Figure 4 plots how $\pi$ is estimated in NuerFS as training progresses. It is
observed that $\pi$ starting from about $0.26$ tends to converge to the true
value $0.2$ as the number of epochs increases. On the other hands, GLM and
dGLM considerably overestimate $\pi$ to compensate the non-linearity,
resulting in large prediction error.

Figure 4: Estimation of $\pi$ in the NeurFS model with respect to the number
of epochs. $x$-axis represents the number of epochs and $y$-axis represents
the estimated $\pi$ at each epoch.

#### Severity model

Let us focus on the severity model defined in (13). We train a neural network
with $2$ hidden layers and $25$ neurons in each layer to estimate
$S^{*}(\cdot)$. See Figure 2 for the whole neural network structure. The
optimization problem (11) for the severity model is solved using AMSGrad for
$50$ epochs with the same setting done in the frequency model except for the
initial learning rate $0.02$.

Figure 5 plots the true function, $\exp(S(\cdot))$, and the fitted curves,
$\exp(S^{*}(\cdot))$, under the three different models. The NeurFS model
generates the highly accurate estimated function, but GLM(and dGLM)
overestimates(and underestimates) $S(\cdot)$ when $x_{1}$ and $x_{2}$ are
close to $1$. Table 1 summarizes simulation results and their accuracy
measures, showing that our model is superior to GLM and dGLM. The proposed
model achieves the lowest MAE and RMSE and recovers accurate distribution
parameters, $\phi$ and $\gamma$, with $0.85\%$ and $0.44\%$ relative errors,
respectively.

Figure 5: Fitted curves for the mean function $\exp(S)$ of the frequency model
defined in (13)

We also plot the mean function of the average severity, $\mu(\cdot,k)$, with
$k=1,2,3$ and the estimated curves to investigate the dependence structure in
Figure 6. It is observed that the estimated mean functions of the average
severity in our model and dGLM are positively correlated to $N$. However, the
mean function in GLM does not change with respect to $N$ due to the
independence assumption between frequency and severity.

Figure 6: Fitted curves for $\mu(\cdot,n)$ with $n=1,2,3$.

#### Aggregate Loss

Predicting the aggregate loss is primary purpose of frequency-severity models.
The performance of our model in predicting the aggregate loss is compared with
GLM and dGLM in Table 1. Our model significantly reduces prediction error for
the aggregate loss. These remarkable results demonstrate that the NeurFS model
is effectively able to extract complex relations in response and explanatory
variables and dependence structure between frequency and severity
distributions. Furthermore, we find that the impact of model misspecification
can cause substantial prediction errors and cannot be negligible when we deal
with nonlinear and dependent frequency-severity data.

##  5 Application to insurance claims

This section is devoted to applying our method to real insurance claims data
and to interpreting the results of the model. This empirical study will show
the feasibility and effectiveness of the NeurFS model in practice.

###  5.1 Data description

We consider the auto-insurance claim data from “freMTPL2freq” and
“freMTPL2sev” in the R package “CASdatasets”, which is also used in Su and Bai
[2020]. After deleting missing data, we have 668,897 records and each record
includes 9 explanatory variables $\textbf{x}_{i}\in\mathbb{R}^{9}$ as well as
the number of claims $n_{i}$, the average severity $\overline{y}_{i}$ and
exposure $t_{i}$. See Table 2 for the description of explanatory variables,
which is available at “http://cas.uqam.ca/pub/R/web/CASdatasets-manual.pdf”.

Variable | Type | Description
---|---|---
VehPower | C | The power of the car (6 classes)
VehAge | C |  The vehicle age in years. ((0,1], (1,4], (4,10], (10, $\infty$))
DrivAge | C | The driver age in years
([18,21], (21,25], (25,35],(35,45], (45,55],(55,70], (70,$\infty$))
VehBrand | C | The car brand (14 classes)
VehGas | C | The car gas (diesel or regular)
Region | C | The policy region in France (22 classes)
Area | C | The density value of the city community where
the car driver lives in (6 classes)
Density | N | the density of inhabitants in the city the driver of the car lives in
BonusMalus | N | Bonus/malus, between 50 and 350
: $<100$ means bonus, $>100$ means malus in France.
Table 2: Explanatory variables and their description.

The dataset is split into $70\%$ samples as training set to train models and
$30\%$ samples as test set to evaluate the trained model. For neural networks,
one-hot encoding is used for converting categorical variables and max-min
normalization is used to scale numeric variables. In addition, the intercept
is included for GLM and dGLM.

###  5.2 Result analysis

We use a fully connected feed forward network with 2 hidden layers and 100
neurons on each layers for both frequency and severity models. i.e, $L=3$,
$N_{0}=9$, $N_{1}=N_{2}=100$ and $N_{3}=1$ (See the notations in subsection
3.1). Furthermore, we employ dropout (Srivastava et al. [2015]) and batch
normalization (Ioffe and Szegedy [2015]) to prevent our models from
overfitting and boost the training speed. Other hyperparameters including the
number of epochs, learning rate, learning decay schedule and so on can be
found in Table 3.

| Frequency model | Severity model
---|---|---
The number of epochs | 100 | 100
Optimizer | AMSGrad | AMSGrad
Learning rate | 0.001 | 0.005
Batch size | 256 | 128
Dropout rate | 0.25 | 0
Batch Normalization | True | False
Activation function | ELU | ELU
Early stopping - patience | 5 | 5
Early stopping - learning decay factor | 0.5 | 0.5
Table 3: Hyperparameters for neural networks

In order to compare the accuracy of models, we measure the mean absolute error
in predicting the frequency and the average severity as follows:

$MAE=\frac{1}{N_{t}}\sum_{i=1}^{N_{t}}|n_{i}-\widehat{n}_{i}|\quad\mbox{for
frequency model}$

and

$MAE=\frac{1}{N_{t}}\sum_{i=1}^{N_{t}}|\overline{y}_{i}-\widehat{\overline{y}}_{i}|\quad\mbox{for
severity model}$

where $N_{t}$ is the number of test set, $n_{i}$ is the observed number of
claims, $\widehat{n}_{i}$ is the model’s prediction of $n_{i}$,
$\overline{y}_{i}$ is the observed average severity and
$\widehat{\overline{y}}_{i}$ is the model’s prediction of $\overline{y}_{i}$.
That is, for the NeurFS model, $\widehat{n}_{i}$ and
$\widehat{\overline{y}}_{i}$ are written as
$\exp\\{\lambda(\textbf{x}_{i},t_{i})\\}$ and
$\exp\\{\mu(\textbf{x}_{i},n_{i};\gamma)\\}$, respectively. In addition, we
report the root mean squared error with respect to the actual observations.
Table 4 shows that NeurFS provides the lowest performance measures for
predicting both frequency and severity, confirming that our model outperforms
GLM and dGLM.

| GLM | dGLM | Neural FS
---|---|---|---
Frequency
MAE | 0.0720 | 0.0720 | $\bm{0.0702}$
RMSE | 0.2047 | 0.2047 | $\bm{0.1998}$
Severity model
MAE | 2023 | 1566 | $\bm{1309}$
RMSE | 2919 | 2510 | $\bm{2267}$
Table 4: MAE and RMSE for the prediction errors of GLM, dGLM and the NeurFS
model in frequency and average severity.

Predicting the aggregate loss is particularly important for an insurer because
the insurer that charges a mispriced premium is prone to adverse selection in
a competitive market (Cohen and Siegelman [2010] and Frees et al. [2014]).
However, when it comes to assessing models in predicting the aggregate loss, a
special care needs to be taken because the distribution of the aggregate loss
is not only zero-inflated, but also heavily skewed. The obstacles make the MAE
or RMSE to be not appropriate for evaluating the accuracy of the predicted
aggregate loss. Alternatively, this paper adopts the ordered Lorenz curve and
the associated Gini index, proposed by Frees et al. [2011b], as a statistical
measure for examining predictions of the aggregate loss distribution.

Define the relative premium $R_{i}$ as $\frac{P_{i}}{B_{i}}$ where $B_{i}$ is
the predicted aggregate loss of a base model and $P_{i}$ is the predicted
aggregate loss of a competing model. Then, the ordered premium distribution,
taking values between $[0,1]$, is given by

$\widehat{F}_{P}(s)=\frac{\sum_{i=1}^{N_{t}}P_{i}{\bf I}(R_{i}\leq
s)}{\sum_{i=1}^{N_{t}}P_{i}}$

and the ordered loss distribution, taking values between $[0,1]$, is

$\widehat{F}_{P}(s)=\frac{\sum_{i=1}^{N_{t}}L_{i}{\bf I}(R_{i}\leq
s)}{\sum_{i=1}^{N_{t}}L_{i}}$

where $\bf I(\cdot)$ is the indicator function. The graph of
$(\widehat{F}_{P}(s),\widehat{F}_{P}(s))$ is called the ordered Lorenz curve
and the $45$-degree line is known as the line of equality. The Gini index is
computed by twice the area between the ordered Lorenz curve and the line of
equality. A larger Gini index implies that the competing model is superior to
the base model in identifying risk groups. We compute the Gini indices by
successively selecting the base model and the competing model among GLM, dGLM
and NFS.

Figure 7 shows the ordered Lorenz curves of all possible pairs of models. It
is clear that the area between the ordered Lorenz curve and the line of
equality is larger when the competing model is NeurFS. On the other hand, if
NeurFS is chosen as the base model, the areas between the ordered Lorenz curve
and the line of equality are close to zero for GLM and dGLM as the competing
model, meaning that GLM and dGLM does not improve the model performance
against NeurFS. We record the corresponding Gini indices in Table 5. These
results suggest that our model that takes account of nonlinearity and
dependency provides effective risk segmentation, and outranks the existing
methods.

| Alternative score
---|---
Base score | GLM | dGLM | Neural FS
GLM | - | -0.1037 | $\bm{0.1787}$
dGLM | 0.1557 | - | $\bm{0.2207}$
Neural FS | $\bm{0.0240}$ | $0.0066$ | -
Table 5: Gini indices for the auto insurance claim data

Figure 7: Ordered Lorenz curves for insurance claims data with GLM, dGLM and
NeurFS models.

###  5.3 Model interpretation

We introduce an excellent tool for explaining outputs of neural networks based
on Shapely value. The Shapley value proposed by Shapley [1953] is one of the
most important object in cooperative game theory to explain marginal
contributions of players. Given a set of players $\mathcal{N}$, a cooperative
game is characterized by a finite set $\mathcal{S}\subseteq\mathcal{N}$ and a
characteristic function $v:\mathcal{P}(\mathcal{S})\rightarrow\mathbb{R}$ with
$v(\emptyset)=0$ where $\mathcal{P}(\mathcal{S})$ is the set of all possible
subset of $\mathcal{S}$. Simply speaking, $v$ maps to the sum of payoffs that
can be obtained from given players. Then, the Shapely value is given by

$\phi_{i}=\sum_{\mathcal{S}\subseteq\mathcal{N}\backslash\\{i\\}}\frac{|\mathcal{S}|!(|\mathcal{N}|-|\mathcal{S}|-1)!}{|\mathcal{N}|!}(v(\mathcal{S}\cup\\{i\\})-v(\mathcal{S})).$

Lundberg and Lee [2017] proposed Shapley additive explanation(SHAP) to
interpret a prediction of highly complex models such as neural networks,
Gradient Boost, XGBoost, by combining local explanation models of Riberiro et
al. [2016] and Shapley values. In their seminal work, it is shown that SHAP
provides the unique additive feature importance values that satisfy three
desirable properties: local accuracy, missingness and consistency. Also, it is
proven that SHAP values are more consistent with human intuition than other
methods such as the dependence plot, the feature importance and saliency map
in the literature. We refer the interested readers to Lundberg and Lee [2017]
and Adebayo et al. [2018].

In the SHAP, each player is considered as each input variable in the model and
$v(\cdot)$ is the output of the model of interest. SHAP specifies the
explanation model $g$ defined as

$g(x^{\prime})=\phi_{0}+\sum_{i=1}^{|N|}\phi_{i}$ (14)

where $\phi_{i}\in\mathbb{R}$ is the importance of a input variable $i$ and
$\phi_{0}$, so called the base value, is the prediction value without any
features. That is, SHAP value $\phi_{i}$ represents the contribution of each
feature to the prediction. We build two SHAP models for $F^{*}(\cdot)$ and
$S^{*}(\cdot)$.

First of all, we investigate the global importance of each feature, which is
computed by the mean absolute value of that feature over given samples. Figure
8 shows the global feature importance of the 9 explanatory variables for the
frequency and severity models. One can observe that significant variables of
the frequency and severity models are completely different. For example,
BonusMalus is the most influential variable for the frequency model. It is
quite natural and intuitive since BonusMalus, also known as the no-claim
discount, is directly determined by the driver’s claim history. More
importantly, it is indirectly verified that a driver with a low BonusMalus is
likely to avoid a car accident due to the driver’s skill or defensive driving,
which giving the justification of a BonusMalus system. Next, it is well
recognized that the age of the driver is an important variable in the
frequency model. We can see that the impact of Driver’s age is second most
influential for explaining the predicted frequency. LogDensity, Area and
Region are the most three relevant variables for the severity model. This
result is consistent with a general consensus that the risk of fatal injury in
rural areas is less than that in urban areas (Castro-Nuno and Arevalo-Quijada
[2018]). On the other hand, BonusMalus is the least important factor. This is
because severity is a quantity conditional on the occurrence of the accident.

Figure 8: Global feature importance plots for frequency (Left) and severity
(Right).

We now perform an individual risk analysis by looking SHAP values for each
observation’s prediction. Imagine that one sample is randomly chosen. We then
compute SHAP values for the frequency model. The model’s output is
approximately $0.08$. Figure 9 illustrates how SHAP values for each feature
contribute to the model’s output. The base value $\phi_{0}$ defined in (14) is
$0.03946$, which is the predicted frequency without any information on
features of the observation. To reach the model’s prediction $0.08$, SHAP
values of all features are added. In this example, the fact that the driver
lives in Area6 increases the frequency rate. In contrast, the predicted value
is negatively associated with BonusMalus of the driver.

Figure 9: Individual risk profile for the frequency model

##  6 Conclusion

NeurFS offers several advantages over the existing methods in the literature.
We provide explicit formulas for the mean and variance of the aggregate loss
when the model takes account of the nonlinearity and the dependency. A wide
range of distributions studied in the literature can be utilized for marginal
distributions without additional effort. We also propose an efficient
estimation process with detailed neural network architectures. To interpret
the results of NeurFS, we apply the SHAP which is a more consistent and
accurate tool for explaining machine learning models. Numerical experiments
including the artificial example and the real dataset confirm the
effectiveness and accuracy of our model. It is also worth noting that
application of NeurFS is not limited to the insurance claim prediction.

This scalable but transparent model can be an important alternative to the
traditional GLM in insurance business where the large volume of data is
available, but strict regulations are imposed. In particular, our research
output will contribute to maintaining the competitiveness of an insurer by
predicting insurance claims accurately, and avoiding adverse selection.
Lastly, the proposed model can be applied to analyzing operations losses,
health care expenditures or credit assessment data where a joint distribution
of frequency and severity is involved.

Frequency distributions
---
Models | $P(N=n)$ | $\bm{\alpha}$ | $E[N]$ | $Var[N]$
Binomial$(m,\lambda)$ | $\binom{m}{n}(\frac{\lambda}{m})^{n}(1-\frac{\lambda}{m})^{n-k}$ | $m$ | $\lambda$ | $\lambda(1-\frac{\lambda}{m})$
Poisson$(\lambda)$ | $\frac{e^{-\lambda}\lambda^{n}}{n!}$ | $\varnothing$ | $\lambda$ | $\lambda$
ZIP$(\pi,\lambda)$ | $\begin{array}[]{ll}&\pi+(1-\pi)e^{-\lambda},\quad n=0\\\ &(1-\pi)\frac{e^{-\lambda}\lambda^{n}}{n!},\quad n=1,2,\cdots.\end{array}$ | $\\{\pi\\}$ | $(1-\pi)\lambda$ | $(1-\pi)(\lambda+\pi\lambda^{2})$
Severity distributions
Models | $f_{\overline{Y}|N}(y)$ | $\phi$ | $E[Y]$ | $V(\mu)$
Gamma$(\mu,\nu)$ | $\frac{1}{y\Gamma(\nu)}\big{(}\frac{y\nu}{\mu}\big{)}^{\nu}e^{-\frac{\nu}{\mu}y}$ | $\frac{1}{\nu}$ | $\mu$ | $\mu^{2}$
Inverse Gaussian$(\mu,\sigma^{2})$ | $\frac{1}{\sqrt{2\mu y^{3}\sigma^{2}}}\exp\big{(}-\frac{1}{2y}\frac{(y-\mu)^{2}}{\mu^{2}\sigma^{2}}\big{)}$ | $\sigma^{2}$ | $\mu$ | $\mu^{3}$
Normal$(\mu,\sigma^{2})$ | $\frac{1}{\sqrt{2\pi}\sigma}\exp\big{(}-\frac{1}{2\sigma^{2}}(y-\mu)^{2}\big{)}$ | $\sigma^{2}$ | $\mu$ | $1$

Table 6: Frequency Severity distributions

Frequency distributions
---
Models | $M_{N}(t)$ | $M_{N}^{(1)}(t)$ | $M_{N}^{(2)}(t)$
Binomial | $\bigg{(}1-\frac{\lambda}{m}+\frac{\lambda}{m}e^{t}\bigg{)}^{m}$ | $\lambda\bigg{(}1-\frac{\lambda}{m}+\frac{\lambda}{m}e^{t}\bigg{)}^{m-1}$ | $\frac{\lambda^{2}(m-1)}{m}\bigg{(}1-\frac{\lambda}{m}+\frac{\lambda}{m}e^{t}\bigg{)}^{m-2}$
Poisson | $e^{\lambda(e^{t}-1)}$ | $\lambda e^{\lambda(e^{t}-1)+t}$ | $\lambda(\lambda e^{t}+1)e^{\lambda(e^{t}-1)+t}$
ZIP | $\pi+(1-\pi)e^{\lambda(e^{t}-1)}$ | $(1-\pi)\lambda e^{\lambda(e^{t}-1)+t}$ | $(1-\pi)\lambda(\lambda e^{t}+1)e^{\lambda(e^{t}-1)+t}$

Table 7: MGF for frequency distributions

## References

  * Adebayo et al. [2018] J. Adebayo, J. Gilmer, M. Muelly, I. Goodfellow, M. Hardt, and B. Kim.  Sanity checks for saliency maps.  _Working paper_ , 2018.
  * Barron [1993] A. Barron.  Universal approximation bounds for superpositions of a sigmoidal function.  _IEEE Transactions on Information theory_ , 39:930–945, 1993.
  * Brechmann et al. [2014] E. Brechmann, C. Czado, and S. Paterlini.  Flexible dependence modeling of operational risk losses and its impact on total capital requirements.  _Journal of Banking & Finance_, 40:271–285, 2014.
  * Castro-Nuno and Arevalo-Quijada [2018] M. Castro-Nuno and M. Teresa Arevalo-Quijada.  Assessing urban road safety through multidimensional indexes: Application of multicriteria decision making analysis to rank the spanish provinces.  _Transport policy_ , 68, 2018.
  * Chavez-Demoulin et al. [2006] V. Chavez-Demoulin, P. Embrechts, and J. Nešlehová.  Quantitative models for operational risk: extremes, dependence and aggregation.  _Journal of Banking & Finance_, 30(10), 2006.
  * Clevert et al. [2016] D-A. Clevert, T. Unterthiner, and S. Hochreiter.  Fast and accurate deep network learning by exponential linear units (elus).  _International Conference on Learning Representations(ICLR)_ , 2016\.
  * Cohen and Siegelman [2010] A. Cohen and P. Siegelman.  Testing for adverse selection in insurance markets.  _Journal of Risk and insurance_ , 77(1), 2010.
  * Czado et al. [2012] C. Czado, R. Kastenmeier, E. C. Brechmann, and A. Min.  A mixed copula model for insurance claims and claim sizes.  _Scandinavian Actuarial Journal_ , 4, 2012.
  * Dahen and Dionne [2010] H. Dahen and G. Dionne.  Scaling models for the severity and frequency of external operational loss data.  _Journal of Banking & Finance_, 34(7), 2010.
  * Frees [2014] E. W. Frees.  _Ch6. Frequency and severity models in Predictie Modeling Applications in Actuarial Science_.  Cambridge University Press, 2014.
  * Frees et al. [2011a] E. W. Frees, J. Gao, and M. A. Rosenberg.  Predicting the frequency and amount of health care expenditures.  _North American Actuarial Journal_ , 15(3):377–392, 2011a.
  * Frees et al. [2011b] E. W. Frees, G. Meyers, and A. D. Cummings.  Summarizing insurance scores using a gini index.  _Journal of the American Statistical Association_ , 106(495):1085–1098, 2011b.
  * Frees et al. [2014] E. W. Frees, G. Meyers, and A. D. Cummings.  Insurance ratemaking and a gini index.  _Journal of Risk and Insurance_ , 81(2):335–366, 2014.
  * Garrido et al. [2016] J. Garrido, C. Genest, and J. Schulz.  Generalized linear models for dependent frequency and severity of insurance claims.  _Insurance: Mathematics and Economics_ , 70:205–215, 2016\.
  * Gschlö$\beta$l and Czado [2007] S. Gschlö$\beta$l and C. Czado.  Spatial modelling of claim frequency and claim size in non-life insurance.  _Scandinavian Actuarial Journal_ , 3:202–225, 2007.
  * Guelman [2012] L. Guelman.  Graient boosting trees for auto insurance loss cost modeling and prediction.  _Expert Systems with Applications_ , 39:3659–3667, 2012\.
  * He et al. [2015] K. He, X. Zhang, S. Ren, and J. Sun.  Delving deep into rectifiers: Surpassing human-level performance on imagenet classification.  _Proceedings of the IEEE international conference on computer vision(ICCV)_ , 2015.
  * Hornik [1991] K. Hornik.  Approximation capabilities of multilayer feedforward networks.  _Neural Networks_ , 4:251–257, 1991.
  * Hornik et al. [1989] K. Hornik, M. Stinchcombe, and H. White.  Multilayer feedforward networks are universal approximators.  _Neural Networks_ , 2:359–366, 1989.
  * Huang and Meng [2019] Y. Huang and S. Meng.  Automobile insurance classification ratemaking based on telematics driving data.  _Decision Support Systems_ , 127, 2019.
  * Ioffe and Szegedy [2015] S. Ioffe and C. Szegedy.  Batch normalization: Accelerating deep network training by reducing internal covariate shift.  _International conference on machine learning_ , pages 448–456, 2015\.
  * Kelly and Nielson [2006] M. Kelly and N. Nielson.  Age as a variable in insurance pricing and risk classification.  _The Geneva Papers on Risk and Insurance-Issues and Practice_ , 31(2):212–232, 2006.
  * Kingma and Ba [2015] D. Kingma and J. Ba.  Adam: A method for stochastic optimization.  _International Conference on Learning Representations(ICLR)_ , 2015\.
  * Kramer et al. [2013] N. Kramer, E. Brechmann, D. Silversrini, and C. Czado.  Total loss estimation using copula-based regression models.  _Insurance: Mathematics and Economics_ , 53:829–839, 2013\.
  * Lundberg and Lee [2017] S. Lundberg and S. Lee.  A unified approach to interpreting model predictions.  _Advances in neural information processing systems(NeurIPS)_ , 2017\.
  * Pan and Vivek [2016] X. Pan and S. Vivek.  Expressiveness of rectifier networks.  _International Conference on Machine Learning(ICML)_ , pages 2427–2435, 2016.
  * Reddi et al. [2018] S. Reddi, S. Kale, and S. Kumar.  On the convergence of ADAM and beyond.  _International Conference on Learning Representations(ICLR)_ , 2018\.
  * Riberiro et al. [2016] M. T. Riberiro, S. Singh, and C. Guestrin.  Why should i trust you?: Explaining the predictions of any classifier.  _22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining_ , pages 1135–1144, 2016.
  * Shapley [1953] L. Shapley.  A value for n-person games.  _Contributions to the Theory of Games_ , 2:307–317, 1953\.
  * Shi et al. [2015] P. Shi, X. Feng, and A. Ivantsova.  Dependent frequency-severity modeling of insurance claims.  _Insurance: Mathematics and Economics_ , 64:417–428, 2015\.
  * Srivastava et al. [2015] N. Srivastava, G. Hinton, A. Krizhevsky, I. Sutskever, and R. Salakhutdinov.  Dropout: a simple way to prevent neural networks from overfitting.  _The journal of machine learning research_ , 15:1929–1958, 2015.
  * Su and Bai [2020] X. Su and M. Bai.  Stochastic gradient boosting frequency-severity model of insurance claims.  _PloS one_ , 15, 2020.
  * Yang et al. [2018] Y. Yang, W. Qian, and H. Zou.  Insurance premium prediction via gradient tree-boosted Tweedie compound Poisson models.  _Journal of Business & Economic Statistics_, 36(3):456–470, 2018.
  * Zhou et al. [2020] H. Zhou, W. Qian, and Y. Yang.  Tweedie gradient boosting for extremely unbalanced zero-inflated data.  _Communications in Statistics-Simulation and Computation_ , online:1–23, 2020.