Multivariate Frequency-Severity Regression Models in Insurance

: In insurance and related industries including healthcare, it is common to have several outcome measures that the analyst wishes to understand using explanatory variables. For example, in automobile insurance, an accident may result in payments for damage to one’s own vehicle, damage to another party’s vehicle, or personal injury. It is also common to be interested in the frequency of accidents in addition to the severity of the claim amounts. This paper synthesizes and extends the literature on multivariate frequency-severity regression modeling with a focus on insurance industry applications. Regression models for understanding the distribution of each outcome continue to be developed yet there now exists a solid body of literature for the marginal outcomes. This paper contributes to this body of literature by focusing on the use of a copula for modeling the dependence among these outcomes; a major advantage of this tool is that it preserves the body of work established for marginal models. We illustrate this approach using data from the Wisconsin Local Government Property Insurance Fund. This fund offers insurance protection for (i) property; (ii) motor vehicle; and (iii) contractors’ equipment claims. In addition to several claim types and frequency-severity components, outcomes can be further categorized by time and space, requiring complex dependency modeling. We ﬁnd signiﬁcant dependencies for these data; speciﬁcally, we ﬁnd that dependencies among lines are stronger than the dependencies between the frequency and average severity within each line.


Introduction and Motivation
Many insurance data sets feature information about how often claims arise, the frequency, in addition to the claim size, the severity.Observable responses can include: • N, the number of claims (events), • y k , k = 1, ..., N, the amount of each claim (loss), and By convention, the set {y j } is empty when N = 0.
Importance of Modeling Frequency.The aggregate claim amount S is the key element for an insurer's balance sheet, as it represents the amount of money paid on claims.So, why do insurance companies regularly track the frequency of claims as well as the claim amounts?As in an earlier review [1], we can segment these reasons into four categories: (i) features of contracts; (ii) policyholder behavior and risk mitigation; (iii) databases that insurers maintain; and (iv) regulatory requirements.
1. Contractually, it is common for insurers to impose deductibles and policy limits on a per occurrence and on a per contract basis.Knowing only the aggregate claim amount for each policy limits any insights one can get into the impact of these contract features.2. Covariates that help explain insurance outcomes can differ dramatically between frequency and severity.For example, in healthcare, the decision to utilize healthcare by individuals (the frequency) is related primarily to personal characteristics whereas the cost per user (the severity) may be more related to characteristics of the healthcare provider (such as the physician).
Covariates may also be used to represent risk mitigation activities whose impact varies by frequency and severity.For example, in fire insurance, lightning rods help to prevent an accident (frequency) whereas fire extinguishers help to reduce the impact of damage (severity).3.Many insurers keep data files that suggest developing separate frequency and severity models.
For example, insurers maintain a "policyholder" file that is established when a policy is written.
A separate file, often known as the "claims" file, records details of the claim against the insurer, including the amount.These separate databases facilitate separate modeling of frequency and severity.4. Insurance is a closely monitored industry sector.Regulators routinely require the reporting of claims numbers as well as amounts.Moreover, insurers often utilize different administrative systems for handling small, frequently occurring, reimbursable losses, e.g., prescription drugs, versus rare occurrence, high impact events, e.g., inland marine.Every insurance claim means that the insurer incurs additional expenses suggesting that claims frequency is an important determinant of expenses.
Importance of Including Covariates.In this work, we assume that the interest is in the joint modeling of frequency and severity of claims.In actuarial science, there is a long history of studying frequency, severity and the aggregate claim for homogeneous portfolios; that is, identically and independently distributed realizations of random variables.See any introductory actuarial text, such as [2], for an introduction to this rich literature.
In contrast, the focus of this review is to assume that explanatory variables (covariates, predictors) are available to the analyst.Historically, this additional information has been available from a policyholder's application form, where various characteristics of the policyholder were supplied to the insurer.For example, in motor vehicle insurance, classic rating variables include the age and sex of the driver, type of the vehicle, region in which the vehicle was driven, and so forth.The current industry trend is towards taking advantage of "big data", with attempts being made to capture additional information about policyholders not available from traditional underwriting sources.An important example is the inclusion of personal credit scores, developed and used in the industry to assess the quality of personal loans, that turn out to also be important predictors of motor vehicle claims experience.Moreover, many insurers are now experimenting with global positioning systems combined with wireless communication to yield real-time policyholder usage data and much more.Through such systems, they gather micro data such as the time of day that the car is driven, sudden changes in acceleration, and so forth.This foray into detailed information is known as "telematics".See, for example, [3] for further discussion.
Importance of Multivariate Modeling.To summarize reasons for examining insurance outcomes on a multivariate basis, we utilize an earlier review in [4].In that paper, frequencies were restricted to binary outcomes, corresponding to a claim or no claim, known as "two-part" modeling.In contrast, this paper describes more general frequency modeling, although the motivation for examining multivariate outcomes are similar.Analysts and managers gain useful insights by studying the joint behavior of insurance risks, i.e., a multivariate approach:

•
For some products, insurers must track payments separately by component to meet contractual obligations.For example, in motor vehicle coverage, deductibles and limits depend on the coverage type, e.g., bodily injury, damage to one's own vehicle, or damage to another party; is natural for the insurer to track claims by coverage type.

•
For other products, there may be no contractual reasons to decompose an obligation by components and yet the insurer does so to help better understand the overall risk.For example, many insurers interested in pricing homeowners insurance are now decomposing the risk by "peril", or cause of loss.Homeowners is typically sold as an all-risk policy, which covers all causes of loss except those specifically excluded.By decomposing losses into homogenous categories of risk, actuaries seek to get a better understanding of the determinants of each component, resulting in a better overall predictor of losses.

•
It is natural to follow the experience of a policyholder over time, resulting in a vector of observations for each policyholder.This special case of multivariate analysis is known as "panel data", see, for example, [5].

•
In the same fashion, policy experience can be organized through other hierarchies.For example, it is common to organize experience geographically and analyze spatial relationships.

•
Multivariate models in insurance need not be restricted to only insurance losses.For example, a study of term and whole life insurance ownership is in [6].As an example in customer retention, both [7,8] advocate for putting the customer at the center of the analysis, meaning that we need to think about the several products that a customer owns simultaneously.
An insurer has a collection of multivariate risks and the interest is managing the distribution of outcomes.Typically, insurers have a collection of tools that can then be used for portfolio management including deductibles, coinsurance, policy limits, renewal underwriting, and reinsurance arrangements.Although pricing of risks can often focus on the mean, with allowances for expenses, profit, and "risk loadings", understanding capital requirements and firm solvency requires understanding of the portfolio distribution.For this purpose, it is important to treat risks as multivariate in order to get an accurate picture of their dependencies.

