Simple models for multivariate regular variations and the Hüsler-Reiss Pareto distribution
Zhen Wai Olivier Ho Clement Dombry
stat.ME






Simple models for multivariate regular variations and the Hüsler-Reiss Pareto distribution
    Zhen Wai Olivier Ho
[Univ. Bourgogne Franche-Comté, Laboratoire de Mathématiques de Besançon, UMR 6623, 16 route de Gray, 25030 Besancon, France. E-mail: zhen_wai_olivier.ho@univ-fcomte.fr] 
Clément Dombry 
[Univ. Bourgogne Franche-Comté, Laboratoire de Mathématiques de Besançon, UMR 6623, 16 route de Gray, 25030 Besancon, France. E-mail: clement.dombry@univ-fcomte.fr] 

    
=========================================================================================================================================================================================================================================================================================================================================================================================



 We revisit multivariate extreme value theory modeling by emphasizing multivariate regular variations and the multivariate Breiman Lemma. This allows us to recover in a simple framework the most popular multivariate extreme value distributions, such as the logistic, negative logistic, Dirichlet, extremal-t and Hüsler-Reiss models. In a second part of the paper, we focus on the Hüsler-Reiss Pareto model and its surprising  exponential family property. After a thorough study of this exponential family structure, we focus on maximum likelihood estimation. We also consider the generalized Hüsler-Reiss Pareto model with different tail indices and  a likelihood ratio test for discriminating constant tail index  versus varying tail indices.





§ INTRODUCTION

Regular variation is a fundamental notion in extreme value theory that was widely popularized  by Resnick  <cit.>. As a simple illustration of the importance a regular variation in univariate extreme value theory, consider an independent and identically distributed (i.i.d.) sequence X,X_1,X_2,… of positive random variables with cumulative distribution F. For n≥ 1, let a_n=F^←(1-1/n) be the quantile of order 1-1/n of F. Then the following statements are equivalent:


  i) the tail function 1- F is regularly varying at infinity with index -α<0, i.e.

    lim_u→∞1-F(ux)/1-F(u)=x^-α,   x>0;


  ii) the rescaled maximum a_n^-1max(X_1,…,X_n)  converge in distribution as n→∞ to a standard α-Fréchet distribution, i.e.

    lim_n→∞ℙ[a_n^-1max(X_1,…,X_n)≤  x ] =exp(-x^-α),   x>0;


  iii) the rescaled exceedance u^-1X of X given X>u converge in distribution as u→∞ to a standard α-Pareto distribution; i.e.

    lim_u→∞ℙ[ u^-1X>x| X>u]=x^-α,   x>1.


  iv) the sample point process {X_i/a_n, 1≤ i≤ n} converge to a Poisson point process on (0,∞) with intensity α x^-α-1dx.

The equivalence i)-ii) dates back to Gnedenko <cit.>, the equivalence ii)-iii) is due to Balkema and de Haan <cit.> and the equivalence i)-iv) can be found in Resnick <cit.>. As will be reviewed in Section <ref>, a similar result holds in the multivariate setting and multivariate regular variations is crucial in multivariate extreme value theory. 


Historically, multivariate extreme value theory has been developped by considerations on the asymptotic behaviour of i.i.d. random vectors.  Key early contributions are the papers by Tiago de Oliveira
<cit.>, Sibuya <cit.>, de Haan and Resnick <cit.>, Deheuvels <cit.>. The general structure of multivariate extreme value distribution has been characterized by de Haan and Resnick <cit.> in terms of the so-called spectral representation. Domain of attractions have been characterized by Deheuvels <cit.> that pointed out the convergence of the dependence structure to an extreme value copula. Since then a rich literature has emerged on modeling or statistical aspects of the theory, of which a nice recent review from the copula viewpoint is provided by  Gudendorf and Segers <cit.>. 

More recent developments focus on exceedances over high threshold in a multivariate setting and the so called multivariate generalized Pareto distributions. Seminal papers in that direction are Coles and Tawn <cit.> and Rootzen and Tajvidi <cit.>. Further recent development on modeling and statistical aspects include Rootzeén et al. <cit.> and Kiriliouk et al. <cit.>.

In this framework, the motivations of the present paper are twofold. In a first part corresponding to Section 2, we revisit  multivariate extreme value theory models and put the emphasis on regular variations and the limiting homogeneous measure. More precisely, a multivariate extension of the celebrated Breiman Lemma due to Davis and Mikosch <cit.> allows us to construct a regularly varying random vectors as a product of a heavy tailed random variable (thought as a radial component) and a sufficiently integrable random vector (thought as a spectral component). The limiting homogeneous measure is easily characterized and, for specific choice of the spectral component, we recover standard parametric models from multivariate extreme value theory such as the Hüsler-Reiss <cit.>, extremal-t <cit.>, logistic, negative logistic  or Dirichlet models <cit.>. We believe putting the emphasis on the exponent measure is important since it is the fundamental notion that unifies maxima, exceedances or point processes approaches in extreme value theory.  On the other hand, from the copula point of view, the multivariate Breiman Lemma provides a general framework for deriving extreme value copula models closely related to the results by Nikoloulopoulas et al. <cit.> or Belzile and  Nešlehová <cit.>. 

The second part of the paper corresponds to Sections 3 and 4 and proposes a thorough study of the so-called Hüsler-Reiss Pareto model, that is the exceedance Pareto model associated with the max-stable Hüsler-Reiss model <cit.>. The exceedances of the related  Brown-Resnick spatial model were considered recently by Wadsworth and Tawn <cit.> who proposed inference via censored maximum likelihood, see also Kiriliouk et al. <cit.>. Here, we focus on the finite-dimensional multivariate Husler-Reiss Pareto model and notice that it has a simple exponential family structure (see Bandorf-Nielsen <cit.>), that seems to have been overlooked in the literature. We propose in Section 3 an extensive study of this exponential family structure and consider also  maximum likelihood inference as well as perfect simulation.  We extend these results in Section 4 where we consider the non-standard Husler-Reiss Pareto model that incorporates different tail parameters for the different margins. Maximum likelihood estimators are shown again to be asymptotically normal and an alternating optimization procedure is considered. To conclude, we propose a maximum likelihood ratio test for testing the equality of the different marginal tail parameters.

Vector notations: we denote by ·_∞ the max-norm on ℝ^d and by · an arbitrary norm, 1_d=(1,…,1) is the vecor with all components equal to 1. Operations on vector are usually meant componentwise. The componentwise maximum of vector is denoted  max(x_1,x_2)=x_1∨ x_2, the comparison of vectors x_1≤ x_2 is  meant componentwise so that x_1≰x_2 means that some components of x_1 is larger than the corresponding component of x_2. For x∈[0,∞)^d, we note [0,x] the cube [0,x_1]×⋯×[0,x_d] and [0,x]^c=[0,∞)^d∖[0,x].





§ A SIMPLE MODEL FOR MULTIVARIATE REGULAR VARIATION




 §.§ Preliminaries on multivariate regular variation

Following Hult and Lindskog <cit.>, we define multivariate regular variation in terms of M_0-convergence in ℝ^d rather than vague convergence in [-∞,∞]^d∖{0}. This is completely equivalent in the multivariate setting but M_0-convergence can be more easily generalized to a metric space. 


Consider the space M_0(ℝ^d) of Borel measures μ on ℝ^d∖{0} that assigns finite mass on sets bounded away from 0, that is μ(ℝ^d∖ O) is finite for all O open neighborhood of 0. A sequence μ_n∈ M_0(ℝ^d) is said to converge to μ∈ M_0(ℝ^d), noted μ_nM_0→μ, if ∫ fdμ_n→∫ fdμ
for all bounded continuous function f that vanishes on a neighborhood of 0. 

A random vector X on ℝ^d is called regularly varying with sequence a_n→+∞ if 

    nℙ(X/a_n∈·) M_0→Λ
 
with a non-zero limit measure Λ∈ M_0(ℝ^d). Necessarily, there exists α>0, called the tail index of X, such that the limit measure is homogeneous of order α, i.e.

    Λ(uA)=u^-αΛ(A)    u>0, A⊂ℝ^d∖{0}.

Furthermore, the sequence (a_n) is  regularly varying at infinity with index 1/α and a possible choice for the normalizing sequence a_n is 

    a_n=inf{x>0; ℙ(X_∞≤ x)≥ 1-1/n },   n≥ 1.

Due to its importance in multivariate extreme value theory, we emphasize here the case of random vectors with non negative components and regular variations on [0,∞)^d. In this simple case, regular variation can be characterized by the convergence of the tail function, see Hult and Lindskog <cit.>: we have equivalence 


  i) the random variable X is regularly varying on [0,∞)^d with limit measure Λ, that is

    nℙ(a_n^-1X∈·)M_0⟶Λ(·),   n→∞;


  i') the tail function 1-F(x) is regularly varying with  limit function V(x)=Λ([0,x]^c), i.e.

    lim_u→+∞1-F(ux)/1-F(u1_d)=V(x),   x∈ [0,∞)^d∖{0};



Paralleling the univariate extreme value theory and the equivalence i)-iv) mentioned in the introduction, we consider a sequence X,X_1,X_2,… of non negative random vectors with cumulative distribution F on [0,∞)^d and we assume for convenience ℙ(X=0)=0.  The following statements are known to be equivalent, see e.g. the monograph by Resnick <cit.> or Coles <cit.>:


  i) (regular variation) the random variable X is regularly varying on [0,∞)^d with α-homogeneous limit measure Λ;

  ii) (componentwise maxima) the rescaled componentwise maximum  a_n^-1max(X_1,…,X_n)  converge in distribution as n→∞ to a jointly α-Fréchet random vector with exponent function V(x)=Λ([0,x]^c), i.e.

    lim_n→∞ℙ[a_n^-1max(X_1,…,X_n)≤  x ] =exp(-V(x)),   x∈ [0,∞)^d∖{0};


  iii) (excess above threshold) the rescaled exceedance u^-1X  given that  some component of X exceeds u>0  converges in distribution as u→∞ to an α-Pareto random vector, i.e.

    lim_u→∞ℙ[ u^-1X ≰x| X≰u1_d]=V(x∨ 1_d)/V(1_d),   x∈ [0,∞)^d∖ [0,1]^d;


  iv) (sample point process) the sample point process {a_n^-1X_i, 1≤ i≤ n} converges in distribution to a Poisson point process on [0,∞)^d∖{0} with intensity Λ.





 §.§ A multivariate version of Breiman Lemma

Before considering its multivariate extension, let us recall the celebrated Breiman Lemma (see Breiman <cit.>). 

[Breiman lemma]
Let R and Z be independent non negative random variables satisfying either of the following conditions:


  i) the tail function 1-F of R is regularly varying at infinity with index -α<0 and 𝔼[Z^α+ε]<∞ for some ε>0;

  ii) 1-F(x)∼ Cx^-α as x→∞ for some C>0 and 𝔼[Z^α]<∞.

Then, the product RZ is regularly varying with index α and

    ℙ(RZ>x) ∼𝔼[Z^α]ℙ(R>x)   x→∞.



The following multivariate extension of Breiman Lemma follows the line of Davis and Mikosch <cit.> .
 
