Multivariate Credibility in Bonus-Malus Systems Distinguishing between Different Types of Claims

Abstract

In the classical bonus-malus system the premium assigned to each policyholder is based only on the number of claims made without having into account the claims size. Thus, a policyholder who has declared a claim that results in a relatively small loss is penalised to the same extent as one who has declared a more expensive claim. Of course, this is seen unfair by many policyholders. In this paper, we study the factors that affect the number of claims in car insurance by using a trivariate discrete distribution. This approach allows us to discern between three types of claims depending wether the claims are above, between or below certain thresholds. Therefore, this model implements the two fundamental random variables in this scenario, the number of claims as well as the amount associated with them. In addition, we introduce a trivariate prior distribution conjugated with this discrete distribution that produce credibility bonus-malus premiums that satisfy appropriate traditional transition rules. A practical example based on real data is shown to examine the differences with respect to the premiums obtained under the traditional system of tarification.

1. Introduction

In an attempt of reducing the economic and casualty losses, the bonus-malus systems (BMS) have been introduced in the actuarial community. BMS is a pricing system mainly used in Europe in vehicle insurance. In this systems, the insured may have his/her premium discounted or penalized based on his/her own experience of claims. Actuarial literature about this topic is extensive see for example Lemaire; (1985, 1995); Boucher et al. (2007); Mert and Saykan (2005); Sarabia et al. (2004); Denuit et al. (2009), among other papers. Different methodologies have been used to determine the fair premium that policyholders must pay for the different classes in which the system is configured. Among these methods, the most popular ones are Bayesian methods and discrete Markov chains.

In Bonus-Malus systems, the premium is usually computed by using only the random variable number of claims. However, as not all the events produce the same individual claim amount Then, as different claims produce different claim sizes, it would be sensible to develop BMS based on both the number of claims and the corresponding severity (see Gómez-Déniz (2016)). In addition to this, if the severity is not included in the bounus-malus premium (BMP), the independence between the claims number and severity is implicitly assumed (see Lemaire (2004)). In this regard, several papers have discussed the question of implementing both variables in the BMS. See, for example, Frangos and Vrontos (2001); Gómez-Déniz et al. (2014) and the recent paper Gómez-Déniz (2016).

In this article we develop a trivariate model that discriminates between three types of claims. Since this distinction is made with respect to claims that are above, between or below certain thresholds, the resulting model takes into account the two fundamental random variables in this scenario, the number of claims as well as the corresponding claim size. Furthermore, we incorporate a trivariate prior distribution which is conjugated with the latter discrete trivariate model. As a result, we obtain credibility BMP’s which satisfy desirable transition rules. We present an example consisting of real data corresponding to an Australian portfolio of automobile insurance claims. Our findings reveal that the BMP’s computed by using the methodology proposed in this work (unlike those ones derived under the traditional Poison-Gamma model) does not modify the discounts make in the absence of claims. However, the methodology used in this paper is different to the recent developments in multivariate credibility models that can be found in the recent actuarial literature (see for example Frees (2003); Bühlmann and Gisler (2005)). Multidimensional credibility models was considered by Englund et al. (1999). In that paper the authors assumed that each dimension of the risk parameter represented one cover from the business. However, they only used frequency information in the credibility approach. Similarly, Thuring (2011) investigated the effect of assuming that one out of two insurance products is inactive when estimating the latent risk profile. Moreover Thuring et al. (2012) used a multivariate credibility model that allows the practitioner to consider the positive correlation in customer behaviour between different financial products and estimate the customer specific risk profiles for a specific product not owned by the customer. Again, this approach uses only two quantities, the a priori expected number of events and the observed number of events.

The Bayesian methodology has been used in actuarial science since the mid-twentieth century and it has proved to be a useful tool for the calculation of insurance premiums. It generally consists of accepting that each policy or insured is represented by a risk parameter that is unknown but random with a certain probability distribution (in the insurance portfolio), called a priori distribution or structure function. This way of proceeding is even more useful in the BMS scenario since the premium obeys certain transition rules that classify the policyholders as a bonus or malus. In other words, it lower the premiums to be paid for the bonuses and it increase the premiumt for the malus. The fundamental Bayesian tool here is simply the Bayes’s Theorem, so by dividing the a posteriori mean of the parameter by the a priori mean, when the net premium principle or the quadratic loss function are used, see Denuit et al. (2009); Lemaire (1985, 1995, 2004); Sarabia et al. (2004), among others, we obtain an estimator of the risk/s parameter/s that will indirectly divide the insured into good risks and bad risks. In addition, the empirical results illustrated in Frangos and Vrontos (2001); Mert and Saykan (2005); Gómez-Déniz et al. (2014) show that many auto insurance portfolios present a positive correlation between these two random variables, and therefore, the assumption of some kind of dependence between them should be considered in the calculation of BMP’s. The pioneering research in this field was that of Picard (1976) (see also Lemaire (1995), chp. 13), who divided the claims into two types: small, those ones that are below a limit value, say $\psi$, and large, above $\psi$. Then, as this assumption did not produce a good fit, the author proposed to distinguish between accidents that caused property damage and those ones that caused personal injury.