Dependence and Contagion.
We have seen in the above discussion that dependencies arise naturally when modeling insurance data.As a first approximation, we typically think about risks in a portfolio as being independent from one another and rely upon risk pooling to diversify portfolio risk.However, in some cases, risks share common elements such as an epidemic in a population, a natural disaster such as a hurricane that affects many policyholders simultaneously, or an interest rate environment shared by policies with investment elements.These common (pandemic) elements, often known as "contagion", induce dependencies that can affect a portfolio's distribution significantly.
Thus, one approach is to model risks as univariate outcomes but to incorporate dependencies through unobserved "latent" risk factors that are common to risks within a portfolio.This approach is viable in some applications of interest.However, one can also incorporate contagion effects into a more general multivariate approach that we adopt this view in this paper.We will also consider situations where data are available to identify models and so we will be able to use the data to guide our decisions when formulating dependence models.
Modeling dependencies is important for many reasons.These include: 1. Dependencies may impact the statistical significance of parameter estimates.
2. When we examine the distribution of one variable conditional on another, dependencies are important.3.For prediction, the degree of dependency affects the degree of reliability of our predictions.4. Insurers want to construct products that do not expose them to extreme variation.They want to understand the distribution of a product that has many identifiable components; to understand the distribution of the overall product, one strategy is to describe the distribution of each product and a relationship among the distributions.
A recent review paper [3] provides additional discussion.
Plan for the Paper.The following is a plan to introduce readers further to the topic.Section 2 gives a brief overview of univariate models, that is, regression models with a single outcome for a response.This section sets the tone and notation for the rest of the paper.Section 3 provides an overview of multivariate modeling, focusing on the "copula" regression approach described here.This section discusses continuous, discrete, and mixed (Tweedie) outcomes.For our regression applications, the focus is mainly on a family of copulas known as "elliptical", because of their flexibility of modeling pairwise dependence and wide usages in multivariate analysis.Section 3 also summarizes a modeling strategy for the empirical approach of copula regression.Section 4 reviews other recent work on multivariate frequency-severity model and describes the benefits of diversification, particularly important in an insurance context.To illustrate our ideas and approach, Sections 5 and 6 provide our analysis using data from the Wisconsin Local Government Property Insurance Fund.Section 7 concludes with a few closing remarks.

Univariate Foundations
For notation, define N for the random number of claims, S for the aggregate claim amount, and S = S/N for the average claim amount (defined to be 0 when N = 0).To model these outcomes, we use a collection of covariates x, some of which may be useful for frequency modeling whereas others will be useful for severity modeling.The dependent variables N, S, and S as well as covariates x vary by the risk i = 1, . . ., n.For each risk, we also are interested in multivariate outcomes indexed by j = 1, . . ., p. So, for example, N i = (N i1 , . . ., N ip ) represents the vector of p claim outcomes from the ith risk.
This section summarizes modeling approaches for a single outcome (p = 1).A more detailed review can be found in [1].

Frequency-Severity
For modeling the joint outcome (N, S) (or equivalently, (N, S)), it is customary to first condition on the frequency and then modeling the severity.Suppressing the {i} subscript, we decompose the distribution of the dependent variables as: joint = frequency × conditional severity, where f(N, S) denotes the joint distribution of (N, S).Through this decomposition, we do not require independence of the frequency and severity components.
There are many ways to model dependence when considering the joint distribution f(N, S) in Equation (1).For example, one may use a latent variable that affects both frequency N and loss amounts S, thus inducing a positive association.Copulas are another tool used regularly by actuaries to model non-linear associations and will be described in subsequent Section 4. The conditional probability framework is a natural method of allowing for potential dependencies and provides a good starting platform for empirical work.

Modeling Frequency Using GLMs
It has become routine for actuarial analysts to model the frequency N i based on covariates x i using generalized linear models, GLMs, cf., [9].For binary outcomes, logit and probit forms are most commonly used, cf., [10].For count outcomes, one begins with a Poisson or negative binomial distribution.Moreover, to handle the excessive number of zeros relative to that implied by these distributions, analysts routinely examine zero-inflated models, as described in [11].
A strength of GLMs relative to other non-linear models is that one can express the mean as a simple function of linear combinations of the covariates.In insurance, it is common to use a "logarithmic link" for this function and so express the mean as µ i = E N i = exp(x i β), where β is a vector of parameters associated with the covariates.This function is used because it yields desirable parameter interpretations, seems to fit data reasonably well, and ties well with other approaches traditionally used in actuarial ratemaking applications [12].
It is also common to identify one of the covariates as an "exposure" that is used to calibrate the size of a potential outcome variable.In frequency modeling, the mean is assumed to vary proportionally with E i , for exposure.To incorporate exposures, we specify one of the explanatory variables to be ln E i and restrict the corresponding regression coefficient to be 1; this term is known as an offset.With this convention, we have ln Since there are inflated numbers of 0 s and 1 s in our data, a "zero-one-inflated" model is introduced.As an extension of the zero-inflated method, a zero-one-inflated model employs two generating processes.The first process is governed by a multinomial distribution that generates structural zeros and ones.The second process is governed by a Poisson or negative binomial distribution that generates counts, some of which may be zero or one.
Denote the latent variable in the first process as I i , i = 1, . . ., n, which follows a multinomial distribution with possible values 0, 1 and 2 with corresponding probabilities Here, N i is frequency.
Here, P i may be a Poisson or negative binomial distribution.With this, the probability mass function of A logit specification is used to parameterize the probabilities for the latent variable I i .Denote the covariates associated with I i as z i .A logit specification is used to parameterize the probabilities for the latent variable I i .Using level 2 as a reference, the specification is Correspondingly, Maximum likelihood estimation is used to fit the parameters.