Let R be a non negative random variable and Z an independent d-dimensional random vector. Assume either of the following conditions is satisfied:


  i) the tail function 1-F of R is regularly varying at infinity with index -α<0 and 𝔼[Z^α+ε]<∞ for some ε>0;

  ii) 1-F(x)∼ Cx^-α as x→∞ for some C>0 and 𝔼[Z^α]<∞.

Then the product X=RZ defines a regularly varying random vector on [-∞,∞]^d∖{0} with index α. More precisely, : 

    nℙ(a_n^-1X ∈·)M_0⟶Λ(·),  ,

where a_n is the quantile of order 1-1/n of R and the limit measure Λ  is homogeneous of order α and  given by

    Λ(A)=∫_0^∞ℙ(uZ ∈ A)α u^-α -1du,   A⊂ℝ^d∖{0} .

Moreover, in the case when Z is non-negative, Λ is supported by [0,∞)^d∖{0} and we have

    V(x):=Λ([0,x]^c)=𝔼[⋁_i=1^d ( Z_i/x_i)^α],   x∈ [0,+∞)∖{0}.



For example, this applies directly to the multivariate Student distribution with ν degrees of freedom that is the product of an inverse χ^2(ν) distribution (with heavy tail of order ν/2) and an independent multivariate Gaussian distributions (with moments of all orders).  See Nikololoupolos et al. <cit.> and Section <ref> below.



Consider an arbitrary norm · on ℝ^d-1 and denote by 𝒮^d-1 the unit sphere. For x>0 and B⊂𝒮^d-1 Borel, define

    A={z∈ℝ^d : z>x,  z/z∈ B }.

We have, as n→∞, 

    nℙ(a_n^-1X ∈ A)    = nℙ( RZ> a_n x, Z/Z∈ B )=nℙ( RZ1_{ Z/Z∈ B } >a_n x )
       ∼ n 𝔼(Z^α1_{Z/Z∈ B}) ℙ(R>a_n x) ∼𝔼(Z^α1_{Z/Z∈ B}) x^-α n ℙ(R>a_n)
       → x^-α𝔼(Z^α1_{Z/Z∈ B}).

We have used here the univariate Breiman Lemma <ref> to go from the first to the second line and then the fact that R has a regularly varying tail with index α>0 and that nℙ(R>a_n)→ 1. Using the fact that the sets of the form (<ref>) form a convergence determining class (Hult and Linskog <cit.>), we deduce from Equation (<ref>) the M_0-convergence nℙ(X/a_n∈·)M_0⟶Λ(·),
where the limit measure Λ is characterized by

    Λ(A)=x^-α𝔼(Z^α1_{Z/Z∈ B})

for all set A of the form (<ref>). We then check that Λ admits the integral representation (<ref>). Computing the right hand side of (<ref>) with A given by (<ref>), we get

    ∫_0^∞ℙ(uZ ∈ A)α u^-α -1du
           =∫_0^∞𝔼(1_{uZ>x, Z/Z∈ B }) α u^-α -1du
       =𝔼(1_{Z/Z∈ B}∫_0^∞ 1_{ u>x/Z}α u^-α -1du)
       =x^-α𝔼(Z^α 1_{Z/Z∈ B})=Λ(A).

Since the sets A of the form (<ref>) form a determining class,  the integral representation (<ref>) holds for all A⊂ℝ^d∖{0} Borel.
We can then check directly that Λ is homogeneous: for v>0, 

    Λ (vA)   =∫_0^∞ℙ(uZ ∈ vA) α u^-α-1du
    =∫_0^∞ℙ(v^-1uZ ∈ A) α u^-α-1du
       =v^-α∫_0^∞ℙ(uZ ∈ A) α u^-α-1du=v^-αΛ(A),

where we used the change of variable u'=u/v on the second line.

Finally, when Z is supported by [0,∞)^d,  Equation (<ref>) implies that Λ is supported by [0,∞)^d∖{0} and the   tail function V is computed as follows:

    V(x)   :=Λ([0,x]^c)
    = ∫_[0,∞)^dℙ(uZ∉ [0,x])α u^-α-1du
       =∫_[0,∞)^dℙ(u> Z_ix_i )α u^-α-1du
       =∫_[0,∞)^dℙ(u> min_1≤ i≤ dx_i/Z_i)α u^-α-1du
    = 𝔼[ ( min_1≤ i≤ dx_i/Z_i)^-α]
       =𝔼[⋁_i=1^d ( Z_i/x_i)^α].




If Z has a density f_Z, then Λ is absolutely continuous with respect to the Lebesgue measure and its Radon-Nikodym derivative is given by

    λ(z)=∫_0^∞ f_Z(z/u) α u^-d-α-1du.

and is homogeneous of order -d-α, that is

    λ(vz)=v^-d-αλ(z),   v>0, z∈ℝ^d∖{0}.




If Z has a density f_Z, the measure Λ writes

    Λ(A)   =∫_0^∞ℙ(uZ∈ A) α u^-α -1du 
     =∫_0^∞∫_ℝ^d1_{ uz ∈ A}f_Z(z) dz u^-α-1du
       =∫_0^∞∫_A f_Z(z/u)α u^-α -d -1dz du= ∫_A λ(z) dz,

where we use the change of variable z'=uz and Fubini Theorem. Furthermore, with the change of variable u'=u/v,

    λ (vz)= ∫_0^∞ f_Z(vz/u) α u^-d-α-1du
    = v^-d-α∫_0^∞ f_Z(z/u) α u^-d-α-1du=v^-d-αλ(z).





 §.§ A copula point of view

When focusing on the dependence structure, Proposition <ref> can be rephrased in terms of copulas (we refer to Joe <cit.> for a background on copulas and Gudendorf and Segers for extreme value copulas <cit.>). Following Krupskii et al. <cit.>, we consider here the simple common factor model 
    X=α E 1_d+ Y

with α>0, E  exponentially distributed and, independently, Y  a  d-dimensional random vector such that 𝔼[e^α Y_i]<∞, i=1,…,d. The different component of X share the common factor E that  introduces dependence in the extremes, because the components of Y are lighter tailed. Since the exponential distribution has a density, all the components X_i=α E+Y_i have a continuous distribution. Sklar Theorem entails that the copula C_X pertaining to X is uniquely defined by 

    C_X(u_1,…,u_d)=F_X(F_X_i^←(u_1),…,F_X_d^←(u_d)),    (u_1,…,u_d)∈[0,1]^d,

where F_X denotes the multivariate cumulative distribution of X and F_X_i^← the quantile function of component X_i.


Consider the copula C_X associated to the random vector X defined by (<ref>).
Then 

    C_X^n(u_1^1/n,…,u_d^1/n) → C_V(u_1,…,u_d),   (u_1,…,u_d)∈[0,1]^d,

where

    C_V(u_1,…,u_d)=exp(-V(σ_1(-log u_1)^1/α,…,σ_d(-log u_d)^1/α) )

with 

    σ_i^α=𝔼[e^α Y_i]      V(x)=𝔼[∨_i=1^d e^α Y_i/x_i^α].

In words, C_X belongs to the domain of attraction of the extreme value copula C_V.

Here, we use the fact that exp(α E) has an α-Pareto distribution but, in view of the proof and the multivariate Breiman Lemma, the result holds  as soon as exp(α E) has an heavy tail with index α and (α+ε)Y_i has a finite exponential moment for i=1,…,d.


By Proposition <ref>, e^X=e^α Ee^Y is regularly varying with exponent function V and hence, the normalized maximum of n independent copies of X converge to an α-Fréchet vector with distribution function e^-V(x). On the other hand, since the exponential transformation operates separately on each component, e^X has copula C_X and  the normalized  maximum of n i.i.d. copies has copula C_X^n(u_1^1/n,…,u_d^1/n). It remains to note that C_V is the copula associated with the limiting 
α-Fréchet vector, where the i-th margin as  shape parameter α and scale parameter σ.
The fact that convergence of pointwise maxima implies convergence of the copula is justified in Deheuvels <cit.>.




 §.§ Examples

In this section, we apply Proposition <ref> and consider various models for the various random vector Z. For these models, we provide an explicit expression for the limit measure Λ that characterizes the regular variation of the product X=RZ. Our computations rely on the general form of the density λ expressed in Proposition <ref> and technical computations.



  §.§.§ Gaussian case

The following result states a regular variation result in connection with the extremal-t model, see Nikoloulopoulas et al. <cit.>.

In the framework of the multivariate Breiman's lemma, if Z ∼𝒩( 0, Σ), then the limit measure Λ has density

    λ(z)   =α/(2π)^d/2 |Σ|^1/2Γ( α+d/2) (z^t Σ^-1z/2)^-(α+d)/2,   z∈ℝ^d∖{0}.



 Starting from Eq. (<ref>) and introducing the Gaussian density, we get

    λ(z)   = ∫_0^∞ f_Z(z/u)α u^-α-d-1du
       =∫_0^∞1/√(2π)^d |Σ|^1/2exp{-1/2u^2 z^tΣ^-1z}α u^-α-d-1du.

The change of variable v=1/u in the integral yields

    λ(z)   =(2π)^-d/2 |Σ|^-1/2α∫_0^∞exp{-v^2/2 z^tΣ^-1z } u^α+d-1du
       =(2π)^-d/2 |Σ|^-1/2α2/z^t Σ^-1z∫_0^∞z^tΣ z/2exp{-v^2/2 z^tΣ^-1z } u^α+d-1du
       =(2π)^-d/2 |Σ|^-1/2α2/z^t Σ^-1z𝔼[ X^α+d-2]

where X has a Weibull distribution with shape parameter equal to 2 and scale parameter equal to √(2/(z^t Σ^-1z)). We deduce 

    𝔼[ X^α+d-2]=(z^t Σ^-1z/2)^-(α+d-2)/2Γ( α+d/2)

and we obtain the claimed formula for λ(z). 




  §.§.§ Log-normal case

The case of log-normal spectral functions is connected with the Hüsler-Resii model <cit.>, see also Wadsworth and Tawn <cit.>.

In the framework of the multivariate Breiman's lemma, if Z ∼ℒ𝒩(m,Σ) with Σ definite positive, then the limit measure Λ has density

    λ(z)= C  exp{-1/2logz^t Q logz + L logz }∏_i=1^d z_i^-1,   z∈(0,∞)^d,

where 

    C= α/(2π)^(d-1)/2|Σ|^-1/2√(1_d^t Σ^-11_d)exp{-1/2 m^t Σ^-1 m +1/2(m^tΣ^-11_d+α)^2/1_d^t Σ^-11_d},
       Q=Σ^-1 -Σ^-11_d 1_d^t  Σ^-1/1_d^t Σ^-11_d, 
       l=( m^t -α+m^tΣ^-11_d/1_d^t Σ^-11_d1_d^t) Σ^-1.

and 

    V(x)=C (2π)^(d-1)/2/α∑_i=1^d x_i^-α |Q_-i|^-1/2exp{1/2l^T_-iQ_-i^-1l_-i}Φ_d-1(logx_-i/x_i; Q_-i^-1l_-i, Q_-i^-1).



 Starting from Eq. (<ref>) and introducing the log-normal Gaussian density, we get

    λ(z)   =∫_0^∞ f_Z( z/u) α u^-α-d-1du
       =∫_0^∞∏_i=1^d z_i^-1α |Σ|^-1/2 (2π)^-d/2exp{-1/2(log(z) - log(u) 1_d-m)^tΣ^-1(log(z) - log(u) 1_d-m) }
         u^-α -1du

