Oskar Laverny
email: oskar.laverny@univ-lyon1.fr

Esterina Masiello
email: esterina.masiello@univ-lyon1.fr

Véronique Maume-Deschamps
email: veronique.maume@univ-lyon1.fr

Didier Rullière
email: didier.rulliere@emse.fr

Estimation of multivariate generalized gamma convolutions through Laguerre expansions

Keywords: Primary 62H12, 60E07; secondary 60E10 Multivariate Generalized Gamma Convolutions, Estimation, Thorin's measure, Laguerre's basis, Infinite Divisibility

come   L'archive ouverte pluridisciplinaire

Introduction

Olof Thorin introduced the Generalized Gamma Convolutions (GGC) class of distributions in 1977 as a tool to study the infinite divisibility of Pareto [START_REF] Thorin | On the Infinite Divisibility of the Pareto Distribution[END_REF] and log-Normal [START_REF] Thorin | On the Infinite Divisibility of the Lognormal Distribution[END_REF] distributions. The quite natural and pragmatic question of identifying if some distributions are infinitely divisible or not turned out to be theoretically fruitful, and notably gave rise to the concept of hyperbolic complete monotone functions. Later, L. Bondesson extended these concepts in a book [START_REF] Bondesson | Generalized Gamma Convolutions and Related Classes of Distributions and Densities[END_REF] which is still today a good reference on the subject. An introduction to more recent literature about GGC can be found in the two surveys [START_REF] James | Generalized Gamma Convolutions, Dirichlet Means, Thorin Measures, with Explicit Examples[END_REF][START_REF] Bondesson | A Class of Probability Distributions That Is Closed with Respect to Addition as Well as Multiplication of Independent Random Variables[END_REF], on the probabilistic side of the problem. This class contains a lot of known and commonly used distributions, including heavy-tailed distributions such as Pareto, log-Normal and α-stable distributions. This versatility makes the class appealing for the statistical practitioners working in various fields such as life and non-life insurance, reinsurance, anomaly or fraud detection, floods analysis and meteorology, etc. Other classes of general approximating distributions with high-order features can be considered. Mixtures of Erlangs [START_REF] Cossette | Moment-Based Approximation with Mixed Erlang Distributions[END_REF][START_REF] Lee | Modeling Dependent Risks with Multivariate Erlang Mixtures[END_REF][START_REF] Tijms | A First Course in Stochastic Models[END_REF], for example, have the property of being dense in L 2 [START_REF] Tijms | A First Course in Stochastic Models[END_REF] and allow for fast k-MLE estimation algorithms [START_REF] Nielsen | k-MLE: A Fast Algorithm for Learning Statistical Mixture Models[END_REF]. On the other hand, they lack crucial closure properties and interpretability that the generalized gamma convolutions have. Surprisingly, on the statistical side, very little work can be found about estimation of generalized gamma convolutions. Only one projection procedure, that is a procedure that projects generalized gamma convolutions onto convolutions of a finite number of gammas, was published recently in [START_REF] Miles | Risk Aggregation: A General Approach via the Class of Generalized Gamma Convolutions[END_REF][START_REF] Furman | On Log-Normal Convolutions: An Analytical-Numerical Method With Applications to Economic Capital Determination[END_REF], Although the resulting convergence result is stunning, this procedure cannot handle cases where the incoming density is not already inside the class. This typically happens for any empirical distribution, due to sampling noise, even if the sampling density is in the class. Furthermore, we show that the theoretical background of the estimation procedure in [START_REF] Miles | Risk Aggregation: A General Approach via the Class of Generalized Gamma Convolutions[END_REF][START_REF] Furman | On Log-Normal Convolutions: An Analytical-Numerical Method With Applications to Economic Capital Determination[END_REF] is inherently univariate, and no direct extension to a multivariate case is possible.

The multivariate analogue is still an active research field, and, as far as we know, has never been considered as a statistical tool in the literature. In applications, however, statisticians might deal with marginals that are in this class. A way to handle dependence structure in this interpretable framework seems therefore appealing as it would allow studying functions of the random vector, like the sum of components, under dependence assumptions. Therefore, an estimation procedure for multivariate generalized gamma convolutions that takes into account the dependence structure could be useful. The probabilistic background is already developed: based on an original idea from [START_REF] Cherian | A Bivariate Correlated Gamma-Type Distribution Function[END_REF], Bondesson introduces a bivariate extension of the Thorin class, and gives some properties about these distributions. Later, [START_REF] Bondesson | On Univariate and Bivariate Generalized Gamma Convolutions[END_REF] extends this concept to the multivariate case. More recently, [START_REF] Pérez-Abreu | A Class of Infinitely Divisible Multivariate and Matrix Gamma Distributions and Cone-Valued Generalised Gamma Convolutions[END_REF][START_REF] Pérez-Abreu | Infinitely Divisible Multivariate and Matrix Gamma Distributions[END_REF] propose another equivalent framework allowing extension on any cone, with application to matrix gamma convolutions. On the other hand, no estimation procedure nor explicit density exist. We discuss this multivariate class and show that the produced dependence structures are highly flexible: non-exchangeability and tail dependence are possible.

Regarding tools, the literature gives two series expansions for the density of univariate gamma convolutions: Mathai [START_REF] Mathai | Storage Capacity of a Dam with Gamma Type Inputs[END_REF] proposes a density based on Kummer's confluent geometric functions and Moschopoulos [START_REF] Moschopoulos | The Distribution of the Sum of Independent Gamma Random Variables[END_REF] refined it as a convergent gamma series, with an explicit truncation error. Both these densities have the problem of being based on the smallest scale parameter, and are not well conditioned when the smallest scale is too small. Unfortunately, this is typically the case when approximating a log-Normal or a Pareto distribution by a finite convolution of gammas. Therefore, no stable procedure for the density computation is available for the entirety of the parameter range. We are not aware of any density estimation available in the literature for the multivariate case.

We consider here the problem of estimation from samples of multivariate generalized gamma convolutions. By the introduction of a specific Laguerre basis, which was already considered for non-parametric deconvolution problems by [START_REF] Dussap | Anisotropic Multivariate Deconvolution Using Projection on the Laguerre Basis[END_REF][START_REF] Mabon | Adaptive Deconvolution on the Non-Negative Real Line: Adaptive Deconvolution on R+[END_REF][START_REF] Benhaddou | Anisotropic Functional Laplace Deconvolution[END_REF], we manage to obtain a convergent series for the density. Classical deconvolution problems usually consider only one source of signal and one source of noise [START_REF] Meister | Deconvolution Problems in Nonparametric Statistics[END_REF][START_REF] Mendel | Maximum-Likelihood Deconvolution. Signal Processing and Digital Filtering[END_REF][START_REF] Zhang | Optimizing Laguerre Expansion Based Deconvolution Methods for Analysing Bi-Exponential Fluorescence Lifetime Images[END_REF][START_REF] Sprott | Estimating the Parameters of a Convolution by Maximum Likelihood[END_REF]]. Although we have here a finite number of signals to be estimated, the Laguerre expansion is still a useful tool. We show that this expansion is quite natural for our problem, as it includes and generalizes Moschopoulos's univariate series. Using this basis, we bridge the gaps and provide a new stable algorithm for density computation in both the univariate and multivariate cases, as well as estimation procedures which handle both clean data (the density is given as a formal function) and dirty data (such as an empirical dataset in 64 bits precision). The resulting algorithms are implemented in the Julia package ThorinDistributions.jl, available on GitHub1 under a MIT License and archived in Zenodo [29].

After fixing notations, Section 2 covers some definitions, properties, and algorithms to set up the stage. Section 3 considers the Laguerre expansion of densities in the multivariate Thorin class, provides some error bounds and develops an estimation procedure for these distributions. In Section 4, we investigate the numerical results of our estimation procedure on several examples. Section 5 concludes and gives leads for further work.

Background

As we deal with a lot of multivariate objects, we start by fixing some notations.

We use bold letters such as a to designate generally indexable objects, e.g., vectors or matrices, and corresponding indexed and possibly unbolded letters, such as a i , designate values in these objects. We use Cartesian indexing. For example, if we consider a row-major matrix r, then r i denotes the i th row of that matrix, r i,j the j th value of that row and if k = (i, j), then r k = r i,j .

We denote |x| = x 1 + . . . + x d the sum of components of a vector x. The product of factorials of component of an integer vector k is denoted k! = k 1 !...k d !, and we set

x k = x k1 1 . . . x k d d .
We denote derivatives of a multivariate function f by f Inequalities, sums, and products between vectors of same shape are always intended componentwise, and standard broadcasting between objects of different shapes applies. For example, 1 + 2x, 1 + 2x and (1 + 2x) x∈x all denote the same object, which has same shape as x.