The rest of this paper is organised as follows. Section 2 describes the basic distributional assumptions formulated for the numbers of either type of claims. Section 3 discusses a trivariate conjugate, with respect to the discrete model, prior distribution. Next, Bayesian BMP’s are derived and written as credibility formula. Numerical applications are displayed in Section 4. Finally, Section 5 concludes the article.

2. Basic Model

There is a lot of criticism about assuming that the number of claims in an auto insurance portfolio can follow a Poisson distribution due to the fact that for this distribution the dispersion index (ratio between the variance and the mean) is one, when in the auto insurance portfolios have been empirically proven to be a value slightly higher than the unit. However, as an initial starting point and facilitating the methodology that we will develop, we will assume that, in effect, the number of claims has a Poisson distribution with parameter $\theta >0$ and probability function given by

$$\begin{array}{c}\hfill f\left(x\right|\theta )=\frac{1}{x!}{\theta }^{x}exp(-\theta ),x=0,1,\dots \end{array}$$

(1)

When the ith policyholder makes a claim $x_i$, it has associated a certain size, say $y_i\ge 0$. It is our interest now to distinguish between different types of claims (three in our case). For that reason, we include two new random variables that give rise to the consideration of three separate sub-events as follows. Let us to consider $Z_i^0$, $Z_i^1$ and $Z_i^2$, the following random variables