The change of variable v=log (u)
yields

    λ(z)   =α |Σ|^-1/2(2π)^-d/2∏_i=1^d z_i^-1∫_-∞^∞exp{P(v) }dv

with 

    P(v)   = -1/2(log(z)- v 1_d -m)^tΣ^-1(log(z)- v 1_d -m)-α v
       = -1/21_d^tΣ^-11_d v^2 + (logz^tΣ^-11_d- m^tΣ^-11_d -α) v -1/2logz^t Σ^-1logz -1/2m^tΣ^-1m  + logz^t Σ^-1 m
       =-1/2C_1 v^2+ C_2v+ C_3.

Recognizing a Gaussian integral, we get with X∼𝒩(0,C_1^-1),

    ∫_-∞^∞exp{ P(v)}dv =√(2π/C_1)e^C_3𝔼[exp{C_2 X}]
        =√(2π/C_1)e^C_3e^C_2^2/2C_1.

We deduce the claimed formula for λ(z) after some straightforward simplifications.




  §.§.§ Independent Fréchet case

The case of independent spectral components is related to the logistic model  <cit.>.

Suppose Z=(Z_1,…, Z_d) with Z_i ∼Frechet(λ_i,β) independent with β>α. Then, the limit measure Λ in multivariate Breiman's lemma has density 

    λ(z)=αβ^d-1Γ(d-α/β)  ∏_i=1^d z_i^-β-1/λ_i^-β(∑_i=1^d (z_i/λ_i)^-β)^(α+1)/β-d

with Γ the Gamma function and 

    V(x):=Λ([0,x]^c)=Γ(1-α/β)(∑_i=1^d(x_i/λ_i)^-β)^α/β.


 Starting from Eq. (<ref>) and introducing the product Fréchet density yields

    λ(z)   =∫_0^∞ f_Z( z/u) α u^-α-d-1du
       =∫_0^∞∏_i=1^d (z_i^-β-1 u^β+1βλ_i^βexp{ -(z_i/λ_iu)^-β}) α u^-α-d-1du 
       =αβ^d ∏_i=1^d z_i^-β-1/λ_i^-β∫_0^∞ u^-α+β d-1exp{ -u^β∑_i=1^d z_i/λ_i^-β}du.

The change of variable v=u^β∑_i=1^d( z_i/λ_i)^-β in the integral gives

    λ(z)=αβ^d ∏_i=1^d z_i^-β-1/λ_i^-β1/β(∑_i=1^d (z_i/λ_i)^-β)^(α+1)/β-d∫_0^∞ e^-v v^d -α/β -1dv

The last integral is the definition of the Gamma function Γ(d-α/β).
Proposition <ref> gives

    V(x)   =𝔼[⋁_i=1^d (Z_i/x_i)^α]
       =∫_0^∞ℙ( ⋁_i=1^d (Z_i/x_i)^α>x) dx
       =∫_0^∞ 1-∏_i=1^d ℙ((Z_i   /x_i)^α≤ x) dx.

Introducing the Fréchet density function yields

    V(x)=∫_0^∞ 1-exp(-x^-β/α∑_i=1^d(x_i/λ_i)^-β)dx.

The change of variable y=x(∑_i=1^d (x_i/λ_i)^-β)^-α/β gives

    V(x)=(∑_i=1^d(x_i/λ_i)^-β)^α/β∫_0^∞ 1-exp(-y^-β/α)dy

The last integral correspond to the expectation of a Fréchet(1,β/α) and therefore, assuming β>α,  we have the result.




  §.§.§ Independent Weibull case

The case of independent spectral components is related to the negative logistic model  <cit.>.

Suppose Z=(Z_1,…, Z_d) with Z_i ∼Weibull(λ_i, β) independent with α>β. Then the limit measure Λ in multivariate Breiman's Lemma has density

    λ(z)=αβ^d-1Γ(d+α/β)(∑_i=1^d ( z_i/λ_i)^β)^-(α+1)/β-d∏_i=1^d  z_i^β-1/λ_i^β

and 

    V(x):=Λ([0,x]^c)=Γ(1+α/β)∑_∅ J ⊂{ 1,⋯,d} (-1)^|J|+1(∑_j∈ J(x_j/λ_j)^β)^-α/β.


 Starting from Eq. (<ref>) and introducing the product Weibull density yields

    λ(z)   =∫_0^∞ f_Z( z/u) α u^-α-d-1du
       =∫_0^∞∏_i=1^d ( β/λ_i(z_i/uλ_i)^β-1exp{-(z_i/u λ_i)^β}) α u^-α-d-1du
       =αβ^d∏_i=1^d (1/λ_i(z_i/λ_i)^β-1) ∫_0^∞exp{-u^-β(∑_i=1^d(z_i/λ_i)^β) } u^-α-β d-1du.

The change of variable v=u^-β( ∑_i=1^d ( z_i/λ_i)^β) in the integral gives

    λ(z)=αβ^d-1(∑_i=1^d ( z_i/λ_i)^β)^-(α+1)/β-d∏_i=1^d  z_i^β-1/λ_i^β∫_0^∞ e^-vv^α/β+d-1dv.

Proposition <ref> yields

    V(x)   =𝔼[ ⋁_i=1^d (Z_i/x_i)^α]
       =∫_0^∞ 1-∏_i=1^d ℙ( (Z_i/x_i)^α≤ x) dx.

Introducing the Weibull density function yields

    V(x)   =∫_0^∞ 1-∏_i=1^d ( 1 -exp(-x^β/α(x_i/λ_i)^β))dx
       =∫_0^∞∑_∅ J⊂{ 1,⋯,d} (-1)^|J|+1exp(-x^β/α∑_j ∈ J(x_j/λ_j)^β)dx.

The change of variable y=x(∑_j∈ J(x_j/λ_j)^β)^α/β yields

    V(x)=∑_∅ J ⊂{ 1,⋯,d} (-1)^|J|+1(∑_j∈ J(x_j/λ_j)^β)^-α/β∫_0^∞exp(-y^β/α)dy.

The last integral correspond to the expectation of a Weibull(1,β/α).




  §.§.§ Independent Gamma case

This last example is related to the max-stable model with Dirichlet spectral density. β_i≡ 1, the restriction of λ on the simplex is proportional to the Dirichlet density.

Suppose Z=(Z_1,…, Z_d) with Z_i ∼Γ(θ_i, β_i) independent. Then

    λ(z)   =αΓ(α+ ∑_i=1^d θ_i)(∑_i=1^d β_i z_i )^-∑_i=1^d θ_i -α∏_i=1^d  (β_i^θ_i z_i^θ_i-1/Γ(θ_i))



    
    λ(z)   =∫_0^∞∏_i=1^d(β_i^θ_i/Γ(θ_i)(z_i/u)^θ_i-1 e^-β_i z_i/u) α u^-α-1-ddu
       =α∏_i=1^d  (β_i^θ_i z_i^θ_i-1/Γ(θ_i)) ∫_0^∞ u^-∑_i=1^d θ_i -α -1exp{-u^-1∑_i=1^d β_i z_i }du.

    Setting v= u^-1∑_i=1^d β_i z_i, we obtain
    
    λ(z)   = α(∑_i=1^d β_i z_i )^-∑_i=1^d θ_i -α∏_i=1^d  (β_i^θ_i z_i^θ_i-1/Γ(θ_i)) ∫_0^∞ e^-v v^∑_i=1^d θ_i +α-1du.





 §.§ Non standard regular variations

 Following Resnick <cit.>, non-standard multivariate regular variations correspond to different tail index for the different components. Proposition <ref> has a simple extension to this case.
        
        Let R be a non negative heavy-tailed random variable with index 1, α=(α_1, ⋯,α_d)∈ (0,∞)^d and  Z=(Z_1,⋯,Z_d) a d-dimensional random vector such that 𝔼| Z_i|^α_i + ε<∞ for some ε>0. Then the product X=(R^1/α_1Z_1,⋯,R^1/α_dZ_d)=R^1/αZ satisfies
        
    nℙ(a_n^-1/αX∈·)M_0⟶Λ(·)

        where a_n is the quantile of order 1-n^-1 of R and the limit measure Λ satisfies
        
    Λ(A)=∫_0^∞ℙ( u ^-1/αZ ∈ A ) du,   A ⊂ℝ^d∖{0} measurable.

        
        
        Proposition <ref> for X̃=RZ^α=(RZ_1^α_1,…,RZ_d^α_d) yields the regular variations for X̃. Then, the change of variable  X=X̃^1/α together with the continuous mapping theorem for M_0-convergence <cit.> imply the non-standard regular variations stated in Proposition <ref> . 
        
       



§ THE HÜSLER-REISS PARETO MODEL




 §.§ Definition and transformation properties


Motivated by Proposition <ref>, we introduce the family of Hüsler-Reiss Pareto distributions and study their properties. The main reason why we focus on that particular class is that it enjoys an exponential family property, see Bandorff-Nielsen <cit.>.


Let d≥ 2, a=(a_1,…,a_d)∈(0,∞)^d,  Q∈ℝ^d× d a  symmetric positive semi-definite matrix such that Ker Q=vect(1_d) and l∈ℝ^d such that  l^T1_d<0. The Hüsler-Reiss Pareto model on [0,∞)^d∖[0,a] with parameters (Q,l) is defined by the density

    f_a(z;Q,l)=1/C_a(Q,l)exp(-1/2log z^T Q log z+l^Tlog z )(∏_i=1^d z_i^-1)1_{z≰ a},    z∈ (0,∞)^d,

with C_a(Q,l) the normalization constant. We call α=-l^T1_d>0 the exponent of the Pareto distribution f_a(z;Q,l).

We write Z⇝HRPar_a(Q,l) for a  random vector Z with density f_a(z;Q,l).
 
    
  
The Hüsler-Reiss Pareto model is closely connected with the exponent measure λ obtained in Proposition 3.2. Indeed, the parameters Q and l introduced there satisfy the constraint stated in Definition <ref>.  The symmetric semi-definite positive matrix Q satisfies 

    Q1_d=(Σ^-1 -Σ^-11_d 1_d^T Σ^-1/1_d^T Σ^-11_d)1_d=0

and, for all vector x∈ℝ^d∖{0} such that x^TΣ^-11_d=0, we have x^TQx>0 whence we deduce Ker Q=vect(1_d). As for l, we check readily

    l^T1_d=( m^T -α+m^TΣ^-11_d/1_d^T Σ^-11_d1_d^T) Σ^-11_d=-α<0.

Conversely, for all (Q,l) as in Definition <ref>, there exist  (non unique) 
Σ∈ℝ^d× d and m∈ℝ^d such that Equations (<ref>) and (<ref>) are satisfied.


In dimension d=2, the model parameters are

    Q=([  c -c; -c  c ])     l=([ l_1; l_2 ])    q>0, l_1+l_2<0.

The exponent is α=-(l_1+l_2)>0 and 

    f_a(z;Q,l)=1/C_a(Q,l)exp(-c/2(log z_1-log z_2)^2+l_1log z_1+l_2log z_2 )1/z_1z_21_{z≰ a}.




Interestingly, Hüsler-Reiss Pareto distributions inherit from log-normal distributions a stability property under scale and power transformations.