Modeling Severity
Modeling Severity Using GLMs.For insurance analysts, one strength of the GLM approach is that the same set of routines can be used for continuous as well as discrete outcomes.For severities, it is common to use a gamma or inverse Gaussian distribution, often with a logarithmic link (primarily for parameter interpretability).One strength of the linear exponential family that forms the basis of GLMs is that a sample average of outcomes comes from the same distribution as the outcomes.Specifically, suppose that we have m independent variables from the same distribution with location parameter θ and scale parameter φ.Then, the sample average comes from the same distributional family with location parameter θ and scale parameter φ/m.This result is helpful as insurance analysts regularly face grouped data as well as individual data.For example, [1] provides a demonstration of this basic property.
To illustrate, in the aggregate claims model, if individual losses have a gamma distribution with mean µ i and scale parameter φ, then, conditional on observing N i losses, the average aggregate loss Si has a gamma distribution with mean µ i and scale parameter φ/N i .
Modeling Severity Using GB2.The GLM is the workhorse for industry analysts interested in analyzing the severity of claims.Naturally, because of the importance of claims severity, a number of alternative approaches have been explored, cf., [13] for an introduction.In this review, we focus on a specific alternative, using a distribution family known as the "generalized beta of the second kind", or GB2, for short.
A random variable with a GB2 distribution can be written as where the constant C 1 = (α 1 /α 2 ) σ , G 1 and G 2 are independent gamma random variables with scale parameter 1 and shape parameters α 1 and α 2 , respectively.Further, the random variable F has an F-distribution with degrees of freedom 2α 1 and 2α 2 , and the random variable Z has a beta distribution with parameters α 1 and α 2 .
Thus, the GB2 family has four parameters (α 1 , α 2 , µ and σ), where µ is the location parameter.Including limiting distributions, the GB2 encompasses the "generalized gamma" (by allowing α 2 → ∞) and hence the exponential, Weibull, and so forth.It also encompasses the "Burr Type 12" (by allowing α 1 = 1), as well as other families of interest, including the Pareto distributions.The GB2 is a flexible distribution that accommodates positive or negative skewness, as well as heavy tails.See, for example, [2] for an introduction to these distributions.
For incorporating covariates, it is straightforward to show that the regression function is of the form where the constant C 2 can be calculated with other (non-location) model parameters.Under the most commonly used way of parametrization for the GB2, where µ is associated with covariates, if −α 1 < σ < α 2 , then we have Thus, one can interpret the regression coefficients in terms of a proportional change.That is, In principle, one could allow for any distribution parameter to be a function of the covariates.However, following this principle would lead to a large number of parameters; this typically yields computational difficulties as well as problems of interpretations, [14].In this paper, µ is used to incorporate covariates.An alternative parametrization, as described in Appendix A.1, is introduced as an extension of the GLM framework.

Tweedie Model
Frequency-severity modeling is widely used in insurance applications.However, for simplicity, it is also common to use only the aggregate loss S as a dependent variable in a regression.Because the distribution of S typically contains a positive mass at zero representing no claims, and a continuous component for positive values representing the amount of a claim, a widely used mixture is the Tweedie (1984) distribution.The Tweedie distribution is defined as a Poisson sum of gamma random variables.Specifically, suppose that N has a Poisson distribution with mean λ, representing the number of claims.Let y j be an i.i.d.sequence, independent of N, with each y j having a gamma distribution with parameters α and β, representing the amount of a claim.Note, β is standard notation for this parameter used in loss-model textbooks, and the reader should understand it is different from the bold-faced β, as the latter is a symbol we will use for the coefficients corresponding to explanatory variables.Then, S = y 1 + . . .+ y N is a Poisson sum of gammas.
To understand the mixture aspect of the Tweedie distribution, first note that it is straightforward to compute the probability of zero claims as Pr(S = 0) = Pr(N = 0) = e −λ .The distribution function can be computed using conditional expectations, Because the sum of i.i.d.gammas is a gamma, S n = y 1 + . . .+ y n (not S) has a gamma distribution with parameters nα and β.For y > 0, the density of the Tweedie distribution is From this, straight-forward calculations show that the Tweedie distribution is a member of the linear exponential family.Now, define a new set of parameters µ, φ, P through the relations Easy calculations show that E S = µ and Var S = φµ P , where 1 < P < 2. The Tweedie distribution can also be viewed as a choice that is intermediate between the Poisson and the gamma distributions.
In the basic form of the Tweedie regression model, the scale (or dispersion) parameter φ is constant.However, if one begins with the frequency-severity structure, calculations show that φ depends on the risk characteristics (i, cf., [1]).Because of this and the varying dispersion (heteroscedasticity) displayed by many data sets, researchers have devised ways of accommodating and/or estimating this structure.The most common way is the so-called "double GLM" procedure proposed in [15] that models the dispersion as a known function of a linear combination of covariates (as for the mean, hence the name "double GLM").

