On a Multiplicative Multivariate Gamma Distribution with Applications in Insurance

One way to formulate a multivariate probability distribution with dependent univariate margins distributed gamma is by using the closure under convolutions property. This direction yields an additive background risk model, and it has been very well-studied. An alternative way to accomplish the same task is via an application of the Bernstein–Widder theorem with respect to a shifted inverse Beta probability density function. This way, which leads to an arguably equally popular multiplicative background risk model (MBRM), has been by far less investigated. In this paper, we reintroduce the multiplicative multivariate gamma (MMG) distribution in the most general form, and we explore its various properties thoroughly. Specifically, we study the links to the MBRM, employ the machinery of divided differences to derive the distribution of the aggregate risk random variable explicitly, look into the corresponding copula function and the measures of nonlinear correlation associated with it, and, last but not least, determine the measures of maximal tail dependence. Our main message is that the MMG distribution is (1) very intuitive and easy to communicate, (2) remarkably tractable, and (3) possesses rich dependence and tail dependence characteristics. Hence, the MMG distribution should be given serious considerations when modelling dependent risks.


Introduction
Let X be a collection of actuarial risks, that is let it contain random variables (r.v.'s) X : Ω → R defined on the probability space (Ω, F , P) and interpreted as the financial risks an insurer is exposed to. Often, for applications in insurance, actuaries would consider those X ∈ X , whose distributions are supported on the non-negative real half-line, have positive skewness, and allow for a certain degree of heavy-tailness. One such distribution, which has been of prominent importance in insurance applications, is gamma. We refer to Hürlimann (2001), Dornheim and Brazauskas (2007), Furman et al. (2018), and Zhou et al. (2018) for applications in solvency assessment, loss reserving, and aggregate risk approximation, respectively. Furthermore, let γ ∈ R + and σ ∈ R + denote, correspondingly, the shape and scale parameters, then the r.v. X is said to be distributed gamma, succinctly X ∼ Ga(γ, σ), if it has the probability density function (p.d.f.) where Γ(·) stands for the complete gamma function. The popularity of the r.v.'s distributed gamma in insurance applications is not surprising: the p.d.f.'s of the (aggregate) insurance losses have as for all (x 1 , . . . , x n ) ∈ R n + , which inconveniently takes distinct forms for each of the n! orderings of x 1 , . . . , x n . Remark 1. The r.v.'s Y 1 , . . . , Y n and Y n+1 are often interpreted as, respectively, the specific and systematic risk factors. The systematic risk factor, Y n+1 , has also been referred to as the background risk (Gollier and Pratt 1996), and so the distribution of the r.v. X = (X 1 , . . . , X n ) can be associated with an Additive Background Risk Model with risk components distributed gamma (G-ABRM). Succinctly, for γ = (γ 1 , . . . , γ n ) , we write X ∼ Ga + n (γ, γ n+1 , σ), where γ n+1 serves as the dependence parameter.
An alternative way to link the specific risk factors and the systematic (or background) risk factor is with the help of multiplication. Namely, in order to formulate a Multiplicative Background Risk Model with the risk components distributed gamma (G-MBRM), we must find a sequence of (n + 1) independent r.v.'s Z 1 , . . . , Z n , Z n+1 , say, such that X = (X 1 , . . . , X n ) = (Z 1 Z n+1 , . . . , Z n Z n+1 ) results in the coordinates of the r.v. X being distributed gamma. One solution of this exercise, which is of pivotal importance for this paper, can be found in Feller (1968) (also, Albrecher et al. 2011;Sarabia et al. 2018). We organize the rest of the paper as follows: in Section 2, we explore the basic distributional properties of-what we call-the multiplicative multivariate gamma (MMG) distribution. Then, in Sections 3 and 4, respectively, we discuss in detail and elucidate with examples of actuarial interest the aggregation and (tail) dependence properties of the MMG distribution. Section 5 concludes the paper. All proofs are relegated to Appendix A to facilitate the reading.

Definition and Basic Properties
Multivariate distributions lay the very foundation of the successful (insurance) risk measurement-and thus of the consequent risk management-processes. However, the toolbox of the available stochastic dependencies that can be used to link stand-alone risk components into risk portfolios is somewhat overwhelming. Indeed, there are infinitely many ways to formulate the joint distribution of two dependent risk r.v.'s, whereas there is a single way only to write this distribution under the assumption of independence. The case of the multivariate distributions with the margins distributed gamma is of course not an exception (e.g., Kotz et al. 2000).
Nevertheless, real applications impose significant constraints on the model choice. Namely, practitioners often opt for those multivariate distributions that: (i) admit meaningful and relevant interpretations; (ii) allow for an adequate fit to the modelled data, be it in the 'tail', in the 'body', and/or in the dependence; and (iii) can be readily implemented. We feel that the multivariate distribution with the univariate margins distributed gamma that we put forward next (also, Albrecher et al. 2011;Sarabia et al. 2018) is exactly such.
Formally, let E λ and Λ denote, respectively, an exponentially distributed r.v. with the rate parameter λ ∈ R + and an arbitrarily distributed r.v. with the range Λ ∈ R + ; assume that the r.v.'s E Λ and Λ are independent. In addition, let ' * ' represent the mixture operator (e.g., Feller 1968;Su and Furman 2017a), such that, for ' We note in passing that the just-mentioned mixture operator is referred to as 'randomization' in Feller (1968), and is closely related-via the Bernstein-Widder theorem-to the notion of the Laplace transform of the p.d.f. of Λ. More specifically, if f Λ and L{ f Λ } denote, correspondingly, the p.d.f. of Λ and its Laplace transform, which is Recall that in this paper we are interested in formulating a multivariate distribution with the univariate margins distributed gamma and a dependence. To this end, we assume that the r.v. Λ is distributed as a special shifted inverse Beta, succinctly Λ ∼ IB(γ), with the p.d.f.
where γ ∈ (0, 1) is the shape parameter. In our context, the choice of p.d.f. (4) is unique, which readily follows from the Bernstein-Widder theorem. The next few facts are used frequently later on in the paper, and are hence formulated as a lemma. In the following, the k-th order derivative of the Laplace transform is denoted by ψ (k) , k ∈ N := {1, 2, . . .}, also R 0,+ := [0, ∞).
Remark 2. Let E j := E 1,j , j = 1, . . . , n denote independent copies of a r.v. distributed exponentially with unit scale, then the joint distribution of the r.v. X = (X 1 , . . . , X n ) in Definition 1 admits the following multiplicative background risk model representation (see, Asimit et al. 2016;Frank et al. 2006, for the corresponding economic implication and application) Above, the r.v. Λ can be interpreted as the systematic risk factor that endangers every risk component of the portfolio X = (X 1 , . . . , X n ) in Equation (6). The Monte Carlo simulation of Equation (6) is immediate.
for all (x 1 , . . . , x n ) ∈ R n 0,+ . (iii) The p.d.f. that corresponds to d.d.f. (7) is, for all (x 1 , . . . , x n ) ∈ R n + , The following facts are immediate from Theorem 1: (i) the distribution of X ∼ Ga × n (γ, σ) is marginally closed, namely, X j ∼ Ga(γ, σ j ), j = 1, . . . , n; (ii) the mathematical expectation of the j-th coordinate is E[X j ] = γσ j ; and (iii) the variance of the j-th coordinate is Var[X j ] = γσ 2 j . We further explore some less obvious properties of the MMG/G-MBRM and note with satisfaction that the risk portfolios with the joint distributions within this class are often more tractable than the portfolios having stochastically independent risk components distributed gamma. At the outset, we note in passing that the MMG distribution put forward in Definition 1 is a non-exchangeable generalization of the multivariate distributions having univariate margins distributed gamma that are discussed in Albrecher et al. (2011);Sarabia et al. (2018). As such, the MMG distribution requires a more delicate treatment when deriving the results below, which hinge crucially on the stochastic characteristics of the univariate margins of the r.v. X ∼ Ga × n (γ, σ). We look into the minima and maxima r.v.'s first; both are of evident importance in insurance. To this end, denote by X min = min i=1,...,n X i ∼ F min and by X max = max i=1,...,n X i ∼ F max the minima and maxima r.v.'s. Then we have -unlike in the independent case -that the coordinates of the r.v. X = (X 1 , . . . , X n ) in Definition 1 are closed under minima.
Another r.v. of pivotal interest in insurance is the aggregate risk r.v. denoted by X + = X 1 + · · · + X n ; in addition, let X + ∼ F + . It is well known that, if X 1 , . . . , X n are mutually independent and distributed gamma with arbitrary parameters, then F + admits an infinite sum representation (Moschopoulos 1985;Provost 1989). We further show that for X ∼ G × n (γ, σ) and when all the scale parameters are distinct, then F + is noticeably more elegant. The derivation of F + in the general case-i.e., for arbitrary (possibly equal) scale parameters-is more cumbersome and is presented in Section 3. Let We often write w i omitting the vector of scale parameters σ for the simplicity of notation.
and assume that all the scale parameters are distinct, which is The last result in this section provides an expression for the higher-order (product) moments of the r.v. X ∼ G × n (γ, σ). We employ a special form of this expression later on in Section 4 to derive the formula for the Pearson index of linear correlation.
We conclude the discussion in the present section by noticing that joint p.d.f. (8) can be used to estimate the parameters of the MMG distribution via the (numerical) maximum likelihood approach, whereas expression (11) is of interest if the moment-based estimation is being pursued.

Aggregation Properties of the Multiplicative Multivariate Gamma Distribution
One of the key paradigms in the modern enterprise risk management requires that all risks are treated on a holistic basis. As a result, risk aggregation is of fundamental importance for the effective conglomerate-wide risk management, risk-sensitive supervision, and a great variety of other business decision making processes. In the context of the MMG distribution, when all scale parameters are distinct, the decumulative distribution function of the aggregate risk r.v. is given by (10). The situation with arbitrary (possibly equal) scales is more involved. We further show that, within the MMG/G-MBRM class, that is for X ∼ Ga × n (γ, σ) with arbitrary scale parameters and γ ∈ (0, 1), the d.d.f. of the aggregate risk r.v. X + admits a finite sum representation. To this end, we employ Risks 2018, 6, 79 6 of 20 the well-studied machinery of divided differences (e.g., Milne-Thomson 2000, for a comprehensive treatment). The rest of the section is divided into two: theoretical considerations and applications.

Theoretical Considerations
We remind at the outset that the divided differences, denoted by ω(y 1 , . . . , y m ), on a grid ∆ = {y 1 , . . . , y m } for a function ω : R → R can be written as (e.g., Milne-Thomson 2000) Denote Then the following corollary is merely a rearrangement of Equation (10).
Obviously, Equation (12) does not yield sensible results when some of the scale parameters of the r.v. X ∼ Ga × n (γ, σ) are equal. To circumvent this inconvenience, we formulate and prove the following lemma.
The next assertion establishes the distribution of the aggregate risk r.v. with arbitrary scale parameters. Its proof follows by rearranging d.d.f. (10) using the divided differences operator and consequently evoking Lemma 2.
where γ ∈ (0, 1) and σ = (σ 1 , . . . , σ n ) with arbitrary coordinates in the latter vector of parameters. Let σ = (σ 1 , . . . , σ 1 n 1 , . . . , σ m , . . . , σ m n m ) for m ∈ N and n 1 + · · · + n m = n, then, for x ∈ R 0,+ , the d.d.f. of X + admits the following finite sum form: Remark 3. A close look at Theorem 4 reveals that the distribution of the r.v. X + can be considered a finite mixture of the r.v.'s distributed Erlang with stochastic scale parameters. To see this, first note that in which e i,h i denotes the r.v. distributed Erlang with the shape parameter h i + 1 and the random scale parameter By setting However, these weights are not necessarily positive. For an example, consider the bivariate case with n = m = 2 and n 1 = n 2 = 1. A simple calculation yields Therefore, depending on the values of σ 1 and σ 2 , one of the weights must be negative.

Applications
Herein we confine the discussion to the individual and collective risk models. In this respect, recall that we call the r.v. S n = X 1 + · · · + X n , n ∈ N the individual risk model, where we let the severity r.v.'s X j , j = 1, . . . , n be possibly non-homogeneous. In the collective risk model case, for N ∈ Z 0,+ := {0, 1, 2, . . .} denoting the frequency r.v., we are interested in exploring the distribution of the random sum S N = X 1 + · · · + X N . In the context of the MMG/G-MBRM, we have and P.d.f.'s-rather than d.d.f.'s-often play an important role in the individual/collective risk model contexts. Therefore, the p.d.f.'s of S n and S N engage us in the rest of this section. We start with the p.d.f. of the former r.v. in the following proposition. Recall that p i,h i are given in (14), and ψ (k) denotes the k-th order derivative of the Laplace transform.
We further derive the p.d.f. of the r.v. S N .
We conclude this section by specializing the p.d.f. of the r.v. S N reported in Proposition 3 for particular choices of the frequency r.v. In actuarial science, some popular choices of the r.v. N are, e.g., the Poisson, negative binomial, and logarithmic (e.g., Klugman et al. 2012). Below, we first remind the reader in passing the probability mass functions (p.m.f.'s) of the just-mentioned r.v.'s, and we then present the p.d.f.'s of the aggregate r.v.'s within the framework of the corresponding collective risk models.
Let Φ 1 and Φ 3 , respectively, denote the two-variable confluent hypergeometric series of the first and third kind (see, e.g., Srivastava and Karlsson 1985), that is with a 1 , a 2 , a 3 ∈ R + , for x, y ∈ R. The following corollary follows readily.

Dependence Properties of the Multiplicative Multivariate Gamma Distribution
At first sight, the dependence structure that underlies the MMG distribution-that is d.d.f. (7)is not as versatile as the one behind the additive counterpart of Mathai and Moschopoulos (1991). This is because the Pearson correlation, ρ, for the former class of distributions does not attain every value in the interval [0, 1], whereas it does so in the context of the latter class of distributions (e.g., Das et al. 2007;Furman 2017a, 2017b, for a similar constraint in the context of default risk). More formally, we have the following proposition, the proof of which is a direct application of Theorem 3 and is thus omitted.
In the rest of this section, we show that the just-mentioned seeming shortcoming should in fact be attributed to the Pearson index of correlation, ρ, itself, rather than to the dependence structure of the MMG distribution. As hitherto, we divide our observations herein into two subsections.

Theoretical Considerations
At the outset, we observe that the dependence structure that underlies the MMG/G-MBRM is not linear in the-common-background r.v. Λ. Therefore, the machinery of copulas lands itself very naturally to exploring the relevant dependence properties. The next theorem states the copula function (e.g., Joe 1997) of X ∼ Ga × n (γ, σ).

Remark 5. Copula function
is a legitimate completely monotonic function-known as the Archimedean generator-and φ −1 is its inverse (e.g., McNeil and Nešlehová 2009). The MMG copula therefore enriches the encompassing toolbox of the distinct Archimedean copulas available to researchers and practitioners.
We have already mentioned at the end of Section 2 that the maximum likelihood approach can be used in order to numerically estimate the parameters of the MMG distribution. An alternative way to estimate the parameters is via the two-step copula approach. That is, we first fit the MMG copula to the pseudo uniform samples based on the empirical c.d.f.'s of X i , i = 1, . . . , n and estimate the γ parameter (e.g., Embrechts and Hofert 2013;Genest et al. 2011, and references therein), and we then estimate the σ i parameters based on the univariate marginal distributions assuming that the γ parameter is known. Given the cumbersome form of p.d.f. (8), the copula-based approach is computationally simpler.
Besides the just-mentioned statistical inference, a useful contribution of copulas to the vast literature of multivariate modelling is that they have given rise to a number of indices of dependence that circumvent the known fallacies of the Pearson ρ. Such indices of dependence are, e.g., the Kendall τ and Spearman ρ S measures of rank correlation, and we derive these two in the next subsection in the context of the MMG copula function C γ . In the rest of this subsection, we build up the theoretical groundwork necessary for exploring the tail dependence of C γ . As tail dependence represents the co-movement of extreme risks, it is of particular importance in the era following the financial crisis of 2007-2009. We note in passing that, since the majority of the existing methods for quantifying tail dependence mainly aim at random pairs, we specialize the discussion in this part of the present report to the bivariate case only.
Recently, an argument has been put forward that tail dependence measures (21) and (22) may underestimate the extent of the tail dependence inherent in a copula. More specifically, Furman et al. (2015) claim and elucidate with numerous examples that as measures (21) and (22) are computed along the main diagonal (u, u), u ∈ [0, 1], their values are not necessarily maximal when alternative paths in [0, 1] 2 are considered. This motivated the following definitions of the admissible paths and the paths of maximal dependence in ibid.
Then the path (ϕ(u), u 2 /ϕ(u)) 0≤u≤1 is admissible whenever the function ϕ is admissible. In addition, we denote by A the set of all admissible functions ϕ.
Obviously, the function ϕ 0 (u) = u is admissible and yields the representation of the diagonal path that serves as a building block for classical indices (21) and (22). For the Archimedean class of copulas, the following property of the maximal dependence path holds. The verification of the condition stated in Lemma 3 below is not trivial, and is carried out for the MMG copula C γ in Theorem 6.
The next lemma on a L'Hospital type rule for monotonicity, plays an importantly auxiliary role when deriving the maximal dependence path for C γ .
Our last result in this subsection implies that measures of tail dependence (21) and (22) are in fact maximal in the context of the MMG copula C γ . Theorem 6. The maximal dependence path of the copula function C γ in (19) is diagonal.
Risks 2018, 6, 79 13 of 20 For a 1 , . . . , a q+1 all positive, and these are the cases of interest in the present report. The radius of convergence of the series is the open disk |z| < 1. On the boundary |z| = 1, the series converges absolutely if d = b 1 + · · · + b q − a 1 − · · · − a q+1 > 0, and it converges except at z = 1 if 0 ≥ d > −1.
Proposition 5. For the copula C γ , the Kendall τ rank correlation is given by and the Spearman ρ S rank correlation is given by Proposition 6. Assume that X ∼ Ga × n (γ, σ) has copula C γ , the lower maximal tail dependence of C γ is The upper maximal tail dependence of C γ is λ U (C γ ) = 2 − 2 γ , and χ U (C γ ) = 1.
Proposition 6 readily implies-recall to this end that the copula C γ is in fact a survival copula (by construction)-that the coordinates of X ∼ Ga × n (γ, σ) are asymptotically dependent in the lower tail, but independent in the upper tail. Speaking bluntly, this means that X is more likely to take smaller values simultaneously, but less likely to form a cluster of large values. This conforms to the already made intuitive observation that the copula C γ can serve as a reflected variant of the well-studied Clayton copula.

Conclusions
In the present report, we have systematically studied a class of multivariate multiplicative gamma distributions. We have demonstrated that the MMG distribution admits a very meaningful background risk model representation, where the interdependencies among risks are implied by a common systematic risk factor. Moreover, we have shown that the MMG distribution enjoys a remarkable level of analytical tractability, that is, the risk r.v.'s distributed MMG are straightforward to simulate, easy to aggregate and take maxima, closed under minima, and have attractive dependence and tail dependence characteristics. In view of the above, we think that the potential applications of the MMG distribution in actuarial science are vast, and we hope to draw the attention of the community to this class of distributions. In fact, reduced forms of the proposed MMG distribution have been recently heuristically adopted in the actuarial literature to model a variety of dependent insurance risks (e.g., Sarabia et al. 2018, also, Albrecher et al. 2011).
This completes the proof of the lemma.
Proof of Theorem 1. The d.d.f.'s of the r.v.'s X and X follow immediately from Lemma 1, Statement (i) and Chapter 4 in Joe (1997). The joint p.d.f. follows from Lemma 1, Statement (iii) since This completes the proof of the theorem.
Proof of Theorem 2. The closure under the minima operation is trivial by evoking Theorem 1, Statement (ii). The distribution of the r.v. X max follows immediately (e.g., Corollary 2.2 in Su and Furman 2017a, for a similar result in the context of a multivariate Pareto distribution). This completes the proof of the theorem.
Proof of Theorem 3. We immediately have where ' (1) = ' holds due to the moments' formula in the case of the exponentially distributed r.v.'s (see, e.g., Klugman et al. 2012). The proof is then completed by evoking Lemma 1, Statement (ii).
Proof of Proposition 2. The proof of the proposition follows from Remark 3 that reports the mixture representation. Namely, This completes the proof.

Proof of Proposition 3.
We have the following string of equations, for all x ∈ R + , This completes the proof of the proposition.
We next turn to study the p.d.f. of C γ . By definition, we readily obtain where f denotes the p.d.f. of Ga(γ, 1). This completes the proof of the theorem.
We now turn to study the upper tail dependence of C γ . Note that the mixture r.v. Λ has d.d.f. F Λ ∈ RV −γ that varies regularly at infinity with order −γ (Bingham et al. 1987). The expressions for λ U and χ U are readily obtained by evoking Corollary 3.3 in Su and Hua (2017). This completes the proof of this proposition.