Let Z⇝HRPar_a(Q,l).


  (i) For all u∈(0,∞)^d, uZ⇝HRPar_ua(Q,l+Qlog u).

  (ii) For all β>0, Z^β⇝HRPar_a^β(β^-2Q,β^-1l). In particular, if Z has exponent α,  Z^β has exponent α/β.



The change of variable z̃=uz implies

    ℙ(uZ ∈ A)=∫_A f_a(z̃/u;Q,l)∏_i=1^du_i^-1dz̃.
    
Simple computations show that 

    f_a(z/u;Q,l)∏_i=1^du_i^-1
         =1/C_a(Q,l)exp(-1/2(log z-log u)^T Q (log z-log u) +l^T (log z -log u))(∏_i=1^d z_i^-1)1_{z ≰ ua}
          = C_ua(Q,l+Qlog u)/exp(1/2log u^T Q log u +l^T log u)C_a(Q,l)f_ua(z;Q,l+Qlog u)
    
This proves (i) as well as the equality

    C_ua(Q,l+Qlog u)=exp(1/2log u^T Q log u +l^T log u)C_a(Q,l).

Using a similar reasoning, the change of variable z̃=z^β yields

    ℙ(Z^β∈ A)=∫_A f_a(z̃^1/β)∏_i=1^d 1/βz̃_i^1/β-1dz̃

and simple computations show that

    f_a(z^1/β)∏_i=1^d 1/β z_i^1/β-1
          =1/C_a(Q,l) β^dexp(-1/2log z^T β^-1 Q β^-1log z + l^T β^-1log z)(∏_i=1^d z_i^-1)1_{z≰ a^β}
          = C_a^β(β^-2,β^-1l)/C_a(Q,l) β^d f_a^β(z; β^-2Q, β^-1l)

This implies (ii) as well as the equality

    C_a^β(β^-2Q,β^-1l)=β^d C_a(Q,l).




As a consequence of Proposition <ref>, Hüsler-Reiss Pareto vectors with a=1_d and α=1 are particularly important, especially for simulation. Indeed, the random vector Z⇝HRPar_a(Q,l) with exponent α=-l^T1_d satisfies Zd=aZ^1/α where the random vector Z⇝HRPar_1_d(α^-2Q,α^-1(l-Qlog a)) takes values in [0,∞)^d∖[0,1_d] and  has exponent 1.


The following equalities on the normalizing constant seen in the proof of Proposition (<ref>) are worth noting:

    C_ua(Q,l+Qlog u) = exp(1/2log u^T Qlog u+l^Tlog u )C_a(Q,l)

and

    C_a^β(β^-2Q,β^-1l) = β^dC_a(Q,l).

As a consequence, we will often assume without loss of generality that a=1_d. The general case a∈(0,∞)^d follows with the relation

    C_a(Q,l) = exp(-1/2log a^T Qlog a+l^Tlog a )C_1_d(Q,l-Qlog a).







 §.§ Exponential family properties

An important property of the Hüsler-Reiss Pareto distributions introduced above is to form an exponential family. Let E be an euclidean space with dot product ⟨·,·⟩. A parametric family of densities (f(z;θ))_θ∈Θ with Θ⊂ E is  a canonical exponential family if it can be written in the form

    f(z;θ)=1/C(θ)e^⟨θ,T(z)⟩h(z),   z∈ℝ^d,

where T:  ℝ^d→ E is the natural sufficient statistic. The exponential family is called a full exponential family if 

    Θ={t∈ E : ∫_ℝ^de^⟨ t,T(z)⟩ h(z)  <∞}
 
is not contained in a strict subspace of E. For a detailed account on exponential family, the reader should refer to Barndorff-Nielsen <cit.>.

Our main result in this section is the following Theorem.

Consider the d(d+1)/2-dimensional Euclidean space

    E={(A,b)∈ℝ^d× d×ℝ^d : A^T=A, A1_d=0}

with inner product 

    ⟨(A,a),(A',a')⟩=∑_1≤ i,j≤ dA_i,jA'_i,j+∑_1≤ k≤ da_ka'_k.

Define

    Θ={(Q,l)∈ E : Q , Ker Q=vect(1_d), l^T1_d<0}.

For all fixed a∈(0,∞)^d, the Hüsler-Reiss Pareto distributions (f_a(z;θ))_θ∈Θ form a full canonical exponential family with parameter θ=(Q,l)∈Θ and sufficient statistic

    T(z)= ( -1/2(log z-log z)(log z-log z)^T,log z),

with log z=d^-1(1_d^T log z) 1_d.


 Without loss of generality, let a=1_d. Consider the intensity function

    λ̃(z)=exp(-1/2log z^T Q log z+l^Tlog z )(∏_i=1^d z_i^-1),   z∈ (0,∞)^d,

The symmetric matrix Q can be diagonalized in an orthonormal basis Q=UΔ U^T with Δ=diag(λ_1,…,λ_d) and U orthonormal.
Thanks to the condition Q1_d=0, we can suppose λ_1=0 and the first column of U is equal to U_1=1_d/√(d). Denote by Δ_-1 (resp. v_-1) the matrix Δ (resp. vector v) with its first row and column removed (resp. first component removed), Ũ the d×( d-1) matrix obtained by removing the first column of U. The change of variable log z=Uv gives

    ∫_(0,∞)^d 1_{z≰ a}λ̃(z) dz=∫_ℝ^dexp(-1/2v_-1^T Δ_-1 v_-1 +l^TŨv_-1+l^TU_1v_1) 1_v ∈ Adv

where A equals

    A= { v ∈ℝ^d : v_1 1_d ≰ -Ũv_-1}={v∈ℝ^d:v_1>a(v_-1); a(v_-1)=min_i -∑_j=1^d-1Ũ_ijv_j+1}.
 
By Fubini theorem,

    ∫_(0,∞)^d 1_{z≰ a}λ̃(z)dz=∫_ℝ^d-1exp(-1/2v_-1^T Δ_-1v_-1 + l^T Ũv_-1) ∫_a(v_-1)^∞exp(l^T U_1 v_1) dv_1 dv_-1.

The inner integral with respect to v_1 converges if and only if l^TU_1<0 and then

    ∫_(0,∞)^d 1_{z≰ a}λ̃(z)dz=∫_ℝ^d-1exp( -1/2v_-1^TΔ_-1 v_-1 +l^TŨv_-1+ l^T U_1 a(v_-1)) dv_-1

is finite if and only if Δ_-1 is positive definite. This proves that the integral converge if and only if (Q,l)∈Θ and that the exponential family is full.


In a general exponential model (<ref>), the logarithm of the normalisation constant C(θ) is related to the cumulant generating function of the natural statistics T by the relation

    log𝔼_θ[e^⟨ t,T(z)⟩]= log C(θ+t)-log C(θ),  θ,θ+t∈Θ.

If θ is an interior point of Θ, this implies 

    𝔼_θ[T(z)]=∂log C/∂θ(θ) 
        Var_θ[T(z)]=∂^2 log C/∂θ∂θ^T(θ).3.6.1

The computation of the normalization constant C(θ) is hence particularly important.

In the Hüsler-Reiss Pareto model described in Theorem <ref>, we have

    C_a(Q,l)= (2π)^(d-1)/21/α∑_i=1^d a_i^-αdet(Q_-i)^-1/2exp{1/2l_-i^TQ_-i^-1l_-i}Φ_d-1( loga_-i/a_i;Q_-i^-1l_-i,Q_-i^-1),

where α=-1_d^Tl, the notation l_-i (resp. a_-i) denotes the vector l (resp. a) with its ith component removed, Q_-i  the matrix Q with its ith column and row removed and Φ_d(z;m,Σ) denotes the cumulative distributive function at z of a d-dimensional multivariate Gaussian distribution with mean m and covariance Σ.

The expression for C_a(Q,l) was first established by Huser and Davison <cit.>. We provide here a direct proof that will be needed for further reference (proof of Proposition <ref>).

  With λ̃ the function defined by Equation (<ref>),
the normalization constant C_a(Q,l) is given by
		
    C_a(Q,l)=		∫_[0,a]^cλ̃(z) dz.

Since [0,a]^c=∪_i=1^d A_i with

    A_i={ z∈ℝ^d  :  z_i>a_i , z_-i/z_i≤  a_-i/a_i},    i=1,…,d,

we have
		
    C_a(Q,l)=		∑_i=1^d ∫_A_iλ̃(z) dz.

Using the homogeneity relation (<ref>) with x=z_i, we get

    λ̃(z)=λ̃(z_i z/z_i)=z_i^-d-αλ̃(z/z_i).

Since the ith component of z/z_i is equal to 1, we have also

    λ̃(z/z_i)=exp(-1/2logz̃_-i^T Q_-ilogz̃_-i+l_-i^Tlogz̃_-i)∏_j≠ i^d z̃_j^-1,  z̃_-i=z_-i/z_i.

These relations imply

    ∫_A_iλ̃(z) dz 
    
    				=   ∫_(0,∞)^d 1_{z_i>a_i, z_-i/z_i≤ a_-i/a_i}  z_i^-d-αλ̃(z/z_i) dz
    
    				=   ∫_(0,∞)^d 1_{z_i>a_i, z̃_-i≤ a_-i/a_i} z_i^-d-αexp(-1/2logz̃_-i^T Q_-ilogz̃_-i+l_-i^Tlogz̃_-i)∏_j≠ i^d z̃_j^-1 dz
    
    				=   ∫_a_i^∞∫_[0,a_-i/a_i] z_i^-α-1exp(-1/2logz̃_-i^T Q_-ilogz̃_-i+l_-i^Tlogz̃_-i)(∏_j≠ i^d z̃_j^-1) dz_idz̃_-i
    
    				=   1/αa_i^-α∫_[0,a_-i/a_i] z_i^-α-1exp(-1/2logz̃_-i^T Q_-ilogz̃_-i+l_-i^Tlogz̃_-i)(∏_j≠ i^d z̃_j^-1) dz̃_-i

		where we have used  for the third inequality the change of variable z→ (z_i,z̃_-i). In the last integral with respect to z̃_-i, we recognize a  log-normal density (up to a multiplicative factor), so that

    ∫_[0,a_-i/a_i]exp(-1/2logz̃_-i^T Q_-ilogz̃_-i+l_-i^Tlogz̃_-i)(∏_j≠ i^d z̃_j^-1) dz̃_-i
    
    =    (2π)^(d-1)/2det(Q_-i)^-1/2exp{1/2l_-i^TQ_-i^-1l_-i}Φ_d-1(log(a_-i/a_i);Q_-i^-1l_-i,Q_-i^-1).

The result follows:

    C_a(Q,l)   =∑_i=1^d ∫_A_iλ̃(z) dz
       = (2π)^(d-1)/21/α∑_i=1^d a_i^-αdet(Q_-i)^-1/2exp{1/2l_-i^TQ_-i^-1l_-i}Φ_d-1(log(a_-i/a_i);Q_-i^-1l_-i,Q_-i^-1).

    
		
		



Let Z⇝HRPar_a(Q,l) with exponent α=-l^T1_d>0. Then,


  (i) for all u=(u_1,…,u_d) such that ∑_i=1^d u_i<α, we have

    𝔼[∏_i=1^d Z_i^u_i]=C_a(Q,l+u)/C_a(Q,l).


  (ii) The expectation and covariance matrix of log Z are given by

    𝔼[ log Z_i]=∂log C_a/∂ l_i(Q,l)   i=1,…,d,