(k) (x) = ∂ ∂ k 1 x1 . . . ∂ ∂ k d x d (
To characterize the distribution of random vectors, we might use several functions. For a (real-valued) random vector

X of dimension d ≥ 1, the cumulative distribution function (cdf) of X is F (x) = P (X ≤ x). When they exist, the probability density function (pdf) of X, f (x) = ∂ ∂x F (x), its Laplace transform (lt) L(t) = e -tx f (x) dx, its moment generating function (mgf) M (t) = L(-t)
or its cumulant generating function (cgf) K(t) = ln(M (t)) can also be used to characterize the distribution. All these functions might be indexed by the random object to avoid ambiguity if needed.

We denote by L 2 (R d + ) the set of functions that are Lebesgue square-integrable onto R d + , and by L 2 (R + , w) the set of square-integrable functions with respect to a weight function w.

The following subsections set the stage of our analysis by reviewing common tools from the literature.

Gamma convolutions classes

A random variable X is said to be Gamma distributed with shape α ∈ R * + , scale s ∈ R * + and rate 1 s , which we denote X ∼ G 1,1 (α, s), if it admits the cumulant generating function (cgf):

K(t) = ln E e tX = -α ln (1 -st) .
We denote by G 1,1 this class of distributions. Recall that the cumulant generating function of the sum of independent random variables is the sum of their cumulant generating functions. Therefore, gamma distributions are clearly infinitely divisible with gamma distributed pieces (same scale and smaller shape). We now consider independent convolutions of these distributions:

Definition 2.1.1 (G 1,n and G 1 ). Let X be a random variable with support R + . (i) X is a gamma convolution with shapes α ∈ R n + and scales s ∈ R n + , denoted X ∼ G 1,n (α, s), if it has cgf: K(t) = - n i=1 α i ln (1 -s i t) .
We denote this class of distributions G 1,n .

(ii) X is an (untranslated) generalized gamma convolutions with Thorin measure ν, denoted X ∼ G 1 (ν), if its cgf writes:

K(t) = - ∞ 0 ln (1 -st) ν(ds),
where ν is such that:

[1,+∞) | log(s)| ν(ds) < ∞ and (0,1) s ν(ds) < ∞.
We denote by G 1 this class of distributions.

The generalized gamma convolutions class G 1 , which is commonly called the one dimensional Thorin class, has received a lot of interest in the literature. These limits of convolutions are not dense in the set of positive continuous random variables, but they contain many interesting distributions, including the log-Normal and Pareto distributions as historically important cases [START_REF] Thorin | On the Infinite Divisibility of the Pareto Distribution[END_REF][START_REF] Thorin | On the Infinite Divisibility of the Lognormal Distribution[END_REF]. They have several interesting properties: noteworthy is the fact that the G 1 class is closed with respect to independent sums of random variables, by construction, but also, as a more recent result shows, by independent products of random variables [START_REF] Bondesson | A Class of Probability Distributions That Is Closed with Respect to Addition as Well as Multiplication of Independent Random Variables[END_REF]. In fact, G 1 can be defined as the smallest class that is closed through convolution and that contains gamma distributions. See [START_REF] Bondesson | Generalized Gamma Convolutions and Related Classes of Distributions and Densities[END_REF] for an extensive study of these distributions2 .

Remark that if the Thorin measure

ν of X ∼ G 1 (ν) has a finite discrete support of cardinal n, X ∼ G 1,n (α, s) for α, s vectors such that ν = i α i δ si (where δ x
is the Dirac measure in x).

Mathai [START_REF] Mathai | Storage Capacity of a Dam with Gamma Type Inputs[END_REF] provides a series expression for the density of G 1,n random variables based on Kummer's confluent hypergeometric functions. Later, Moschopoulos [START_REF] Moschopoulos | The Distribution of the Sum of Independent Gamma Random Variables[END_REF] refines the result by providing the following gamma series:

Property 2.1.2 (Moschopoulos expansion). For X ∼ G 1,n (α, s), denoting without loss of generality s 1 = min s, the density f of X is given by the following series:

f (x) = C ∞ k=0 δ k f ρ+k, 1 s 1 (x)
where:

• f a,b is the density of a G 1,1 (a, b) distribution,
• ρ = |α| and C = s ρ 1 s -α , • δ 0 = 0 and for k ≥ 0,

δ k+1 = 1 k + 1 k+1 i=1 iγ i δ k+1-i , where γ k+1 = n i=1 α i (1 -s1 si ) k+1 k + 1 .
Furthermore, for S = 1 -s1 max s , we have that |δ k+1 | ≤ S k+1 Γ(ρ+k) Γ(ρ)(k+1)! , ensuring the uniform convergence of the series.

The dependence of this expansion onto the smallest scale parameter s 1 has a major drawback. When the smallest scale is close to zero, the corresponding truncated series is quite unstable. A simple show-case is the simulation of random numbers from a G 1,2 10, 10 -3 , 1, 10 -3 distribution: the Moschopoulos implementation from the R package coga [START_REF] Hu | Density and Distribution Evaluation for Convolution of Independent Gamma Variables[END_REF] gives a density that evaluates to 0 on all simulated random numbers, which is obviously wrong. Mathai's method, implemented in the same package, produces the same result. Sadly, as later implementation will show, parameters that correspond to approximations of useful distributions such as log-Normal, Weibull or Pareto usually have very small scales and trigger the same instability.

A likelihood approach to fit distributions in the Thorin class is therefore not practical, as there is no stable density, which might explain why there are currently no estimation procedures in the literature.

In [START_REF] Bondesson | Generalized Gamma Convolutions and Related Classes of Distributions and Densities[END_REF][START_REF] Bondesson | On Univariate and Bivariate Generalized Gamma Convolutions[END_REF], Bondesson defined a class of multivariate convolutions of gamma distributions based on the following idea from Cherian [START_REF] Cherian | A Bivariate Correlated Gamma-Type Distribution Function[END_REF]. Suppose that Z = (Z i ) i∈{0,1,2} are three independent gamma random variables with (respective) shapes α i and scales s i such that s 0 = 1 w.l.o.g. Then the cgf of the vector (Z 0 + Z 1 , Z 0 + Z 2 ) writes:

K(t) = -α 1 ln (1 -s 1 t 1 ) -α 2 ln (1 -s 2 t 2 ) -α 0 ln (1 -s 1 t 1 + s 2 t 2 ) .

This construction can be extended to what

K(t) = n j=1 -α j ln (1 -s j , t ) .
We denote by G d,n this class of distributions.

(ii) X is an (untranslated) multivariate generalized gamma convolutions with Thorin measure ν, denoted X ∼ G d (ν), if it is a weak limit of d-variates gamma convolutions. Its cgf writes:

K(t) = -ln (1 -s, t ) ν(ds),
for a positive measure ν of R d + , the Thorin measure, with suitable integration conditions (see [START_REF] Pérez-Abreu | A Class of Infinitely Divisible Multivariate and Matrix Gamma Distributions and Cone-Valued Generalised Gamma Convolutions[END_REF]). We denote by G d this class of distributions.

Note that X ∼ G d,n (α, s) if and only if there exists a vector Y of n independent unit scale gamma random variables with shapes α 1 , ..., α n such that X = s Y . The i th marginal X i has distribution G 1,n (α, s i ). We can interpret this as a decomposition of the random vector X into an additive risk factor model with gamma distributed factors. Note that we allow scales to be zero, and some factors might therefore appear only in some marginals. As in the univariate case, we close the class by taking convolutional limits. For the analysis of these distributions, we refer to [START_REF] Bondesson | Generalized Gamma Convolutions and Related Classes of Distributions and Densities[END_REF][START_REF] Bondesson | On Univariate and Bivariate Generalized Gamma Convolutions[END_REF], but also [START_REF] Pérez-Abreu | A Class of Infinitely Divisible Multivariate and Matrix Gamma Distributions and Cone-Valued Generalised Gamma Convolutions[END_REF] which uses a slightly different but equivalent parametrization, allowing a generalization to other cones than R d + .

We are here interested in the construction of distributions estimators in G d , with particular interest for estimators in G d,n . Indeed, G d,n models have a structure that allows fast simulation, and that provides meaningful insight about the dependence structure, since a G d,n distribution follows an additive risk factor model. On the other hand, distributions in G d with a non-atomic Thorin measure are hard to sample (see [START_REF] Bondesson | On Simulation from Infinitely Divisible Distributions[END_REF][START_REF] Ridout | Generating Random Numbers from a Distribution Specified by Its Laplace Transform[END_REF] for potential solutions to this problem). They have no known closed-form density or distribution function, and even the cumulant generating function requires integration to be evaluated.

Before diving into statistical considerations and estimation of these distributions in Section 3, we give here some properties about the dependence structure that can be achieved in the G d class. The following properties give more insights about the relation between the dependence structure of the random vector and the dependence structure of its Thorin measure:

Property 2.1.4 (The support of the Thorin measure). Let X ∼ G d (ν), and let Sup(ν) be the support of ν. Then we have:

(i) The marginals of X are mutually independent if and only if

Sup(ν) ⊂ S ⊥ ⊥ = d j=1 j-1 × i=1 {0} × R * + × d × i=j+1 {0} .
In this case, ν(s) =

Sup(ν) ⊂ S c = r ∈ R d + , r r = c .
In this case, the multivariate Thorin measure and its marginals have all the same total mass.

(iii) The r.v. X has a Lebesgue-square-integrable density if D ≥ d, where D denotes the minimum integer such that there exists constants c 1 , ..., c D such that Sup(ν)

⊂ d i=1
S ci . On the other hand if D < d, X is a singular random vector with a D-dimensional support (this is obviously the case when ν is atomic with less than d atoms).

Proof. For (i) and (ii), consider that X ∼ G d,n (α, s). Then there exists a matrix Y of gamma random variables Y i,j ∼ G 1,1 (α j , s i,j ) such that

X i = Y i,1 + . . . + Y i,n ,
where row vectors Y i,. are independent and column vectors Y .,j have two separate groups of marginals: a first comonotonous group (with gamma marginals) corresponding to scales s i,j that are strictly positive and a second group (with identically 0 marginals) corresponding to null scales. Note that one of these two groups might be empty. From this expansion, we can easily deduce the structure of maximal supports for ν that leads to independent or comonotonous random vectors X. Consider non-atomic Thorin measures as limits of the previous case.

For (iii), the problem of regularity of X and the existence of a density in L 2 , simply consider the rank of the (possibly non-square) matrix s such that X = sZ where Z j ∼ G 1,1 (α j , 1) are independent. If D = d = n, then X is an invertible linear transformation of an independent, d-dimensional random vector, and therefore has a density in L 2 that we can express easily. If D = d < n, then X is the convolution of the previous case with something else, and therefore also has a density in L 2 . Same goes for D > d.

Hence, the only comonotonous structure that is allowed inside the G d class is X i = c i X 0 for X 0 ∈ G 1 and c > 0. Between the independent and the comonotonous cases, there is a wide range of achievable dependence structures. Furthermore, since the class G d is closed w.r.t. convex convolution, every (positive) value of dependence measures, such as Kendall tau or Spearman rho, are achievable. Some examples in Section 4 also exhibit highly asymmetrical dependence structure and/or tail dependency.

For the study of the relation between the dependence structure of X and the dependence structure of ν, the Thorin measure proposed in Example 2.1.5 is of interest.

Example 2.1.5 (A curious distribution). Let ν be a one-dimensional measure defined through ν(dx) = 1 x∈[0,1] dx. Then G 1 (ν) has moment generating function:

M (t) = e 1-(t-1) t ln(1-t) Moreover, if for all n ∈ N, G j,n ∼ G 1,1 1 n , j n+1
are independent, then:

lim n→+∞ n j=1 G j,n ∼ G 1 (ν).
Proof. The second result comes from a discretization of the Thorin measure, from continuously uniform on [0, 1] to uniform on j n+1 , j ∈ {1, ..., n} . For the first result, note that ν satisfies the integration constraints given by Definition 2.1.1. Then, simply compute the cumulant generating function as:

K(-t) = - ∞ 0 ln (1 + rt) 1 r∈[0,1] dr = - 1 0 ln (1 + rt) dr = - 1 t 1+t 1 ln(x) dx = 1 - 1 + t t ln(1 + t).
This random variable has the nice property that its Thorin measure is uniform on [0,1]. It could be simulated more efficiently than by summing gammas through [START_REF] Bondesson | On Simulation from Infinitely Divisible Distributions[END_REF][START_REF] Ridout | Generating Random Numbers from a Distribution Specified by Its Laplace Transform[END_REF] if we wanted. However, its use is limited by the following negative result:

Remark 2.1.6 (The no-bijection result). Let X ∼ G d (ν) have marginal cdf F 1 , .
.., F d , and denote by G the cdf associated to the distribution of Example 2.1.5. Then, assuming that ν is absolutely continuous, the random vector

Y = G -1 (F i (X i )) i∈1,...,d has uniform marginal Thorin measures, but is not G d (c ν )-distributed, where c ν is the copula of ν.
There is in fact a bijection between the copula of the random vector and the copula of its Thorin measure, but conditionally on the marginal distributions. Therefore, an estimation scheme that separates the dependence structure from the marginals is not possible.

Before discussing the estimation of these distributions, we still need to introduce in the next subsection some specific integrals, which correspond to derivatives of the moment generating function and the cumulant generating function of a random vector.

Shifted moments and shifted cumulants

We formally introduce the shifted moments and the cumulants of a random vector, and present a known and useful result that maps moments to cumulants and vice-versa.

For X a random vector of dimension d, we recall the following notations for, respectively, the moment generating function M and the cumulant generating function K of X:

∀t ∈ C d , M (t) = E e t,X
and K(t) = ln M (t).

We now introduce specific notations for the derivatives of these two functions.

Definition 2.2.1 (Shifted moments and cumulants). For X a random vector of dimension d, for i ∈ N d , define the t-shifted i th moments and cumulants as the i th derivatives of (respectively) M and K, taken in t:

µ i,t = M (i) (t) and κ i,t = K (i) (t).
They correspond to the moment and cumulants of an exponentially tilted version of the random vector X, also called the Escher transform [START_REF] Escher | On the Probability Function in the Collective Theory of Risk[END_REF] of X.

Although standard Monte-Carlo estimators for shifted moments are unbiased and easy to compute, shifted cumulants are a little harder to estimate from i.i.d. samples of a random vector. We refer to [START_REF] Smith | A Recursive Formulation of the Old Problem of Obtaining Moments from Cumulants and Vice Versa[END_REF][START_REF] Smith | A Tutorial on Multivariate k-Statistics and Their Computation[END_REF] for some unbiased estimators of multivariate cumulants, and [START_REF] Krutto | Empirical Cumulant Function Based Parameter Estimation in Stable Laws[END_REF] for an application to the estimation of stable laws. To switch from moments to cumulants, Property 2.2.2 gives a bijection between the two: Property 2.2.2 (Bijection between moments and cumulants). For a given shift t ∈ C d and m ∈ N d , the sets (µ i,t ) i≤m and (κ i,t ) i≤m can express one from one another through:

µ i,t = e κ0,t π∈Π(i) B∈π κ i(B),t κ i,t = π∈Π(i) (|π| -1)!(-1) |π|-1 (µ 0,t ) -|π| B∈π µ i(B),t
where Π(i) = P(B(i)) is the set of all partitions of the multi-set

B(i) = d j=1 ij p=1
{j}, a partition π in Π(i) is constituted of multi-sets B composed of integers, and i(B) = (#{occurences of j in B}) j∈{1,...,d} .

Proof. This is a direct application of the multivariate Faà Di Bruno's formula [START_REF] Constantine | A Multivariate Faa Di Bruno Formula with Applications[END_REF][START_REF] Mishkov | Generalization of the Formula of Faa Di Bruno for a Composite Function with a Vector Argument[END_REF][START_REF] Nardo | A New Algorithm for Computing the Multivariate Faà Di Bruno's Formula[END_REF] on M = exp •K in one hand and K = ln •M in the other hand.

Remark that this property can be expressed in several other ways: most analysis in the literature uses multivariate Bell Polynomials [START_REF] Bernardini | Multidimensional Bell Polynomials of Higher Order[END_REF][START_REF] Cvijović | New Identities for the Partial Bell Polynomials[END_REF][START_REF] Kim | On Central Complete and Incomplete Bell Polynomials I[END_REF][START_REF] Melman | Realization of a Method for Calculating Bell Polynomials Based on Compositae of Generating Functions[END_REF] for this task. For computational purposes, there also exists some practical functions that give coefficients of these polynomials [START_REF] Withers | Expressions for Moments and Cumulants in Terms of Each Other Using Bell Polynomials[END_REF][START_REF] Hardy | Combinatorics of Partial Derivatives[END_REF], even for d > 1. However, in all generality, the number of coefficients to compute and store is exponential with m making this method quite unpractical, as the Bell polynomials are combinatorial in nature. To switch from cumulants to moments, the good property of the exponential function of being its own derivative allows a recursive approach, which is described in full details in [START_REF] Miatto | Recursive Multivariate Derivatives of e f (X1[END_REF]. We adapted the notation and simplified the indexes conventions from [START_REF] Miatto | Recursive Multivariate Derivatives of e f (X1[END_REF], which allowed, after adding a slack variable to handle the dimensionality index, for an even faster implementation that the one from [START_REF] Miatto | Recursive Multivariate Derivatives of e f (X1[END_REF]. The expression in Algorithm 1 below uses derivatives of the cgf and mgf as it is the case that we need, but the derivative of the exponential of any function could be computed through exactly the same algorithm, knowing the derivatives of the exponent.

Algorithm 1: Recursive computation of µ from κ.

Input: m ∈ N d , Shifted Cumulants (κ k,t ) k≤m Result: Shifted Moments (µ k,t ) k≤m Set µ 0,t = exp (κ 0,t ) Set µ k,t = 0 for all k = 0 foreach k : 0 = k ≤ m do Set d as the index of the first k i that is non-zero. Set p = k Set p d = p d -1 foreach l : l ≤ p do µ k,t += (µ l,t ) (κ k-l,t ) p l end end Return µ
Note that the main loop of Algorithm 1 must be done in the right order such that each µ k is computed before being used. This recursive formulation has stunning results and is at the moment the fastest way this computation can be done 3 . Unfortunately, no equivalently fast procedure can be found the other way around, as the derivatives of the log function are more complicated.

We now dive into the estimation of multivariate gamma convolutions. We highlight in the next subsection the univariate projection procedure from [START_REF] Furman | On Log-Normal Convolutions: An Analytical-Numerical Method With Applications to Economic Capital Determination[END_REF][START_REF] Miles | Risk Aggregation: A General Approach via the Class of Generalized Gamma Convolutions[END_REF].

Projection from G 1 to G 1,n

Suppose that we have a density in G 1 . Then, for a certain n ∈ N, the procedure provided by Miles, Furman & Kuznetsov in [START_REF] Miles | Risk Aggregation: A General Approach via the Class of Generalized Gamma Convolutions[END_REF][START_REF] Furman | On Log-Normal Convolutions: An Analytical-Numerical Method With Applications to Economic Capital Determination[END_REF], hereafter denoted MFK, gives expression for the shapes α and scales s of a G 1,n (α, s) that is exponentially close to the desired density. We will not describe the algorithm here, we refer to the original article for the details as the exposition is quite long, but we will comment on the results and caveats and use it as a comparison basis for later experiments.

A deep result about Padé approximation of Stieltjes functions supports the MFK procedure. Indeed, the first derivative of the cumulant generating function of a random variable in G 1 writes (using the rate notation r = 1 s and the corresponding push-forward ν r of the Thorin measure ν):

K (1) (-t) = 1 r + t ν r (dr).
When ν r is discrete with n atoms, this expression is a one-dimensional Stieltjes function. This Stieltjes function can be approximated from (µ i ) i≤2n by [m-1/m] standard Padé approximants, (see [START_REF] Derevyagin | On the Convergence of Padé Approximations for Generalized Nevanlinna Functions[END_REF]), whose denominator has zeros that are all simple, real, and negative. Applying a partial fraction decomposition to the Padé approximation is a direct way to estimate the rates and shapes of the gammas, since the obtained form for K (1) , and therefore K(t), matches the form given in Definition 2.1.1 for a G 1,n model.

The MFK algorithm has two major drawbacks. The first one is that the partial fraction decomposition of the Padé approximation, central in the procedure, is typical to univariate Stieltjes functions. Since the first derivative of the cgf is a Stieltjes function only in G 1 , the method cannot be used to approximate a random variable that is not in the G 1 class. Since equivalent results do not exist for multivariate Stieltjes functions, the method is restricted to one-dimensional problems. The second one is that precision of computations needs to be at 1024 bits or higher for the algorithm to produce a result at all.

We do have an explanation for each of these drawbacks. In all generality, the Padé approximation of the Stieltjes function is a special case of a more general problem, a moment problem for the Thorin measure. Indeed, remark that:

κ k,t = (k -1)! s k (1 + st) k ν(ds),
and use the push-forward x = s 1+st , mapping R + to the one-dimensional tsimplex ∆ 1 (t) = [0, 1 t ]. You obtain a moment problem of the form:

ξ k = ∆1(t)
x k ξ(dx) , ∀k ≤ m,

where ξ = (ξ k ) k≤m = κ k,t (k-1)! k≤m
is the vector of moments that the measure ξ needs to fit.

We can then obtain a solution with n atoms for the measures ξ and ν by solving this moment problem using classical tools, see, e.g., [START_REF] Schmudgen | The Moment Problem[END_REF] for an extensive discussion of potential algorithms. When the vector of cumulants κ comes from a density in the G 1 class, the vector of moments ξ of the measure ξ is inside the moment cone of positive measure on ∆ 1 (t). We therefore obtain a proper measure ξ matching these moments. On the other hand, if ξ is not a sequence of moments of some measure on the simplex, which happens as soon as the random variable is not in G 1 , then the moment problem becomes unsolvable. In this case, the approximation of the Stieltjes function fails to provide correct (within boundaries) atoms and weights for the measure ξ, and moreover for the Thorin measure.

Experiments from [START_REF] Miles | Risk Aggregation: A General Approach via the Class of Generalized Gamma Convolutions[END_REF][START_REF] Furman | On Log-Normal Convolutions: An Analytical-Numerical Method With Applications to Economic Capital Determination[END_REF] suggests that, if the density used to compute the shifted moments µ is in G 1 , the integrations of moments µ and the derivation of κ and ξ from µ need to be carried out with at least 300 digits of accuracy (about 1024 bits) so that ξ is in the desired cone of moments. From an empirical dataset, even if the true distribution lies in G 1 , evaluation of µ through Monte-Carlo with enough precision is therefore impossible.

Last but not least, the results about the complete factorization of Padé approximants denominator of a Stieltjes function, which simplifies the moment problem, are not true for d > 1. The approach could nevertheless be carried out, as the corresponding moment problem still has a solution if the used cumulants are true cumulants from a G d distribution, which is an even stronger precision requirement as the dimension increases. This possibility is discussed in the next section.

Construction of an estimator in G d,n

We attempt a multivariate extension of the MFK algorithm. In a multivariate setting, the Padé approximants of the cumulant generating function exists but has a denominator that is not known to be completely factorizable. Thus, we cannot take a partial fraction expansion to obtain the cumulant generating function in its canonical form to match the coefficients. Nevertheless, through the moment problem interpretation of MFK algorithm developed in Subsection 2.3, we can still express and try to solve the corresponding multivariate moment problem.

For X ∼ G d (ν), consider K the cumulant generating function of X, given in Definition 2.1.3 as:

K(t) = -ln (1 -s, t ) ν(ds).
Assuming |i| ≥ 1 and t ≤ 0, we have that

∂ i ∂ i t ln (1 -s, t ) = -s i (|i| -1)! (1 -s, t ) -|i| ,
and, therefore,

κ i,t = K (i) (t) = - ∂ i ∂ i t ln (1 -s, t ) ν(ds) = s i (|i| -1)! (-1) |i|-1 (1 -s, t ) -|i| ν(ds) = (|i| -1)! s i (1 -s, t ) |i| ν(ds). (1)
Note that t ≤ 0 is inside the region of convergence of the function, and use now the (bijective, continuous) push-forward:

x = s 1 -s, t ⇐⇒ s = x 1 + x, t , which goes from R d + to the simplex ∆ d (-t) = x ∈ R d + ,
x, -t ≤ 1 , yielding an equivalent problem of finding the measure ξ that solves the more standard moment problem:

∀ i ≤ m, ξ i = ∆ d (-t) x i ξ (dx) , (2)
where we denoted ξ = κ i,t (|i|-1)! i≤m the moments that ξ needs to match according to Equation [START_REF] Arnold | Möbius Transformations Revealed[END_REF].

Finding an n-atomic measure ξ solution to these equations is equivalent to finding the parametrization of a G d,n distribution that fulfills the cumulant constraints. The corresponding Generalized Moment Problem is known as a hard problem in the literature, but can be solved through semi-definite positive relaxations, following, e.g., [START_REF] Helton | A Semidefinite Approach for Truncated K-Moment Problems[END_REF][START_REF] Nie | The A-Truncated K-Moment Problem[END_REF][START_REF] Henrion | Projection Methods in Conic Optimization[END_REF]. Sadly, these algorithms and the corresponding literature focus on the first case, where ξ are, indeed, moments of a measure supported on the simplex ∆ d (-t), which is the case only if X ∈ G d,n and if ξ are computed with enough precision.

If the knowledge about the distribution of X is given through empirical observations, with sampling noise, or if the true distribution is not in the G d,n class, these moments will not belong to the moment cone, and the moment problem will have no solution. By projecting the moments onto the moment cone, we could force the known algorithms to provide a solution, but the projection onto this cone is not an easy task (see [START_REF] Di Dio | The Multidimensional Truncated Moment Problem: The Moment Cone[END_REF][START_REF] Di Dio | The Multidimensional Truncated Moment Problem: Atoms, Determinacy, and Core Variety[END_REF] for details on moment cones).

The optimization of an n-atomic Thorin measure from samples of the random vector X therefore rises several questions, both in univariate and multivariate cases. If we empirically compute cumulants, we will stand outside the moment cone and the moment problem will have no solution. Then, should we consider Equation (2) as a multi-objective optimization problem ? Do we seek a Pareto front, or can we consider that some moments are less important than others ? If we do, how do we weight the several objectives ?

To produce a meaningful loss to minimize, we will weight these moment conditions. For this purpose, we leverage a certain Laguerre basis of L 2 (R d + ) to expand the density of random vectors in G d,n , and to construct a coherent loss for this, potentially impossible to solve, moment problem.

The tensorized Laguerre basis

The standard Laguerre polynomials [START_REF] Mustapha | Generalized Laguerre Expansions of Multivariate Probability Densities with Moments[END_REF] form an orthogonal basis of the set L 2 (R + , e -x ) of square integrable functions with respect to the weight function x → e -x . From these polynomials, we can extract the following orthonormal basis of the set L 2 (R d + ) of functions that are Lebesgue square-integrable on R d + :

Definition 3.1.1 (Laguerre function). [START_REF] Comte | Adaptive Laguerre Density Estimation for Mixed Poisson Models[END_REF][START_REF] Mabon | Adaptive Deconvolution on the Non-Negative Real Line: Adaptive Deconvolution on R+[END_REF][START_REF] Dussap | Anisotropic Multivariate Deconvolution Using Projection on the Laguerre Basis[END_REF] For all k ∈ N d , define the Laguerre functions ϕ k with support R d + by:

ϕ k (x) = d i=1 ϕ ki (x i ) where for k ∈ N, ϕ k (x) = √ 2e -x k =0 k (-2x) ! .
For every function f ∈ L 2 (R d + ), for every k ∈ N d , coefficients of f in the basis are denoted by:

a k (f ) = R d + f (x)ϕ k (x) dx,
and we have, since the basis is orthonormal, that f (x) =

k∈N d a k (f )ϕ k (x) and that f 2 2 = k∈N d a 2 k (f ).
Furthermore, for X a random vector with density f ∈ L 2 R d + and shifted moments µ, we have by simple plug-in of f and ϕ k into the expression of a k that:

a k (f ) = √ 2 d ≤k k (-2) | | ! µ ,-1 . (3)
We now denote a k (α, s) the k th Laguerre coefficient of the G d,n (α, s) distribution.

Property 3.1.2 (X ∼ G d,1 (α, s)). For X ∼ G d,1 (α, s), we have a bijection between the d + 1 first Laguerre coefficients and the d + 1 parameters, given through the following equations:

a k (α, s) = √ 2 d ≤k k (-2s) ! Γ (α + | |) Γ (α) (1 + |s|) -α-| | ,
on one hand, and

α = 1 c 2 - 1 ln (c 1 ) W 0 ln (c 1 ) c 2 e ln(c 1 ) c 2 s i = a 0 -a 1(i) 2αa 0 for all i ∈ 1, .., d
reciprocally, where:

• 1(i) is the d-variate vector with j th component equal to 1 i=j , for all j.

• c 1 = a 0 √ 2 -d and c 2 = d 2 -1 2a0 d i=1 a 1(i)
• W 0 is the zeroth branch of the Lambert function W defined by:

y = W (x) ⇐⇒ x = ye y .
Proof. For the first part, note that there exists X 0 ∼ G 1,1 (α, 1) such that X = s, (X 0 , ..., X 0 ) . Therefore, denoting by M the mgf of X 0 ,

µ k,-1 = E X k e -1,X = s k E X |k| 0 e -|s|X0 = s k M (|k|) (-|s|) = s k Γ (α + |k|) Γ (α) (1 + |s|) -α-|k| .
Therefore, the Laguerre coefficients a k (α, s) of X are given by:

a k (α, s) = √ 2 d ≤k k (-2s) ! Γ (α + | |) Γ (α) (1 + |s|) -α-| | .
Then, for the reverse assertion, consider the expression of the first coefficients:

a 0 = √ 2 d (1 -|s|) and a 1(i) = a 0 1 -2α s i 1 + |s| .
Use the simplex transformation x = s 1+|s| to solve for s as a function of α in the second equation, and inject in the first. We obtain an equation of the form ae x +bx+c which is solved by the zeroth branch of the Lambert W function. This example is quite peculiar since by Property 2.1.4 the random vector is singular. Note that the edge case ∃i ∈ 1, .., d, r i = 0, in which the i th marginal is identically zero, is included in the previous result.

We do not know if the same kind of inversion could be carried out analytically for the coefficients of G d,n densities, since Laguerre coefficients are not as trivially expressed in all generality. Furthermore, this kind of parametric inversion would not work for empirical datasets or densities outside the class, and is therefore not our focus here. However, we can efficiently compute coefficients from parameters of a G d,n random vector through Algorithm 2, based on Algorithm 1:

x i = si 1+|si| for all i ∈ 1, ..., n. Let κ 0 = - n i=1 α i ln (1 -|x i |) Let a 0 = µ 0 = exp (κ 0 ) foreach 0 = k ≤ m do Let a k = µ k = 0 Let d be the index of the first k i that is non-zero. Let p = k Set p d = p d -1 Let κ k = (|k| -1)! n i=1 α i x k i according to Eq. (1) foreach l ≤ p do Set µ k += (µ l ) (κ k-l ) p l according to Algorithm 1 Set a k += µ l k l (-2) |l| l!
according to Eq. ( 3) end

Set a k += µ k (-2) |k| k! end a = √ 2 d a

Return a

The complexity of Algorithm 2 is quadratic in m. Note that, as in Algorithm 1, computations need to be performed in the right order. Sometimes, the Laguerre coefficients a = (a k ) k≤m overflow the Float64 limits, but the implementation we provide in the Julia package ThorinDistributions.jl4 , archived in Zenodo [29], ensures that the computations do not overflow by using multiple precision arithmetic when needed. Furthermore, as in Algorithm 1, the use of Miatto's fast recursion gives this algorithm a good efficiency, even for reasonably large d, n, m.

By reordering the coefficients and leveraging generalizations of Laguerre polynomials, we have the following link with Moschopoulos's density:

Remark 3.1.3 (Generalized Laguerre basis and Moschopoulos density). For a random vector X ∼ G d,n (α, s) with Laguerre coefficients a, the density f could be rewritten as a gamma series:

f (x) = k∈N d b k f k (x), where f k (x) = x k e |x| k! and b k = √ 2 d (-2) |k| ≥k k a .
We note that f k is the density of a vector with independent G 1,1 (k i , 1) marginals. Laguerre's polynomials can be generalized by adding a parameter ρ ≥ 0, which we took here to be 0. This alteration would reflect in the series by yielding densities of independent random vectors with G 1,1 (k i +ρ, 1) marginals instead. This hyperparameter could be varied to produce different concurrent models, and could even be chosen dimension wise. Noteworthy is the fact that, in the one-dimensional case, if we chose ρ to be the total mass of the Thorin measure, we retrieve Moschopoulos series from Property 2.1.2 as a Laguerre expansion.

Sadly, even if mixtures of gammas are easier to fit by k-MLE procedure [START_REF] Nielsen | k-MLE: A Fast Algorithm for Learning Statistical Mixture Models[END_REF][START_REF] Schwander | Fast Learning of Gamma Mixture Models with K-MLE[END_REF], we have no way of identifying the subspace of coefficients that would match a true convolution of gammas: as [START_REF] Bondesson | Generalized Gamma Convolutions and Related Classes of Distributions and Densities[END_REF] noted in the univariate case, generalized gamma convolutions can be expressed as mixtures of gammas, but there is no simple reverse mapping. Last but not least, if the expansion through generalized Laguerre polynomials with a parameter corresponding to the total mass of the Thorin measure would converge faster, we would have no way of estimating this parameter beforehand from data, and for many useful cases it would be infinite (log-Normals among others).

However, under simple technical assumptions on the parameters of a G d,n random vector, Laguerre coefficients are exponentially decreasing. Property 3.1.8 will provide such a bound. Before exposing the main result, we start by a few technical statements. Lemma 3.1.4 gathers some useful properties about a particular univariate Möbius transform.

Lemma 3.1.4 (Property of a Möbius transform). Let h be the Möbius transform:

h(t) = t + 1 t -1 .
For a ∈ C and b ∈ R + , we denote the disc with center a and radius b by:

D(a, b) = {t ∈ C, |t -a| < b} .
Also denote c(t) = t 2 +1 t 2 -1 and r(t) = 2t t 2 -1 , and finally denote C -and C + the left half and right half of the complex plane. The following holds:

(i) h is its own inverse, i.e., h(h(t)) = t, and h has a simple pole at t = 1 which has limit +∞ if Re(t) > 1 and -∞ if Re(t) < 1, t → 1.

(ii) h(1/t) = -h(t) for all t ∈ C and |h(t)| > 1 for all Re(t) > 0,

Im(t) < +∞ (iii) h(D(0, 1)) = C - (iv) If b < 1, h(D(0, b)) = D(c(b), r(b)) ⊂ C - (v) If b > 1, C \ h(D(0, b)) = D(c(b), r(b)) ⊂ C +
Proof. Point (i) and (ii) are obvious. To prove the 3 remaining points, we use the fact that h is a Möbius transform, see [START_REF] Arnold | Möbius Transformations Revealed[END_REF], and hence maps generalized circles to generalized circles. Furthermore, h maps reals to reals, and we obtain c(b) and r(b) from the expression of h(b) and h(-b). Last, by checking if the disk with center c(b) and radius r(b) contains h(0) = -1, we conclude.

Definition 3.1.5 provides a specific regularity assumption on a parametrization of generalized gamma convolution that is a little cumbersome to describe but will be helpful later.

(i) G d,n (α, s) is well-behaved. (ii) ∃ ε > 0 such that G d,n (α, s) is ε-well-behaved.
(iii) |α| > 1 and ∀I ⊆ {1, . . . , n} such that i∈I α i > i / ∈I α i , Ker(s I ) = {0}.

(iv) ν(R d + ) > 1 and there exists no c 1 , . . . , c p linearly independent, p < d, such that

ν p i=1 S ci > ν(R d + ) 2 ,
where S c was defined in Property 2.1.4.

Proof. First, (i) ⇐⇒ (ii) is just the definition of well-behavior. Furthermore, (iv) is equivalent to (iii): a subset of scales s I such that i∈I α i > i / ∈I α i is a subset of scale that has the majority of the weight of the Thorin measure assigned to it, and Ker(s I ) = 0 if and only if at least d linearly independent directions are inside s I . We therefore show only (ii) ⇐⇒ (iii).

By Rouché-Capelli's Theorem, see [START_REF] Suetin | Linear algebra and geometry[END_REF], any complex solution (if there are some) of the linear system s I t = 1 is unique if and only if Ker(s I ) = {0}. If it exists, call this solution t * ∈ C d . Now, since for a given i ∈ I, s i,j are all non-negative reals and one of them must not be zero, t * must have at least one component t * j that has a positive real part. Furthermore, this component is clearly of bounded modulus since the solution is unique and all s i,j are finite.

examples).

As soon as the total mass of their Thorin measure is greater than one, the following gamma convolutions are wellbehaved:

(i) Any univariate gamma convolution (the best ε is given in Definition 3.1.5).

(ii) Any independent random vector in G d,n .

(iii) Any invertible linear transformation of an independent random vector in G d,n .

(iv) Any finite convolution of well-behaved gamma convolutions.

Proof. The three first cases are deduced directly from Definition 3.1.5 and Property 3.1.6. For point (iv), recall that for two gamma convolutions with Thorin measures ν 1 and ν 2 , the Thorin measure of the sum is:

ν(s) = ν 1 (s) + ν 2 (s). Now,
|a k | ≤ B(d, ε )(1 + ε ) -|k| .
Proof. The sketch of the proof is as follows: we construct a multivariate complex function R that has the Laguerre coefficients as its Taylor coefficients, and we show by a singularity analysis that this function is analytic around the origin, in a polydisc with polyradius 1 + ε. We conclude by applying Cauchy's inequality.

As the one-dimensional case allows for more detailed computations, we start by working for all d, and split later into the two cases d = 1 and d > 1.

We follow the path of [START_REF] Kabardov | Acceleration of the Convergence of the Laguerre Series in the Problem of Inverting the Laplace Transform[END_REF] to express the Laguerre coefficients (a k ) k∈N d as the Taylor coefficients of a certain (d-variate, complex) function R. From the computation of the moment generating function M (t) of the Laguerre expanded density f (x) = k∈N d a k ϕ k (x), we have:

M (t) = R d + e t,x f (x) dx = √ 2 d k∈N d a k d i=1 ki j=0 k i j (-2) j j! R+ e (ti-1)xi x j i dx i = √ 2 d k∈N d a k d i=1 ki j=0 k i j (-2) j (1 -t i ) j+1 = √ 2 d (1 -t) -1 k∈N d a k 1 - 2 1 -t k = √ 2 d (1 -t) -1 k∈N d a k t + 1 t -1 k ,
which implies that:

k∈N d a k t + 1 t -1 k = √ 2 -d (1 -t) 1 M (t).
Now denote h(t) = (h(t)) t∈t , so that h applies the Möbius transform h(t) = t+1 t-1 from Lemma 3.1.4 componentwise. Using the substitutions y = h(t) ⇐⇒ t = h(y) (since h is its own inverse), which implies that (1 -t)

1 = 2 d (1 -y) -1
, we define the function R as:

R(y) := k∈N d a k y k = √ 2 -d M (h(y)) (1 -h(y)) 1 = √ 2 d M (h(y)) (1 -y) -1 .
We now study the regularity of the function R, which is tightly related to the regularity of M . Recall first the expression of the moment generating function for our random vector:

M (t) = n i=1 (1 -s i , t ) -αi
Denote by V R the singular variety of R, i.e., the set of points where R is not analytic. From the singularities of y → (1 -y) -1 and from those of M , we have:

V R ⊆ {y : ∃j, y j = 1} n i=1 {y : s i , h(y) = 1} .
Consider first the fact that, since the model is ε-w.b., |α| > 1. Hence, M (h(y)) dominates (1 -y) -1 when y goes to 1 (as h(y) goes to infinity) and therefore R is analytic at 1. We therefore only need to consider the singularities of M (h(y)), and we can further restrict V R :

V R ⊆ n i=1
{y : s i , h(y) = 1} .

We now split the argument into the cases d = 1 and d > 1, to show that the ε-w.b. condition is equivalent to the analyticity of R on the polydisc with polyradius 1 + ε.

When d = 1, all singularities of R are such that sh(y) = 1 for a scale s ∈ s. In other words, the singularities are at the points y = h( 1 s ) for each scale s ∈ s. Recall that -h(x) = h(x -1 ), and that h is (strictly) decreasing before and after its vertical asymptote at 1. Three cases arise:

• If s ∈]0, 1[, then |y| = |h 1 s | = h 1 s > 1 + ε ⇐⇒ 1 s < h(1 + ε) ⇐⇒ s > h(1 + ε) -1 = ε 2+ε . • If s ∈]1, +∞[, then |y| = |h 1 s | = h(s) > 1 + ε ⇐⇒ s < h(1 + ε) = 2+ε ε . • If s = 1, h 1 s = ±∞ and there is no singularity. Hence, min y∈V R |y| > 1 + ε ⇐⇒ ε 2 + ε < s < 2 + ε ε ,
which corresponds (with |α| > 1) to the ε-w.b. condition.

We now turn ourselves to the multivariate case. When d > 1, the singularities of t → M (t) are the zeros of the function:

G : t → n i=1 (1 -s i , t ) αi ,
By the definition of a ε-w.b model, for any t * a zero of G, we have that:

∃i ∈ 1, ..., d, |h(t * i )| > 1 + ε.
Hence, the singularities of R, which are images by h of the zeros of G, must have at least one dimension that is outside the centered polydisc with radius 1 + ε. Hence, when d > 1, R is also analytic on the centered polydisc with polyradius 1 + ε.

Now, for any d, by Cauchy's inequality, see Theorem 2.2.7 in [START_REF] Hormander | An Introduction to Complex Analysis in Several Variables[END_REF], taking a ε strictly between 0 and ε we have the wanted finite upper bound on the coefficients:

B(d, ε ) = sup (α,s):G d,n (α,s) is ε-w.b. ∀i,|yi|<1+ε R(y) < ∞.
Since we can now compute efficiently the Laguerre coefficients of densities in G d,n , with an insurance that they would decrease fast if the model is wellbehaved, we propose in the next subsection to finally discuss an estimation algorithm based on an Integrated Square Error loss that leverages these coefficients.

An incomplete but however consistent empirical loss

Expressing the density of G d,n random vectors into the basis from Definition 3.1.1 allows us to truncate the basis and effectively compute an approximated density. From this, by evaluating empirical Laguerre coefficients from data, we will fit the parameters to data by minimizing the L 2 distance between coefficients (since the basis is orthonormal). For parameters α, s, denote by f α,s the density of the G d,n (α, s) distribution, and f the density underlying the observations x 1 , ..., x N .

We would like to minimize the integrated square error:

f -f α,s 2 2 = k∈N d (a k (f ) -a k (α, s)) 2 .
Estimating a k (f ) = f (x)ϕ k (x) dx by a simple Monte-Carlo plug-in

a k = 1 N N i=1 ϕ k (x i ),
we could compute an approximation of this loss. However, we cannot compute all couples of coefficients a k and a k (α, s) for all k ∈ N d , so that we are forced to truncate the basis.

Definition 3.2.1 (G d,n parameters estimator). For a dataset x = (x 1 , ..., x N ), we define parameters estimator of a well-behaved G d,n distribution that matches the observations as:

( α, ŝ) = arg min G d,n (α,s) w.b. L m (x, α, s),
where the approximated truncated loss L m is given by:

L m (x, α, s) = k≤m ( a k -a k (α, s)) 2 ,
and where the empirical Laguerre coefficients are given by:

a k = 1 N N i=1 ϕ k (x i ).
The loss L m (x, α, s) can be efficiently computed through Algorithm 2. It will be 0 if the first coefficients of the Laguerre expansion of f and of our estimator match perfectly, i.e., assuming f ∈ G d,n and m big enough, if we found the right Thorin measure. Note that if the goal was given theoretically, like in a projection from G d to G d,n , we could compute a k by formal integration with high precision instead of Monte-Carlo. We will discuss this option in the next section. However, even in this case, the error that comes from the truncation of the basis cannot be computed.

To show the consistency of this loss, we use in the following the results from Property 3.1.8.

Property 3.2.2 (Consistency). Consider that f ∼ G d,n (α 0 , s 0 ) is well-behaved,
and denote x a dataset generated from this density with N observations. Fix m ∈ N d . Suppose that we have a global minimizer

( α, ŝ) = arg min G d,n (α,s) w.b. L m (x, α, s),
which depends on the threshold m and on the N observations x. Then:

f -f α,ŝ 2 2 a.s ----→ N →∞ m→∞ 0.
Proof. To show the result, we start by expressing the loss in the Laguerre basis, and we split the basis on indices smaller and greater than m. For the sake of simplicity, we will omit f in a k (f ) and denote a k = a k (f ).

f -f α,ŝ 2 2 = k∈N d (a k -a k ( α, ŝ)) 2 = k≤m (a k -a k ( α, ŝ)) 2 A + k>m (a k -a k ( α, ŝ)) 2 B
We first show that A a.s.

----→ N →∞ 0 for any fixed m. For A, we can treat each term in the summation independently as the sum is finite. A can be then expanded further by plugging-in the Monte-Carlo estimator

a k = 1 N N i=1
ϕ k (x i ) and expanding the square:

A = k≤m (a k -a k ) 2 A1 + 2 k≤m (a k -a k ) ( a k -a k ( α, ŝ)) A2 + k≤m ( a k -a k ( α, ŝ)) 2 Lm(x, α,ŝ) . Now, since a k a.s.
----→ N →∞ a k by the strong law of large numbers, A

----→

N →∞ 0.
Furthermore, since we restricted the optimization to well-behaved estimators, each a k ( α, ŝ) is bounded by Property 3.

1.8. Now | a k | ≤ √ 2 d because |ϕ k (x)| ≤ √ 2 d whatever x ∈ R d + . Therefore, for any k, ( a k -a k ( α, ŝ)) is bounded and (a k -a k ) ----→ N →∞ 0, which makes A 2 a.s. ----→ N →∞ 0.
Last but not least, by definition of ( α, ŝ) as global minimizers,

L m (x, α, ŝ) ≤ L m (x, α 0 , s 0 ) a.s. ----→ N →∞ 0,
since Monte-Carlo estimators are not biased, and (α 0 , s 0 ) is inside the range of optimizations. Therefore, we can conclude that A a.s.

----→ N →∞ 0 for any m fixed. Now consider the remainder B. We show that B a.s.

----→ m→∞ 0 uniformly in N . Since both the true model and the approximation are well-behaved, there exists ε > ε > 0 so that they are both ε-w.b. By Property 3.1.8 we have then,

|a k | ≤ B(d, ε )(1 + ε ) -|k| and |a k ( α, ŝ)| ≤ B(d, ε )(1 + ε ) -|k| .
Therefore:

B = k>m (a k -a k ( α, ŝ)) 2 ≤ 4B(d, ε ) 2 k>m (1 + ε ) 2 -|k| a.s.
----→ m→∞ 0, uniformly in N , which concludes the argument.

Although the consistency result restricts the parameter space to only wellbehaved Thorin measures, in practice we run our optimization procedures without this restriction, and simply check at the end that the resulting estimators are well-behaved. When d = 1, we can always fix ε to make the estimated gamma convolution ε-w.b. if | α| > 1. In practice, numerical computations have troubles producing collinear scales, so every estimated multivariate gamma convolution is usually well-behaved as soon as the total mass of the produced Thorin measure is greater than 1. Recall that univariate distributions like log-Normal and Pareto have an infinite total mass for the Thorin measure.

We now have a reliable loss to minimize and an efficient algorithm to compute it, allowing us to estimate multivariate gamma convolutions from data. To test the approach, we propose in the next section to discuss some potential applications.

Investigation of performance

In this section, we show the results of our algorithm on several examples, both in univariate and multivariate settings.

Note that, through the recursive Faà di Bruno's formula in Algorithm 2, we produce Laguerre coefficients as very high degree multinomials into α, µ 0 and ξ (the simplex projection of the scales). The loss will therefore have a myriad of local minima, making it highly non-convex and forcing us to use global minimization routines. We found the Particle Swarm [START_REF] Koyuncu | A PSO Based Approach: Scout Particle Swarm Algorithm for Continuous Global Optimization Problems[END_REF][START_REF] Li | Experiment of Standard Particle Swarm Optimization 2011 in[END_REF] global optimization routine to work particularly well on our test sets. All the implementation was gathered in the provided Julia package, along with the code corresponding to this investigation.

We start by the projection of known densities onto the G 1,n class, differentiating the cases when the known density is inside G 1 or not. Then we discuss the estimation of G 1,n models from empirical data, with heavy tailed simulated data, and we finally look at the estimation of G d,n models from multivariate empirical data, with a particular emphasis on the properties of the dependence structure underlying the data.

Projection from a known density

The two first examples are a log-Normal and a Weibull, that have the particularity of being respectively inside and outside G 1 . We will also take a look at Although MFK has a very strong convergence result, it seems that for n = 2 gammas, our estimator performs better with the same information (theoretical shifted moments given with 1024 bits precision), and that for n > 2, the performance is overall the same. Indeed, if we look at estimated parameters, we see that in most cases the two algorithms found approximations that are close to each other. We reported the shapes and scales that both algorithms produced for n = 2, 3, 4 and 5 in Table 1. Although the produced estimations are close to each other on this example, our method has several clear benefits: the full density does not need to belong to G 1 , and we do not require shifted moments up to 300 digits precision. In the next subsection, we describe the same experiment on a Weibull density that is known to make MFK simply fail, even with accurate enough computations.

4.1.2. Projection from densities outside G 1 to G 1,n
Let X be a Weibull r.v. with shape k > 0. Its pdf can then be written as:

f (x) = kx k-1 e -x k 1 x≥0 .
When k > 1, this distribution does not belong to G 1 . In [START_REF] Miles | Risk Aggregation: A General Approach via the Class of Generalized Gamma Convolutions[END_REF][START_REF] Furman | On Log-Normal Convolutions: An Analytical-Numerical Method With Applications to Economic Capital Determination[END_REF], authors report negative shapes for the MFK approximation in G 1,2 of the Weibull with shape k = 3 2 , which is clearly a failure. This failure is occurring because the truncated moment problem that MFK is solving has no solution, as we detailed in Section 2.3, since the (normalized) cumulants ξ are not inside the moment cone.

Instead, the Laguerre basis is an orthonormal basis of L 2 (R + ), which contains this density. The projection procedure detailed at the top of Subsection 4.1 still works correctly: it produces meaningful shapes and scales that are positive, whatever the number of gammas we ask for. Figure 2 displays a Violin plot of KS distances for several numbers of gammas: In abscissa is n ∈ {2, 3, 4, 5, 10, 20, 30, 40}, the number of gammas used in the approximation. Only the Laguerre approximation results are presented.

On Figure 2, we see that the KS distance is decreasing for an increasing number of gammas n, when n is small, but reaches a minimum quite quickly. Indeed, the distance is lower-bounded since the target is not in the class. Although none of the proposed approximations passes the KS test, we will revisit this example later from an estimation point of view to show that the projection is still an acceptable approximation of the input Weibull distribution.

We also note that the KS distance is not stable for increasing n: the algorithm has troubles to fix parameters to 0 if he does not need more precision. An additional penalty term in this direction might be a good solution to overcome this behavior.

Through these two first examples, we showed the projection possibilities of our algorithm, and compared them with MFK, the only existing method in the literature. We saw that for input densities inside the G 1 class our method performs as good as MFK. Moreover, our method still works properly and produces meaningful results for examples outside the class, contrary to MFK which does not give any results. We will now talk about estimation from data, a problem that has never been solved for gamma convolutions, and showcase how our algorithm can find a good representation for a given dataset, whatever the dimension of the dataset.

In the following, we will treat empirical datasets given in standard IEE 745 float64 precision.

Estimation from univariate empirical data

Before coming back in more details to the Weibull case, we wanted to perform some tests with empirical data coming from simulations of two heavy-tailed distributions: log-Normal and Pareto.

Heavy-tailed examples

We consider first a dataset of 100 000 samples from a log-Normal distribution with the same parameter values as in Section 4.1, LN (µ = 0, σ = 0.83), in 64 bits precision.

From this sample, we compute empirical Laguerre coefficients, and then we optimize the set of shapes and scales of the estimated density to minimize the truncated L 2 loss between them and the ones produced through Algorithm 2. Following [START_REF] Furman | On Log-Normal Convolutions: An Analytical-Numerical Method With Applications to Economic Capital Determination[END_REF][START_REF] Miles | Risk Aggregation: A General Approach via the Class of Generalized Gamma Convolutions[END_REF], we choose the size of the basis as being m = 2n + 1, such that if we remove the first moment that is irrelevant, we have a number of moments equal to the number of parameters5 . We ran the experiment for different number of gammas n = 10, 20 and 40. The corresponding results are respectively available in Figures 3, 4 and5. We see on these graphs that the approximation is good enough to fool the Kolmogorov-Smirnov test. The absolute pointwise error between the Laguerre approximation of the density of the gamma convolution and the true log-Normal density is lower than 0.005 for n = 10, 20, and twice smaller for n = 40. A positive but surprising result is that, although the distance between the Laguerre approximation and the true density does not reduce with increasing m, the distance between the gamma approximation and the true density does.

The KS test setup is as follows: we simulate 250 resamples of 100 000-sized dataset from the estimated density, and compare it with a one-sample exact KS test to the true density (and not to the samples that the estimator was computed from). For all values of n, we observe close-to-uniform distribution of the p-values, which is a good result. We also treated several datasets simulated from Pareto's distribution, with exactly the same experiment. Pareto's distribution has a shape parameter k > 0 that influences the thickness of the tail: the distribution has a variance if and only if k > 2, an expectation if and only if k > 1 and the density belongs to L 2 if and only if k > 1/2. Note that these inequalities are strict.

The fact that the distribution has no variance or no expectation is not really a problem for our procedure, since we use only shifted moments that all Pareto distributions have, but our convergence result requires a square-integrable density. Therefore, we propose to check several values of the shape parameter k, covering all cases, while varying the number n of gammas. Figure 6 summarizes the results by only showing quantile-quantile plots of the estimated distributions against the starting Pareto distribution. All the models follow the same procedure as for the previous log-Normal estimators. Note that even quantile-quantile plots between two simulations from the same Pareto distribution have a tendency to deviate in the extremes (due to the fat tail), which is why we excluded the last points. In each case, we also ran the same Kolmogorov-Smirnov testing procedure as for the log-Normal, and obtained very good uniformity of the p-values, except for k = 0.25 when n is too small. We see on the graph that increasing the number of gammas usually produces better results, but also that the heavier the tail (smaller shapes) the harder it is to get a good fit. Note that the KS distance is not the best distance to take into account the tails, which could explain why the distances from Figure 1 did not decrease for high values of n.

We see on these graphs that the method is accurate enough for reasonably high quantiles, and very good in the body of the distribution. The k = 0.25 case is not in the theoretical boundaries for the convergence result to take place, and from sample to sample we have to wait until at least 20 gammas to have an acceptable fit.

An example outside the class

We revisit the Weibull(k = 1.5) case of Section 4.1, a distribution that is not inside the class.

We ran our algorithm on 100 000 samples from a Weibull with shape k = 3 2 , with successive numbers of gammas of 2, 10, 20 and 40. The corresponding graphs are shown in Figure 7. Here, we show only the difference on the density estimation and the quantile-quantile plot as one notable thing is that the Laguerre expansion of the density gets better and better as we increase the size of the basis (jointly with the number of gammas), but the estimator does not: the difference to the density is always the same, and always makes Kolmogorov-Smirnov reject the proposed estimator, although quantile-quantile plots are all OK. On these estimations, several things should be noted:

• The KS tests reject the hypothesis that the distributions are the same whatever the number of gammas, which was expected.

• The pointwise difference on the density of the estimator does not match the error of the Laguerre expansion (which is getting smaller as the basis grows) but stays the same whatever the size of the basis and the number of gammas. This error represents some kind of distance between the Weibull distribution and the class.

• The quantile-quantile plot is on the other hand getting better and better as we add more gammas.

These three points show that our procedure works even for data that are far from the G 1 class, finding some orthogonal projection, relative to the basis, onto the class. This is a very good point for the procedure.

In the next examples, we will showcase some multivariate uses of the procedure, with a particular emphasis on the dependence structure of the simulated data and of the produced estimators.

Estimation from multivariate empirical data

Now that we extensively discussed the univariate case, we turn ourselves to the multivariate capabilities of our algorithm. As no known algorithm estimates multivariate gamma convolutions, we are not able to compare our results to other procedures. However, we can still check quantile-quantile plots for marginals and verify that the shape of the dependence structure matches the inputted one.

As use-cases of our models might vary, we propose two different examples: a multivariate log-Normal, with a Gaussian copula that induces no tail dependency, and a (survival) Clayton copula, including a strong tail dependency. Note that we reduced the exposition to estimation from empirical data, but formal integration to obtain (µ k,-1 ) k≤m from a formal density at any desired computational precision is also a perfectly fine use case. However, we do not have a clever way of choosing the parameter m yet. We found that setting all m i 's equal to n provide enough evaluation points to obtain a good fit on our simple examples.

A Gaussian example

We do not know if the d-variate log-Normal distribution, defined as the componentwise exponentiation of a d-variate Gaussian random vector, is inside G d . However, we can still run our procedure and observe the resulting approximations.

A first result is given in Figure 8, that corresponds to a G 2,10 estimation taken from 100 000 data points simulated from a M LN (µ = 0, σ 1 = σ 2 = 1, ρ = 0.5) with two standard log-Normal margins and a Gaussian copula with ρ = 0.5. Figure 8 proposes quantile-quantile plots for the marginals, and a simple Gaussian kernel density estimator for the pseudo-observations on the original and estimated models. We see on Figure 8 that the dependence structure was correctly reproduced by the approximation as a convolution of 10 (comonotonous) gamma random vectors, but marginal tails are not estimated correctly. Thankfully, pushing the number of multivariate gammas to n = 20 solves this issue, as shown in Figure 9: On Figure 9, we see that the dependence structure is still very good, and that the problem of the marginal tails has been solved by adding a few more gammas.

Note that since we have 100 000 points in the quantile-quantile plots, one point out of the line corresponds to a deviation at the 0.00001-quantiles, which are of little importance in many use cases.

Overall, marginal quantile-quantile plots are correct when the number of gammas is big enough, and the Gaussian dependence structure is well fitted. The Gaussian copula is a copula that exhibits no tail dependency. We already know that our model can exhibit asymmetric behaviors, but the possibility of tail dependency property is also an important feature to have. Therefore, in the next example we run the model on datasets that have such a feature.

A copula with tail dependency

To convince ourselves that tail dependency is a property that our estimator can exhibit, we will try to fit data whose dependence structure is given through a (survival) Clayton copula. This copula is a symmetric copula that exhibits upper tail dependency. We took here a parameter θ = 7, yielding a very high upper tail dependence.

In Figure 10, you can see the results corresponding to data sampled from a Clayton copula, a Pareto(k = 2.5) marginal and a LN (0, 0.83) marginal. The conditions of the experiments are identical to the previous subsection ones. We note that marginal upper tails are not perfectly fitted: this could probably be overcome by taking a greater number of gammas. However, we see that the shape of the estimated copula is not perfectly symmetric: there is no constraint on the model to produce a symmetric dependence structure, and therefore producing one is hard. The small misfit around the (0.25, 0.02) quantile is probably due to this: marginals are not on the same scale, and therefore the Laguerre coefficients with k 1 k 2 are probably on a different scale than those with k 2 k 1 . The goodness of fit is therefore not symmetric.

Conclusion

Although the class of gamma convolutions has very appealing properties for the practitioner, the estimation of gamma convolutions from empirical data had not been considered in the literature. The closest algorithm we could find is the recent MFK algorithm, that we described in Subsection 2.3, which projects densities from the G 1 class to the G 1,n class. Unfortunately, extending this algorithm for general estimation purposes is hard if not impossible.

Using a well-chosen Laguerre basis, we constructed a series expansion for the density of a distribution in G d,n . Coupled with a simple L 2 loss for the density, we were able to project densities from G d , d > 1 onto G d,n , which was not possible, but also to project any density onto G d,n , which was not possible either. Finally, through the same loss, we can estimate G d,n models for any reasonable d, n, which was also not treated in the literature.

Although the convergence result is not as strong as MFK's one, the precision needs and computational speeds are far better. We saw that given the same information, shifted moments of the random variable with high precision, the performance of the two algorithms is overall the same, except maybe for small numbers of gammas where we beat MFK. For both algorithms, Figure 1 also shows that the model's goodness of fit seems to be decreasing for large numbers of gammas, and an adaptive estimator might be a suitable way out.

Maybe a better loss could be found for the estimation of these gamma convolutions. Furthermore, automatic penalization of the model could be done through the number of gammas, but also through sparsity of the Thorin measure, according to Property 2.1.4, yielding more independent factors when possible. Such a modification will also provide a clearer view at the additive risk factor structure.

Last but not least, many features are still to be discovered about multivariate Thorin classes. The estimation of multiplicative structures could yield interesting results. Bondesson showed that the product of random variables in G 1 is also in G 1 , but we still do not have a corresponding result for the multivariate case.

  x) and scalar products of vectors by x, y = x 1 y 1 + . . . + x d y d . Finally, we denote x y = i xi yi the product of binomial coefficients.

ν

  j (s j )1 s k =0 ∀ k =j where ν j denotes the Thorin measure of the j th marginal, and the total masses of ν, ν 1 , ..., ν d are such that ν(R d + ) = ν(S ⊥ ⊥ ) = d i=1 ν i (R + ). (ii) The marginals of X are comonotonous if and only if there exists a constant c ∈ R d + , c = 1 such that

Algorithm 2 :

 2 Laguerre coefficients of G d,n (α, s) random vectors. Input: Shapes α ∈ R n + , scales s ∈ M n,d (R + ), and truncation threshold m ∈ N d Result: Laguerre coefficients a = (a k ) k≤m of the G d,n (α, s) density Side-product: -1-shifted cumulants κ and moments µ op to order m of the G d,n (α, s) distribution. Compute the simplex version of the scales

  An equivalent statement when d = 1 is simply that |α| > 1 and s ∈] ε 2+ε , 2+ε ε [ n If there exists ε > 0 such that the model is ε-w.b., we say that the model is well-behaved, in short w.b.Of course, any ε 1 -w.b. gamma convolution is ε 2 -w.b. for any 0 < ε 2 < ε 1 . Property 3.1.6 gives equivalent statements for the well-behavior of a gamma convolution, hopefully clarifying the requirements from Definition 3.1.5. Property 3.1.6 (Well-behaved gamma convolution). For a gamma convolution G d,n (α, s) with Thorin measure ν, the following statements are equivalent:

Finally, the requirementExample 3 . 1 . 7 (

 317 |h(t * j )| > 1+ε is equivalent, through Lemma 3.1.4, point (v), to |t * j -c(1 + ε)| < r(1 + ε) (where c and r were defined in Lemma 3.1.4). Since the disk D(c(1 + ε), r(1 + ε)) is in C + , and tends to the whole C + as ε goes to 0 by Lemma 3.1.4, point (iv), there always exist a ε > 0 so that a bounded t * j with positive real part is inside this disk. The proof of Property 3.1.6 gives a way of finding the greatest ε that makes a well-behaved gamma convolution ε-w.b. for the general multivariate case. In a certain sense, the ε-w.b. property is a quantification of the statements from Property 2.1.4: the Thorin measure must not assign half of its mass or more to any subset that would correspond to a non absolutely continuous gamma convolution. A surprise comes in the next examples, where we show the unexpected fact that the w.b. subclass of gamma convolutions is closed w.r.t finite convolutions. Simple w.b.

Fig 1 .

 1 Fig 1. Violin densities of resampled KS distances (smaller is better). Blue: MFK projection, Orange: Laguerre Projection. In abscissa is n ∈ {2, 3, 4, 5, 10, 20, 30, 40}, the number of gammas used in the approximation.

Fig 2 .

 2 Fig 2. Violin densities of resampled KS distances for approximation of a Weibull(k = 1.5).In abscissa is n ∈ {2, 3, 4, 5, 10, 20, 30, 40}, the number of gammas used in the approximation. Only the Laguerre approximation results are presented.

Fig 3 .

 3 Fig 3. Log-Normal results with 10 gammas: Upper left: the comparison of the density approximations. Lower left: the difference of to the true density. Upper right: a quantile-quantile plot of the approximation against the true distribution. Lower right: p-values of one-sample KS tests of simulations from the estimator against the true distribution.

Fig 4 .

 4 Fig 4. Log-Normal results with 20 gammas. Same legend as Figure 3.

Fig 5 .

 5 Fig 5. Log-Normal results with 40 gammas. Same legend as Figure 3.

Fig 6 .

 6 Fig 6. Quantile-quantile plots for Pareto experiments, with N = 1000 samples (log-scaled).We only show the 995 first points: 5 points are excluded in the tail for clarity. Each row corresponds to a shape parameter for the Pareto, and each column to a number of gammas.

Fig 7 .

 7 Fig 7. Weibull results with 2, 10, 20 and 40 gammas, respectively on each row. Same legend as Figure 3.

Fig 8 .

 8 Fig 8. Multivariate log-Normal results with 10 gammas. Top left and bottom right: quantilequantile plots for, respectively, the first and second marginals. Bottom left and top right: levels plots of pseudo-observations from, respectively, simulations from the estimated model and the dataset.

Fig 9 .

 9 Fig 9. Multivariate log-Normal results with 20 gammas. Same legend as Figure 8.

Fig 10 .

 10 Fig 10. Results with 20 gammas. Same legend as Figure 8.

  ∈I α i , denoting s I = (s i ) i∈I , for any solution t * of the linear system of equations s I t = 1, at least for one j ∈ 1, ..., d, |h(t *

Definition 3.1.5 (ε-well-behaved gamma convolution). A gamma convolution G d,n (α, s) is said to be ε-well-behaved, in short ε-w.b., if |α| > 1 and ∀I ⊆ {1, ..., n} such that i∈I α i > i / j )| > 1 + ε.

  according to the reading of Property 3.1.6, point (iv), no set of p < d rays S c1 , ..., S cp concentrate more than half of ν 1 weights and more than half of ν 2 weights. Hence, no set of p < d rays can concentrate more than half of ν weights. Therefore, the class of w.b. gamma convolutions is closed w.r.t finite convolutions.Property 3.1.8 (Exponential decay of Laguerre coefficients). For any ε, ε s.t. ε > ε > 0 and any dimension d, there exists a finite positive constant B(d, ε ) such that the Laguerre coefficients (a k ) k∈N d of any d-variate ε-wellbehaved gamma convolution verify:

	An interesting fact is that any well-behaved gamma convolution has exponen-
	tially decreasing Laguerre coefficients, as the following property shows:

Table 1

 1 Estimated parameters from both algorithms on the log-Normal example

		MFK		Laguerre
		α	ŝ	α	ŝ
	n=2	0.5688	1.5622	0.5458	1.6283
		2.4596	0.1941	2.4539	0.1999
	n=3	0.2195	2.4902	0.2070	2.5781
		0.8934	0.6718	0.8919	0.6875
		2.6081	0.0972	2.6071	0.0987
	n=4	0.0919	3.4076	0.0844	3.5307
		0.4596	1.2208	0.4555	1.2513
		1.0042	0.3737	1.0063	0.3792
		2.6956	0.0588	2.6957	0.0594
	n=5	0.0398	4.3466	0.0346	4.5447
		0.2569	1.7963	0.2492	1.8576
		0.5708	0.7234	0.5721	0.7394
		1.0576	0.2399	1.0609	0.2428
		2.7574	0.0396	2.7582	0.0399

ThorinDistributions.jl: https://github.com/lrnv/ThorinDistributions.jl. O. Laverny et al./Estimation of multivariate generalized gamma convolutions

Note that in[START_REF] Bondesson | Generalized Gamma Convolutions and Related Classes of Distributions and Densities[END_REF], distributions were parametrized through rates β = 1 s instead. We chose here to parametrize by scales, as it simplifies some of our notations.

The complexity is not exponentially increasing in the size of the arrays, on the contrary to naive implementations of Faà Di Bruno's formula, see[START_REF] Miatto | Recursive Multivariate Derivatives of e f (X1[END_REF].

ThorinDistributions.jl: https://github.com/lrnv/ThorinDistributions.jl.

Since Julia starts counting at 1 and not 0, we actually have m moments and not m + 1 as we did in the previous theoretical analysis.

the Pareto case, which is inside G 1 , but might be outside L 2 (R + ).

Only one projection algorithm is available in the literature that projects G 1 densities onto G 1,n , finding a set of scales and shapes that matches the input distribution: the MFK procedure that we already discussed in Section 2.3.

We want to compare the two algorithm, MFK and our, on fair grounds. Therefore, we slightly modify our algorithm to use Laguerre coefficients based on the theoretical integrals µ k,-1 , as MFK. The choice of the input distributions is also guided by [START_REF] Miles | Risk Aggregation: A General Approach via the Class of Generalized Gamma Convolutions[END_REF][START_REF] Furman | On Log-Normal Convolutions: An Analytical-Numerical Method With Applications to Economic Capital Determination[END_REF]. The experiment is as follows:

• We compute shifted moments (µ k,-1 ) k≤2n of the distribution through direct integration of the density, at 1024 bits of precision (to match the 300 digits that MFK requires).

• We use MFK to compute an approximation in G 1,n based on these moments.

• Through Eq 3, we compute Laguerre coefficients based on the same shifted moments µ k,-1 , and we minimize the loss given through Algorithm 2 to obtain our approximation.

• We compute Kolmogorov-Smirnov (one-sample, exact) distances between N = 10 000 simulations from the estimated G 1,n model to the theoretical distribution, on B = 100 different simulations.