Copula Regression
Copulas have been applied with GLMs in the biomedical literature since the mid-1990s ( [16][17][18]).In the actuarial literature, the t-copula and the Gaussian copula with GLMs as marginal distributions were used to develop credibility predictions in [19].In more general cases, [20] provides a detailed introduction of copula regression that focuses on the Gaussian copula and [21] surveys copula regression applications.
Introducing Copulas.Specifically, a copula is a multivariate distribution function with uniform marginals.Let U 1 , . . ., U p be p uniform random variables on (0, 1).Their joint distribution function is a copula.
Of course, we seek to use copulas in applications that are based on more than just uniformly distributed data.Thus, consider arbitrary marginal distribution functions F 1 (y 1 ), . . ., F p (y p ).Then, we can define a multivariate distribution function using the copula such that F(y 1 , . . . ,y p ) = C(F 1 (y 1 ), . . ., F p (y p )). ( If outcomes are continuous, then we can differentiate the distribution functions and write the density function as f(y 1 , . . . ,y p ) = c(F 1 (y 1 ), . . ., F p (y p )) where f j is the density of the marginal distribution F j and c is the copula density function.
It is easy to check from the construction in Equation ( 3) that F(•) is a multivariate distribution function.Sklar established the converse in [22].He showed that any multivariate distribution function F(•) can be written in the form of Equation ( 3), that is, using a copula representation.Sklar also showed that, if the marginal distributions are continuous, then there is a unique copula representation.See, for example, the introductory book to copulas [23,24] for an introduction to copulas from an insurance perspective and [25] for a comprehensive modern treatment.

Regression with
Copulas.In a regression context, we assume that there are covariates x associated with outcomes y = (y 1 , . . ., y p ) .In a parametric context, we can incorporate covariates by allowing them to be functions of the distributional parameters.
Specifically, we assume that there are n independent risks and p outcomes for each risk i, i = 1, . . ., n.For this section, consider an outcome y i = (y i1 , . . ., y ip ) and K × 1 vector of covariates x i , where K is the number of covariates.The marginal distribution of y ij is a function of x ij , β j and θ j .Here, x ij is a K j × 1 vector of explanatory variables for risk i and outcome type j, a subset of x i , and β j is a K j × 1 vector of marginal parameters to be estimated.The systematic component x ij β j determines the location parameter.The vector θ j summarizes additional parameters of the marginal distribution that determine the scale and shape.Let F ij = F(y ij ; x ij , β j , θ j ) denote the marginal distribution function.
This describes a classical approach to regression modeling, treating explanatory variances/ covariates as non-stochastic ("fixed") variables.An alternative is to think of the covariates themselves as random and perform statistical inference conditional on them.Some advantages of this alternative approach are that one can model the time-changing behavior of covariates, as in [26], or investigate non-parametric alternatives, as in [27].These represent excellent future steps in copula regression modeling that are not addressed further in this article.

Multivariate Severity
For continuous severity outcomes, we may consider the density function f ij = f(y ij ; x ij , β j , θ j ) associated with the distribution function F ij and c the copula density function with parameter vector α.With this, using Equation ( 4), the log-likelihood function of the ith risk is written as where β = (β 1 , . . ., β p ) and θ = (θ 1 , . . ., θ p ) are collections of parameters over the p outcomes.This is a fully parametric set-up; the usual maximum likelihood techniques enjoy certain optimality properties and are the preferred estimation method.
If we consider only a single outcome, say y i1 , then the associated log-likelihood is ln f i1 .Thus, the set of outcomes y 11 , . . ., y n1 allows for the usual "root-n" consistent estimator of β 1 , and similarly for the other outcomes y ij , j = 2 . . ., p.By considering each outcome in isolation of the others, we can get desirable estimators of the regression coefficients β j , j = 1, . . ., p.These provide excellent starting values to calculate the fully efficient maximum likelihood estimators using the log-likelihood from Equation ( 5).Joe coined the phrase "inference for margins", sometimes known by the acronym IFM in [28], to describe this approach to estimation.
In the same way, one can consider any pair of outcomes, (y ij , y ik ) for j = k.This permits consistent estimation of the marginal regression parameters as well as the association parameters between the jth and kth outcomes.As with the IFM, this technique provides excellent starting values of a fully efficient maximum likelihood estimation recursion.Moreover, they provide the basis for an alternative estimation method known as "composite likelihood", cf., [29] or [30], for a description in a copula regression context.

Multivariate Frequency
If outcomes are discrete, then one can take differences of the distribution function in Equation ( 3) to write the probability mass function as Here, u j,1 = F j (y j −) and u j,2 = F j (y j ) are the left-and right-hand limits of F j at y j , respectively.For example, when p = 2, we have It is straightforward in principle to estimate parameters using Equation ( 6) and standard maximum likelihood theory.
In practice, two caveats should be mentioned.The first is that the result of [22] only guarantees that the copula is unique over the range of the outcomes, a point emphasized with several interesting examples in [31].In a regression context, this non-identifiability is less likely to be a concern, as noted in [25,29,30,32].Moreover, the latter reference emphasizes that the Gaussian copula with binary data has been used for decades by researchers as this is just another form for the commonly used multivariate probit.
The second issue is computational.As can be seen in Equation ( 6), likelihood inference involves the computation of multidimensional rectangle probabilities.The review article [30] describes several variations of maximum likelihood that can be useful as the dimension p increases, see also [33].
As in [34], the composite likelihood method is used for computation.For large values of p, the pair (also known as "vine") copula approach described in [35] for discrete outcomes seems to be a promising approach.

Multivariate Tweedie
As emphasized in [25] (p.226), in copula regression it is possible to have outcomes that are combinations of continuous, discrete, and mixture distributions.One case of special interest in insurance modeling is the multivariate Tweedie, where each marginal distribution is a Tweedie and the margins are joined by a copula.Specifically, Shi considers different types of insurance coverages with Tweedie margins in [36] .
To illustrate the general principles, consider the bivariate case (p = 2).Suppressing the i index and covariate notation, the joint distribution is where ∂ j C denotes the partial derivative of copula with respect to jth component.
See [36] for additional details of this estimation where he also described a double GLM approach to accommodate varying dispersion parameters.

Association Structures and Elliptical Copulas
First consider an outline of the evolution of multivariate regression modeling.
1.The multivariate normal (Gaussian) distribution has provided a foundation for multivariate data analysis, including regression.By permitting a flexible structure for the mean, one can readily incorporate complex mean structures including high order polynomials, categorical variables, interactions, semi-parametric additive structures, and so forth.Moreover, the variance structure readily permits incorporating time series patterns in panel data, variance components in longitudinal data, spatial patterns, and so forth.One way to get a feel for the breadth of variance structures readily accommodated is to examine options in standard statistical software packages such as PROC Mixed in [37] (for example, the TYPE switch in the RANDOM statement permits the choice of over 40 variance patterns).2. In many applications, appropriately modeling the mean and second moment structure (variances and covariances) suffices.However, for other applications, it is important to recognize the underlying outcome distribution and this is where copulas come into play.As we have seen, copulas are available for any distribution function and thus readily accommodate binary, count, and long-tail distributions that cannot be adequately approximated with a normal distribution.Moreover, marginal distributions need not be the same, e.g., the first outcome may be a count Poisson distribution and the second may be a long-tail gamma.3. Pair copulas (cf., [25]) may well represent the next step in the evolution of regression modeling.
A copula imposes the same dependence structure on all p outcomes whereas a pair copula has the flexibility to allow the dependence structure itself to vary in a disciplined way.This is done by focusing on the relationship between pairs of outcomes and examining conditional structures to form the dependence of the entire vector of outcomes.This approach is useful for high dimensional outcomes (where p is large), an important developing area of statistics.This represents an excellent future step in copula regression modeling that is not addressed further in this article.
As described in [25], there is a host of copulas available depending on the interests of the analyst and the scientific purpose of the investigation.Considerations for the choice of a copula may include computational convenience, interpretability of coefficients, a latent structure for interpretability, and a wide range of dependence, allowing both positive and negative associations.
For our applications of regression modeling, we typically begin with the elliptical copula family.This family is based on the family of elliptical distributions that includes the multivariate normal and t-distributions, see more in [38].
This family has most of the desirable traits that one would seek in a copula family.From our perspective, the most important feature is that it permits the same family association matrices found in the multivariate Gaussian distribution.This not only allows the analyst to investigate a wide degree of association patterns, but also allows estimation to be accomplished in a familiar way using the same structure as in the Gaussian family, e.g., [37].
For example, if the ith risk evolves over time, we might use a familiar time series model to represent associations, e.g., such as an autoregressive of order 1 (AR1).See [19] for an actuarial application.For a more complex example, suppose that y i = (y i1 , y i2 , y i3 ) represents three types of expenses for the ith company observed over 4 time periods.Then, we might use the following dependence structure as in [39].This is a commonly used specification in models of several time series in econometrics where σ jk represents a cross-sectional association between y ij and y ik and Σ jk represents cross-associations with time lags.

Assessing Dependence
Dependence can be assessed at all stages of the model fitting process: • Copula identification begins after marginal models have been fit.Then, use the "Cox-Snell" residuals from these models to check for association.Create simple correlation statistics (Spearman, polychoric) as well as plots (pp and tail dependence plots) to look for dependence structures and identify a parametric copula.

•
After a model identification, estimate the model and examine how well the model fits.Examine the residuals to search for additional patterns using, for example, correlation statistics and t-plot (for elliptical copulas).Examine the statistical significance of fitted association parameters to seek a simpler fit that captures the important tendencies of the data.

•
Compare the fitted model to alternatives.Use overall goodness of fit statistics for comparisons, including AIC and BIC, as well as cross-validation techniques.For nested models, compare via the likelihood ratio test and use Vuong's procedure for comparing non-nested alternative specifications.

•
Compare the models based on a held-out sample.Use statistical measures and economically meaningful alternatives such as the Gini statistic.
In the first and second step of this process, a variety of hypothesis tests and graphical methods can be employed to identify the specific type of copula (e.g., Archmedian, elliptical, extreme value, and so forth) that corresponds to the given data.Researchers have developed a graphical tool called the Kendall plot, or the K-plot for short, to detect dependence.See [40].To determine whether a joint distribution corresponds to an Archimedean copula or a specific extreme-value copula, goodness-of-fit tests developed by [41][42][43][44] can be helpful.The reader may also reference [25] for a comprehensive coverage of the various assessment methods for dependence.

Frequency-Severity Modeling Strategy
Multivariate frequency-severity modeling strategies are a subset of the usual regression and copula identification and inference strategies.In absence of a compelling theory to suggest the appropriate covariates and predictors (which is generally the case for insurance applications), the modeling strategy consists of model identification, estimation, and inference.Typically, this is done in a recursive fashion where one sets aside a random portion of the data for identification and estimation (the "training" sample), and one proceeds to validate and conduct inference on another portion (the "test" sample).See, for example, [45] for a description of this and many procedures for variable selection, mainly in a cross-sectional regression context.
You can think about the identification and estimation procedures as three components in a copula regression model: 1. Fit the mean structure.Historically, this is the most important aspect.One can apply robust standard error procedures to get consistent and approximately normally distributed coefficients, assuming a correct mean structure.2. Fit the variance structure with a selected distribution.In GLMs, the choice of the distribution dictates the variance structure that can be over-ruled with a separately specified variance, e.g., a "double GLM." 3. Fit the dependence structure with a choice of copula.
For frequency-severity modeling, there are two mean and variance structures to work with, one for the frequency and one for the severity.

Identification and Estimation
Although the estimation of parametric copulas is fairly established, the literature on identification of copulas is still in the early stage of development.As described in Section 3.2, maximum likelihood is the usual choice with an inference for margins and/or composite likelihood approach for starting values of the iterations.A description of composite likelihood in the context of copula modeling can be found in [25].As noted here, composite likelihood may be particularly useful for multivariate discrete data when univariate margins have common parameters (p.233, [25]).Another variation of maximum likelihood in copula regression is the "maximization by parts" method, as described in [20] and utilized in [46].In the context of copula regressions, the idea behind this is to split the likelihood into two pieces, an easier part corresponding to the marginals and a more difficult part corresponding to the copula.The estimation routine takes advantage of these differing levels of difficulty in the calculations.
Identification of copula models used in regression typically starts with residuals from marginal fits.For severities, the idea is to estimate a parametric fit to the marginal distribution, such as normal or gamma regression.Then, one applies the distribution function (that depends on covariates) to the observation.Using notation, we can write this as F i (y i ) = ˆ i .This is known as the "probability integral transformation."If the model is correctly specified, then the ˆ i has a uniform (0,1) distribution.This is an idea that dates back to works by Cox and Snell in [47] and so these are often known as "Cox-Snell" residuals.For copula identification, it is recommended in [25] to take an inverse normal distribution transform (i.e., Φ −1 ( ˆ i ), for a standard normal distribution function Φ) to produce "normal scores." Because of the discreteness with frequencies, these residuals are not uniformly distributed even if the model is correctly specified.In this case, one can "jitter" the residuals.Specifically, define a modified distribution function F i (y, λ) = Pr(Y i < y) + λ Pr(Y i = y) and let V be a uniform random number that is independent of Y i .Then, we can define the jittered residual to be F i (y i , V) = ˜ i .If the model is correctly specified, then jittered residuals have a uniform (0,1) distribution, cf., [48].
Compared to classical residuals, residuals from probability integral transforms have less ability to guide model development-we can only tell if the marginal models are approximately correct.The main advantage of this residual is that it is applicable to all (parametric) distributions.If you are working with a distribution that supports other definitions of residuals, then these are likely to be more useful because they may tell you how to improve your model specification, not whether or not it is approximately correct.If the marginal model fit is adequate, then we can think of the residuals as approximate realizations from a uniform distribution and use standard techniques from copula theory to identify a copula.We refer to [25] for a summary of this literature.

Model Validation
After identification and estimation, it is customary to compare a number of alternative models based on the training and on the test samples.For the training sample, the "in-sample" comparisons are typically based on the significance of the coefficients, overall goodness of fit measures (including information criteria such as AIC and BIC), cross-validation, as well as likelihood ratio comparisons for nested models.
For comparisons among non-nested parametric models, it is now common in the literature to cite a statistic due to Vuong [49].For this statistic, one calculates the contribution to the logarithmic likelihood such as in Equation ( 5) for two models, say, l (1) i and l (2) and m is the size of the validation sample.To assess the significance of this difference, one can apply approximate normality with approximate standard errors given as SD D / √ m where 2 .In a copula context, see (p. 257, [25]) for a detailed description of this procedure, where sample size adjustments similar to those used in AIC and BIC are also introduced.
Comparison among models using test data, or "out-of-sample" comparisons are also important in insurance because many of these models are used for predictive purposes such as setting rates for new customers.Out-of-sample measures compare held-out observations to those predicted by the model.Traditionally, absolute values and squared differences have been used to summarize differences between these two.However, for many insurance data sets, there are large masses at zero, meaning that these traditional metrics are less helpful.To address this problem, a newer measure is developed in [50] that they call the "Gini index."In this context, the Gini index is twice the average covariance between the predicted outcome and the rank of the predictor.In order to compare models, Theorem 5 of [50] provides standard errors for the difference of two Gini indices.

Frequency Severity Dependency Models
In traditional models of insurance data, the claim frequency is assumed to be independent of claim severity.We emphasize in Appendix A.4 that the average severity may depend on frequency, even when this classical assumption holds.
One way of modeling the dependence is through the conditioning argument developed in Section 2.1.An advantage of this approach is that the frequency can be used as a covariate to model the average severity.See [51] for a healthcare application of this approach.For another application, a Bayesian approach for modeling claim frequency and size was proposed in [52], with both covariates as well as spatial random effects taken into account.The frequency was incorporated into the severity model as covariate.In addition, they checked both individual and average claim modeling and found the results were similar in their application.
As an alternative approach, copulas are widely used for frequency severity dependence modeling.In [46], Czado et al. fit Gaussian copula on Poisson frequency and gamma severity and used an optimization by parts method from [53] to do the estimation.They derived the conditional distribution of frequency given severity.In [54], the distribution of policy loss is derived without the independence assumption between frequency and severity.They also showed that the ignoring of dependence can lead to underestimation of loss.A Vuong test was adopted to select the copula.
To see how the copula approach works, recall that S represents average severity of claims and N denotes frequency.Using a copula, we can express the likelihood as With this, This yields the following expression for the likelihood For another approach, Shi et al. also built a dependence model between frequency and severity in [55].They used an extra indicator variable for occurrence of claim to deal with the zero-inflated part, and built a dependence model between frequency and severity conditional on positive claim.The two approaches described previously were compared; one approach using frequency as a covariate for the severity model, and the other using copulas.They used a zero-truncated negative binomial for positive frequency and the GG model for severity.In [56], a mixed copula regression based on GGS copula (see [25] for an explanation of this copula) was applied on a medical expenditure panel survey (MEPS) dataset.In this way, the negative tail dependence between frequency and average severity can be captured.
Brechmann et al. applied the idea of the dependence between frequency and severity to the modeling of losses from operational risks in [57].For each risk class, they considered the dependence between aggregate loss and the presence of loss.Another application of this methodology in operational risk aggregation can be found in [58].Li et al. focused on two dependence models; one for the dependence of frequencies across different business lines, and another for the aggregate losses.They applied the method on Chinese banking data and found significant difference between these two methods.

LGPIF Case Study
We demonstrate the multivariate frequency severity modeling approach using a data set from the Wisconsin Local Government Property Insurance Fund (LGPIF).The LGPIF was established to provide property insurance for local government entities that include counties, cities, towns, villages, school districts, fire departments, and other miscellaneous entities, and is administered by the Wisconsin Office of the Insurance Commissioner.Properties covered under this fund include government buildings, vehicles, and equipment.For example, a county entity may need coverage for its snow plowing trucks, in addition to its building and contents [59].These data provide a good example of a typical multi-line insurance company encountered in practice.More details about the project may be found at the Local Government Property Insurance Fund project website [60].

Data / Problem Description
The data consist of six coverage groups; building and content (BC), contractor's equipment (IM), comprehensive new (PN), comprehensive old (PO), collision new (CN), collision old (CO) coverage.The data are longitudinal, and Tables 1 and 2

BC Building and Contents
This coverage provides insurance for buildings and the properties within.
In case the policyholder has purchased a rider, claims in this group may reflect additional amounts covered under endorsements.

IM
Contractor's Equipment IM, an abbreviation for "inland marine" is used as the coverage code for equipments coverage, which originally belong to contractors.

C Collision
This provides coverage for impact of a vehicle with an object, impact of vehicle with an attached vehicle, or overturn of a vehicle.

P Comprehensive
Direct and accidental loss or damage to motor vehicle, including breakage of glass, loss caused by missiles, falling objects, fire, theft, explosion, earthquake, windstorm, hail, water, flood, malicious mischief or vandalism, riot or civil common, or colliding with a bird or animal.

N New
This code is used as an indication that the coverage is for vehicles of current model year, or 1∼2 years prior to the current model year.
O Old This code is used as an indication that the coverage is for vehicles three or more years prior to the current model year.
From Table 3, there are collision and comprehensive coverages, each for new and old vehicles of the entity.Hence, an entity can potentially have collision coverage for new vehicles (CN), collision coverage for old vehicles (CO), comprehensive coverage for new vehicles (PN), and comprehensive coverage for old vehicles (PO).Hence, in our analysis, we consider these sub-coverages as individual lines of businesses, and work with six separate lines, including building and contents (BC), and contractor's equipment (IM) as separate lines also.
Preliminary dependence measures for discrete claim frequencies and continuous average severities can be obtained using polychoric and polyserial correlations.These dependence measures both assume latent normal variables, whose values fall within the cut-points of the discrete variables.The polychoric correlation is the inferred latent correlation between two ordered categorical variables; the polyserial correlation is the inferred latent correlation between a continuous variable and an ordered categorical variable, cf.[25].
Table 4 shows the polychoric correlation among the frequencies of the six coverage groups.Note that these dependencies in Table 4 are measured before controlling for the effects of explanatory variables on the frequencies.As Table 4 shows, there is evidence of correlation across different lines, however these cross-sectional dependencies may be due to correlations in the exposure amounts or, in other words, the sizes of the entities.The dependence between frequencies and average claim severities is often of interest to modelers.In Appendix A.4 we show that average severity may depend on frequency, even when the classical assumption, independence of frequency and individual severities, holds.Our data are consistent with this result.The diagonal entries of Table 5 show the polyserial correlations between the frequency and severity of each coverage group.According to Table 5, the observed correlation between frequency and severity is small.For the CN line, a positive correlation can be observed although very small (0.032, while the other correlations between frequency and severity are negative).Again, these numbers only provide a rough idea of the dependency.Table 6 shows the Spearman correlation between the average severities, for those observations with at least one positive claim.The correlation among the severities of new and old car comprehensive coverage is high.In summary, these summary statistics show that there are potentially interesting dependencies among the response variables.

Explanatory Variables
Table 7 shows the number of observations available in the data set, for years 2006-2010.Explanatory variables used are summarized in Table 8.The marginal analyses for each line are performed on the subset for which the coverage amounts shown in Table 7 are positive.

Marginal Model Fitting-Zero/One Frequency, GB2 Severity
For each coverage type, a frequency-severity model is fit marginally.

BC (Building and Contents) Frequency Modeling
In the frequency part, we fit several commonly employed count models: Poisson, negative binomial (NB), zero-inflated Poisson (zeroinfPoisson), zero-inflated negative binomial (zeroinflNB).Our data not only exhibit a large mass at 0, as with many other insurance claims data, but also an inflated number of 1 s.For BC, there are 997 policies with 1 claim.This can be compared to the expected number under zero-inflated Poisson, 754, and under the zero-inflated negative binomial, 791.(See Table 9 for details).These zero-inflated models underestimate the point mass at 1 due to the shrinkage to 0. Thus, alternative "zero-one-inflated" models are introduced in Section 2.2.
Table 9 shows the expected count for each frequency value under different models and the empirical values from the data.A Poisson distribution underestimates the zero proportions while zero-inflated and negative binomial models underestimate the proportion of 1 s.The zero-one inflated models do provide the best fits for simultaneously estimating the probability of a zero and a one.
Chi-square goodness of fit statistics can be used to compare different models.Table 10 shows the result.It is calculated depending on Table 9.The zero-one-inflated negative binomial is significantly better than other methods.

BC (Building and Contents) Severity Modeling
In the average severity part, the most commonly used distribution, gamma, is fit and compared with the GB2 model.To do the goodness of fit test, the quantiles of normal Cox-Snell residuals are compared with normal quantiles.
Figure 1 shows the residual plot of severity fitted with gamma and GB2.Clearly, the gamma does not fit well especially in the tail part.

Building and Contents Model Summary
Table 11 shows the coefficients for the fitted marginal models.Here, coefficients of GB2, NB and the zero-one-inflated parts are provided.

Marginal Models for Other Lines
Appendix A.2 provides the model selection and marginal model results for lines other than building and contents.

Copula Identification and Fitting
Dependence is fit at two levels.The first is between frequency and average severity within each line.The second is among different lines.

Frequency Severity Dependence
Vuong's test, as described in Section 3.7.2, is used for copula selection.Specifically, we consider two models M (1) and M (2) , in our example, M (1) is Gaussian copula while M (2) is t copula.Let ∆ 12 be the difference in divergence from models M (1) and M (2) .When the true density is g, this can be written as A large sample 95% confidence interval for ∆ 12 , D ± 1.96 × n −1/2 SD D , is provided in Table 12.Table 12 shows the comparison of Gaussian copula against t copula with commonly used degrees of freedom for frequency and severity dependence in BC line.An interval completely below 0 indicates that copula 1 is significantly better than copula 2. Thus, the Gaussian copula is preferred.Maximum likelihood estimation with the full multivariate likelihood, which estimates parameters in marginal and copula models simultaneously, is fit here.Table 13 shows parameters of BC line with the full likelihood method.Here the marginal dispersion parameters are fixed from marginal models.By comparing Tables 11 and 13, it can be seen that the coefficients are close.As pointed out in [25], inference functions for margins, with the results in Table 11, is efficient and can provide a good starting point for the full likelihood method, as in Table 13.For other lines, the results of the full likelihood method are summarized in Table 14.As described in Appendix A.2, the other lines use a negative binomial model for claim frequencies, not the 0-1 inflated model introduced in Section 2.2.For the CO line severity, 1  σ is fitted for the purpose of computation.Model selection and marginal model results can be found in Appendix A.2. Table 13 shows significantly strong negative association between frequency and average severity for the building and contents (BC) line.In contrast, the results are mixed for other lines.Table 14 shows no significant relationships for the CO and IM lines, mild negative relationships for the PN and PO lines, and a strong positive relationship for the CN line.For the BC and CN lines, these results are consistent with the polyserial correlations in Table 5, calculated without covariates.

Dependence between Different Lines
The second level of dependence lies between different lines.In this section, the dependence model for frequencies, severities and aggregate loss with Tweedie margins, as in Section 3.4, are fit.Here, we use marginal results from the inference functions for margins method.In principle, full likelihood can be used.As mentioned previously in this section, in our case, the results of inference functions for margins are close to full likelihood estimation.
Tables 15 and 16 show the dependence parameters of copula models for frequencies and severities, respectively.A Gaussian copula is applied and the composite likelihood method is used for computation.Comparing Tables 4 and 15, it can be seen that frequency dependence parameters decrease substantially.This is due to controlling for the effects of explanatory variables.In contrast, comparing Tables 6 and 16, there appears to be little change in the dependence parameters.This may be due to the smaller impact that explanatory variables have on the severity modeling when compared to frequency modeling.

Out-of-Sample Validation
For out-of-sample validation, the coefficient estimates from the marginals and the dependence parameters are used to obtain the predicted claims with the held-out 2011 data.We compare the independent case, dependent frequency-severity model, and the dependent pure premium approach.

Spearman Correlation
The out of sample validation is performed on the 2011 held-out data, where there are 1098 observations.The claim scores for the pure premium approach are obtained using the conditional mean of the Tweedie distribution for each policyholder.For the frequency-severity approach, the conditional mean for the zero-one-inflated negative binomial distribution is multiplied to the first moment of the GB2 severity distribution for the policyholder.Claim scores for the dependent pure premium approach and the dependent frequency-severity approach are computed using a Monte Carlo simulation of the normal copula.
We first consider the nonparametric Spearman correlation between model predictions and the held-out claims.Four models are considered: the frequency-severity and pure premium (Tweedie) model, assuming independence among lines, and assuming a Gaussian copula among lines.As can be seen from Table 18, the predicted claims are about the same whether dependence is considered or not.The interesting question is how much improvement the zero-one-inflated negative binomial model, and the long-tail distribution (GB2) marginals bring.We observe that the long-tail nature of the severity distribution sometimes results in a large predicted claim.We found that prediction of the mean, using the first moment, can be numerically sensitive.Figure 2 shows a plot of the predicted claims against the out-of-sample claims, for the independent pure premium approach.Figure 3 shows the dependent frequency-severity approach.(In these plots, the claim scores for each line is simulated from the frequency-severity model with dependence, using a Monte Carlo approach with B = 50,000 samples from the normal copula.The model with 01-NB and GB2 marginals show clear improvement for the BC line, in particular for the upper tail prediction.For other lines such as CO, the GB2 marginal results in miss-scaling).

Gini Index
To further validate our results, we use the Gini index to measure the satisfaction of the fund manager with each score.The Gini index is calculated using relativities computed with the actual premium collected by the LGPIF in 2011 as the denominator, with the scores predicted by each model as numerator.This means we are looking for improvements over the original premium scores used by the LGPIF.We expect the Gini index to be higher with the frequency-severity approach, as the fit for the upper tail is better.Figure 4 compares the independent Tweedie approach, and the dependent frequency-severity approach.For the dependent frequency-severity approach, a random 6-dimensional vector is sampled from the normal copula, and the quantiles are converted to frequencies and severities.
For the BC line scores calculated using the Tweedie model, we obtain −1.74% Gini index, meaning this model does not improve the existing premium scores used by the fund.Note that in [59], where the interest is more in the regularization problem, a constant premium is used as denominator for assessing the relativity.Here, the denominator used is the original premiums, which means in order for the index to be positive, there must be an improvement over the original premiums.The dependent frequency-severity scores with B = 50,000, normal copula, and zero-one-inflated NB and GB2 margins results in a Gini index of 22.77%, meaning a clear improvement from the original premium scores.As a side note, the Spearman correlations are: original BC premiums 42.59%, Tweedie model 40.97%, and Frequency-severity model 43.52%, with the out-of-sample claims.Also the reader may observe from Table 18 that the improvement is mostly due to the better marginal model fit, instead of the dependence modeling.

Concluding Remarks
This study shows standard procedures for dependence modeling for a multiple lines insurance company.We have demonstrated that sophisticated marginal models may improve claim score calculations, and further demonstrated how multivariate claim distributions may be estimated using copulas.Our study also verifies that dependence modeling has little influence on the claim scores, but rather is a potentially useful tool for assessing risk measures of liabilities when losses are dependent on one another.A potentially interesting study would be to analyze the difference in risk measures associated with the correlation in liabilities carried by a multiple lines insurer, with and without dependence modeling.We leave this study for future work.
An interesting question is how to predict claims that are large relative to the rest of the distribution.For example, when simulating the GB2 distribution, we regularly generated large predicted values (that our software converted into "infinite values").This suggests that more sophisticated consideration of the upper limits (possibly due to policy limits of the LGPIF) may be necessary to model the claim severities using long-tail distributions.This is actually the log linear model used with errors following the log-F distribution.Thus µ + σ(log α 1 − log α 2 ) can be used as the location parameter associated with covariates.
As a special case of GB2, the location parameter of GG can be derived based on GB2.The density of GG(a, b, α 1 ) is GG(y; a, b, α 1 ) = a Γ(α 1 )y (y/b) aα 1 e −(y/b) a .
When a = 1, the GG distribution becomes the gamma distribution with shape parameter α 1 and scale parameter b. log(b) + log(α 1 ) is the location parameter, which is the log-mean of the gamma distribution, and is hence consistent with the GLM framework.

A.2.1. IM (Contractor's Equipment)
The property fund uses IM as a symbol to denote contractor's equipment, and we follow this notation.Figure A1 shows the residual plot of severity fitted with gamma and GB2 in the IM line.Based on the plot, the GB2 is chosen.Table A1 shows the expected count for each frequency value under different models and empirical values from the data.The proportion of 1 s for the IM line is not high, and hence most models were able to capture this.The property fund uses P for comprehensive, and N to denote new vehicles.Hence PN would mean comprehensive coverage for new vehicles.Figure A2 shows the residual plot of severity, fitted with gamma and GB2 for PN.Based on the plot, the GB2 is chosen for the PN line.Table A3 shows the expected count for each frequency value under different models, and the empirical values from the data.The property fund uses symbol O to denote old, hence PO would be comprehensive coverage for old vehicles.Figure A3 shows the residual plot of severity fitted with gamma and GB2, for the PO line.Table A5 shows the expected count for each frequency value under different models and empirical values from the data.Figure A5 shows the residual plot of severity fitted with Gamma and GB2 for the CO line.GB2 is preferred.Note, here 1  σ instead of σ is fitted for computational stability.Table A9 shows the expected count for each frequency value under different models, and the empirical values from the data.Table A10 shows the goodness of fit tests result.It was calculated using the results in Table A9.The parsimonious model, negative binomial, is selected.Figure A6 shows the cdf plot of jittered aggregate losses, as described in Section 3.7.For most lines, the plots do not show a uniform trend.This tells us that the Tweedie model may not be ideal for such cases.

A.4. Dependence of Frequency and Severity
To motivate this problem, let us think about a classical "aggregate loss" model insurance.In this model, we have a count random variable N representing the number of claims of a policyholder.The random variable S is said to have a "compound distribution."Its moments are readably computable in terms of the frequency and severity moments.Using the law of iterated expectations, we have A.4.2.Average Severity These are calculations basic to the actuarial curriculum.Less common is an expression for the average severity, defined as S = S/N.Note that when N = 0, we define S = 0. Also, use the notation p n = Pr(N = n).Again, using the rule of iterated expectations, we have However, the case for the average severity and frequency is not so clear.
Cov( S, N) = E(S) − E( S)E(N) Thus, S and N are uncorrelated when p 0 = 0. Are S and N independent when p 0 = 0? Basic calculations show that this is not the case.To this end, define the n-fold convolution of the Y random variables, So, S and N provide a nice example of two random variables that are uncorrelated but not independent.
provide summary statistics for the frequencies and severities of claims within the in-sample years 2006 to 2010, and the validation sample 2011.

Figure 1 .
Figure 1.QQ Plot for Residuals of Gamma and GB2 Distribution for BC.

Figure 2 .Figure 3 .
Figure 2. Out-of-Sample Validation for Independent Tweedie.(In these plots, the conditional mean for each policyholder is plotted against the claims.)

Figure A1 .
Figure A1.QQ Plot for Residuals of Gamma and GB2 Distribution for contractor's equipment (IM).

Figure A2 .
Figure A2.QQ Plot for Residuals of Gamma and GB2 Distribution for comprehensive new (PN).

Figure A3 .
Figure A3.QQ Plot for Residuals of Gamma and GB2 Distribution for comprehensive old (PO).

Figure A5 .
Figure A5.QQ Plot for Residuals of Gamma and GB2 Distribution for collision old (CO).
The positive claim amounts are denoted as Y 1 , Y 2 , . . .which are independent and identically distributed.During a specified time interval, the policyholder incurs N claims Y 1 , . . ., Y N , that sum to aggregate loss.In the classical model, the claims frequency distribution N is assumed independent of the amount distribution Y.A.4.1.Moments

Table 3
describes each coverage group.Automobile coverage is subdivided into four subcategories, which correspond to combinations for collision versus comprehensive and for new versus old cars.

Table 3 .
Description of Coverage Groups.

Table 4 .
Polychoric Correlation among Frequencies of Claims.

Table 7 .
Number of Observations.

Table 8 .
Summary of Explanatory Variables.

Table 9 .
Comparison between Empirical Values and ExpectedValues for the building and contents (BC) Line.

Table 10 .
Goodness of Fit Statistics for BC Line.

Table 11 .
Coefficients of Marginal Models for BC Line.

Table 12 .
Vuong Test of Copulas for BC Frequency and Severity Dependence.

Table 13 .
Coefficients of Total Likelihood for BC Line.

Table 14 .
Coefficients of Total Likelihood for Other Lines.

Table 15 .
Dependence Parameters for Frequency.

Table 16 .
Dependence Parameters for Severity.

Table 17
shows the result of dependence parameters for different lines with Tweedie margins.The coefficients of marginal models are in Appendix A.3.

Table A1 .
Comparison between Empirical Values and Expected Values for IM line.

Table A2 .
Goodness of Fit Statistics for IM Line.

Table A3 .
Comparison between Empirical Values and Expected Values for PN Line.

Table A4
shows goodness of fit tests result.It was calculated using the results in TableA3.The simpler model, negative binomial, is preferred.

Table A4 .
Goodness of Fit Statistics for PN Line.

Table A5 .
Comparison between Empirical Values and Expected Values for PO Line.TableA6shows goodness of fit tests result.It was calculated using the results in TableA5.Negative binomial model is selected, based on the test results.

Table A8 .
Goodness of Fit Statistics for CN Line.

Table A9 .
Comparison between Empirical Values and Expected Values for CO Line.

Table A10 .
Goodness of Fit Statistics for CO Line.TableA11shows marginal coefficients for each line.

Table A11 .
Marginal Coefficients of Tweedie Model.