and

    Cov(log Z_i,log Z_j)=∂^2 log C_a/∂ l_i∂ l_j(Q,l)   i,j=1,…,d.


  (iii) Moreover, the expectation and covariance matrix of log Z satisfies

    𝔼[(log z-log z)(log z-log z)^T]=∂log C_a/∂ Q(Q,l)




For θ=(Q,l) ∈Θ, we have for all u=(u_1,…,u_d) such that ∑_i=1^d u_i<α

    θ + (0,u) ∈Θ

by definition of α. Using equality (<ref>) with t=(0,u), we have

    log𝔼_θ[e^⟨ u,log z⟩]=log C(Q,L+u)-log C(Q,L),

taking to the exponential and developing the product implies (i).

The results (ii) and (iii) are straightforwards applications of (<ref>).

 In dimension d=2 with a=(1,1) and the same notations as in Example (<ref>), we have

    C(Q,l)=√(2π)/α√(c){ e^l_1^2/2cΦ(-l_1/√(c)) +  e^l_2^2/2cΦ(-l_2/√(c))}.

The first order partial derivatives of log C are equal to

    ∂log C/∂ l_1   =-1/l_1+l_2+ cl_1Φ(-l_1/√(c))-φ(-l_1/√(c))/√(c)/Φ(-l_1/√(c)) +  e^c(l_2^2-l_1^2)/2Φ(-l_2/√(c))
    ∂log C/∂ l_2   = -1/l_1+l_2+ cl_2Φ(-l_2/√(c))-φ(-l_2/√(c))/√(c)/Φ(-l_2/√(c)) +  e^c(l_1^2-l_2^2)/2Φ(-l_1/√(c))
    ∂log C/∂ c   = -1/2c +1/2l_1Φ(-l_1/√(c))-l_1c^-3/2ϕ(-l_1/√(c))/Φ(-l_1/√(c)) +  e^c(l_2^2-l_1^2)/2Φ(-l_2/√(c))
       +1/2l_2Φ(-l_2/√(c))-l_2c^-3/2ϕ(-l_2/√(c))/Φ(-l_2/√(c)) +  e^c(l_1^2-l_2^2)/2Φ(-l_1/√(c))
 
This formulas provides respectively the expectations 𝔼[log Z_1], 𝔼[log Z_2] and -1/8𝔼[(log Z_1 -log Z_2)^2].
Formulas for the general case a=(a_1,a_2) can be deduced using Proposition <ref>.






 §.§ Simulation of HR-Pareto random vectors

We now consider the simulation of an Hüsler-Reiss Pareto random vector Z⇝HRPar_a(Q,l).
Thanks to the transformation property (<ref>), we focus on the case a=1_d.
In the following proposition, we denote by S={x∈ (0,∞)^d : x_∞=1} the unit sphere and 
we use S=∪_i=1^d S_i with S_i={x∈ S :  x_i=1}.

Let Z⇝HRPar_1_d(Q,l) with exponent α>0. 
Then R=Z and Θ=Z/Z  are independent and such that


  - R is a Pareto(α)-distributed real random variable, i.e. ℙ(R>r)=r^-α, r>1; 

  - Θ is a random vector on S satisfying, for i=1,…,d,

    ℙ(Θ∈ S_i)= det(Q_-i)^-1/2exp{1/2l_-i^TQ_-i^-1l_-i}Φ_d-1( 0;Q_-i^-1l_-i,Q_-i^-1)/∑_j=1^d det(Q_-j)^-1/2exp{1/2l_-j^TQ_-j^-1l_-j}Φ_d-1(0;Q_-j^-1l_-j,Q_-j^-1)

and, given Θ∈ S_i, Θ_i=1 and

    ℒ(Θ_-i|Θ∈ S_i)=ℒ(exp(G_i)| G_i≤ 0)       G_i ⇝𝒩_d-1(Q_-i^-1l_-i,Q_-i^-1).
 



 The proof is mostly a reinterpretation of the computations from the  proof of Proposition <ref>. The density of Z⇝HRPar_1_d(Q,l)  is given by

    f_1_d(z;Q,l)=1/C_1_d(Q,l)1_{z>1}λ̃(z),   z∈ (0,∞)^d,

with λ̃ the function defined by Equation (<ref>). From the proof of Proposition <ref>, we have 

    C_1_d(Q,l)=∑_i=1^d ∫_A_iλ̃(z)dz

with  

    A_i={z∈ (0,∞)^d : z>1, z/z∈ S_i},   i=1,…,d.

The expression  for A_i here is slightly different but equivalent since a=1_d. Consequently, we get

    ℙ(Θ∈ S_i)=∫_(0,∞)^d 1_{z/z∈ S_i }f_1_d(z;Q,l)dz=1/C_1_d(Q,l)∫_A_iλ̃(z)dz

which yields Equation (<ref>) in view of Propositio <ref> and its proof.

On the other hand, when Z∈ A_i or equivalently Θ∈ S_i, we have R=Z=Z_i whence the change of variable z→ (z_i,z̃_-i) in Equation (<ref>) provides exactly the joint distribution of (R,Θ_-i). This amounts to be the product  of  α-Pareto and  log-normal distributions, proving the  independence of R and Θ and the form of their distribution. 


In order to simulate the Gaussian random variable G_i conditioned on G_i≤ 0, we propose a recursive sampling procedure. Let i∈{1,…,d} be fixed and denote by G_i,j  the components of G_i. We first set G_i,i=0 and  J={1,…,d}∖{i} the set of indices to sample. For j∈ J, the conditional distribution of  G_i,j given the already sampled components G_i,J^c has a Gaussian distribution with mean and variance
 
    m_i,j=(Q_J,J^-1(l_J -Q_J,J^c G_i,J^c))_j    σ^2_i,j=(Q_J,J^-1)_j,j

 subject to the constraint G_i,j≤ 0.  By the inversion method, we can sample from G_i,j as

    G_i,j=m_i,j+σ_i,jΦ^-1(Φ(-m_i,j/σ_i,j)U_j),   U_j⇝Unif([0,1]),

where Φ denotes the standard normal cumulative distribution function.
Then, we replace J by J∖{j} and repeat the procedure with the next component to sample until J is empty.



Note that the above computations are closely related to the distribution of extremal functions in the conditional sampling procedure of the Brown-Resnick max-stable process, see Dombry et al. <cit.>. Based on Proposition <ref> and the above recursive scheme, Algorithm <ref> describes a simulation procedure for Hüsler-Reiss Pareto random vectors.































 §.§ Maximum likelihood inference

The exponential family property of the Hüsler-Reiss Pareto distributions makes maximum likelihood inference particularly convenient. We always suppose the threshold  a∈(0,∞)^d to be known and estimate the parameter θ=(Q,l)∈Θ from observations  z^(1),…,z^(n)∈ (0,∞)^d∖[0,a]. In the Hüsler-Reiss Pareto model, the log-likelihood of the sample writes, for θ=(Q,l)∈Θ,
    
    L_n(θ;z^(1),⋯,z^(n))   =1/n∑_i=1^n log f_a(z^(i);Q,l)
       = ⟨  (Q,l),T_n ⟩ - log C_a(Q,l) + cst

    where T_n is the sufficient statistic defined by
    
    T_n=(-1/2n∑_i=1^n(log z^(i)- log z^(i))( log z^(i)-log z^(i))^T ,1/n∑_i=1^nlog z^(i))

    and the constant term cst does not  depend on the parameter θ=(Q,l). 
    Using the classical theory of maximum likelihood estimation for exponential families, we obtain the following result, regarding existence, uniqueness and asymptotic normality of the maximum likelihood estimator
    
    θ̂_n=_θ∈ΘL_n(θ;z^(1),⋯,z^(n)).


 Let a∈(0,∞)^d and n≥ 1. 

  (i) (existence and uniqueness) For observations z^(1),…,z^(n)∈ [0,a]^c, the log-likelihood (Q,l)↦ L_n(Q,l;z^(1),⋯,z^(n)) is strictly concave on Θ.
 and a maximum likelihood estimator θ̂_n exists if and only the sample covariance matrix 
 
    V_n=1/n∑_i=1^n log z^(i)log z^(i)T-(1/n∑_i=1^n log z^(i))(1/n∑_i=1^n log z^(i))^T
 is conditionally definite positive in the sense that v^T V_n v >0 for all v ∈ℝ^d∖{0_d} such that v^T 1_d=0.
If it exists, the maximum likelihood θ̂_n^mle is the unique solution of the score equation

    ∂log C_a/∂θ (θ)=T_n,  θ∈Θ.


  (ii) (asymptotic normality) Let θ=(Q,l)∈Θ and assume Z^(1),…,Z^(n) are generated from the distribution HRPar_a(Q,l). Then, for  n≥ d-1, there exists almost surely a unique maximum likelihood estimator θ̂_n^mle which is asymptotically normal and efficient, that is  
     
    √(n)(θ̂_n^mle -θ) d⟶𝒩(0, I(θ)^-1),   n→∞,

where I(θ) is the Fisher Information matrix given by

    I(θ)=- ∂^2 log C_a/∂θ∂θ^T (θ).





In statement i), if  1_d^Tlog z^(1),…,1_d^Tlog z^(n) are not all equal for i=1,⋯,n then the condition V_n conditionally definite positive is equivalent to V_n definite positive.


The proof of Theorem <ref> relies on the following Lemma.

Recall the definition (<ref>)  of  the sufficient statistic T(z). Then, the closed convex hull of the set

    S={T(z) ; z∈ (0,∞)^d, z≰ 1_d}

is equal to 

    C={ (Q,l)∈ E:Q≼ -1/2(l-l)(l-l)^T },

where Q_1≼ Q_2 means that the symmetric matrix Q_2-Q_1 is semi-definite positive. 


The change of variable u=log z shows that 

    S={(-1/2(u-u)(u-u)^T, u ) , u≰ 0}

where u=d^-1(1_d^T log z)1_d.
It is easily shown that C is closed, convex and contains  S, so that  conv(S)⊂conv(C)=C. We consider now the reverse inclusion. Consider U, U^(1),U^(2),… i.i.d. with mean l, variance Σ and such that U≰ 0 a.s. The random element

    S_n=1/n∑_i=1^n (-1/2( U^(i)- U^(i))(U^(i)- U^(i))^T,U^(i)).

belongs to conv(S) and, by the  law of large numbers, 

    S_n a.s.⟶S_∞= (-1/2𝔼((U-U)(U-U)^T),𝔼(U)),   n →∞,
 
so that S_∞∈conv(S). We prove below that for all (Q,l) ∈ C, one can choose Σ such that S_∞=(Q,l)∈conv(S), proving the reverse inclusion C⊂conv(S). Using U=d^-11_d 1_d^T U, we deduce

    𝔼((U-U)(U-U)^T )   =𝔼((U-d^-11_d1_d^T U)(U-d^-11_d1_d^T U)^T )
       =𝔼((I-1/d1_d1_d^T )UU^T(I-1/d1_d1_d^T )^T )
    
           =(I-1/d1_d1_d^T)(Σ+ll^T)(I-1/d1_d1_d^T)^T

