Bivariate Mixed Poisson and Normal Generalised Linear Models with Sarmanov Dependence—An Application to Model Claim Frequency and Optimal Transformed Average Severity

Abstract

The aim of this paper is to introduce dependence between the claim frequency and the average severity of a policyholder or of an insurance portfolio using a bivariate Sarmanov distribution, that allows to join variables of different types and with different distributions, thus being a good candidate for modeling the dependence between the two previously mentioned random variables. To model the claim frequency, a generalized linear model based on a mixed Poisson distribution -like for example, the Negative Binomial (NB), usually works. However, finding a distribution for the claim severity is not that easy. In practice, the Lognormal distribution fits well in many cases. Since the natural logarithm of a Lognormal variable is Normal distributed, this relation is generalised using the Box-Cox transformation to model the average claim severity. Therefore, we propose a bivariate Sarmanov model having as marginals a Negative Binomial and a Normal Generalized Linear Models (GLMs), also depending on the parameters of the Box-Cox transformation. We apply this model to the analysis of the frequency-severity bivariate distribution associated to a pay-as-you-drive motor insurance portfolio with explanatory telematic variables.

1. Introduction

Calculating premiums is a fundamental task for an insurance company. To this purpose, a simple procedure consists of considering the aggregate claims as the product of the random variable (r.v.) number of claims and of the r.v. average cost of these claims, then of fitting appropriate distributions to these two random variables; if, moreover, the premium is evaluated for a given policyholder, some of its characteristics are often included in the calculation as covariates in Generalized Linear Models (GLMs) used for both variables. Then the pure premium is obtained as the product of the means of the number of claims and of the average claim cost, procedure that implies assuming that the two r.v.s are uncorrelated; however, there are evidences in the literature showing how, in practice, the correlation between both variables is not zero (see References [1,2,3,4] for some illustrations using GLM). Therefore, in this paper, we shall assume that the two just mentioned r.v.s are dependent and jointly distributed using a bivariate Sarmanov distribution with GLM marginals.

Because it provides a flexible structure in joining different types of marginals, the bivariate Sarmanov distribution has recently found its place in actuarial studies like—modeling continuous claim costs [5], modeling discrete claim frequencies [1,6], modeling dependence between discrete frequency and continuous claims [7], ruin probability calculation [8,9] or modeling telematic variables [10].