$$\begin{array}{ccc}{Z}_i^0=\left\{\begin{array}{cc}1,\hfill & \hfill y_i\le {\varphi }_1,\\ 0,\hfill & \hfill \mathrm{otherwise},\end{array}\right.\hfill & Z_i^1=\left\{\begin{array}{cc}1,\hfill & \hfill {\varphi }_1<y_i\le {\varphi }_2,\\ 0,\hfill & \hfill \mathrm{otherwise},\end{array}\right.\hfill & Z_i^2=\left\{\begin{array}{cc}1,\hfill & \hfill {\varphi }_2<y_i,\\ 0,\hfill & \hfill \mathrm{otherwise},\end{array}\right.\hfill \end{array}$$

where ${\varphi }_1$ and ${\varphi }_2,\in {\mathbb{R}}^{+}$ with ${\varphi }_2>{\varphi }_1$.

The $Z_i^j$’s, $j=0,1,2$ are modelled as independent and identically distributed random variables with the following Bernoulli distributions:

$$f(z_i^j|p_i)=\left\{\begin{array}{ccc}{p}_i,\hfill & \mathrm{if}& \hfill z_i^j=1,\\ 1-p_i,\hfill & \mathrm{if}& \hfill z_i^j=0,\end{array}\right.$$

where $0<p_i<1$. Observe that these assumptions imply that $E\left(Z_i\right)=p_i$, $(i=0,1,2)$. Since in practise majority of policyholders in the portfolio does not make a claim and those ones that declare claims with large claim size are sparse, we will also assume also $p_0\ge p_1\ge p_2$ with ${\sum }_{j=0}^2{p}_j=1$.1

We now assume that $Z_1={\sum }_{i=1}^x{Z}_i^1$ is the total number of claims with a claim size between ${\varphi }_1>0$ and ${\varphi }_2>0$ and $Z_2={\sum }_{i=1}^{x-z_1}{Z}_i^2$ is the total number of claims with a claims size larger than ${\varphi }_2$. Thus, if the $Z_i^j$, $i=1,\dots ,x$, $j=1,2$, are assumed to be mutually independent, then the conditional probability function of $Z_1$, given that $X=x$, is binomial with parameters x and $p_1$ and the conditional probability function of $Z_2$, given that $X=x$ and $Z_1=z_1$ is also binomial with parameters $x-z_1$ and $p_2$. That is,

$$\begin{array}{ccc}\hfill f\left(z_1\right|x,p_1)& =& \left(\genfrac{}{}{0pt}{}{x}{z_1}\right)p_1^{z_1}{(1-p_1)}^{x-z_1},z_1=0,1,\dots ,x,\hfill \end{array}$$

(2)

$$\begin{array}{ccc}\hfill f\left(z_2\right|x,z_1,p_2)& =& \left(\genfrac{}{}{0pt}{}{x-z_1}{z_2}\right)p_2^{z_2}{(1-p_2)}^{x-z_1-z_2},z_2=0,1,\dots ,x-z_1.\hfill \end{array}$$

(3)

Thus, the conditional mean and variance are linear. They are given by $E\left(Z_i\right|x,(i-1)z_{i-1},p_i)=(x-(i-1)z_{i-1})p_i$ and $var\left(Z_i\right|x,(i-1)z_{i-1},p_i)=(x-(i-1)z_{i-1})p_i(1-p_i)$, $i=1,2$.

Now, by conditioning it is easy to get the joint probability function of the random variable $(X,Z_1,Z_2)$ which results,

$$\begin{array}{c}\hfill f(x,z_1,z_2|\theta ,p_1,p_2)=\frac{{\left(\theta q_1{q}_2\right)}^{x}exp(-\theta )}{z_1!z_2!(x-z_1-z_2)!}{\left(\frac{p_1}{q_1{q}_2}\right)}^{z_1}{\left(\frac{p_2}{q_2}\right)}^{z_2},\end{array}$$

(4)

for $x=0,1,\dots$, $z_1=0,1,\dots ,x$, $z_2=0,1,\dots ,x-z_1$, being $q_i=1-p_i$, $i=1,2$. Observe that the distribution depends on three parameters every one related with the three types of claims.

Straightforward algebra provides moments and the cross moment,

$$\begin{array}{ccc}\hfill E\left(X\right|\theta )& =& \theta ,\hfill \end{array}$$

(5)

$$\begin{array}{ccc}\hfill E\left(Z_1\right|\theta ,p_1)& =& \theta p_1,\hfill \end{array}$$

(6)

$$\begin{array}{ccc}\hfill E\left(Z_2\right|\theta ,p_1,p_2)& =& \theta p_2{q}_1,\hfill \end{array}$$

(7)

$$\begin{array}{ccc}\hfill E\left(X{Z}_1{Z}_2\right|\theta ,p_1,p_2)& =& (2+\theta ){\theta }^2{p}_1{p}_2{q}_1.\hfill \end{array}$$

(8)

Estimation

Given a sample $(\tilde{x},{\tilde{z}}_1,{\tilde{z}}_2)$, where t is the sample size, estimation of the parameters $\theta$, $p_1$ and $p_2$ via maximum likelihood method are easily obtained. They result $\widehat{\theta }=\overline{x}$, ${\widehat{p}}_1={\overline{z}}_1/\overline{x}$ and ${\widehat{p}}_2={\overline{z}}_2/(\overline{x}-{\overline{z}}_1)$, where $\overline{x}$, ${\overline{z}}_1$ and ${\overline{z}}_2$ are the sample mean of the three random variables, respectively. These estimators coincide with the moment estimators obtained by using (5)–(7). The Fisher information matrix is given by

$$\mathcal{J}(\widehat{\theta },{\widehat{p}}_1,{\widehat{p}}_2)=\mathrm{diag}\left[\frac{t}{\widehat{\theta }},\frac{t\widehat{\theta }}{{\widehat{p}}_1(1-{\widehat{p}}_1)},\frac{t\widehat{\theta }(1-{\widehat{p}}_1)}{{\widehat{p}}_2(1-{\widehat{p}}_2)}\right],$$

from which the asymptotic variance-covariance matrix of $(\widehat{\theta },{\widehat{p}}_1,{\widehat{p}}_2)$ is obtained by inverting this information matrix. The score equations used to estimate the parameters and the Fisher’s information matrix are provided in the Appendix A.

3. Contemplating Heterogeneity

Let us assume now that the model includes a certain level of heterogeneity and it allows parameters $\theta$, $p_1$ and $p_2$ to vary among insureds in the portfolio. In this regard, we suppose that the parameter $\theta$ follows a gamma prior distribution (structure function) with a shape parameter $\alpha >0$, a scale parameter $\beta >0$ and a probability density function given by

$$\begin{array}{c}\hfill {\pi }_1\left(\theta \right)=\frac{{\beta }^{\alpha }}{\Gamma \left(\alpha \right)}{\theta }^{\alpha -1}exp(-\beta \theta ),\theta >0,\end{array}$$

with mean and variance are given by $E\left(\theta \right)=\alpha /\beta$ and $var\left(\theta \right)=\alpha /{\beta }^2$, respectively.

The $p_i$ parameters are assumed to follow a beta prior distribution with parameters ${\alpha }_i>0$ and ${\beta }_i>0$, $i=1,2$. That is, the probability density function of $p_i$ are given by

$$\begin{array}{c}\hfill {\pi }_i\left(p_i\right)=\frac{1}{B({\alpha }_i,{\beta }_i)}{p}_i^{{\alpha }_i-1}{q}_i^{{\beta }_i-1},0<p_i<1,\end{array}$$

(9)

respectively where $B(a,b)$ represents the beta function given by $B(a,b)=\Gamma \left(a\right)\Gamma \left(b\right)/\Gamma (a+b)$ and $\Gamma (·)$ is the Euler gamma function. The mean and variance of these prior distributions, given by (9) for $i=1,2$, are provided by $E\left(p_i\right)={\alpha }_i/({\alpha }_i+{\beta }_i)$ and $var\left(p_i\right)={\alpha }_i{\beta }_i/({({\alpha }_i+{\beta }_i)}^2({\alpha }_i+{\beta }_i+1))$, $i=1,2$, respectively.

The flexibility of the beta distribution allows it to take up different shapes depending on the values of its two parameters. The choice of these prior distributions obeys to the fact that they are conjugate with the likelihoods (see Heilmann (1989); Denuit et al. (2009); Klugman et al. (2008); among others). For that reason, we will assume independence between the random variables $\theta$, $p_1$ and $p_2$ by taking $\pi (\theta ,p_1,p_2)={\pi }_1\left(\theta \right){\pi }_2\left(p_1\right){\pi }_3\left(p_2\right)$.

Given a sample $(\tilde{x},{\tilde{z}}_1,{\tilde{z}}_2)$, where t is the sample size, the posterior distribution of $\vartheta =(\theta ,p_1,p_2)$ given the sample information is computed according to Bayes’s theorem and is proportional to the product of the likelihood and the prior distribution. Thus we find that the likelihood function is proportional to

$$\begin{array}{c}\hfill \mathcal{L}\left((\tilde{x},{\tilde{z}}_1,{\tilde{z}}_2)\right|\vartheta )\propto {\left(\theta q_1{q}_2\right)}^{t\overline{x}}{\left(\frac{p_1}{q_1{q}_2}\right)}^{t{\overline{z}}_1}{\left(\frac{p_2}{q_2}\right)}^{t{\overline{z}}_2}exp(-t\theta )\end{array}$$

(10)

and the prior distribution is proportional to

$$\begin{array}{c}\hfill \pi \left(\vartheta \right)\propto {\theta }^{\alpha -1}exp(-\beta \theta )\prod _{i=1}^2{p}_i^{{\alpha }_i-1}{q}_i^{{\beta }_i-1}.\end{array}$$

Thus, the posterior distribution is conjugated with respect to the likelihood (10) and it is described by

$$\begin{array}{c}\hfill {\pi }^{*}\left(\vartheta \right|(\tilde{x},{\tilde{z}}_1,{\tilde{z}}_2))\propto {\theta }^{\alpha +t\overline{x}-1}exp\left(-(\beta +t)\theta \right)\prod _{i=1}^2{p}_i^{{\alpha }_i+t{\overline{z}}_i-1}{q}_i^{{\beta }_i+t(\overline{x}-{\overline{z}}_1-(i-1){\overline{z}}_2)-1},\end{array}$$

where the constant of proportionality does not depend on $\theta$, $p_1$ and $p_2$. Here, $\overline{x}=(1/t){\sum }_{i=1}^t{\tilde{x}}_i$, ${\overline{z}}_i=(1/t){\sum }_{i=1}^t{\tilde{z}}_i$, $i=1,2$, are the sample means of X, $Z_1$ and $Z_2$, respectively.

Therefore, the posterior distribution is the product of a gamma and two beta distributions, with the updated parameters given in

Table 1. Updated parameters of the posterior distribution.

Parameter	Updated Parameter
$\alpha$	${\alpha }^{*}=\alpha +t\overline{x}$
$\beta$	${\beta }^{*}=\beta +t$
${\alpha }_1$	${\alpha }_1^{*}={\alpha }_1+t{\overline{z}}_1$
${\beta }_1$	${\beta }_1^{*}={\beta }_1+t(\overline{x}-{\overline{z}}_1)$
${\alpha }_2$	${\alpha }_2^{*}={\alpha }_2+t{\overline{z}}_2$
${\beta }_2$	${\beta }_2^{*}={\beta }_2+t(\overline{x}-{\overline{z}}_1-{\overline{z}}_2)$

.

In the numerical applications Section later we will adopt an empirical Bayes approach where the parameters of the prior distributions can be estimated from the data (see Robbins (1964); Casella (1985)). In order to do this, we need the marginal (unconditional) distribution of $(X,Z_1,Z_2)$, that can be easily obtained by compounding. Due that the variables are separated the integration process is simple. Thus, the unconditional distribution results,

$$\begin{array}{cc}\hfill f(x,z_1,z_2)& ={\int }_0^{\infty }{\int }_0^1{\int }_0^{1}f(x,z_1,z_2|\vartheta )\pi \left(\vartheta \right)d\vartheta \hfill \\ & =\frac{1}{z_1!z_2!(x-z_1-z_2)!}\mathcal{NB}\left(\alpha ,\frac{\beta }{1+\beta }\right)\mathcal{BB}\left({\alpha }_1+z_1,{\beta }_1+x-z_1\right)\hfill \\ & \times \mathcal{BB}\left({\alpha }_2+z_2,{\beta }_2+x-z_1-z_2\right),\hfill \end{array}$$

(11)

where $\mathcal{NB}$ represents the negative binomial distribution and $\mathcal{BB}$ the compound binomial-beta distribution. Some algebra provides that the probability function (12) of this trivariate unconditional model can be rewritten in a compact form as,

$$\begin{array}{c}\hfill f(x,z_1,z_2)=\frac{{\beta }^{\alpha }\left(\genfrac{}{}{0pt}{}{x+\alpha -1}{x}\right)\left(\genfrac{}{}{0pt}{}{z_1+{\alpha }_1-1}{z_1}\right)\left(\genfrac{}{}{0pt}{}{z_2+{\alpha }_2-1}{z_2}\right)\left(\genfrac{}{}{0pt}{}{x-z_1+{\beta }_1-1}{x-z_1}\right)\left(\genfrac{}{}{0pt}{}{x-z_1-z_2+{\beta }_2-1}{x-z_1-z_2}\right)}{{(1+\beta )}^{\alpha +x}\left(\genfrac{}{}{0pt}{}{x+{\alpha }_1+{\beta }_1-1}{x}\right)\left(\genfrac{}{}{0pt}{}{x-z_1+{\alpha }_2+{\beta }_2-1}{x-z_1}\right)}.\end{array}$$

(12)

The unconditional means are obtained using (5)–(8) by compounding and given by

$$\begin{array}{ccc}\hfill E\left(X\right)& =& \frac{\alpha }{\beta },\hfill \\ \hfill E\left(Z_1\right)& =& \frac{\alpha }{\beta }\frac{{\alpha }_1}{{\alpha }_1+{\beta }_1},\hfill \\ \hfill E\left(Z_2\right)& =& \frac{\alpha }{\beta }\frac{{\alpha }_2}{{\alpha }_2+{\beta }_2}\frac{{\beta }_1}{{\alpha }_1+{\beta }_1},\hfill \\ \hfill E\left(X{Z}_1{Z}_2\right)& =& \frac{\alpha {\alpha }_1{\alpha }_2{\beta }_1(\alpha +1)(\alpha +2\beta +2)}{{\beta }^3\left({\alpha }_1+{\beta }_1\right)\left({\alpha }_2+{\beta }_2\right)\left({\alpha }_1+{\beta }_1+1\right)}.\hfill \end{array}$$

The Premiums

Premiums can be derived by following the ideas displayed in Gómez-Déniz (2016). Let

$$\begin{array}{c}\hfill g(x,z_1,z_2)=p_{z_2}{z}_2+p_{z_1}{z}_1+p_x(x-z_1-z_2),\end{array}$$

(13)

be an appropriate function of the number of claims with claim size below ${\varphi }_1>0$, between ${\varphi }_1$ and ${\varphi }_2>0$ and above ${\varphi }_2$, where $p_x$, $p_{z_1}$ and $p_{z_2}$ are appropriate weights assigned to the number of claims with different types of size. It is also reasonable to assume that $p_{z_2}>p_{z_1}>p_x$.

Now, by using the net premium principle, i.e., the squared-error loss function, and simple algebra, we obtain the risk premium,

$$\begin{array}{c}\hfill P\left(\vartheta \right)=\left(q_1{p}_{z_2}{p}_2+p_1{p}_{z_1}+q_1{q}_2{p}_x\right)\theta .\end{array}$$

(14)

Observe that if $p_{z_1}=p_{z_2}=p_x=1$, then the risk premium in (14) is simply $P\left(\theta \right)=\theta$, that is, the risk premium obtained under the traditional model (net premium principle). The premium depends only on the number of claims, irrespective of their size.

This is the fair premium to be charged to a policyholder if $\theta$, $p_1$ and $p_2$ were known. However, these quantities are unobservable in practice and then, the risk premium is a theoretical one which cannot be determined exactly and it must be estimated from the data. On the other hand, the a priori premium is obtained for a policyholder about whom nothing is known, i.e., the average premium for all possible risk premiums.

We now obtain the a priori (collective) premium, as follows:

$$\begin{array}{c}\hfill P={\int }_0^{\infty }{\int }_0^1{\int }_0^{1}P\left(\vartheta \right)\pi \left(\vartheta \right)d\vartheta =\frac{\alpha }{\beta }\frac{p_{z_1}{\alpha }_1({\alpha }_2+{\beta }_2)+{\beta }_1(p_{z_2}{\alpha }_2+p_x{\beta }_2)}{({\alpha }_1+{\beta }_1)({\alpha }_2+{\beta }_2)}.\end{array}$$

(15)

Again, by inserting $p_{z_1}=p_{z_2}=p_x=1$ in (15) we obtain the collective premium computed under the traditional model. That is, $P=\alpha /\beta$.

The Bayesian premium $P^{*}(t,x,z_1,z_2)$, which is no reproduced here, is derived from (15) by interchanging the parameters $\alpha$, $\beta$, ${\alpha }_i$ and ${\beta }_i$$(i=1,2)$ with the updated parameters by using the expressions displayed in Table 1.

Note that $P^{*}(0,0,0,0)=P$. That is, the Bayesian premium coincides with the a priori premium when no information is available. Furthermore, the expression of the Bayesian premium can be written as

$$\begin{array}{c}\hfill P^{*}(x,z_1,z_2,t)=\gamma (x,z_1,t)P+[1-\gamma (x,z_1,t)]h(x,z_1,z_2,t),\end{array}$$

where

$$\begin{array}{ccc}\hfill \gamma (x,z_1,t)& =& \frac{{\alpha }^{*}({\alpha }_1+{\beta }_1)({\alpha }_2+{\beta }_2)}{{\beta }^{*}({\alpha }_1^{*}+{\beta }_1^{*})({\alpha }_2^{*}+{\beta }_2^{*})},\hfill \\ \hfill h(x,z_1,z_2,t)& =& \frac{{\alpha }^{*}\left[{\mathcal{H}}_1+{\mathcal{H}}_2+(x-z_1)(p_{z_2}{\alpha }_2+p_x{\beta }_2)\right]}{{\beta }^{*}({\alpha }_1^{*}+{\beta }_1^{*})({\alpha }_2^{*}+{\beta }_2^{*})-{\alpha }^{*}({\alpha }_1+{\beta }_1)({\alpha }_2+{\beta }_2)},\hfill \end{array}$$

(16)

where

$$\begin{array}{ccc}\hfill {\mathcal{H}}_1& =& p_{z_1}{\alpha }_1^{*}(x-z_1)+p_{z_1}{z}_1({\alpha }_2+{\beta }_2),\hfill \\ \hfill {\mathcal{H}}_2& =& {\beta }_1^{*}\left[p_{z_2}{z}_2+p_x(x-z_1-z_2)\right].\hfill \end{array}$$

Additionally, $\gamma (x,z_1,z_2,t)$ can also be written as

$$\begin{array}{c}\hfill \gamma (x,z_1,z_2,t)=Z\left(t\right)\frac{\alpha }{\beta }\frac{({\alpha }_1+{\beta }_1)({\alpha }_2+{\beta }_2)}{({\alpha }_1^{*}+{\beta }_1^{*})({\alpha }_2^{*}+{\beta }_2^{*})}+\left[1-Z\left(t\right)\right]\frac{({\alpha }_1+{\beta }_1)({\alpha }_2+{\beta }_2)\overline{x}}{({\alpha }_1^{*}+{\beta }_1^{*})({\alpha }_2^{*}+{\beta }_2^{*})},\end{array}$$

where

$$\begin{array}{c}\hfill Z\left(t\right)=\frac{\beta }{\beta +t}=\frac{\kappa }{t+\kappa }\in [0,1],\end{array}$$

with $\kappa =E\left[var\right(X\left|\theta \right)]/var[E\left(X\right|\theta \left)\right]$, coincides with the classical credibility factor usually appearing in this setting in actuarial science (see Bühlmann (1967); Jewell (1974); Gómez-Déniz (2008), among others for details). Now it is simple to see that:

• When $t\to 0$, $Z\left(0\right)\to 1$, $\gamma (x,z_1,z_2,0)\to \frac{\alpha }{\beta }$ and therefore $P^{*}(x,z_1,z_2,0)\to \frac{\alpha }{\beta }P$. Then, the premium is based only in the prior information about the risk. Therefore, the case is the one in which experience is ignored and external information is used as the sole basis for the process of ratemaking.
• When $t\to \infty$, $Z(\infty )\to 0$, $\gamma (x,z_1,z_2,\infty )\to 0$ and therefore $P^{*}(x,z_1,z_2,\infty )\to h(x,z_1,z_2,\infty )$. Then, the premium is based only in the sample information.

The Bayesian BMP can now be obtained by considering the rate (see Lemaire (1995); Gómez-Déniz et al. (2002); among others).

$$\begin{array}{c}\hfill P^{**}(x,z_1,z_2,t)=\frac{P^{*}(x,z_1,z_2,t)}{P^{*}(0,0,0,0)}=\frac{P^{*}(x,z_1,z_2,t)}{P},\end{array}$$

(17)

which ensures that the initial premium, i.e., when $(x,z_1,z_2,t)=(0,0,0,0)$, the rate in (17) is precisely P and the rates achieved for the first year is given by $P^{*}(x,z_1,z_2,1)$, for the second year $P^{*}(x,z_1,z_2,2)$, etc.

4. Numerical Applications

Now, in order to compute the premiums based on the models introduced in this paper, we will examine a dataset that include information based on one-year vehicle insurance policies taken out in 2004 or 2005. This dataset is available on the website of the Faculty of Business and Economics, Macquarie University (Sydney, Australia) (see also De Jong and Heller (2008)). The total portfolio contains 67,856 policies of which 4624 have at least one claim. With respect to the number of claims, the minimum and maximum are 0 and 4 respectively. The mean is 0.072 and standard deviation is 0.278. On the other hand, regarding the claim size, the minimum and maximum are 0 and 55,922.10 respectively. The mean is 137.27 and the standard deviation is 1056.30. This latter measure is very large for the size of the claims, therefore it means that a premium based only on the mean claim size is not adequate for calculating the bonus-malus premiums. Due to this portfolio only includes the aggregate value of the claims severity, a simulation analysis was completed to randomly determine the exact value that corresponds to all claims. Then, we proceed to allocate the claims that correspond to each interval, i.e., $0–$\$500$, $500–$1000 and >$1000. Thus, we are assuming that ${\varphi }_1=500$ and ${\varphi }_2=1000$. This simulation analysis was carried by using Mathematica software package. We have taken the integer part of the individual claim amount, this does not seem very relevant in our analysis. It is important to mention that due to RandomChoice function, the partition of the aggregate claim amount is different every time the program is run.

Empirical values and fitted values by using the discrete trivariate distribution, Fitted (1), and the mixing model (12), Fitted (2), are illustrated in

Table 2. Empirical and fitted data for the basic model (1) and mixture model (2).

$\mathbf{(}\mathit{x}\mathbf{,}{\mathit{z}}_{\mathbf{1}}\mathbf{,}{\mathit{z}}_{\mathbf{2}}\mathbf{)}$	Empirical	Fitted (1)	Fitted (2)	$\mathbf{(}\mathit{x}\mathbf{,}{\mathit{z}}_{\mathbf{1}}\mathbf{,}{\mathit{z}}_{\mathbf{2}}\mathbf{)}$	Empirical	Fitted (1)	Fitted (2)
$(0,0,0)$	63,232	63,094.30	63,233.20	$(3,0,0)$	0	0.29	1.49
$(1,0,0)$	1840	1921.02	1812.24	$(3,1,0)$	5	1.05	4.75
$(1,1,0)$	2084	2257.62	2128.08	$(3,0,1)$	0	0.19	0.44
$(1,0,1)$	409	411.91	387.96	$(3,0,2)$	0	0.04	0.22
$(2,0,0)$	31	29.24	51.83	$(3,1,1)$	3	0.45	1.28
$(2,1,0)$	134	68.73	113.61	$(3,2,1)$	0	0.26	1.11
$(2,0,1)$	7	12.54	13.98	$(3,1,2)$	0	0.05	0.51
$(2,1,1)$	16	14.74	24.32	$(3,3,0)$	3	0.48	2.04
$(2,2,0)$	79	40.40	66.82	$(3,0,3)$	0	0.00	0.10
$(2,0,2)$	4	1.34	5.60	$(4,·,·)$	4	0.35	0.89

. These sample values were taken from the results obtained after dividing the claims by using the simulation scheme mentioned above. The estimated parameter values (the standard errors appear in brackets) are shown in

Table 3. Parameters estimated and standard errors (SE) for the basic and mixture model without including covariates.

Basic Model			Mixture Model
Parameter	Estimate	SE	Parameter	Estimate	SE
$\widehat{\theta }$	0.072	0.001	$\widehat{\alpha }$	1.157	0.121
${\widehat{p}}_1$	0.492	0.007	$\widehat{\beta }$	15.903	1.681
${\widehat{p}}_2$	0.176	0.007	${\widehat{\alpha }}_1$	575.261	2.779
			${\widehat{\beta }}_1$	594.757	2.961
			${\widehat{\alpha }}_2$	0.365	0.268
			${\widehat{\beta }}_2$	1.705	1.257
${\chi }^2$	137.06			7.18
p-value	0.00			0.007
df	4			1
AIC	45,129.80			45,027.70
CAIC	45,160.20			45,088.50

. We also show the value obtained for two measures of model selection: Akaike’s information criterion (AIC) and the consistent Akaike information criterion (CAIC). See Akaike (1974); Bozdogan (1987) for details. The goodness of fit was determined by standard Pearson’s chi squared test statistics with the following grouping procedure: the outermost classes were consolidated to produce theoretical class sizes of 5 or larger. It is observable that the fit to data is reasonably good for the mixture model and not very promising for the basic model. For the mixture model, the maximum likelihood estimates were obtained by directly maximizing the log-likelihood surface.

The Proposed Premiums

Table 4. BMP’s for claims when there are x claims, $z_1$ with a claim size between ${\psi }_1$ and ${\psi }_2$, $z_2$ claims with a size larger than ${\psi }_2$ and $x-z_1-z_2$ claims with a claim size smaller than ${\psi }_1$ with $p_x=0.25$, $p_y=0.50$ and $p_z=0.75$.

$(\mathit{x},{\mathit{z}}_1,{\mathit{z}}_2)$	t
	0	1	2	3	4	5
$(0,0,0)$	1.000	0.940	0.888	0.841	0.799	0.760
$(1,0,0)$	1.798	1.692	1.597	1.513	1.437	1.368
$(1,1,0)$	1.864	1.754	1.656	1.568	1.489	1.418
$(1,0,1)$	2.168	2.040	1.926	1.824	1.732	1.649
$(2,0,0)$	2.583	2.430	2.295	2.173	2.064	1.965
$(2,1,0)$	2.633	2.477	2.339	2.215	2.104	2.003
$(2,1,1)$	3.174	2.986	2.819	2.670	2.536	2.414
$(2,2,0)$	2.729	2.568	2.424	2.296	2.180	2.076
$(2,2,1)$	3.530	3.321	3.135	2.969	2.820	2.685

illustrates the relativities (Bayesian BMP’s) obtained by applying (17) and the parameter estimates displayed in Table 3. It is noticeable that the structure of this table is built in a similar way the one derived in traditional BMS. Namely, at the beginning of the system the relativity is set equal to 1.000; then this relativity decreases within the year in the absence of claims, and it increases when claims are declared. Nevertheless, for $x>1$ the system now discerns whether the number of claims corresponds to those below the size ${\varphi }_1$, between ${\varphi }_1$ and ${\varphi }_2$ and above ${\varphi }_2$. For the sake of comparison, the reader is referred to Gómez-Déniz (2016) where the Bayesian bonus-malus premiums calculated for the Poisson-Gamma model under the net premium principle (see also Dionne and Vanasse (1989)) and those ones computed by using expression (20) in Gómez-Déniz (2016). It is observable that the bonus-malus premiums are the same for the bonus class $(x=0)$ and different for the rest of the malus $(x\ge 1)$ classes with respect to the first of the models mentioned above. In this regard, it can be distinguished now between claims with a severity below ${\varphi }_1$, between ${\varphi }_1$ and ${\varphi }_2$ and above ${\varphi }_2$. In this sense if we consider Table 3 in Gómez-Déniz (2016), for example for the case $(t,x,z)=(1,0,1)$ (i.e., at the end of the first year the policyholder has declared one claim with size higher than $500) the Bayesian bonus-malus premium is 1.800. However, under the scheme introduced in Table 4, the Bayesian BMP to be paid is 1.754 if the claim size belong to the interval $500–$1000, i.e., $(t,x,z_1,z_2)=(1,1,1,0)$. On the other hand, when the claim amount is >$1000, i.e., $(t,x,z_1,z_2)=(1,1,0,1)$ the Bayesian BMP is 2.040. Therefore, the premiums as shown in Table 4, may be larger or smaller than those ones shown in Table 3 in Gómez-Déniz (2016). This methodology would ensure the financial viability of the company.

5. Final Comments and Future Research

In this paper, a simple model that distinguishes, among three types of claims in bonus-malus settings has been introduced. This distinction is based on discriminating between those claims with associated amount below a threshold, between two values of thresholds or greater than a certain threshold. This methodology presented is based on the use of a trivariate distribution (not common in any statistical scenario) that depends on parameters that in turn are considered as random variables that follow certain a priori probability distributions. As a consequence, it is possible to express the bonus-malus premium based on the net premium principle (quadratic error loss function) as a credibility formula that writes the premium as a convex combination of sample information and a priori information. The bonuses of the premium obtained are undoubtedly fairer than those ones computed by using the classical methodology that does not discern between different types of claims. We shall conclude with an interesting comment made by a referee with respect to the range of values that parameter $p_2$ can take on. In this work we have assumed that $p_2\in (0,1)$, however, in the third section when the probabilities are randomized it could be sensible to suppose that $p_2$ is dependent on the value of $p_1$. In this sense, we believe it should be more realistic to consider $\pi \left(\vartheta \right)=\pi \left(\theta \right)\pi (p_1,p_2)$, where the latter factor is a bivariate distribution that assumes dependency between $p_1$ and $p_2$. This could be subject of future research. In addition, it would be interesting to examine how the premiums behave when both the distribution based on the classical model and the marginal model derived after including heterogeneity are normalized to implement a generalized linear model. This is likely to refine the premiums according to the individual factors of each insured.