It is proved in Lemma <ref> that the  linear operator on the space of symmetric d× d matrices  defined by

    P: M ↦(I-1/d1_d1_d^T)M(I-1/d1_d1_d^T)^T

is the orthogonal projection on the linear subspace {M: M1_d=0}. 
Therefore, for Σ such that
P(Σ+ll^T)=Q, we have

    S_n ⟶ (P(Σ+ll^T),l)=(Q,l).

In particular, since we can take Σ among all symmetric positive semi-definite matrix, the choice Σ=-2Q-P(ll^T) which is positive by definition of C leads to the result. Therefore C⊂conv(S). 



We assume here without loss of generality that a=1_d. The cumulant transform θ∈ (Q,l)∈Θ↦log C_a(Q,l) is a strictly convex function. Therefore the log-likelihood L_n is strictly concave as a difference of a linear function and a strictly convex function. The general theory for exponential families (see e.g. Bandorff-Nielsen <cit.>) ensures that the maximum likelihood estimator exists if an only if the sufficient statistic T_n belongs to the interior of the closed convex hull of the support of T, that is  T_n ∈int(conv(S))=int(C) with S and C defined in Lemma <ref>. In this case, Theorem 9.13 in Barndorff-Nielsen <cit.> implies that the maximum likelihood estimator is unique and solves the score equation (<ref>). So in order to prove statement (i), it remains  to prove that T_n ∈int(conv(S)) if and only if V_n is conditionally definite positive. Note that, by Lemma <ref>,

    int (conv(S))=int(C)={ (Q,l)∈ E :  Q≺ -1/2(l-l)(l-l)^⊥ on vect(1_d)^T },

where Q_1≺ Q_2  on vect(1_d)^T means that  v^T(Q_2-Q_1)v>0 for all ℝ^d∖{0_d} such that v^T 1_d=0. For such v and for (Q,l)=T_n, we have

    v^T(-Q -1/2(l-l)(l-l)^T)v
    
        =   v^T( 1/2n∑_i=1^n(log z^(i))( log z^(i))^T -1/2(1/n∑_i=1^nlog z^(i))(1/n∑_i=1^nlog z^(i))^T)v≥ 0
    
        =    v^T V_n v

whence we deduce that T_n∈int(conv(S)) if and only if  V_n is conditionally positive. 

Statement ii) follows directly from the general theory of exponential families since the Hüsler-Reiss distributions form a full rank exponential family (see e.g. Van der Vaart <cit.>).






















