A Priori Ratemaking Selection Using Multivariate Regression Models Allowing Different Coverages in Auto Insurance

: A comprehensive auto insurance policy usually provides the broadest protection for the most common events for which the policyholder would ﬁle a claim. On the other hand, some insurers offer extended third-party car insurance to adapt to the personal needs of every policyholder. The extra coverage includes cover against ﬁre, natural hazards, theft, windscreen repair, and legal expenses, among some other coverages that apply to speciﬁc events that may cause damage to the insured’s vehicle. In this paper, a multivariate distribution, based on a conditional speciﬁcation, is proposed to account for different numbers of claims for different coverages. Then, the premium is computed for each type of coverage separately rather than for the total claims number. Closed-form expressions are given for moments and cross-moments, parameter estimates, and for a priori premiums when different premiums principles are considered. In addition, the severity of claims can be incorporated into this multivariate model to derive multivariate claims’ severity distributions. The model is extended by developing a zero-inﬂated version. Regression models for both multivariate families are derived. These models are used to ﬁt a real auto insurance portfolio that includes ﬁve types of coverage. Our ﬁndings show that some speciﬁc covariates are statistically signiﬁcant in some coverages, yet they are not so for others.


Introduction
In the automobile insurance sector, it is natural to calculate the a priori premium taking into account the number of claims and individual characteristics of each insured, such as gender, age, years of validity of the policy, etc. This procedure to compute the a priori premium is usually completed via parametric models rather than using the ordinary regression model, which can predict values of the number of claims even if negative. For this purpose, parametric models based on the use of the Poisson, negative binomial, and Poisson-inverse Gaussian distributions, among others, are the standard models considered in the univariate case. As of today, most insurance companies distinguish, apart from the total number of claims, individualized claims for different coverages, such as windscreen claims, thefts and fire claims, etc. So far, most actuarial models aim to differentiate only between two types of coverage when computing an appropriate premium based on different coverage. Perhaps one of the reasons for that is due to the lack of models capable of describing more than two coverages. The most often considered approach to tackle this problem is the one based on the bivariate Poisson distribution (see Bermúdez 2009;Bermúdez and Karlis 2017, among others). See also Gómez-Déniz (2016); Gómez-Déniz and Calderín-Ojeda (2018); Gómez-Déniz and Calderín-Ojeda (2020); Denuit et al. (2009), andFrees (2010) for more details related to this topic. Alternative references for a review of count regression are Cameron and Trivedi (1986); Cameron and Trivedi (1998); Winkelmann (2003), and Boucher et al. (2007).
A copula-based correlated random effects model that accommodates dependence between claim frequency and severity was examined in Oh et al. (2020).
Traditionally, the business associated with insurance consists of selling risk coverage to buyers. In particular, in automobile insurance, the insurer provides financial protection against physical damage or bodily injury resulting from an incident (see Frees et al. 2016). However, it is common today, mainly due to the existing competition, that the insurance companies offer coverage of different claims within the same product not only to gain in competitiveness, but also to benefit from risk diversification and volatility. In this paper, we consider a motor vehicle insurance portfolio with policies observed during some time period that contain, apart from other known factors (gender, age, years of validity of the driver's license, etc.), information about the claims number concerning different coverages that are considered as response variables. This includes windscreen, parking, theft and fire, etc.
Therefore, it is assumed that the insurance company collects information on the claims for these coverages and the total number of claims given by the sum of the claims in all the coverages. Thus, every policyholder generates a sequence of claims numbers for each coverage; one of them is the total claims number, which includes the sum of the coverages' claims. Then, based on a conditional specification, a multivariate model that allows a simple way to describe the use of a finite but sufficiently large number of coverages is proposed. The resulting multivariate discrete distribution obtained enables us to study the dependence structure of a limited number of coverages in automobile insurance and include covariates such as gender, age, etc. We start by using a Poisson model for the random variable total claims number, and then by conditioning, we introduce the remaining variables in a branch architecture structure. Finally, closed-form expressions are given for parameter estimates, and a priori premiums are provided when different premium principles are used.
The purpose of this paper is to introduce a novel methodology based on a multivariate distribution via a conditional specification, proposed to account for different numbers of claims in different coverages and also for the total claims frequency. This approach enable us to examine the dependence structure of a finite number of coverages in motor vehicle insurance and also incorporate heterogeneity in the model through explanatory variables. Then, we use this procedure to calculate premiums based only on the claims frequency. Next, we show that the amount of claims can be incorporated into this multivariate model to derive multivariate claims' severity distributions. For this, we assume that the claims size in the joint coverages follows a multivariate Erlang distribution. As multivariate probability distributions are complex, it is argued that analytical solutions are highly unlikely as compared to those derived under univariate and bivariate cases (see Wiltbank 1983, 1984); nevertheless, in this work, we derive a multivariate model where the total number of claims that affect the portfolio is the result of the interaction of multivariate processes. The main advantage of the modelization presented in this work is that it avoids working with copulas (see, for instance, Balakrishnan and Lai 2009, chp. 1, p. 59). Although the copula approach for modeling multivariate models has been proven to be very useful, it has also been criticized due to the difficulty of choosing an appropriate copula structure and the complication of estimating the parameters that control the dependency. In addition, a multivariate zero-inflated model to account for the excess of common zeros in the empirical distribution is developed. Finally, these two multivariate distributions can be reparameterized to incorporate covariates to determine which factors and explanatory variables have an influence on the mean of the corresponding coverage. As an illustration, in this work, we use the French Motor Personal Line datasets available in the package "CASdatasets" in R, which include five response variables.
Although the modeling proposed here was developed ad hoc for the auto insurance market, it is unquestionable that other insurance lines in general insurance might benefit from it. For example, in home insurance, the whole premium could be split into different coverages such as moisture damage, theft, pipe repairs, locksmiths, and even protection against tenant rent default.
The rest of the paper is structured as follows. Section 2 describes the primary model and some of its properties. Then, premium calculations based on this basic model are discussed. Finally, a multivariate zero-inflated model and multivariate regression procedures are shown. Some methods of estimation are provided in Section 3. Next, a numerical application pertaining to a private motor French insurer is developed in Section 4. Finally, conclusions are drawn in Section 5.

The Branch Architecture Model
Let us consider a portfolio with N observed policies during T periods of time and also assume that the insurance company gathers information on the number of claims related to several types of coverages. Example of these coverages may include windscreens, fire and theft, etc. Therefore, the insurer collects information about these coverages, as well as the total number of claims for each policyholder given by the sum of the claims in all different coverages. For the ith policyholder, we consider the multivariate random variable expressed as the following sequence, N ji = (N ji1 , N ji2 , . . . , N jiT ) of claims numbers for coverage j, with j = 1, 2, . . . , J , assuming that one of them, i.e., the first one, is the total number of claims, which includes the sum of the claims for all types of coverages purchased by this policyholder.
Furthermore, we assume that N 1i , the total number of claims recorded in the auto insurance portfolio, follows a Poisson distribution with mean Θ 1i > 0 for i = 1, 2, . . . , N, where N is the total number of policyholders. Now, let us suppose that the policyholders have purchased some of the types of coverages, such as windscreen protection, fire and theft, parking, etc. That is, once the policyholder has made a claim, this can be of any of these types. Let us denote by Z 1iκ , κ = 1, . . . , N 1i , a random variable associated with the number of claims corresponding to the first type of coverage and policyholder i, resulting from the κth claim of the total claims reported by the ith policyholder assumed to be independent and identically distributed, following also a Poisson distribution with mean Θ 2i > 0. Then, the conditional distribution of N 2i given N 1i = n 1i , N 2i = ∑ N 1i κ=1 Z 1iκ , the total number of claims of this first coverage, among the N 1i total claims is a Poisson distribution with parameter n 1i Θ 2i and the joint distribution of (N 1i , N 2i ) has a probability function given by, for n 1i , n 2i = 0, 1, . . . , and i = 1, 2, . . . , N and with the convention that 0 j = 0 for j = 0 and 1 otherwise. This bivariate distribution appears in Leiter and Hamdan (1973) (see also Cacaoullos and Papageorgiou 1980;Johnson et al. 1996, chp. 37, p. 136 in the context of accident analysis) Let Z 2iκ , κ = 1, . . . , N 1i now be a random variable associated with the total number of claims corresponding to the second type of coverage and policyholder i, resulting from the κth claim of the total claims reported by the ith policyholder assumed to be independent and identically distributed Poisson distribution with mean Θ 3i > 0 and conditionally independent of N 2i . Then, the conditional distribution of N 3i given N 1i = n 1i , N 3i = ∑ N 1i κ=1 Z 2iκ , the total number of claims of this second coverage, among the N 1i total claims is a Poisson distribution with parameter n 1i Θ 3i , and now, the joint distribution of (N 1i , N 2i , N 3i ) has a probability function given by, f (n 1i , n 2i , n 3i |Θ 1i , Θ 2i , Θ 3i ) = Pr(n 1i , n 2i ) Pr(n 3i |(n 1i , n 2i )) = Pr(n 1i , n 2i ) Pr(n 3i |n 1i ) where the hypotheses of conditional independence between the two types of coverages were assumed. Following the same argument, it is easy to see that if we have J types of coverages, then the joint probability function of (N 1i , N 2i , . . . , N J i ) ∈ N J is given by: where Θ i = (Θ 1i , . . . , Θ J i ). For this multivariate distribution, it is allowed that n 1i takes larger or smaller values than n ji ; however, in the proposed model, it is verified that n 1i is larger than or equal to n J i for all J > 1. In this case, it is obvious that Θ 1 must be larger than Θ j , j = 2, . . . , J . The latter statement is confirmed in the numerical application section. This distribution is a multivariate extension of the bivariate one proposed in Leiter and Hamdan (1973) (see also Cacaoullos and Papageorgiou 1980).
The ordinary probability-generating function of (N 1 , . . . , N J ) with the probability mass function (pmf) given in (1) is given by: From here, it is easy to see that the marginal distribution of N 1i is Poisson with parameter Θ 1i , while N ji , j = 1, . . . , J have a Neyman Type A distribution with parameters Θ 1i and Θ ji . Recall that the probability function of the Neyman Type A distribution (see Neyman 1939;Douglas 1955;Kemp 1967;Johnson et al. 2005, chp. 8, among others) is given by: for j = 1, . . . , J . Some computations provide the marginal and cross-moments, which are given by: from which it is simple to see that: and therefore, the model admits only a positive correlation between pairs of random variables. The marginal variances are given by: Observe that, using (5) and (6) together with (3) and (4), the model is equidispersed (variance equal to the mean) for N 1i and overdispersed (variance larger than the mean) for the rest of the coverages.
Finally, the correlation can easily be computed as: One can be interested also in the distribution of N 1i given N ji = n ji . The probabilitygenerating function of this conditional distribution is given by: where Λ j = Θ 1i exp −Θ ji and B n (τ) are the Bell numbers given by: being the Stirling number of the second kind. 1 Now, the conditional mean of N 1i given by N ji = n ji can be written as:

Some Results in Risk Theory
Observe that due to the model construction, we have that ∑ J j=2 Θ ji = 1, i.e., every claim in coverage j is a proportion of the total claims N 1i . Then, if the actuary decides to use the net premium principle, i.e., P(X) = E(X), to compute the premium, then for the ith policyholder and coverage s with s ∈ {2, . . . , J }, the premium results P si = Θ 1i Θ si = P 1i Θ si , where P 1i = Θ 1i is the net premium for the total coverage, that is the sum of the premiums in each of the coverages purchased. A similar result is obtained by using the expected value principle. A catalog of premium principles can be found in Young (2006).
Let us now consider that the actuary decides to use the variance premium principle, i.e., P(X) = E(X) + var(X)/E(X) with E(X) > 0, to calculate the premium. Then, in this case, we obtain that P 1i = 1 + Θ 1i and P ji = 1 + Θ ji P 1i . However, in this case, we have that ∑ J j=2 P ji = J − 1 + P 1i , which is different from P 1i , except for the case in which J = 1 and no coverages exist. However, a model solely based on the number of claims is not realistic. In risk theory, it is common to incorporate the amount associated with each of the claims to build the compound model. That is, the property and/or casualty ratemaking are generally based on a claim frequency distribution and a loss distribution. Due to the complex derivation of this multivariate compound model, the subscript i is removed from the text in the remainder of this section. For this purpose, let us now assume that where E i is the random variable denoting the size or amount of the ith claim, following an exponential distribution with probability density function (pdf) h(y j ) = σ −1 exp(−y j /σ). Furthermore, we assume that E 1 , E 2 , . . . , are independent and identically distributed random variables and also independent of the number of claims N j . It is well known (see for example Rolski et al. 1999) that Y j follows a piecewise distribution with pdf given by g(y j ) = ∑ ∞ n j =1 f N j (n j )h * n j (y j ), y j > 0, and g(0) = f N j (0). Then, by following the methodology given in Lee and Lin (2012), we have that Y Y Y = (Y 1 , . . . , Y J ) follows a multivariate Erlang distribution with scale parameter σ > 0 and shape parameter n j > 0, j = 1, 2, . . . , J . Their marginal distributions are a univariate Erlang mixture.
Then, simple computations provide, Here, I 1 (·) represents the modified Bessel function of the first kind, which admits the following series representation, The distribution for the coverages can be computed by using (2) in the following way, Now, taking into account that: we finally obtain the aggregate claim size pdf for the different coverages given by, for j = 2, . . . , J . Thus, they are also given as a piecewise distribution. For practical purposes, the infinite sum that appears in this expression can be replaced by a finite sum from one to k, where k can take values around one-hundred. From the assumption of the independence between the number of claims and the claims size, we have that: which can be considered as the net premium when both the number and size are considered at the same time.
Finally, we have that: while the covariance (see Lee and Lin 2012) is given by:

Multivariate Zero-Inflated Model
In many automobile insurance portfolios, the claims are rarely observed as compared to the no-claims situation. Univariate and bivariate zero-inflated models have been introduced in the statistical literature in many fields. In the setting of auto insurance, we refer to Boucher et al. (2007) and Frees et al. (2016) for the univariate case and Bermúdez (2009) and Bermúdez and Karlis (2017) for the bivariate case. Multivariate ones are scarce in the general statistical literature. References in the statistical literature are Li et al. (1999) and Liu and Tian (2015). In the actuarial literature, there are no references of models of this nature that go beyond the two variables. However, multivariate zero-truncated models were considered in Zhang et al. (2020).
A multivariate zero-inflated model can be constructed as a mixture of the multivariate distribution given in (1) and a point mass at (0, . . . , 0) ∈ R J in the following way, where 0 ≤ Φ ≤ 1 is an inflation parameter. Obviously, this model reduces to (1) for Φ = 1. Under this model, the marginal means and cross-moments are given by: from which the covariance between pairs of marginal random variables can be obtained. They are given by: Again, if the actuary computes the premium by using the net premium principle, then for each coverage, the premiums are not affected by the inflation parameter Φ. A complete model would allow inflating each coverage with inflation parameters Φ j ; however, they are not included in this work due to the computational cost of estimating a large number of parameters.
The marginal variances are given by: Finally, the correlations are: for j, l = 2, . . . , J , j = l.

A regression Model
For the sake of convenience, the model (1) can be rewritten in a different way to facilitate the implementation of covariates to determine which factors and explanatory variables have an influence on the mean of the corresponding coverage. Then, by equating Θ 1i to µ 1i and Θ ji to µ ji /µ 1i , j = 2, . . . , J , µ 1i = 0, we obtain the normalized joint distribution, which can be expressed as, for n 1i = 0, 1, . . . , n ji = 0, 1, . . . , n 1i , j = 2, . . . , J . The probability function (12) satisfies the condition that the marginal means are given by E(N ji ) = µ ji , j = 1, 2, . . . , J , assuming that µ ji = 0, for all j. Thus, it is suitable for including covariates. Then, to carry out this regression model, we suppose that the observed counts (N 1i , . . . , N J i ) have independent distributions given by (12) with E(N ji ) = µ ji , j = 1, 2, . . . , J . Now, it is assumed that a set of observable covariates useful to subdivide the portfolio into classes of risks with homogeneous characteristics are included in the linear predictor, η ji . To guarantee a positive expected value of the response variables, it is reasonable and common to use a logarithmic link for this function and therefore express the mean as: where x ji = (x ji1 , . . . , x jim ) is a vector of m covariates for the ith observation µ ji and γ j = (γ j1 , . . . , γ jm ) denotes the corresponding vector of regression coefficients to be estimated, which usually includes a constant term. Without loss of generality, it is assumed that for each j = 1, . . . , J , N ji is related to the same set of covariates. In addition, one of the covariates may be identified as an exposure term to calibrate the size of a potential outcome variable by assuming that the mean varies proportionally with the exposure e ji (see Frees 2010;Frees et al. 2016), Similarly, the covariates can be implemented in the multivariate zero-inflated regression model by simply regressing the mean value of the different coverages. It should be pointed out, although it is not considered here, that it could also be assumed that the inflated parameter Φ could depend on certain regressors. This issue seems not to be possible here, and thus, it could be a subject that merits further investigation in future research.

Estimation of the Parameters
In this section, we firstly describe the methodology for the maximum likelihood estimation and derivation of the entries of Fisher's information matrix for the basic model. Next, the same development is illustrated for the associated regression model. Finally, the expression of the log-likelihood function, score equations, and the second derivative of the log-likelihood function with respect to the parameters for the zero-inflated model are exhibited.
In general, the statistical inference for multivariate models is not trivial, and the computational procedure is often expensive (see, for instance, Selch and Scherer 2010). Nevertheless, the estimation procedure for the model proposed here is straightforward. To see this, we first consider the case without covariates. Let us assume that a samplẽ n := {(ñ 11 , . . . ,ñ J 1 ), . . . , (ñ 1n , . . . ,ñ J n )} that includes n independent observations in each one of the J types of coverage is collected. The log-likelihood function is proportional to: After differentiating the latter expression, it is possible to obtain in closed-form the maximum likelihood estimators of the parameters. They are given by: , j = 2, . . . , J , wheren k , k = 1, . . . , J , represents the sample mean, i.e., ∑ n i=1ñ1k /n.
The elements that provide the entries of Fisher's information matrix are as follows: For the regression model, the log-likelihood function contains J × m parameters, and it is proportional to: where µ ji = exp{x ji γ j }, i = 1, 2, . . . , n and j = 1, . . . , J .
Fisher's information matrix is made up of four blocks, as can be seen below: where: where j, h = 2, . . . , J ; l, k = 1, . . . , m, and 0 is the zero matrix with dimension m × m.
For the zero-inflated model, the log-likelihood is proportional to: where n * is the number of zeroes of the random variable N 1i . The normal equations that provide the maximum likelihood estimates are given by, From (13), we obtain that Φ = (n * − n)/(n(exp(−Θ 1 ) − 1)), which can be carried out to (15) to obtain the estimator of Θ 1 , say Θ 1 . This value is carried out now to (13) to obtain the estimator of the inflated parameter, Φ. Finally, from (15), the estimator of Θ j , j = 2, . . . , J is obtained in the closed-form expression given by Θ j = ∑ n i=n * +1 n ji / ∑ n i=n * +1 n 1i . The second partial derivatives are as follows, Observe that the analytic expressions of ∑ n i=n * +1 n 1i and ∑ n i=n * +1 n ji are not feasible. For computational reasons, for large values of n, this is evaluated by ignoring the expectation operator and replacing it by ∑ n i=n * +1 n 1i and ∑ n i=n * +1 n ji . The asymptotic variance-covariance matrix is approximated by inverting the observed information matrix.
When covariates are introduced under the inflated model, we proceed first by replacing in (9) the pmf f by its corresponding (12), where again, µ ji = exp{x ji γ j }, i = 1, 2, . . . , n, and j = 1, . . . , J . In practice, as shown in the numerical applications below, the parameter estimation and computation of standard errors were carried out by the method of maximum likelihood using Mathematica v.12.0. We directly maximized the log-likelihood function by using different maximum search methods available in the Find-Maximum built-in function in the Mathematicasoftware package. This software package also provides at least two methods of obtaining the elements of the Hessian matrix. The first one consists of retrieving them from the Cholesky factors (this package is available on the web upon request). The second one, which is faster, derives them by finite differentiation. Results were also confirmed with WinRATS v.7.0.

Numerical Application
For our empirical analysis, we used the French Motor Personal Line datasets available in the package "CASdatasets" in R. This is a collection of ten datasets that comes from a private motor French insurer. Each dataset includes risk features such as claim amount, risk area, gender of the policyholder, number of claims for different coverages, etc. In particular, we chose the freMPL10 dataset, which includes 32,100 policies for the year 2004. In our study, we considered six response variables, which are shown in Table 1.
Note that the dependent variable Claims for each policyholder comprises the sum of the individual claims in all other variables. The details of the joint claims frequency for all types of coverage and the total number of claims are illustrated in Appendix A (Table A1). Note that the maximum number of claims reported by an insured is six. The number of policyholders who did not report a claim is 12,257 (38.18%), and the number of customers that only declared a claim in any of the coverages is 10,803 (33.65%). Together with all the responses, this dataset includes a set of explanatory variables. Table 2 below describes the factors and explanatory variables used in the investigation. We also considered an offset variable when modeling the claims frequency, exposure, the time exposed to risk during the investigation period.
In Table 3, the parameter estimates and their corresponding p-values are provided for the basic and zero-inflated models without covariates. Some measures of model selection are also provided in the bottom part of this table. For comparisons purposes, we used the multivariate negative binomial distribution (MNB) provided in (Johnson et al. 1996, chp. 36, p. 94) with the pmf given by: where α > 0, 0 < Θ i < 1, i = 1, 2, . . . , k, and ∑ k i=1 Θ i = 1. As can be seen in Table 3, the multivariate Poisson distribution studied in this paper has a better performance than the MNB for this dataset. Furthermore, it is observable that the the zero-inflated model improves the basic one due to the high frequency of zeros. On the other hand, we also tried to fit the two multivariate Poisson distributions provided in Bermúdez and Karlis (2011); however, we were unable to derive the maximum likelihood estimates of this model for this dataset. takes the value 1 if the usage of the vehicle is private + trip to office, 0 otherwise; professional takes the value 1 if the usage of the vehicle is professional, 0 otherwise; driver age the driver age, in years; has km limit takes the value 1 if there is a km limit, 0 otherwise; risk area Unknown risk area between 1 and 13, possibly ordered; bonus takes the value 1 if the numerical value for bonus/malus is less than 100, 0 otherwise; malus takes the value 1 if the numerical value for bonus/malus is larger than 100, 0 otherwise.  Table 4 below exhibits the empirical Pearson's correlation between the different frequencies associated with each response variable for the total number of claims, and each one of the different coverages (first row), the correlation derived computed via the basic model (second row) and zero-inflated model (third row), and that computed by using (7) and (10), respectively, are also shown. It is observable that there exists a weak positive correlation between Claims and the rest of the dependent variables for each one of the coverages, and the empirical values are near the theoretical values. These figures were calculated before incorporating the effect of the explanatory variables for the different coverages. We also calculated the correlation coefficient between the rest of the response variables. Again, there is a weak positive correlation, ranging from 0.0480 between Responsible and Nonresponsible and 0.0035 between Parking and Windscreen. Table 4. Correlation between empirical frequencies associated with the total claims number and each response variables (first row) and the correlation derived by means of the basic model (second row) and zero-inflated model (third row). Empirical marginal distributions and fitted marginals under the basic model (Fit 1) and zero-inflated model (Fit 2) are illustrated below in Table 5 using the estimates computed in the previous section. Note that the total number of observations equals 32,100.

Nonresponsible
We fit the multivariate regression model in (12) and the zero-inflated regression model described in the second section. Parameter estimates and their corresponding p-values are displayed in Tables 6 and 7, respectively. It is observable that for the former regression model, the explanatory variables private 1 and risk area are statistically significant at the 5% significance level. This is also verified by the intercept term of the model, i.e., constant. Furthermore, some other covariates (private 2, profession and has km limit) are significant at the same level for all the responses except for Fire and Theft; similarly, driver age is not significant for the dependent variable Responsible. On the other hand, with respect to the zero-inflated regression model, the explanatory variables risk area and constant are statistically significant at the same significance level for all responses; moreover, the covariate private 1 is not significant for the response variable Parking and has km limit for Fire and Theft. In terms of the four measures of model selection considered, the zero-inflated regression model is preferable over the model (12). Now, we are interested in comparing the six mixed random variables' aggregate claims amount for Claims (Y 1 ) and the different coverages, i.e., Nonresponsible, Responsible, Parking, Windscreen, and Fire and Theft, (Y 2 , . . . , Y 6 ). In order to estimate the scale parameter σ of the exponential distribution, we considered the variable ClaimAmount available in the dataset freMPL10. The estimate of this parameter is σ = 1.96265. The pdf/pmf associated with the mixed random variables is displayed in Figure 1. As expected, the density of the random variable Claims fades away to zero slower than the random variables of the different coverages. Among the different coverages, the Responsible variable is the one that approaches zero faster compared to the other coverages. Table 6. Parameter estimates and p-values for the regression model with the density function (12). Some measures of the model selection are also exhibited.  Figure 1. pdf/pmf associated with the mixed random variables of the aggregate claims amount for Claims (Y 1 ) and the different coverages (Y 2 , . . . , Y 6 ) obtained from the estimated value of the mean of the claims size, σ = 1.96265.

Final Comments and Future Research
It is common today, mainly due to the existing competition, that insurers offer coverage of different claims within the same product not only to gain in competitiveness, but also to benefit from risk diversification and volatility. Up to date, most insurance companies differentiate, apart from the total number of claims, individualized claims for different coverages, such as windscreen claims and thefts and fire claims, among others. Therefore, it seems reasonable to assume that every policyholder generates a sequence of claims numbers for each coverage; one of them is the total claims number, which includes the sum of the claims in all the coverages. In this work, we introduced a new methodology based on a multivariate discrete distribution via conditional specification to explain the claims frequency in different coverages and the total claims number. This procedure allows us to analyze the dependence structure of a finite number of coverages in motor vehicle insurance and also to include heterogeneity via explanatory variables. Closed-form expressions were given for model parameter estimates, and a priori premiums were provided when different premiums principles were used. Numerical applications revealed that specific covariates are statistically significant in some coverages, yet they are not so for others. In this way, it allows us to discern how the different explanatory variables affect each coverage when calculating the corresponding premiums.
The approach introduced in this work avoids the use of copula-based modeling. The latter methodology has been very useful, but at the same time, very criticized in the statistical literature when modeling multivariate data. Although there exists a wide catalog of copulas, it has been mentioned that a weakness of the copula approach is in choosing an appropriate copula structure for the model at hand (Balakrishnan and Lai 2009, chp. 1, p. 59). Furthermore, any copula includes a parameter that controls the dependence structure, and this parameter is sometimes difficult to estimate since it must fall into the admissible support. As explored in the second section of this work, the model depends extremely on the parameter Θ 1 , and for that reason, a more flexible dependence structure based on multivariate subordination is an issue that deserves to be studied. In this regard, using this approach would be interesting to compare this family of distributions with the multivariate regression model based on the multivariate Sarmanov distribution, similar to the models derived in Bolancé and Vernic (2019). This model could be used to explain situations where the policyholder wishes to extend the third-party motor vehicle insurance to account for different coverages that adapt to their personal needs. Alternatively, it could be feasible to implement a multivariate version with elliptical copula-based models to accommodate a wide range of dependence. It is essential to mention that the properties of the copula are not the same as for continuous random variables since the probability of ties in the data is positive. Thus, the estimation cannot be directly carried out, and a continuous extension of integer-valued random variables is needed by using the approach proposed by Denuit and Lamber (2005).
The purpose of the work is not to compare other models, as models of this nature are not known to our knowledge in the actuarial literature. However, the cases with two coverages were discussed via the bivariate Poisson case (see Bermúdez and Karlis 2017) and the case with all the coverages using the multivariate negative binomial distribution in (Johnson et al. 1996, chp. 36, p. 94). Obviously, the fit obtained with the proposed modeling does not seem entirely reasonable (as judged by the chi-squared test statistics, which was not shown in the paper). Then, the model could be improved by using a similar model, but assuming that the total number of claims and all the coverages follow a negative binomial distribution instead. It would be also possible to zero-inflate all the different coverages. This issue could be the subject of future research. Acknowledgments: The authors are grateful to the four anonymous referees for their constructive comments and suggestions, which greatly helped us improve the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest. Table A1 displays the empirical joint claims frequency for all types of coverage and the total number of claims. Table A1. Empirical joint claims frequency corresponding to the response variables Claims (N 1 ), Nonresponsible (N 2 ), Responsible (N 3 ), Parking (N 4 ), Windscreen (N 5 ), and Fire and Theft (N 6