The traditional approach to modeling count data by means of Poisson regression failed due to over or under dispersion of the data. Therefore, in practice, count data models based on mixed Poisson distributions are used for the number of claims (or frequency) (see Reference [11]. As an alternative, generalized Poisson regression models have also been considered in the literature (see Reference [12]). Specifically, for this variable, in this paper a GLM based on the Negative Binomial distribution will be used.

To model the claim cost r.v. (severity), Gamma and Lognormal distributions are the most common choices (see Reference [7] for comparing both cost distributions in a real motor data set). However, it may be the case that neither of these two distributions fits the data well enough. It may even be difficult to find a known distribution for the cost of claims. To overcome this problem, in this paper it is assumed that by applying a Box-Cox transformation to the claim cost r.v., the transformed data follow a Normal distribution. In the insurance context, Harrington [13] proposed the Box-Cox transformation to generalize the log-linear model in order to calculate pure premiums (see also Reference [14]).

The Box-Cox transformation [15] provides a family of power transformations that is useful in linear regression when the normality assumption is violated, its aim being to obtain a transformed r.v. that is Normal distributed. Furthermore, the Gamma and Lognormal distributions, which have been traditionally used to adjust the distribution of claims costs, are particular cases of this family; therefore, if the original r.v. follows one of these two distributions, there is an optimal transformation so that the transformed variable becomes Normal. This idea is generalized for some unknown distribution for which we can find an optimal transformed r.v. that follows a Normal distribution. Using a Box-Cox transformation of the distribution of the variable in the original scale can be considered as a generalization of the Lognormal distribution, allowing a longer and heavier right tail as ${\lambda }_1$ decreases, and shifting the mode as ${\lambda }_2$ increases.

To conclude, in this paper it is considered a bivariate Sarmanov distribution that allows us to join a Normal GLM distribution for the transformed average claim cost r.v. (see Reference [16] for Sarmanov with Normal marginal distributions) with a Negative Binomial GLM distribution for the number of claims (see Reference [6] for Sarmanov with Negative Binomial GLM marginals). This model can be used to obtain the distribution of the total cost of claims based on the collective model, for a policyholder with specific characteristics. The bivariate Sarmanov distribution allows to fit a non-linear dependence between frequency and severity, and to analyze if the riskier profiles implies larger dependency. Furthermore, the marginal severity distribution based on a Box-Cox transformation allows a longer and heavier right tail than classical models like Gamma and Lognormal.

As numerical illustration, a database containing information on a portfolio of policyholders that have contracted an auto insurance which involved the installation of a GPS/inertial device in their vehicle is used. This device provides a source of information (telematic variables) to motor insurers that complements the variables that have traditionally been used in a-priori pricing in auto insurance. The telematic variables provide an innovative way to calculate car insurance pricing (see References [17,18,19,20,21,22] where a similar database is considered and the telematic variables are used to predict accident rate). Furthermore, for each policyholder in the database, the number of claims and their mean cost are known; these variables are right skewed and have excess of zeros, characteristics that are common in insurance auto databases.

Traditionally, tariffication in auto insurance has been based on the total number of claims of the policyholders. Given the large number of zeros, the models commonly used for fitting the claim frequency variable are the NB and the Zero Inflated Poisson [22,23]. Alternatively, some works analyze the number of claims using multivariate models, that is, considering that there are different types of claims—with civil liability, with personal injury, and so forth [1,6,24]. More recently, tariffication based on the collective model has been proposed, that is, the premium is deduced from the distribution of the total cost variable, which in turn is deduced from the bivariate distribution of frequency and severity claims. Two alternative approaches have been proposed—on the one hand, the bivariate model is deduced from the frequency distribution and from the severity distribution conditioned to frequency defined for the average claim size distribution using the number of claims as covariate [2,4]. On the other hand, Czado et al. [25] proposed bivariate models based on copulae (see also References [26] for a more general approach). Furthermore, the bivariate Sarmanov model has also been analyzed in the collective model framework [7].

The rest of the paper is organized as follows—in Section 2.1, we describe the proposed Sarmanov model relating the r.v. number of claims and the r.v. average claim cost; its properties and particular cases for the proposed marginal distributions are presented in Section 2.2; moreover, some characteristics of the inference from margin (IFM) estimation procedure are pointed out in Section 2.3. The application on a real database using telematic variables to predict the total cost of claims of a given insured is discussed in Section 3. Finally, Section 4 concludes. This paper ends with Appendix A and Appendix B containing the proofs and additional results.

2. The Dependence Model and Its Properties

Let $N_i$, $i=1,...,m$, be the counting r.v. representing the number of claims of individual i and let $Y_i$ be the r.v. representing the average claim cost per insured. Then, the resulting aggregate claims of individual i can be represented as

$$S_i=N_i{Y}_i,$$

(1)

where the usual assumption under which this model is considered is that $Y_i$ is independent of $N_i$. We assume that $Y_i>0$ if and only if $N_i>0$; otherwise, both variables take the value zero. Therefore, we define the r.v. ${\tilde{Y}}_i$ representing $Y_i>0$ and, based on it, the expectation and variance of $S_i$ (which are useful to obtain a-priori ratemaking in insurance) are given by:

$$\begin{array}{ccc}\hfill \mathbb{E}{S}_i& =& \mathbb{E}{N}_i\mathbb{E}{\tilde{Y}}_i,\hfill \\ \hfill Var{S}_i& =& \mathbb{E}\left[{\tilde{Y}}_i^2\right]Var{N}_i+{\left(\mathbb{E}{N}_i\right)}^{2}Var{\tilde{Y}}_i.\hfill \end{array}$$

However, in practice, the independence assumption is not always true, in which case the moments of $S_i$ must be deduced from a bivariate distribution associated to the random vector $(N_i,Y_i)$, distribution that takes into account the dependence between the two variables.

Moreover, we assume that for each i, the distributions of both random variables in vector $(N_i,Y_i)$ depend on a set of k quantitative or binary covariates, which are represented by the vector ${\mathbf{X}}_i={\left(X_{i1},...,X_{ik}\right)}^{\prime }$. For simplicity, the covariates are assumed to be common, but they can also be different between the two variables. In practice, we specify the relation as a GLM and define the linear predictor ${\mathbf{X}}_{\mathbf{i}}^{\prime }{\beta }^j$, where ${\beta }^j=\left({\beta }_1^j,...,{\beta }_k^j\right)$, $j\in \{N,Y\}$, are vectors of parameters to be estimated; note that throughout the paper, such a $j\in \{N,Y\}$ should be interpreted as an upper index and not as a power.

2.1. The Model Relating the Counting and Average Claim Cost r.v.s

This dependence model is defined in two parts: the first part is the probability mass function (pmf) of $N_i=0$ and the second part is the conditional bivariate density function that joins the pmf $p_N$ of the discrete r.v. $N_i$ with the probability density function (pdf) $f_{\tilde{Y}}$ of the continuous r.v. ${\tilde{Y}}_i$ corresponding to $Y_i>0$. The model is:

$$f_{Y_i,N_i}\left(y,n|{\mathbf{X}}_{\mathbf{i}}\right)=\left\{\begin{array}{c}{p}_N\left(0|{\mathbf{X}}_i^{\prime }{\beta }^N\right),n=y=0\hfill \\ p_N\left(n|{\mathbf{X}}_i^{\prime }{\beta }^N\right)f_{\tilde{Y}}\left(y|{\mathbf{X}}_i^{\prime }{\beta }^Y\right)\left(1+\omega \psi \left(n|{\mathbf{X}}_i^{\prime }{\beta }^N\right)\varphi \left(y|{\mathbf{X}}_i^{\prime }{\beta }^Y\right)\right),n\ge 1,y>0,\hfill \end{array}\right.$$

(2)

where $\omega$ is the Sarmanov dependence parameter and $\psi (·)$ and $\varphi (·)$ are bounded kernel functions.

In order for (2) to define a proper pdf for all i, we impose the conditions:

$$\begin{array}{ccc}\hfill \sum _{n\ge 1}\psi \left(n|{\mathbf{X}}_i^{\prime }{\beta }^N\right)p_N\left(n|{\mathbf{X}}_i^{\prime }{\beta }^N\right)& =& {\int }_{{\mathbb{R}}^{+}}\varphi \left(y|{\mathbf{X}}_i^{\prime }{\beta }^Y\right)f_{\tilde{Y}}\left(y|{\mathbf{X}}_i^{\prime }{\beta }^Y\right)dy=0,\mathrm{and}\hfill \end{array}$$

(3)

$$\begin{array}{ccc}\hfill 1+\omega \psi \left(n|{\mathbf{X}}_i^{\prime }{\beta }^N\right)\varphi \left(y|{\mathbf{X}}_i^{\prime }{\beta }^Y\right)& \ge & 0,\mathrm{for}\mathrm{all}n\ge 1,y>0.\hfill \end{array}$$

(4)

Then, the kernel functions limits for each i are defined as:

$$\begin{array}{ccc}\hfill m_N\left({\mathbf{X}}_i^{\prime }{\beta }^N\right)& =& \underset{n\ge 1}{inf}\psi \left(n|{\mathbf{X}}_i^{\prime }{\beta }^N\right),M_N\left({\mathbf{X}}_i^{\prime }{\beta }^N\right)=\underset{n\ge 1}{sup}\psi \left(n|{\mathbf{X}}_i^{\prime }{\beta }^N\right),\hfill \\ \hfill m_Y\left({\mathbf{X}}_i^{\prime }{\beta }^Y\right)& =& \underset{y>0}{inf}\varphi \left(y\right|\left({\mathbf{X}}_i^{\prime }{\beta }^Y\right)),M_Y\left({\mathbf{X}}_i^{\prime }{\beta }^Y\right)=\underset{y>0}{sup}\varphi \left(y|{\mathbf{X}}_i^{\prime }{\beta }^Y\right).\hfill \end{array}$$

Taking into account the conditions defined in (3) and (4), limits are also defined for the dependence parameter $\omega$. However, since this parameter does not depend on the linear predictor, new maximum and minimum values are defined as: $m_j^{\star }=\underset{\forall {\mathbf{X}}_i}{max}{m}_j\left({\mathbf{X}}_i^{\prime }{\beta }^j\right)$ and $M_j^{\star }=\underset{\forall {\mathbf{X}}_i}{min}{M}_j\left({\mathbf{X}}_{\mathbf{i}}^{\prime }{\beta }^j\right)$, $j\in \{N,Y\}$, so that the limits of the dependence parameter are:

$$max\left\{-\frac{1}{m_N^{\star }{m}_Y^{\star }},-\frac{1}{M_N^{\star }{M}_Y^{\star }}\right\}\le \omega \le min\left\{-\frac{1}{m_N^{\star }{M}_Y^{\star }},-\frac{1}{M_N^{\star }{m}_Y^{\star }}\right\}.$$

(5)

In the following proposition, the new formulae for the expectation and variance of $S_i$ are presented. To simplify the notation, the functions used in the bivariate density expressed in (2) are rewritten as follows: $p_i\left(n\right)=p_N\left(n|{\mathbf{X}}_i^{\prime }{\beta }^N\right)$, $f_i\left(y\right)=f_{\tilde{Y}}\left(y|{\mathbf{X}}_i^{\prime }{\beta }^Y\right)$, ${\psi }_i\left(n\right)=\psi \left(n|{\mathbf{X}}_i^{\prime }{\beta }^N\right)$ and ${\varphi }_i\left(y\right)=\varphi \left(y|{\mathbf{X}}_i^{\prime }{\beta }^Y\right)$.

Proposition 1.

Under the dependence introduced by the bivariate Sarmanov distribution (2), the expected value and variance of the aggregate claims $S_i$ defined in (1) for individual i are given respectively, by

$$\begin{array}{ccc}\hfill \mathbb{E}{S}_i& =& \mathbb{E}{N}_i\mathbb{E}{\tilde{Y}}_i+\omega \mathbb{E}\left[N_i{\psi }_i\left(N_i\right)\right]\mathbb{E}\left[{\tilde{Y}}_i{\varphi }_i\left({\tilde{Y}}_i\right)\right],\hfill \\ \hfill Var{S}_i& =& \mathbb{E}\left[{\tilde{Y}}_i^2\right]Var{N}_i+{\left(\mathbb{E}{N}_i\right)}^{2}Var{\tilde{Y}}_i-{\omega }^2{\mathbb{E}}^2\left[N_i{\psi }_i\left(N_i\right)\right]{\mathbb{E}}^2\left[{\tilde{Y}}_i{\varphi }_i\left({\tilde{Y}}_i\right)\right]\hfill \\ & & +\omega \left(\mathbb{E}\left[N_i^2{\psi }_i\left(N_i\right)\right]\mathbb{E}\left[{\tilde{Y}}_i^2{\varphi }_i\left({\tilde{Y}}_i\right)\right]-2\mathbb{E}{N}_i\mathbb{E}\left[N_i{\psi }_i\left(N_i\right)\right]\mathbb{E}{\tilde{Y}}_i\mathbb{E}\left[{\tilde{Y}}_i{\varphi }_i\left({\tilde{Y}}_i\right)\right]\right).\hfill \end{array}$$

From Proposition 3 in Bolancé and Vernic [7], the following correlation for each i can be easily deduced:

$$corr\left(Y_i,N_i\right)=\frac{\omega \mathbb{E}\left[N_i{\psi }_i\left(N_i\right)\right]\mathbb{E}\left[{\tilde{Y}}_i{\varphi }_i\left({\tilde{Y}}_i\right)\right]+p_i\left(0\right)\mathbb{E}{N}_i\mathbb{E}{\tilde{Y}}_i}{\sqrt{\left(1-p_i\left(0\right)\right)\left(Var{\tilde{Y}}_i+p_i\left(0\right){\mathbb{E}}^2\left[{\tilde{Y}}_i\right]\right)Var{N}_i}}.$$

(6)

To calculate the moments of $S_i$ and the correlation expressed in (6), we need to define the marginal GLMs and the resulting particular bivariate Sarmanov model with given kernel functions.

We propose to use exponential kernels. For the cost per insured $Y_i$, the kernel function is expressed as ${\varphi }_i\left(y\right)=e^{-\gamma y}-{\mathcal{L}}_{{\tilde{Y}}_i}\left(\gamma \right)$, where ${\mathcal{L}}_{{\tilde{Y}}_i}$ denotes the Laplace transform of ${\tilde{Y}}_i$. For the number of claims with $n\ge 1$ in (2), we let ${\psi }_i\left(n\right)=e^{-\delta n}-k_i,$ and, to find $k_i$, we impose condition (3) as follows

$$\begin{array}{ccc}\hfill \sum _{n\ge 1}{\psi }_i\left(n\right)p_i\left(n\right)& =& \sum _{n\ge 1}\left(e^{-\delta n}-k_i\right)p_i\left(n\right)\hfill \\ & =& \sum _{n\ge 0}{e}^{-\delta n}{p}_i\left(n\right)-p_i\left(0\right)-k_i\left(\sum _{n\ge 0}{p}_i\left(n\right)-p_i\left(0\right)\right)\hfill \\ & =& {\mathcal{L}}_{N_i}\left(\delta \right)-p_i\left(0\right)-k_i\left(1-p\left(0\right)\right)=0.\hfill \end{array}$$

Therefore, $k_i=\frac{{\mathcal{L}}_{N_i}\left(\delta \right)-p_i\left(0\right)}{1-p_i\left(0\right)}$ and ${\psi }_i\left(n\right)=e^{-\delta n}-\frac{{\mathcal{L}}_{N_i}\left(\delta \right)-p_i\left(0\right)}{1-p_i\left(0\right)}$, where ${\mathcal{L}}_{N_i}$ denotes the Laplace transform of $N_i$.

2.2. Marginal Distributions

2.2.1. Counting Distribution

We assume that the distribution of the counting process $N_i$ is the Negative Binomial (NB), where we take ${\mu }_i^N=\mathbb{E}{N}_i=e^{{\mathbf{X}}_i^{\prime }{\beta }^N}$, so that in the GLM specification $ln{\mu }_i^N={\mathbf{X}}_i^{\prime }{\beta }^N$, and the pmf is:

$$p_i\left(n\right)=\frac{\Gamma \left(\alpha +n\right)}{n!\Gamma \left(\alpha \right)}{\left(\frac{\alpha }{\alpha +{\mu }_i^N}\right)}^{\alpha }{\left(\frac{{\mu }_i^N}{\alpha +{\mu }_i^N}\right)}^n,$$

(7)

where $\alpha >0$. The previous model is heteroscedastic and its variance is $Var{N}_i=\frac{{\mu }_i^N\left(\alpha +{\mu }_i^N\right)}{\alpha }$ for each i. Furthermore, the Laplace transform is ${\mathcal{L}}_{N_i}\left(\delta \right)={\left(\frac{\alpha }{\alpha +{\mu }_i^N\left(1-e^{-\delta }\right)}\right)}^{\alpha }$.

To obtain the correlation between the frequency and the severity per policyholder, and the moments of $S_i$ from Proposition 1, we need to calculate $\mathbb{E}\left[N_i^2\right]$, $\mathbb{E}\left[N_i{\psi }_i\left(N_i\right)\right]$ and $\mathbb{E}\left[N_i^2{\psi }_i\left(N_i\right)\right]$. The first one is direct, given that $Var{N}_i=\mathbb{E}\left[N_i^2\right]-{\left(\mathbb{E}{N}_i\right)}^2$, while the other two expectations can be directly deduced from the Proposition 5 of Bolancé and Vernic [7]. The results are:

$$\begin{array}{ccc}\hfill \mathbb{E}\left[N_i^2\right]& =& \frac{{\mu }_i^N\left(\alpha +{\mu }_i^N+\alpha {\mu }_i^N\right)}{\alpha },\hfill \end{array}$$

(8)

$$\begin{array}{ccc}\hfill \mathbb{E}\left[N_i{\psi }_i\left(N_i\right)\right]& =& \frac{{\mu }_i^N{\alpha }^{\alpha +1}}{\left(\alpha +{\mu }_i^N\right){\left[\alpha +{\mu }_i^N\left(1-e^{-\delta }\right)\right]}^{\alpha }}\times \left\{\frac{\alpha +{\mu }_i^N}{\alpha e^{\delta }-{\mu }_i^N\left(1-e^{\delta }\right)}\right.\hfill \\ & -& \left.\frac{{\left(\alpha +{\mu }_i^N\right)}^{\alpha +1}-\left(\alpha +{\mu }_i^N\right){\left(\alpha +{\mu }_i^N\left(1-e^{-\delta }\right)\right)}^{\alpha }}{\alpha \left[{\left(\alpha +{\mu }_i^N\right)}^{\alpha }-{\alpha }^{\alpha }\right]}\right\},\hfill \end{array}$$

(9)

$$\begin{array}{ccc}\hfill \mathbb{E}\left[N_i^2{\psi }_i\left(N_i\right)\right]& =& \frac{{\mu }_i^N{\alpha }^{\alpha +1}}{\left(\alpha +{\mu }_i^N\right){\left[\alpha +{\mu }_i^N\left(1-e^{-\delta }\right)\right]}^{\alpha }}\times \left\{\frac{\alpha {\mu }_i^N\left(\alpha +{\mu }_i^N\right)+{\left(\alpha +{\mu }_i^N\right)}^2{e}^{\delta }}{{\left[\alpha e^{\delta }-{\mu }_i^N\left(1-e^{\delta }\right)\right]}^2}\right.\hfill \\ & -& \left.\left({\mu }_i^N+\alpha \left(1+{\mu }_i^N\right)\right)\frac{{\left(\alpha +{\mu }_i^N\right)}^{\alpha +1}-\left(\alpha +{\mu }_i^N\right){\left(\alpha +{\mu }_i^N\left(1-e^{-\delta }\right)\right)}^{\alpha }}{{\alpha }^2\left[{\left(\alpha +{\mu }_i^N\right)}^{\alpha }-{\alpha }^{\alpha }\right]}\right\}.\hfill \end{array}$$

(10)

2.2.2. Severity Distribution

In what concerns the mean cost per policyholder represented by the r.v. $Y_i$, we assume that its distribution is unknown, but the variable ${\tilde{Y}}_i$ can be normalised using a Box-Cox transformation. The one parameter Box-Cox transformation is given by:

$$T_{\lambda }\left(y\right)=\left\{\begin{array}{c}\frac{y^{\lambda }-1}{\lambda }\mathrm{if}\lambda \ne 0,y>0,\hfill \\ ln\left(y\right)\mathrm{if}\lambda =0,y>0,\hfill \end{array}\right.$$

(11)

while the two parameters Box-Cox transformation is:

$$T_{{\lambda }_1,{\lambda }_2}\left(y\right)=\left\{\begin{array}{c}\frac{{\left(y+{\lambda }_2\right)}^{{\lambda }_1}-1}{{\lambda }_1}\mathrm{if}{\lambda }_1\ne 0,y>-{\lambda }_2,\hfill \\ ln\left(y+{\lambda }_2\right)\mathrm{if}{\lambda }_1=0,y>-{\lambda }_2.\hfill \end{array}\right.$$

(12)

We assume that the two parameters Box-Cox transformation $T_{{\lambda }_1,{\lambda }_2}(·)$ is applied to the r.v. ${\tilde{Y}}_i$ and obtain the truncated normal r.v. $Z_i=T_{{\lambda }_1,{\lambda }_2}\left({\tilde{Y}}_i\right).$ Since the domain of ${\tilde{Y}}_i$ is left bounded (also necessary in order to have a bounded kernel function), when ${\lambda }_1\ge 0,$ the r.v. $Z_i$ is left truncated normal (LTN), $Z_i\sim LTN\left({\mu }_i^Z,{\sigma }^2;a\right)$, while when ${\lambda }_1<0,$ the resulting r.v. is doubly truncated normal (DTN), $Z_i\sim DTN\left({\mu }_i^Z,{\sigma }^2;a,b\right),$ with $a<b\le -{\lambda }_1^{-1}$ as we shall see below.

Then, with $\phi \left(·\right)$ and $\Phi \left(·\right)$ denoting the pdf and, respectively, the cumulative distribution function of the standard normal distribution, the density of ${\tilde{Y}}_i$ is given by

$$f_{{\tilde{Y}}_i}\left(y\right)=\left\{\begin{array}{c}{T}_{{\lambda }_1,{\lambda }_2}^{\prime }\left(y\right)\phi \left(\frac{T_{{\lambda }_1,{\lambda }_2}\left(y\right)-{\mu }_i^Z}{\sigma }\right)\frac{1}{\sigma \overline{\Phi }\left(z_{a,i}\right)},{\lambda }_1\ge 0\hfill \\ T_{{\lambda }_1,{\lambda }_2}^{\prime }\left(y\right)\phi \left(\frac{T_{{\lambda }_1,{\lambda }_2}\left(y\right)-{\mu }_i^Z}{\sigma }\right)\frac{1}{\sigma \left(\Phi \left(z_{b,i}\right)-\Phi \left(z_{a,i}\right)\right)},{\lambda }_1<0,\hfill \end{array}\right.$$

where $\overline{\Phi }\left(z\right)=1-\Phi \left(z\right)$ and $z_{a,i}=\frac{a-{\mu }_i^Z}{\sigma },z_{b,i}=\frac{b-{\mu }_i^Z}{\sigma }$. More precisely:

• If ${\lambda }_1>0,$ the condition $y>-{\lambda }_2$ implies

  $$z=T_{{\lambda }_1,{\lambda }_2}\left(y\right)=\frac{{\left(y+{\lambda }_2\right)}^{{\lambda }_1}-1}{{\lambda }_1}>-\frac{1}{{\lambda }_1}=T_{{\lambda }_1,{\lambda }_2}\left(-{\lambda }_2\right),$$

  hence the left truncation point a must satisfy the condition $a\ge -\frac{1}{{\lambda }_1},$ for which we obtain that the left truncation condition gives

  $$z=T_{{\lambda }_1,{\lambda }_2}\left(y\right)>a\Rightarrow y>{\left(1+a{\lambda }_1\right)}^{\frac{1}{{\lambda }_1}}-{\lambda }_2.$$
  Therefore, the pdf of ${\tilde{Y}}_i$ is

  $$f_{{\tilde{Y}}_i}\left(y\right)=\frac{{\left(y+{\lambda }_2\right)}^{{\lambda }_1-1}}{\sigma \sqrt{2\pi }\overline{\Phi }\left(z_{a,i}\right)}{e}^{-\frac{1}{2{\sigma }^2{\lambda }_1^2}{\left({\left(y+{\lambda }_2\right)}^{{\lambda }_1}-1-{\lambda }_1{\mu }_i^Z\right)}^2},y>{\left(1+a{\lambda }_1\right)}^{\frac{1}{{\lambda }_1}}-{\lambda }_2,a\ge -\frac{1}{{\lambda }_1}.$$

  (13)
• If ${\lambda }_1=0,$ then $y>-{\lambda }_2$ implies

  $$z=T_{{\lambda }_1,{\lambda }_2}\left(y\right)=ln\left(y+{\lambda }_2\right)>T_{{\lambda }_1,{\lambda }_2}\left(-{\lambda }_2\right)=-\infty ,$$

  so that the left truncation point a can be any real value. Then the left truncation condition becomes

  $$z=T_{{\lambda }_1,{\lambda }_2}\left(y\right)>a\Rightarrow y>e^a-{\lambda }_2.$$
  The distribution of ${\tilde{Y}}_i$ becomes the left truncated lognormal $LTLN\left({\mu }_i^Z,{\sigma }^2;a\right)$ having pdf

  $$f_{{\tilde{Y}}_i}\left(y\right)=\frac{1}{\left(y+{\lambda }_2\right)\sigma \sqrt{2\pi }\overline{\Phi }\left(z_{a,i}\right)}{e}^{-\frac{1}{2{\sigma }^2}{\left(ln\left(y+{\lambda }_2\right)-{\mu }_i^Z\right)}^2},y>e^a-{\lambda }_2.$$
• If ${\lambda }_1<0,$ the condition $y>-{\lambda }_2$ implies

  $$z=T_{{\lambda }_1,{\lambda }_2}\left(y\right)=\frac{{\left(y+{\lambda }_2\right)}^{{\lambda }_1}-1}{{\lambda }_1}<-\frac{1}{{\lambda }_1}=T_{{\lambda }_1,{\lambda }_2}\left(-{\lambda }_2\right),$$

  and the left truncation point a must be such that $a<-\frac{1}{{\lambda }_1}.$ Again, the left truncation condition gives

  $$z=T_{{\lambda }_1,{\lambda }_2}\left(y\right)>a\Rightarrow y>{\left(1+a{\lambda }_1\right)}^{\frac{1}{{\lambda }_1}}-{\lambda }_2.$$
  However, note that in this case, when $y\to \infty$ then $z\to -\frac{1}{{\lambda }_1},$ hence the domain of $Z_i$ is also right bounded, yielding the doubly truncated normal distribution $DTN\left({\mu }_i^Z,{\sigma }^2;a,b\right)$ with $b=-{\lambda }_1^{-1}$ (note that we can take $b<$$-{\lambda }_1^{-1}$ as long as $a<b$). We therefore write the pdf of ${\tilde{Y}}_i$ as

  $$f_{{\tilde{Y}}_i}\left(y\right)=\frac{{\left(y+{\lambda }_2\right)}^{{\lambda }_1-1}}{\sigma \sqrt{2\pi }\left(\Phi \left(z_{b,i}\right)-\Phi \left(z_{a,i}\right)\right)}{e}^{-\frac{1}{2{\sigma }^2{\lambda }_1^2}{\left({\left(y+{\lambda }_2\right)}^{{\lambda }_1}-1-{\lambda }_1{\mu }_i^Z\right)}^2},y>{\left(1+a{\lambda }_1\right)}^{\frac{1}{{\lambda }_1}}-{\lambda }_2,a<-\frac{1}{{\lambda }_1}.$$

In order to write the exponential kernel function corresponding to $Z_i,{\varphi }_{Z_i}\left(z\right)=e^{-\gamma z}-{\mathcal{L}}_{Z_i}\left(\gamma \right),$ we note that the Laplace transform of the DTN distribution $DTN\left({\mu }_i^Z,{\sigma }^2;a,b\right)$ is

$${\mathcal{L}}_{Z_i}\left(\gamma \right)=e^{-\gamma {\mu }_i^Z+\frac{{\gamma }^2{\sigma }^2}{2}}\frac{\Phi \left(\frac{b-{\mu }_i^Z+\gamma {\sigma }^2}{\sigma }\right)-\Phi \left(\frac{a-{\mu }_i^Z+\gamma {\sigma }^2}{\sigma }\right)}{\Phi \left(z_{b,i}\right)-\Phi \left(z_{a,i}\right)}.$$

Taking $b=\infty ,$ we obtain the formula for the LTN distribution.

Therefore, the kernel corresponding to ${\tilde{Y}}_i$ becomes

$${\tilde{\varphi }}_i\left(y\right)=\left\{\begin{array}{c}{e}^{-\gamma \frac{{\left(y+{\lambda }_2\right)}^{{\lambda }_1}-1}{{\lambda }_1}}-e^{-\gamma {\mu }_i^Z+\frac{{\gamma }^2{\sigma }^2}{2}}\frac{\overline{\Phi }\left(\frac{a-{\mu }_i^Z+\gamma {\sigma }^2}{\sigma }\right)}{\overline{\Phi }\left(z_{a,i}\right)},{\lambda }_1>0\hfill \\ \frac{1}{{\left(y+{\lambda }_2\right)}^{\gamma }}-e^{-\gamma {\mu }_i^Z+\frac{{\gamma }^2{\sigma }^2}{2}}\frac{\overline{\Phi }\left(\frac{a-{\mu }_i^Z+\gamma {\sigma }^2}{\sigma }\right)}{\overline{\Phi }\left(z_{a,i}\right)},{\lambda }_1=0\hfill \\ e^{-\gamma \frac{{\left(y+{\lambda }_2\right)}^{{\lambda }_1}-1}{{\lambda }_1}}-e^{-\gamma {\mu }_i^Z+\frac{{\gamma }^2{\sigma }^2}{2}}\frac{\Phi \left(\frac{b-{\mu }_i^Z+\gamma {\sigma }^2}{\sigma }\right)-\Phi \left(\frac{a-{\mu }_i^Z+\gamma {\sigma }^2}{\sigma }\right)}{\Phi \left(z_{b,i}\right)-\Phi \left(z_{a,i}\right)},{\lambda }_1<0\hfill \end{array}\right..$$

As discussed above for the r.v. $N_i$, the following quantities are needed: $\mathbb{E}{\tilde{Y}}_i$, $\mathbb{E}\left[{\tilde{Y}}_i^2\right]$, $\mathbb{E}\left[{\tilde{Y}}_i{\tilde{\varphi }}_i\left({\tilde{Y}}_i\right)\right]$, $\mathbb{E}\left[{\tilde{Y}}_i^2{\tilde{\varphi }}_i\left({\tilde{Y}}_i\right)\right]$, which will be separately calculated for ${\lambda }_1\ne 0$ and for ${\lambda }_1=0$.

Case${\lambda }_1\ne 0.$

We have

$$\mathbb{E}{\tilde{Y}}_i=\left\{\begin{array}{c}\frac{1}{\sigma \sqrt{2\pi }\overline{\Phi }\left(z_{a,i}\right)}{\int }_{{\left(1+a{\lambda }_1\right)}^{\frac{1}{{\lambda }_1}}-{\lambda }_2}^{\infty }y{\left(y+{\lambda }_2\right)}^{{\lambda }_1-1}{e}^{-\frac{1}{2{\sigma }^2{\lambda }_1^2}{\left({\left(y+{\lambda }_2\right)}^{{\lambda }_1}-1-{\lambda }_1{\mu }_i^Z\right)}^2}dy,{\lambda }_1>0\\ \frac{1}{\sigma \sqrt{2\pi }\left(\Phi \left(z_{b,i}\right)-\Phi \left(z_{a,i}\right)\right)}{\int }_{{\left(1+a{\lambda }_1\right)}^{\frac{1}{{\lambda }_1}}-{\lambda }_2}^{\infty }y{\left(y+{\lambda }_2\right)}^{{\lambda }_1-1}{e}^{-\frac{1}{2{\sigma }^2{\lambda }_1^2}{\left({\left(y+{\lambda }_2\right)}^{{\lambda }_1}-1-{\lambda }_1{\mu }_i^Z\right)}^2}dy,{\lambda }_1<0\end{array}\right..$$

We change variable $\frac{{\left(y+{\lambda }_2\right)}^{{\lambda }_1}-1}{{\lambda }_1}-{\mu }_i^Z=t\Rightarrow {\left(y+{\lambda }_2\right)}^{{\lambda }_1-1}dy=dt,$ hence

$$\mathbb{E}{\tilde{Y}}_i=\left\{\begin{array}{c}\frac{1}{\sigma \sqrt{2\pi }\overline{\Phi }\left(z_{a,i}\right)}{\int }_{a-{\mu }_i^Z}^{\infty }\left[{\left({\lambda }_1\left(t+{\mu }_i^Z\right)+1\right)}^{\frac{1}{{\lambda }_1}}-{\lambda }_2\right]e^{-\frac{t^2}{2{\sigma }^2}}dt,{\lambda }_1>0\hfill \\ \frac{1}{\sigma \sqrt{2\pi }\left(\Phi \left(z_{b,i}\right)-\Phi \left(z_{a,i}\right)\right)}{\int }_{a-{\mu }_i^Z}^{\infty }\left[{\left({\lambda }_1\left(t+{\mu }_i^Z\right)+1\right)}^{\frac{1}{{\lambda }_1}}-{\lambda }_2\right]e^{-\frac{t^2}{2{\sigma }^2}}dt,{\lambda }_1<0\hfill \end{array}\right..$$

In general, there is no closed type formula for this integral and for similar integrals associated with $\mathbb{E}\left[{\tilde{Y}}_i^2\right]$, $\mathbb{E}\left[{\tilde{Y}}_i{\tilde{\varphi }}_i\left({\tilde{Y}}_i\right)\right]$, $\mathbb{E}\left[{\tilde{Y}}_i^2{\tilde{\varphi }}_i\left({\tilde{Y}}_i\right)\right]$ and they must be evaluated numerically. However, in the following particular case some recursive formulas can be used.

Particular case: ${\lambda }_1=\frac{1}{m},$ where m is a positive integer. The following notation for $k\in \mathbb{N}$, $k>0$ is introduced:

$$\begin{array}{ccc}\hfill A_{m,k}^i& =& \mathbb{E}\left[{\left(\frac{Z_i}{m}+1\right)}^k\right],\hfill \\ \hfill B_{m,k}^i& =& \mathbb{E}\left[{\left(\frac{Z_i}{m}+1\right)}^k{e}^{-\gamma Z_i}\right].\hfill \end{array}$$

Here the upper index i emphasises the connection with individual i. The proof of the following lemma is very easy and we skip it since it can also be obtained as a particular case of a result in Burkardt [27].

Lemma 1.

Let $L_j\left(\alpha \right)$ be the j-th moment of the standard left truncated normal distribution with left truncation point α, that is,

$$L_j\left(\alpha \right)=\frac{1}{\sqrt{2\pi }\overline{\Phi }\left(\alpha \right)}{\int }_{\alpha }^{\infty }{z}^j{e}^{-\frac{z^2}{2}}dz,j\in \mathbb{N}.$$

Then $L_j\left(\alpha \right)$ can be recursively evaluated as

$$\begin{array}{ccc}\hfill L_0\left(\alpha \right)& =& 1,L_1\left(\alpha \right)=\frac{\phi \left(\alpha \right)}{\overline{\Phi }\left(\alpha \right)},\hfill \\ \hfill L_j\left(\alpha \right)& =& {\alpha }^{j-1}\frac{\phi \left(\alpha \right)}{\overline{\Phi }\left(\alpha \right)}+\left(j-1\right)L_{j-2}\left(\alpha \right),j\ge 2.\hfill \end{array}$$

Lemma 2.

With the above notation, for $Z_i\sim LTN\left({\mu }_i^Z,{\sigma }^2;a\right)$ it holds that

$$\begin{array}{ccc}\hfill \mathbb{E}\left[Z_i^r\right]& =& \sum _{j=0}^r\left(\begin{array}{c}r\\ j\end{array}\right){\sigma }^j{\left({\mu }_i^Z\right)}^{r-j}{L}_j\left(\frac{a-{\mu }_i^Z}{\sigma }\right),r\in \mathbb{N},\hfill \\ \hfill A_{m,k}^i& =& \sum _{r=0}^k\left(\begin{array}{c}k\\ r\end{array}\right)\frac{\mathbb{E}\left[Z_i^r\right]}{m^r},\hfill \\ \hfill B_{m,k}^i& =& e^{-\gamma {\mu }_i^Z+\frac{{\gamma }^2{\sigma }^2}{2}}\frac{\overline{\Phi }\left(\frac{a-{\mu }_i^Z}{\sigma }+\sigma \gamma \right)}{\overline{\Phi }\left(\frac{a-{\mu }_i^Z}{\sigma }\right)}\sum _{j=0}^k\left(\begin{array}{c}k\\ j\end{array}\right){\left(\frac{{\mu }_i^Z-{\sigma }^2\gamma }{m}+1\right)}^{k-j}{\left(\frac{\sigma }{m}\right)}^j{L}_j\left(\frac{a-{\mu }_i^Z}{\sigma }+\sigma \gamma \right).\hfill \end{array}$$

The following proposition provides the needed formulas.

Proposition 2.

Let ${\tilde{Y}}_i$ be the r.v. with pdf (13) and ${\lambda }_1=\frac{1}{m}$. Then

$$\begin{array}{ccc}\hfill \mathbb{E}{\tilde{Y}}_i& =& A_{m,m}^i-{\lambda }_2,\hfill \\ \hfill \mathbb{E}\left[{\tilde{Y}}_i^2\right]& =& A_{m,2m}^i-2{\lambda }_2{A}_{m,m}^i+{\lambda }_2^2,\hfill \\ \hfill \mathbb{E}\left[{\tilde{Y}}_i{\tilde{\varphi }}_i\left({\tilde{Y}}_i\right)\right]& =& B_{m,m}^i-A_{m,m}^i{\mathcal{L}}_{Z_i}\left(\gamma \right),\hfill \\ \hfill \mathbb{E}\left[{\tilde{Y}}_i^2{\tilde{\varphi }}_i\left({\tilde{Y}}_i\right)\right]& =& B_{m,2m}^i-A_{m,2m}^i{\mathcal{L}}_{Z_i}\left(\gamma \right)-2{\lambda }_2\left(B_{m,m}^i-{\mathcal{L}}_{Z_i}\left(\gamma \right)A_{m,m}^i\right).\hfill \end{array}$$

Case${\lambda }_1=0.$

Let now ${\lambda }_1=0.$ As noted before, in this case, ${\tilde{Y}}_i$ follows a left truncated Lognormal distribution $LTLN\left({\mu }_i^Z,{\sigma }^2;a\right)$ and the following result holds without further assumption on the parameter ${\lambda }_1$.

Proposition 3.

For ${\tilde{Y}}_i\sim LTLN\left({\mu }_i^Z,{\sigma }^2;a\right)$ and assuming that $\gamma =1$, it holds that

$$\begin{array}{ccc}\hfill \mathbb{E}{\tilde{Y}}_i& =& e^{{\mu }_i^Z+\frac{{\sigma }^2}{2}}\frac{\overline{\Phi }\left(z_{a,i}-\sigma \right)}{\overline{\Phi }\left(z_{a,i}\right)}-{\lambda }_2,\hfill \\ \hfill \mathbb{E}\left[{\tilde{Y}}_i^2\right]& =& e^{2{\mu }_i^Z+2{\sigma }^2}\frac{\overline{\Phi }\left(z_{a,i}-2\sigma \right)}{\overline{\Phi }\left(z_{a,i}\right)}-2{\lambda }_2{e}^{{\mu }_i^Z+\frac{{\sigma }^2}{2}}\frac{\overline{\Phi }\left(z_{a,i}-\sigma \right)}{{\overline{\Phi }}_i\left(z_{a,i}\right)}+{\lambda }_2^2,\hfill \\ \hfill \mathbb{E}\left[{\tilde{Y}}_i{\tilde{\varphi }}_i\left({\tilde{Y}}_i\right)\right]& =& 1-e^{{\sigma }^2}\frac{\overline{\Phi }\left(z_{a,i}+\sigma \right)\overline{\Phi }\left(z_{a,i}-\sigma \right)}{{\overline{\Phi }}^2\left(z_{a,i}\right)},\hfill \\ \hfill \mathbb{E}\left[{\tilde{Y}}_i^2{\tilde{\varphi }}_i\left({\tilde{Y}}_i\right)\right]& =& e^{{\mu }_i^Z+\frac{{\sigma }^2}{2}}\frac{\overline{\Phi }\left(z_{a,i}-\sigma \right)}{\overline{\Phi }\left(z_{a,i}\right)}-e^{{\mu }_i^Z+\frac{5{\sigma }^2}{2}}\frac{\overline{\Phi }\left(z_{a,i}+\sigma \right)\overline{\Phi }\left(z_{a,i}-2\sigma \right)}{{\overline{\Phi }}^2\left(z_{a,i}\right)}\hfill \\ & & +2{\lambda }_2\left(e^{{\sigma }^2}\frac{\overline{\Phi }\left(z_{a,i}+\sigma \right)\overline{\Phi }\left(z_{a,i}-\sigma \right)}{{\overline{\Phi }}^2\left(z_{a,i}\right)}-1\right),\hfill \end{array}$$

where, as before, $z_{a,i}=\frac{a-{\mu }_i^Z}{\sigma }.$

2.3. Parameter Estimation

Since the distribution associated with the r.v. $Y_i$ is unknown, the parameter estimation is based on the two parameters Box-Cox transformed variable $T_{{\lambda }_1,{\lambda }_2}\left({\tilde{Y}}_i\right)=Z_i$, whose distribution is LTN or DTN. Therefore, we can estimate the bivariate Sarmanov with NB and TN marginal distributions and then the results on the original scale can be deduced as we have shown in Section 2.2. The bivariate Sarmanov distribution with NB and LTN or DTN marginals has pdf:

$$f_{Z_i,N_i}\left(z,n|{\mathbf{X}}_{\mathbf{i}}\right)=\left\{\begin{array}{c}{\left(\frac{\alpha }{\alpha +{\mu }_i^N}\right)}^{\alpha },n=z=0\hfill \\ \frac{\Gamma \left(\alpha +n\right)}{n!\Gamma \left(\alpha \right)}{\left(\frac{\alpha }{\alpha +{\mu }_i^N}\right)}^{\alpha }{\left(\frac{{\mu }_i^N}{\alpha +{\mu }_i^N}\right)}^n\frac{\phi \left(\frac{z-{\mu }_i^Z}{\sigma }\right)}{\sigma \overline{\Phi }\left(z_{a,i}\right)}\left(1+\omega {\psi }_i\left(n\right){\varphi }_{Z_i}\left(z\right)\right),n\ge 1,z>0,{\lambda }_1\ge 0\hfill \\ \frac{\Gamma \left(\alpha +n\right)}{n!\Gamma \left(\alpha \right)}{\left(\frac{\alpha }{\alpha +{\mu }_i^N}\right)}^{\alpha }{\left(\frac{{\mu }_i^N}{\alpha +{\mu }_i^N}\right)}^n\frac{\phi \left(\frac{z-{\mu }_i^Z}{\sigma }\right)}{\sigma \left(\Phi \left(z_{b,i}\right)-\Phi \left(z_{a,i}\right)\right)}\left(1+\omega {\psi }_i\left(n\right){\varphi }_{Z_i}\left(z\right)\right),n\ge 1,z>0,{\lambda }_1<0\hfill \end{array}\right.,$$

(14)

where we recall that ${\mu }_i^N=e^{{\mathbf{X}}_i^{\prime }{\beta }^N}$ and ${\mu }_i^Z={\mathbf{X}}_i^{\prime }{\beta }^Z$. In conclusion, we have to estimate the transformation parameters ${\lambda }_1$ and ${\lambda }_2$, the vectors of parameters ${\beta }^N$ and ${\beta }^Z$ associated with the covariate vector, the parameter $\alpha$ of NB marginal distribution, the parameter $\sigma$ of the LTN or DTN marginal distribution and the dependence parameter $\omega$.

Let $(n_i,y_i)$, $i=1,...,m$, be a sample of observed values of frequency and severity per policyholder and let $z_i=T_{{\lambda }_1,{\lambda }_2}\left(y_i\right)$, $i=1,...,m$, be the transformed severity. The logarithm of the likelihood function $l\left(\Theta \right)$ to be maximised with respect to the vector of parameters which is defined as $\Theta =\left({\beta }^N,{\beta }^Z,\alpha ,\sigma ,{\lambda }_1,{\lambda }_2,\omega \right)$ to be estimated, depends on the value of ${\lambda }_1$ as follows: letting $m_0$ be the number of insured with $n_i=0$ and $m_1$ the number of insured with $n_i\ge 1$, then

• If ${\lambda }_1\ne 0$ and ${\lambda }_2>-min\left(y_1,...,y_m\right)$

  $$\begin{array}{ccc}\hfill l\left(\Theta \right)& =& m\alpha ln\alpha -\alpha \sum _{i=1}^{m}ln\left(\alpha +e^{{\mathbf{X}}_i^{\prime }{\beta }^N}\right)+\sum _{i=1}^{m_1}ln\Gamma \left(\alpha +n_i\right)-\sum _{i=1}^{m_1}ln\left(n_i!\right)\hfill \\ & -& m_{1}ln\Gamma \left(\alpha \right)+\sum _{i=1}^{m_1}{n}_i{\mathbf{X}}_i^{\prime }{\beta }^N-\sum _{i=1}^{m_1}{n}_{i}ln\left(\alpha +e^{{\mathbf{X}}_i^{\prime }{\beta }^N}\right)\hfill \\ & +& \sum _{i=1}^{m_1}ln\left[\phi \left(\frac{\frac{{\left(y_i+{\lambda }_2\right)}^{{\lambda }_1}-1}{{\lambda }_1}-{\mathbf{X}}_i^{\prime }{\beta }^Z}{\sigma }\right)\right]-\sum _{i=1}^{m_1}ln\left[\Phi \left(\frac{b-{\mathbf{X}}_i^{\prime }{\beta }^Z}{\sigma }\right)-\Phi \left(\frac{a-{\mathbf{X}}_i^{\prime }{\beta }^Z}{\sigma }\right)\right]\hfill \\ & -& m_{1}ln\sigma +\sum _{i=1}^{m_1}ln\left[1+\omega {\psi }_i\left(n_i\right){\varphi }_{Z_i}\left(\frac{{\left(y_i+{\lambda }_2\right)}^{{\lambda }_1}-1}{{\lambda }_1}\right)\right],\hfill \end{array}$$

  (15)

  where $b=\infty$ if ${\lambda }_1\ge 0$ and $b\le -{\lambda }_1^{-1}$ if ${\lambda }_1<0.$
• If ${\lambda }_1=0$ and ${\lambda }_2>-min\left(y_1,...,y_m\right)$

  $$\begin{array}{ccc}\hfill l\left(\Theta \right)& =& m\alpha ln\alpha -\alpha \sum _{i=1}^{m}ln\left(\alpha +e^{{\mathbf{X}}_i^{\prime }{\beta }^N}\right)+\sum _{i=1}^{m_1}ln\Gamma \left(\alpha +n_i\right)-\sum _{i=1}^{m_1}ln\left(n_i!\right)\hfill \\ & -& m_{1}ln\Gamma \left(\alpha \right)+\sum _{i=1}^{m_1}{n}_i{\mathbf{X}}_i^{\prime }{\beta }^N-\sum _{i=1}^{m_1}{n}_{i}ln\left(\alpha +e^{{\mathbf{X}}_i^{\prime }{\beta }^N}\right)\hfill \\ & +& \sum _{i=1}^{m_1}ln\left[\phi \left(\frac{ln\left(y_i+{\lambda }_2\right)-{\mathbf{X}}_i^{\prime }{\beta }^Z}{\sigma }\right)\right]-\sum _{i=1}^{m_1}ln\left[\overline{\Phi }\left(\frac{a-{\mathbf{X}}_i^{\prime }{\beta }^Z}{\sigma }\right)\right]\hfill \\ & -& m_{1}ln\sigma +\sum _{i=1}^{m_1}ln\left[1+\omega {\psi }_i\left(n_i\right){\varphi }_{Z_i}\left(ln\left(y_i+{\lambda }_2\right)\right)\right].\hfill \end{array}$$

  (16)

To maximise the log-likelihood functions defined in (15) and (16), a procedure similar to the one described and tested in Bolancé and Vernic [7] can be used, which is divided in two phases: the first phase consists in using the inference from margins (IFM) technique to estimate the marginal parameters, and the second phase in using these results as starting values to obtain estimations for all parameters based on the maximization of the full log-likelihood.

The IFM is a two step estimation method which starts from an initial estimation of the marginal parameters, estimation obtained by the maximum likelihood (ML) method for the parameters of the two GLMs—GLM of the frequency variable with NB distribution, and GLM of the severity variable whose distribution is defined by the transformation parameters and by the parameters of the truncated normal distribution with given truncation values a and b. We will take the left truncation value such that $\Phi \left(\frac{a-{\mathbf{X}}_i^{\prime }{\beta }^Z}{\sigma }\right)$ is almost zero; note that since we consider heterogeneity between policyholders, in practice, a good choice will be $a=min(a_1,...,a_m)$, where $a_i=-3\sigma +{\mathbf{X}}_i^{\prime }{\beta }^Z$ (another simple choice for a would be the minimum of the transformed data) (see References [28] for ML estimation of the univariate model based on Box-Cox transformation). However, note that for (15) we must check that $a\ge -\frac{1}{{\lambda }_1}$ when ${\lambda }_1\ge 0$ and that $a<-\frac{1}{{\lambda }_1}$ if ${\lambda }_1<0$. In practice, for ${\lambda }_1<0$, we used the maximum value for the right truncation point $b=-{\lambda }_1^{-1}$, which does not affect the likelihood value because $\Phi \left(\right(b-mu)/sigma)\approx 1$.

Therefore, for the Sarmanov distribution defined in (14), each iteration of the IFM consists of the following two steps:

Step 1

(iteration r) Given the parameters of the marginal distributions obtained at iteration $r-1$, find the estimation ${\widehat{\omega }}^r$ of the dependence parameter within the interval defined in (5), estimation that maximises the log-likelihood in (15) or (16).

Step 2

Given ${\widehat{\omega }}^r$ obtained in Step 1, find new values for the parameters of the marginals that maximise the log-likelihood function in (15) or (16).

In practice, in the numerical analysis that will be presented in Section 3, two assumptions on the Box-Cox transformation parameters are used. On the one hand, as we have described before, we can include the transformation parameters in the optimization procedure. On the other hand, we a-priori select the values of the transformation parameters, which can be calculated by maximizing the log-likelihood of the normal distribution associated with the transformed variable. The results obtained in both cases are practically the same.

3. Numerical Analysis

In this section, a database corresponding to a car insurance portfolio is analyzed, in which some of the variables have been measured via a telematic system. It is interesting to check if, in the collective risk model context, the telematic variables can replace or complement classical a-priori ratemaking variables, as the age, the driving license age, the power of the car, and so forth. The dependent variables are the frequency and mean severity of claims per policyholder. Information on 25,014 policyholders with a car insurance is available, who have contracted a policy that incorporates a GPS device in the vehicle. The data were provided by a Spanish insurer. The information corresponds to all the drivers in the insurer’s database who had a telematics insurance product in 2010. In general, this type of insurance is mainly chosen by young drivers who value the fact that their car can be located in case of an accident. So, the group of drivers is generally composed of new drivers. Our dataset contains all the available drivers; most of them did not report any claim during 2010, but a few reported at least one accident. Only accidents at fault are considered, and the cost of the accident was also collected. The cost of the accident is generally known once damages, medical expenses and bodily injury compensations are paid. The sum of all these costs was included in the dataset in a variable measuring the cost of claim.

Traditionally, the variables used to model the a-priori premium in auto insurance have been classified in policyholder and auto characteristics; furthermore, if there is some experience with the insureds, that is, they are not new in the company, some tariff variables also could be considered (see Reference [11], for some examples, [chapters 12 and 13]). In practice, these variables are those available for the insurance company. In our case, we have a portfolio with new policyholders in the company and for whom we know their age, their driving license age, and if night parking is used; we also know the gender, but this variable is not included in the analysis because official regulations prevent differentiating premiums based on the gender of the insured. About the auto of the policyholder, it is know that all are cars of private use and the available variables are the age of the car and its power. More recently, thanks to GPS devices, telematic variables are available for the insurance company; in general, the variables that are used in the premium calculation are the total of kilometers (km) driven and three percentages: in urban area, at night and over the speed limit (see Reference [22] for an example with similar portfolio).

Table 1. Definition of variables and descriptive statistics: mean, median, standard deviation (STD), minimum (Min), maximum (Max), skewness (Skew) and kurtosis (Kur). The last row shows the linear correlation between the dependent variables and a confidence interval (CI) at a $95\%$ level.

	Description	Mean	Median	STD	Min	Max	Skew	Kur
N	Number of claims per policyholder	0.106	0.000	0.370	0.000	5.000	4.005	21.598
Y	Cost of claims per policyholder *	1444.175	696.720	4141.511	17.750	130870.360	23.365	678.291
$X_1$	Age of the driver	27.565	27.491	3.094	19.849	36.903	−0.059	2.077
$X_2$	Age of driver License	7.174	6.616	3.053	1.810	15.910	1.402	9.959
$X_3$	Age of the vehicle	8.749	7.775	4.174	1.938	20.468	0.770	2.996
$X_4$	Power of the vehicle	97.225	97.000	27.773	12.000	500.000	−1.309	2.713
$X_5$	Night parking (1 = yes, 0 = no)	0.774	1.000	0.418	0.000	1.000	−1.467	9.955
$X_6$	Total of km. driven in logarithm	8.681	8.774	0.698	0.466	10.820	1.049	4.126
$X_7$	% of km. driven at urban area	25.875	22.922	14.357	0.000	100.000	2.382	10.764
$X_8$	% of km. driven over the speed limit	6.332	3.999	6.827	0.000	70.433	1.830	9.512
$X_9$	% of km. driven at night	6.908	5.148	6.351	0.000	100.000	0.497	2.393
$\rho$	Correlation between N and Y (CI)	0.094 (0.052, 0.135)

* Given $N>0$ equal to 2177 policyholders.

shows the main descriptive statistics of the dependent variables and of the covariates used in the bivariate Sarmanov model defined in (14). It can be noticed that the dependent variables have right skewness. Considering this shape and comparing with alternative count data models like the Poisson or Zero Inflated Poisson, the NB GLM defined in (7) proved to be the best option for our data. Focusing on the mean severity, the distributions that have been classically used are the Gamma or the Lognormal (see Reference [7] for an example). However, in our case, the shape of the severity distribution is unknown and has a heavier tail than the Gamma and the Lognormal; in this case, the Box-Cox transformation allows us to work with the normal distribution for the transformed variable. Regarding the explanatory variables, it can be checked that the portfolio is made up of young policyholders with a symmetric age distribution. The variables with the largest skewness, which in our case is positive, are the % of km driven over the speed limit ($X_8$) and the % of km driven at night ($X_9$); these variables are those directly related with the risk exposure and, therefore, with the frequency and the severity of claims.

At the bottom of Table 1, the linear Pearson correlation coefficient and its confidence interval are shown. The value of this coefficient is low, but significant. Given the right skewness shape of both dependent variables, the linear correlation is not the best dependence measure. Alternatively, Reference [25] proposed fitting copula models with GLM marginals, and, using Gaussian copula, they obtain a dependence parameter equal to 0.13, that is, the dependence between frequency and severity is not very strong, but it affects the insurance premium considerably. The Sarmanov model with GLM marginals proposed in this work is an alternative way to copula for modelling dependence; furthermore, it allows heterogeneity dependence between policyholders with different covariate values.

Focusing on the marginal distributions for the frequency and severity variables, the NB and Gamma GLM models (see Reference [2] for example), respectively, are the most used to model the variables number and cost of claims in auto insurance. Alternatively, in some cases, the Zero Inflated Poisson (ZIP) and Lognormal distributions could improve the fit of NB and Gamma, respectively (see References [7,29], for example). In

Table 2. AIC and BIC for Generalized Linear Models (GLMs) for frequency (N) and severity (Y) variables.

	N			Y
	Poisson	ZIP	NB	Gamma	Lognormal	Box-Cox
AIC	17,372	16,861	16,906	36,009	35,312	25,835
BIC	17,454	17,024	16,996	36,071	35,374	25,898

the AIC (Akaike Information Criterion) and BIC (Bayesian Information Criterion) are compared for alternative univariate GLM models considered for marginals, including the Poisson, the ZIP and the NB for frequency, and the Gamma, Lognormal and Box-Cox transformation based model for severity. These results show that ZIP is better than NB if AIC is used, but, if BIC is used, the best model is the NB. For the severity variable, the best model is the one based on the Box-Cox transformation.

Table 3. Parameter estimates (p-values) and goodness of fit statistics for the bivariate Sarmanov models with GLM marginals for number (N) and mean cost (Y) per policyholder.

	${\mathit{\lambda }}_1=$−0.0745,${\mathit{\lambda }}_2=$13.2501
	Model I		Model II		Model III
	$\mathbf{N}$	$\mathbf{Y}$	$\mathbf{N}$	$\mathbf{Y}$	$\mathbf{N}$	$\mathbf{Y}$
Int.	−8.073 (<0.001)	4.990 (<0.001)	−1.794 (<0.001)	5.171 (<0.001)	−8.277 (<0.001)	5.016 (<0.001)
$X_1$	0.006 ( 0.259)	−0.004 ( 0.260)	−0.007 ( 0.227)	−0.005 ( 0.194)
$X_2$	−0.083 (<0.001)	−0.001 ( 0.480)	−0.081 (<0.001)	0.002 ( 0.488)
$X_3$	0.012 ( 0.013)	0.007 ( 0.021)	0.008 ( 0.082)	0.008 ( 0.015)
$X_4$	0.002 ( 0.037)	0.001 ( 0.071)	0.003 ( 0.001)	0.001 ( 0.019)
$X_5$	−0.025 ( 0.321)	−0.009 ( 0.396)	−0.022 ( 0.339	−0.005 ( 0.446)
$X_6$	0.603 (<0.001)	0.012 ( 0.323)			0.603 (<0.001)	0.010 ( 0.344)
$X_7$	0.022 (<0.001)	0.001 ( 0.200)			0.024 (<0.001)	0.001 ( 0.153)
$X_8$	0.007 ( 0.024)	0.004 ( 0.041)			0.006 ( 0.046)	0.005 ( 0.018)
$X_9$	0.007 ( 0.026)	0.004 ( 0.051)			0.009 ( 0.007)	0.004 ( 0.031)
$\alpha$	0.407 (<0.001)		0.339 (<0.001)		0.381 (<0.001)
$\sigma$		0.650 (<0.001)		0.652 (<0.001)		0.652 (<0.001)
$\omega$	110.236 ( 0.036)		91.508 ( 0.065)		129.375 ( 0.019)
$\widehat{l}\left(\Theta \right)$	−10,555.071		−10,706.315		−10,605.962
AIC	21,156.143		21,442.630		21,237.923
BIC	242,786.895		246,265.501		243,957.371

presents the results of the estimated bivariate Sarmanov model used to fit the joint distribution of the frequency and of the Box-Cox transformation of the mean severity of claims for a given insured with the characteristics represented in a given vector of covariates. The transformation parameters are calculated a-priori by maximizing the normal log-likelihood function of the transformed variable; alternatively, the transformation parameters are included in step 2 of the estimation procedure described in Section 2.3, and it is observed that the results do not improve. A negative value for ${\lambda }_1$ is obtained, indicating that the distribution of the mean cost per policyholder has a longer and heavier tail than the Lognormal distribution. The estimated model with all the explanatory variables that are shown in Table 3 is compared with alternative proposals that can be found in the literature (see References [2,7,25]), to estimate collective risk model distributions (results are shown in Appendix B); for the data used in this paper, these models do not improve our proposed model.

Three Sarmanov GLM models are estimated—Model I incorporates all the explanatory variables, Model II only includes classical a-priori characteristics of the policyholder and car, and, finally, Model III only includes telematic variables. Furthermore, Models I and II are estimated without the variable $X_1$ (Age of the driver) given its high correlation with $X_2$ (Age of driver license), and noticed that the signs and p-values of the parameters of the remaining variables change only slightly. The dependence parameter $\omega$ is positive and indicates a significant dependence between frequency and severity that need to be considered in this collective model in order to calculate insurance premiums.

Comparing the three models, we can see that the best one is Model I, this model incorporating non telematic and telematic variables. The effect of the driving experience on the frequency and severity is negative; on the contrary, the effects of the car characteristics are positive in the sense that older car and higher power increase the accident rates. The coefficients associated with the telematic variables are positive, their p-values indicating significant effects on the frequency. For the severity, the significant effects are associated to $X_8$ (Percentage of kilometers driven over the speed limit) and $X_9$ (Percentage of kilometers driven at night).

Focusing on Models II and III, we note that the best fit is obtained with the telematic variables; although, as we have seen in Model I, these variables are complemented by non telematic variables improving the goodness of fit in the full model. In short, the effects of the covariates hardly change if we compare models I, II and III, which implies that multicollinearity hardly affects the results, and that telematic and non-telematic variables complement each other and take into account different characteristics of the insured. The telematic variables are taking into account the exposure to risk that the non-telematic variables do not detect by themselves.

In

Table 4. Pure premium in collective model using parameters of Model I.

	Profile 1	Profile 2	Profile 3	Profile 4
$E\left(S_i\right)$, $\omega =110.2364$	687.7501	659.5172	359.0895	344.4958
$E\left(S_i\right)$, $\omega =0$	675.2375	647.6507	353.0899	338.8105
Difference	12.51264	11.86648	5.999624	5.685362

, we show the mean of $S_i$, that is, the pure premium, for four different insured profiles, using the full sample (with extremes) and with the positive dependence parameter $\omega$ estimated in Model I; this mean $E\left(S_i\right)$, is compared with the one obtained without dependence (i.e., dependence parameter equal to zero); we also calculate the differences between these two pure premiums, with and without frequency and severity dependence. The first and second profiles correspond to policyholders with high risk—they are 25 years old, with 6 years old driver license, with a car with 150 horsepower, 8100 yearly kilometers driven, 50% on urban area, 40% over the speed limit and 20% at night; the difference between these two profiles is that the first one does not use night parking, while the second does. The profiles 3 and 4 are policyholders with lower risk, having the same personal and car characteristics, same total km driven as the profiles 1 and 2, but with only 5% over the speed limit and 10% at night. We observe how the Sarmanov dependence parameter affects the pure premium. The difference between premiums calculated with and without dependence increases with risk, from 5.69 Euros for the lower risk profile (Profile 4) to 12.51 Euros for the higher risk profile (Profile 1).

4. Conclusions

In this paper, dependence between the claim frequency and the average severity of a policyholder was introduced, using a bivariate Sarmanov distribution with Negative Binomial GLM and Normal GLM marginals. The Normal GLM distribution was obtained by applying a Box-Cox transformation with two parameters on the original distribution, motivated by the fact that, in the collective risk model context, the claim cost distribution may not coincide with the traditionally used distributions (Exponential, Gamma, Lognormal, etc.). For some specific values of the transformation parameters, some useful closed type expressions for calculating the moments of the aggregate claims were obtained.

The proposed model was fitted on a sample of policyholders from an auto insurance portfolio. The peculiarity of these data is that they are associated with an insurance which involved the installation of a GPS/inertial device in the vehicle, device that allows the collection of telematic information related to the total kilometers and to the way which they are traveled. We have shown how these variables complement and do not substitute the variables that have traditionally been used in car insurance pricing. Furthermore, we have analyzed the importance of incorporating the dependence between frequency and severity in the collective model for calculating the premium, and how our Sarmanov model allows this dependence to be incorporated in a simple way.

In general, the bivariate Sarmanov distribution shows how the effect of dependency between frequency and severity is different depending on the risk profiles. For our dataset, the results show that riskier profiles have a greater effect of dependency. Furthermore, the Box-Cox transformation improves the distribution fit of the cost of claims variable, given that it allows a longer and heavier right tail than classical models like Gamma and Lognormal.