§ THE GENERALIZED HÜSLER-REISS PARETO MODEL

       
        

 §.§ Definition and transformation properties

        
        Let d≥ 2 and define Θ the  set of all θ=(α,Q,l) such that:
        
            
  - α∈(0,∞)^d,
            
  - Q∈ℝ^d× d is symmetric semi-definite positive and Ker Q= vect(1_d),
            
  - l ∈ℝ^d satisfies l^T 1_d=-1.
         
        For a∈ (0,∞)^d, the generalized Hüsler-Reiss Pareto model on [0,a]^c=[0,∞)^d ∖[0,a] with parameters θ=(α,Q,l) is defined by the density
        
    f_a(z;θ)=1/C_a(θ)exp(-1/2log z^T D_α Q D_αlog z + l^t D_αlog z )(∏_i=1^d z_i^-1)1_{z ≰ a}

         with C_a(θ) the normalization constant and D_α the diagonal matrix with diagonal α.

         We write Z⇝HRPar_a(α,Q,l) for a  random vector Z with density f_a(z;α,Q,l).
        
        For λ>0, the substitution (α,Q,l)↦ (λα,λ^-1/2Q,λ^-1l)  leaves Equation (<ref>) invariant so that the condition l^T 1_d=-1 is meant to ensure that the model is identifiable. In the case  α=α̅1_d with α̅>0,  the generalized Hüsler-Reiss model coincides with the Hüsler-Reiss Pareto model since f_a(z;α,Q,l)=f_a(z;α̅^2Q,α̅l) and α̅ is the tail index.
        
        Similarly as  HR-Pareto distributions, generalized HR-Pareto distributions enjoy a stability property under scale and power transformations.
        
        Let Z ⇝HRPar_a(α,Q,l).
        
            
  (i) For all u∈ (0,∞)^d, uZ⇝HRPar_ua(α,Q,l+QD_αlog u).
            
  (ii) For all β∈ (0,∞)^d, Z^β⇝HRPar_a^β(α/β,Q,l).
        
        
        
        The change of variable z̃=uz implies
        
    ℙ(uZ ∈ A)=∫_A f_a(z̃/u;α,Q,l)∏_i=1^d u_i^-1dz̃.

        Similarly as in  the proof of Proposition (<ref>), we check that
        
    f_a(z/u;α,Q,l)∏_i=1^du_i^-1
          =C_ua(α,Q,l+Qlog u)/exp{1/2log u^T D_α QD_αlog u + l^TD_αlog u }C_a(α,Q,l)f_ua(z;α,Q,l+QD_αlog u)

         whence statement (i) follows.
        The change of variable z̃= z^β implies 
        
    ℙ(Z^β∈ A)=∫_A f_a(z̃^1/β;α,Q,l)∏_i=1^d β_i^-1z̃_i^1/β_i -1dz̃

        and simple computations result in 
        
    f_a(z^1/β;α,Q,l)∏_i=1^d β_i^-1z_i^1/β_i -1
                =C_a^β(α/β,Q,l)/C_a(α,Q,l)∏_i=1^d β_if_a^β(z;α/β,Q,l)

        whence statement (ii) follows.
        
        
        We deduce a simple relation between generalized HR-Pareto distribution and (standard) HR-Pareto distribution.
        
            Let Z⇝HRPar_a(α,Q,l) with α∈ (0,∞)^d. We have Zd=aZ̃^c/α where Z̃⇝HRPar_1_d(Q,l-QD_αlog a) with exponent c>0. 
            Moreover, we have the relationships 
            
                
  (i) C_ua(α,Q,l+Qlog u)=exp{1/2log u^T D_α Q D_αlog u+ l^T D_αlog u}C_a(α,Q,l)
                
  (ii) C_a^β(α/β,Q,l)=C_a(α,Q,l)∏_i=1^d β_i.
            
        
        The following proposition relates the moments of generalized Hüsler-Reiss Pareto model with those of the Hüsler-Reiss Pareto model.
         
        Without loss of generality, assume a=1_d and let Z⇝HRPar(α,Q,l). Then, the expectation and the covariance matrix of log Z are given by
            
    𝔼_α,Q,l[log Z_i]=α_i^-1𝔼_Q,l[log Z_i]   i=1,⋯,d

            and 
            
    Cov_α,Q,l(log Z_i,log Z_j)=α_i^-1α_j^-1Cov_Q,l(log Z_i,log Z_j )

            where 𝔼_Q,l and Cov_Q,l are the expectation and covariance of Hüsler-Reiss Pareto distribution with exponent 1.
        
        
        Proposition <ref> yields 𝔼_α,Q,l[log Z_i]=∫1/α_ilog z_i f_1(z;1_d,Q,l)dz. 
        Similarly, we have 𝔼_α,Q,l[log Z_i log Z_j]=1/α^i α^j𝔼_Q,l[log Z_i log Z_j].
        Thus the result.
        
        
        The family of the generalized Hüsler-Reiss Pareto distributions form a curved exponential family with minimal sufficient statistic T given by
        
    T(z)=( log zlog z^T,log z).

        The associated natural parameter space contain positive definite matrices and the set of parameters of interest (α,Q,l) is included in the boundary of the natural parameter space, making the theory difficult. 
        
        
        

 §.§ Maximum likelihood inference

        We assume without loss of generality that  a=1_d is the known threshold. Based on independent observation  Z^(1), Z^(2),⋯ with distribution  HRPar_1_d(θ_0), θ_0∈θ we define the log-likelihood
        
    L_n(θ;Z^(1),⋯,Z^(n))=∑_i=1^n log f_1_d(Z^(i),θ),  θ∈Θ,

        and consider maximum likelihood estimation. It should be noted that we were not able to apply directly the 'classical' maximum likelihood estimation theory from Lehman <cit.> that uses differentiability properties of the likelihood. Indeed, despite some substantial efforts, we could not prove the relations
        
    ∂^k/∂θ^k∫_(0,∞)^df_1_d(z;θ)dz=∫_(0,∞)^d∂^k/∂θ^kf_1_d(θ,z)dz,   k=1,2,3,

        that are required (assumption M7 in <cit.>). Instead, we use differentiability in quadratic mean and local expansion of the likelihood process as in van der Vaart <cit.>.
       
        
        The statistical model {f_1_d(θ;z), θ∈Θ} is differentiable in quadratic mean. Furthermore, the local likelihood process defined by
        
    L̃_n(h)=L_n(θ_0+h/√(n);Z^(1),⋯,Z^(n)),  θ_0+h/√(n)∈Θ,
       
        satisfies,  uniformly on compact sets, 
        
    L̃_n(h)   =L̃_n(0)+∂L̃_n/∂ h(0)^T h - 1/2 h^T I_θ_0 h + o_p(1),
    ∂L̃_n/∂ h(h)   =∂L̃_n/∂ h(0)-I_θ_0h+o_p(1),
    ∂^2 L̃_n/∂ h ∂ h^T(h)   =-I_θ_0 + o_p(1)

            with  I_θ_0  the Fisher information matrix at θ_0.
            Furthermore,  in Equations (<ref>)-(<ref>), 
            
    ∂L̃_n/∂ h(0)=1/√(n)∑_i=1^n ∂log f_1_d/∂θ(Z^(i),θ_0)⇝𝒩(0,I_θ_0)

            and in  Equation (<ref>), the o_P(1)-term is even uniform on {h≤ n^1/2-ε} for all ε>0.
            
         Differentiability in quadratic mean is proved thanks to Lemma 7.6 in van der Vaart <cit.>. It is easily checked that θ↦√(f_1_d(z;θ)) is continuously differentiable for every z. Then we need to check that, with ℓ(θ,z)=log f_1_d(θ,z), the matrix 
        
    I(θ)=𝔼_θ [∂ℓ/∂θ(θ,Z)∂ℓ/∂θ(θ,Z)^T]

        is well defined and continuous in θ. This follows easily from the fact that the log-likelihood has the specific form
        
    ∂ℓ/∂θ(θ,Z)=⟨ A(θ), T(z)⟩ + B(θ)

        with A(θ), B(θ)  continuous in θ and  T(Z)=(log Zlog Z^T,log Z). Since T(z) has moment of all orders  that depend continuously of θ (this is true for the exponential family HR-Pareto and hence for the generalized HR-Pareto distributions),  I(θ) is well defined and continuous in θ. From <cit.>, we deduce that the model is differentiable in quadratic mean. For further reference, note that by  <cit.>, we have
        
    𝔼_θ[∂ℓ/∂θ(θ,Z)]=0 ,    I(θ)=𝔼[∂ℓ/∂θ(θ,Z)∂ℓ/∂θ(θ,Z)^T ]

        and Equation (<ref>) holds for all fixed h (we don't have uniformity at this point).
        
        We now prove the uniform asymptotic expansion (<ref>). The change of variable θ=θ_0+h/√(n) yields
        
    ∂^2 L̃_n /∂ h∂ h^T (h)=1/n∂^2 L_n /∂θ∂θ^T (θ)=1/n∑_i=1^n ∂^2 ℓ/∂θ∂θ^T(θ;Z^(i)),
 
      so that, by the law of large numbers,
           
    ∂^2 L̃_n /∂ h∂ h^T (0) a.s.⟶𝔼_θ_0[ ∂^2 ℓ/∂θ∂θ^T(θ_0;Z)]:=-J_θ_0,   .

        We don't know at this point that J_θ_0=I_θ_0, this will be proven in a final step. Thanks to Taylor-Lagrange formula,  the second-order derivative increment
      
    ∂^2 L̃_n /∂ h∂ h^T (h)-∂^2 L̃_n /∂ h∂ h^T (0)

     has  norm upper bounded, for h≤ n^1/2-ε, by 
     
    Cn^1/2-εmax_h≤ n^1/2-ε∂^3L̃_n/∂ h^3 (h) 
         =Cn^-1-εmax_θ-θ_0≤ n^-ε∂^3 L_n/∂θ^3 (θ).

    The specific form 
    
    ℓ(θ,Z)=-1/2⟨log z log z^T, D_α Q D_α⟩ +⟨log z, D_α l⟩ -log C_1_d(θ)

    implies that the third order derivative is upper bounded  by 
    
    ∂^3 L_n/∂θ^3(θ,z)≤ C_1+C_2∑_i=1^n log Z^(i)log Z^(i)T

    for some constants C_1,C_2>0 that does not depend on z∈ [0,1_d]^c and θ in a neighborhood of θ_0. We deduce
    
    ∂^2 L̃_n /∂ h∂ h^T (h)-∂^2 L̃_n /∂ h∂ h^T (0)≤ cn^-ε(C_1+C_21/n∑_i=1^n log Z^(i)log Z^(i)^T).

    By the law of large number, the sample mean converge almost surely  so that the right hand side is O_P(n^-ε)=o_P(1) uniformly in h≤ n^1/2-ε. Equations (<ref>) and (<ref>) together imply Equation (<ref>) with J_θ_0 instead of I_θ_0 for the moment. Equations (<ref>) and (<ref>) with J_θ_0 instead of I_θ_0 follow from  (<ref>) by integration with the o_P(1) term uniform on compact set. We have already noticed that differentiability in quadratic mean implies (<ref>) with I_θ_0, so that necessarily the two asymptotic expansion must coincide and J_θ_0=I_θ_0. This proves Equations (<ref>), (<ref>) and (<ref>) in their final form. Finally, in view of (<ref>), the asymptotic normality (<ref>) is a direct consequence of the central limit Theorem.
        
        
        The asymptotic development of the likelihood process stated in Proposition <ref> together with the Argmax Theorem allows us to study the properties of the maximum likelihood estimator (existence, consistency, asymptotic normality). An important argument is that, provided I_θ_0 is definite positive, the asymptotic expansion of the second order differentiate (<ref>) implies that the local likelihood process L̃_n(h) is strictly concave on {h< n^1/2-ε} with high probability. As we will see in the proof below, this entails that with high probability, the likelihood process L_n(θ) as a unique local maximizer in {θ-θ_0< n^-ε} that we define as θ̂^mle_n.
        
        
     
     Let θ_0∈Θ with I_θ_0  definite positive and assume the observations Z^(1), Z^(2),… are independent with distribution HRPar_a(θ_0). Then, there exists a  maximum likelihood estimators θ̂_n^mle that is asymptotically normal and efficient, i.e. 
     
    √(n)(θ̂_n^mle -θ_0) d⟶𝒩(0, I_θ_0^-1)   n→∞.

     
     
         The proof relies on Proposition <ref> and the Argmax theorem (van der Vaart <cit.> Corollary 5.58). Consider the stochastic processes 
        
    M_n(h)=L̃_n(h)-L̃_n(0),  h|≤ n^1/2-ε

        and
        
    M(h)=G h-1/2h^T I_θ_0h
 
        where G is a centered gaussian random vector with variance I_θ_0. 
        Proposition <ref> implies the convergence of M_n to M in distribution in l^∞(K) for all compact K. The limit process M is continuous and has an unique maximizer h given by ĥ=I_θ_0^-1G⇝𝒩(0,I_θ_0^-1). Define the maximizer
        
    ĥ_n=_h≤ n^1/2-ε M_n(h),
 
        where the argmax exists because M_n is continuous on a compact set.
        The argmax theorem implies that provided ĥ_n is tight, ĥ_nd→ĥ as n→∞. 
        
        We now prove the tightness of the sequence ĥ_n, n≥ 1. For all δ>0, there exists R>0 such that 
        
    ℙ(‖ĥ‖≤ R)≥ 1-δ.

        The relation
        
    M(h)=M(ĥ)-1/2(h-ĥ)I_θ_0(h-ĥ)

        implies 
        
    M(ĥ)-max_‖ĥ -h‖≥ 1M(h)≥1/2λ_min

        with λ_min>0  the smallest eigenvalue of I_θ_0. Therefore, with probability at least 1- δ, we have
        
    max_‖ h‖ =R+1M(h)≤ M(ĥ)- 1/2λ_min.

            The convergence in distribution of M_n to M in l^∞(K) with K={h:‖ h‖≤ R+1} implies, for large n,
            
    max_‖ h ‖≤ RM_n(h) -max_‖ h‖ =R+1M_n(h)≥1/4λ_min

            with probability at least 1-2δ. The  convergence (<ref>) together with the positive definiteness of I_θ_0   imply that M_n is strictly concave on {h≤ n^1/2-ε} with probability at least  1-δ for n large. Hence, Equation (<ref>) implies that the maximizer ĥ_n of M_n belongs to  {‖ h‖≤ R+1}. We have proved that for large n,  ℙ(ĥ_n≤ R+1)≥ 1-3δ, establishing  the tightness of ĥ_n. 
            
            Finally, on the event ĥ_n≤ R+1, ĥ_n belongs to the interior of {h≤ n^1/2-ε and is  therefore  a local maximizer of L̃_n such that ∂L̃_n/∂ h(ĥ_n)=0. Then θ̂_n^mle=θ_0+ĥ_n/√(n) is a local maximizer of L_n and such that ∂ L_n/∂θ(θ̂_n^mle)=0, that is a maximum likelihood estimator. Asymptotic normality of θ̂_n is a direct consequence of the convergence of ĥ_n to ĥ since √(n)(θ̂_n -θ_0) = ĥ_n d→ĥ∼𝒩(0,I_θ_0^-1).
     
     
     

 §.§ Optimizing the likelihood

        We have proved in the previous section that, with high probability, the likelihood function L_n is strictly concave on a neighborhood of θ_0 of size n^-ε, ε>0. However, there is no reason why it should be globally convex. We discuss here to issues associated with the likelihood optimization. The first is the initialization of an optimization algorithm and will be addressed thanks to a simple moment estimator that is √(n)-consistent and can serve as a starting point of optimization routines. The second point is how we can take advantage of the biconcavity of the problem: although not globally concave, the log-likelihood is biconcave in the sense that both partial applications α↦ L_n(α,Q,l) and (Q,l)↦ L_n(α,Q,l) are concave. In this context, it is natural to consider alternate convex optimization. 
        
        
         Let θ=(α,Q,l) ∈Θ and assume the observations Z^(1),Z^(2)⋯ independent with distribution HRPar_1_d(θ). For j=1,…,d define 
         
    N_n,j=1/n∑_i=1^n  1_{Z^(i)_j>1}      O_n,j=1/n∑_i=1^n1_{Z^(i)_j>1}log Z_j^(i).

         Then the estimator θ̂_0 =(α̂_0 , Q̂_0 , l̂_0) defined by
         
    α̂_0=(N_n,j/O_n,j)_1≤ j≤ d     
    (Q̂_0, l̂_0)=_Q,l L_n(α̂_0,Q,l)

         is strongly consistent and asymptotically normal.
        
        
        For j=1,…, d, the thresholded marginal Z_j|Z_j>1 are distributed according to a Pareto distribution with parameter α_j , so that 𝔼_θ[log Z_j| Z_j>1]=α_j^-1. Hence, by the law of large numbers
        
    N_n,j/O_n,ja.s⟶ℙ_θ(Z_j>1)/𝔼_θ( 1_{Z_j>1}log Z_j)=( 𝔼_θ[log Z_j| Z_j>1])^-1=α_j

        so that α̂_0 is a consistent estimator for α. 
        
        On the other hand, the vector Z^α has Hüsler-Reiss distribution HRP(Q,l), so that Theorem <ref> suggests the maximum-likelihood estimator
        
    (Q̂,l̂)=_Q,lL_n(α,Q,l)=Ψ(T_n(D_αlog Z^(1),…,D_αlog Z^(n))),

        where Ψ(t̅) denotes the unique solution of the score equation ∂log C/∂θ(Q,l)=t̅. As a general result for full exponential families (see e.g. Bandorff-Nielsen <cit.>), Ψ is a diffeomorphism. Since α is unknown and estimated by α̂_0, we set rather
        
    θ̂_0=(Q̂_0,l̂_0)=Ψ(T_n(D_α̂_0log Z^(1),…,D_α̂_0log Z^(n))).

        Some simple computations show
        
    T_n(D_α̂_0log Z^(1),…,D_α̂_0log Z^(n))= D_N_n/O_n

        where N_n, O_n and N_n/O_n denotes the vectors with components N_n,j, O_n,j and N_n,j/O_n,j respectively, and  
        
    M_n=1/n∑_i=1^d log Z^(i)     
            V_n=1/n∑_i=1^d log Z^(i) (log Z^(i))^T.

        Hence θ̂_0 can be written in the form  
        
    θ̂_0=Θ(N_n,O_n,M_n,V_n)

        with differentiable function Θ. The law of large number ensures the almost sure convergence of N_n,O_n,M_n,V_n as n→∞, whence strong consistency  θ̂_0  a.s.→θ follows. The central limit theorem ensures the asymptotic normality of (N_n,O_n,M_n,V_n), whence the asymptotic normality of θ̂_0 is deduced via the  δ-method (van der Vaart <cit.>).
        
        
        
        
        Let θ_0=(α_0,Q_0,l_0)∈Θ and assume the observations Z^(1),Z^(2),⋯ independent with distributions HRPar(θ_0). 
        Define θ̂_0 as in proposition (<ref>) and 
    V_n={θ∈Θ: ‖θ-θ_0‖< n^1/2-ε}.
         Define θ̂_n^mle as the unique minimizer of the negative log-likelihood on V_n, i.e.
        
    θ̂_n^mle=_θ∈ V_n   -L_n(θ;Z^(1),⋯,Z^(n)).

        Consider the alternating minimization estimators θ̂^(i)=(α̂^(i),Q̂^(i),l̂^(i)) defined by the recursive algorithm
        
    {[                             α̂^(i+1)                                  =_α -L_n(α,Q̂^(i),l̂^(i); Z^(1),⋯,Z^(n));                  (Q̂^(i+1),l̂^(i+1))                                =_Q,l    -L_n(α̂^(i+1),Q,l; Z^(1),⋯,Z^(n)) ].  for i>0

        and initialized with θ̂^(0)=θ̂_0.
        Then, with high probability, the sequence of estimators (θ̂^(i))_i≥ 0 converges almost surely to θ̂_n^mle, i.e.
        
    ℙ(lim_i →∞θ̂^(i)=θ̂_n^mle)→ 1,   as n→∞.

        
        
      The starting point estimator writes
      
    θ̂_0=θ_0+1/√(n)(√(n)(θ̂_0-θ_0)).

      Proposition <ref> and Prohorov's theorem implies that θ̂_0 ∈ V_n with high probability. Assuming the log-likelihood strictly concave on V_n, we show by recurrence that each iterate of the alternating minization algorithm belongs to V_n.
      Define the level set  
      
    ℒ_i={θ: -L_n(θ)≤ -L_n(θ̂^(i))+δ},   i≥ 0,

      where δ>0 is such that  ℒ_i∩∂ V_n=∅.  By convex optimization theory, the intersection between ℒ_i and V_n is a convex set. Let B_1 and B_2 be open balls centered at θ̂^(i) and (α̂^(i+1),Q̂^(i),l̂^(i)) such that B_1 and B_2 are subset of ℒ_i. The biconvex property of -L_n implies that  the convex hull conv(B_1,B_2)  is a subset of ℒ_i. It results that conv(B_1,B_2)⊂ℒ_i ∩ V_n. A similar reasoning concludes that θ̂^(i)∈ V_n for all  i≥ 0 and therefore the alternating minimization estimators θ̂^(i) converge to the unique minimizer in V_n. 
        

        

 §.§ A likelihood ratio test for  α_1=⋯=α_d

        Following the development of generalized pareto models, a natural question that arise when one is given a i.i.d. sample Z^(1),… ,Z^(n) with distribution HRPar(α,Q,l) is whether the Pareto model would be enough to modelize the data.  The following theorem provides a likelihood ratio test for testing α_1=⋯=α_d. 
        
        
        Let θ_0=(α,Q,l) ∈Θ with α=(α_1,…,α_d).
        Let Z^(1),…, Z^(n) be i.i.d. with distribution HRPar(θ_0). Denote by θ̂_n the maximum likelihood estimator in the Generalised Hüssler-Reiss Pareto model and θ̂_0 the maximum likelihood estimation in the Hüssler-Reiss Pareto model and define the likelihood log-ratio by
        
    Δ_n=L_n(θ̂_n)-L_n(θ̂_0).

        Then, under the null hypothesis α_1=⋯=α_d, the distribution of 2Δ_n converge to a chi-squared distribution with d-1 degree of freedom, i.e.
        
    2 (L_n(θ̂_n)-L_n(θ̂_0)) d→χ^2(d-1).

        
        
            Denote by Θ_0 the subset of Θ defined by { (α,Q,l)∈Θ : α_1=⋯,⋯= α_d }. Consider the local log-likelihood process L̃_n and its maximiser ĥ_n on Θ. Likewise, denote by ĥ^0_n the maximiser of L̃_n on Θ_0. 
            We prove below  that 2(L̃_n(ĥ_n)-L̃_n(ĥ_n^0)) d→χ^2(p-1).
            Simple calculations imply that the Taylor expansion of L̃_n at ĥ_n writes 
            
    L̃_n(h)=L̃_n(ĥ_n)-1/2(h-ĥ_n)I_θ_0(h-ĥ_n)+o_p(1)

            where the o_p term is uniform on compact sets containing ĥ_n. Taking a compact K large enough to contain both ĥ_n and ĥ_n^0, we have
            
    2(L̃_n(ĥ_n)-L̃_n(ĥ_n^0))   =min_h∈ K ∩Θ_0 2(L̃_n(ĥ_n)-L̃_n(h))
       =min_h∈ K ∩Θ_0 (h-ĥ _n)I_θ_0(h-ĥ_n)+o_p(1)

            Defining ⟨·,·⟩_I_θ_0 as the inner product induced by I_θ_0, i.e. ⟨ a,b ⟩_I_θ_0=a^t I_θ_0 b, we get
            
    2(L̃_n(ĥ_n)-L̃_n(ĥ_n^0))   =min_h∈ K ∩Θ_0‖ h-ĥ_n‖_I_θ_0^2 +o_p(1)

            The minimum is reached for h the orthogonal projection of ĥ_n into Θ_0 for the ‖.‖_I_θ_0 norm. Thus, we have
            
    2(L̃_n(ĥ_n)-L̃_n(ĥ_n^0))   =∑_i=1^p-1⟨ĥ_n, e_i⟩^2_I_θ_0 +o_p(1)

            where (e_1,⋯,e_p-1) is an orthonormal basis of Θ_0^⊥. Theorem (<ref>)  implies 
    (⟨ĥ_n,e_i⟩_I_θ_0)_1≤ i≤ p-1d→𝒩(0_p-1,I_p-1)

            which in turn results in 
    2(L̃_n(ĥ_n)-L̃_n(ĥ_n^0))d→χ^2(p-1).

        
       
        
       
        
            Let 𝒮_d denote the linear space of symmetric d× d matrices and P the linear operator defined as 
            
    P:𝒮_d→𝒮_d,   Q ↦ (I-1/d1_d1_d^T)Q(I-1/d1_d1_d^T).

            Then P is the orthogonal projection on the linear subspace 𝒮_d^0={S ∈ E: S1_d=0}.
        
        
        Let S ∈ E, we have 
        
    P^2(S)   =(I-1/d1_d1_d^T)^2 S (I-1/d1_d1_d^T)^2
       =(I-2/d1_d1_d^T+1/d^21_d1_d^T1_d1_d^T )S (I-2/d1_d1_d^T+1/d^21_d1_d^T1_d1_d^T )
       =(I-1/d1_d1_d^T) S (I-1/d1_d1_d^T)
       =P(S).

        Therefore P is idempotent.

        For S ∈𝒮_d^0, we have
        
    P(S)   =(I-1/d1_d1_d^T) S (I-1/d1_d1_d^T)
       =S-2/dS1_d1_d^T + 1/d^21_d1_d^TS1_d1_d^T
       =S.

        Therefore P acts  as the identity on 𝒮_d.

        For S ∈ (𝒮_d^0)^⊥, we have 
        
    P(S)   =...=0

        Therefore P is null on (𝒮_d^0)^⊥. This concludes the proof.
        
        
        
        
        





10

BdH74
A. A. Balkema and L. de Haan, Residual life time at great age, Ann.
  Probability 2 (1974), 792–804. 0359049

MR3221776
O. Barndorff-Nielsen, Information and exponential families in statistical
  theory, Wiley Series in Probability and Statistics, John Wiley & Sons,
  Ltd., Chichester, 2014, Reprint of the 1978 original [MR0489333].
  3221776

BN17
L. R. Belzile and J. G. Neslehová, Extremal attractors of liouville
  copulas, arXiv:1704.03377, 2017.

breiman:1965
L. Breiman, On some limit theorems similar to the arc-sin law, Theory of
  Probability and its Applications 10 (1965), no. 2, 323–331.

Coles01
S. Coles, An introduction to statistical modeling of extreme values,
  Springer Series in Statistics, Springer-Verlag London, Ltd., London, 2001.
  1932132

CT91
S. G. Coles and J. A. Tawn, Modelling extreme multivariate events, J.
  Roy. Statist. Soc. Ser. B 53 (1991), no. 2, 377–392. 1108334

DM08
R. A. Davis and T. Mikosch, Extreme value theory for space-time processes
  with heavy-tailed distributions, Stochastic Process. Appl. 118
  (2008), no. 4, 560–584. 2394763

dHR77
L. de Haan and S. I. Resnick, Limit theory for multivariate sample
  extremes, Z. Wahrscheinlichkeitstheorie und Verw. Gebiete 40
  (1977), no. 4, 317–337. 0478290

D78
P. Deheuvels, Caractérisation complête des lois extrêmes
  multivariées et de la convergence aux types extrêmes, Publ. Inst.
  Statist. Univ. Paris 23 (1978), 1–36.

DER13
C. Dombry, F. Éyi-Minko, and M. Ribatet, Conditional simulation of
  max-stable processes, Biometrika 100 (2013), no. 1, 111–124.
  3034327

G43
B. Gnedenko, Sur la distribution limite du terme maximum d'une série
  aléatoire, Ann. of Math. (2) 44 (1943), 423–453. 0008655

GS10
G. Gudendorf and J. Segers, Extreme-value copulas, Copula theory and its
  applications, Lect. Notes Stat. Proc., vol. 198, Springer, Heidelberg, 2010,
  pp. 127–145. 3051266

hult:lindskog:2006
H. Hult and F. Lindskog, Regular variation for measures on metric
  spaces, Publ. Inst. Math. (Beograd) (N.S.) 80(94) (2006), 121–140.
  2281910 (2008g:28016)

HD13
D. Huser and Davison A. C., Composite likelihood estimation for the
  brown–resnick process, Biometrika 100 (2013), no. 1, 511–518.

HR89
J. Hüsler and R.-D. Reiss, Maxima of normal random vectors: between
  independence and complete dependence, Statist. Probab. Lett. 7
  (1989), no. 4, 283–286. 980699

J15
H. Joe, Dependence modeling with copulas, Monographs on Statistics and
  Applied Probability, vol. 134, CRC Press, Boca Raton, FL, 2015. 3328438

KRSW16
A. Kiriliouk, H. Rootzén, J. Segers, and J. Wadsorth, Peaks over
  thresholds modelling with multivariate generalized pareto distributions,
  arXiv:1612.01773, 2016.

L99
E. L. Lehmann, Elements of large-sample theory, Springer Texts in
  Statistics, Springer-Verlag, New York, 1999. 1663158

NJL09
A. K. Nikoloulopoulos, H. Joe, and H. Li, Extreme vale properties of
  multivariate t copulas, Extremes 12 (2009), no. 2, 129–148.

KJLG17
Krupskii P., Joe H., Lee D., and Genton M. G., Extreme value limit of
  convolution of exponential and multivariate normal distribution - link to
  hüsler-reiss distribution, arXiv:??, 2017.

R87
S. I. Resnick, Extreme values, regular variation, and point processes,
  Applied Probability. A Series of the Applied Probability Trust, vol. 4,
  Springer-Verlag, New York, 1987. 900810

R07
, Heavy-tail phenomena, Springer Series in Operations Research
  and Financial Engineering, Springer, New York, 2007, Probabilistic and
  statistical modeling. 2271424

RSW17
H. Rootzén, J. Segers, and J. Wadsorth, Multivariate peaks over
  thresholds models, Extremes (2017), Forthcoming.

RT06
H. Rootzén and N. Tajvidi, Multivariate generalized Pareto
  distributions, Bernoulli 12 (2006), no. 5, 917–930. 2265668

S60
M. Sibuya, Bivariate extreme statistics. I, Ann. Inst. Statist. Math.
  Tokyo 11 (1960), 195–210. 0115241

TdO58
J. Tiago de Oliveira, Extremal distributions, Rev. Fac. Sci. Lisboa,
  Ser. A 7 (1958), 215–227.

vdV98
A. W. van der Vaart, Asymptotic statistics, Cambridge Series in
  Statistical and Probabilistic Mathematics, vol. 3, Cambridge University
  Press, Cambridge, 1998. 1652247 (2000c:62003)

WT14
J. L. Wadsworth and J. A. Tawn, Efficient inference for spatial extreme
  value processes associated to log-Gaussian random functions, Biometrika
  101 (2014), no. 1, 1–15. 3180654


