Stefka Asenova
email: stefka.asenova@uclouvain.be

Gildas Mazo
email: gildas.mazo@inra.fr

Johan Segers
email: jo-han.segers@uclouvain.be

Inference on extremal dependence in the domain of attraction of a structured Hüsler-Reiss distribution motivated by a Markov tree with latent variables

Keywords: multivariate extremes, tail dependence, graphical models, latent variables, Hüsler-Reiss distribution, Markov tree, tail tree, river network

A Markov tree is a probabilistic graphical model for a random vector indexed by the nodes of an undirected tree encoding conditional independence relations between variables. One possible limit distribution of partial maxima of samples from such a Markov tree is a max-stable Hüsler-Reiss distribution whose parameter matrix inherits its structure from the tree, each edge contributing one free dependence parameter. Our central assumption is that, upon marginal standardization, the data-generating distribution is in the max-domain of attraction of the said Hüsler-Reiss distribution, an assumption much weaker than the one that data are generated according to a graphical model. Even if some of the variables are unobservable (latent), we show that the underlying model parameters are still identifiable if and only if every node corresponding to a latent variable has degree at least three. Three estimation procedures, based on the method of moments, maximum composite likelihood, and pairwise extremal coefficients, are proposed for usage on multivariate peaks over thresholds data when some variables are latent. A typical application is a river network in the form of a tree where, on some locations, no data are available. We illustrate the model and the identifiability criterion on a data set of high water levels on the Seine, France, with two latent variables. The structured Hüsler-Reiss distribution is found to fit the observed extremal dependence patterns well. The parameters being identifiable we are able to quantify tail dependence between locations for which there are no data.

Introduction

A major topic in multivariate extreme value theory is the modeling of tail dependence between a finite number of variables. Informally, tail dependence represents the degree of association between the extreme values of these variables. Probabilistic graphical models [START_REF] Lauritzen | Graphical Models[END_REF][START_REF] Koller | Probabilistic Graphical Models: Principles and Techniques[END_REF][START_REF] Wainwright | Graphical models, exponential families, and variational inference[END_REF], are distributions which embody a set of conditional independence relations and have a graph-based representation, according to which the nodes of the graph are associated to the variables and the set of edges encode the conditional independence relations. The intersection of the two fields, extreme value theory and probabilistic graphical models, gives rise to the study of the tail behavior of graphical models.

Consider a river network where the interest is in extreme water levels or water flow in relation to flood risks. Figure 1 illustrates part of the Seine network. The graph fixed by the seven labeled nodes and the river channels between them can be a base for building a model for extremal dependence between the water levels at these sites.

Hydrological data are often used to fit models for multivariate extremes based on graphs. Water flows of the Bavarian Danube are analyzed in [START_REF] Engelke | Graphical Models for Extremes[END_REF]. [START_REF] Lee | Multivariate extreme value copulas with factor and tree dependence structures[END_REF] study water flows of the Fraser river, British Colombia. Precipitation data in the Japanese archipelago is treated in [START_REF] Yu | Modeling spatial extremes via ensemble-of-trees of pairwise copulas[END_REF], where the model is based on a spatial grid viewed as an ensemble of trees. Other extreme-value models involving graphs appear in [START_REF] Einmahl | A continuous updating weighted least squares estimator of tail dependence in high dimensions[END_REF] and [START_REF] Lee | Multivariate extreme value copulas with factor and tree dependence structures[END_REF], who study financial data under different models. The first paper uses max-linear models on a directed acyclic graph (DAG) [START_REF] Gissibl | Max-linear models on directed acyclic graphs[END_REF], and the second one a 1-factor model. [START_REF] Klüppelberg | Max-linear models on infinite graphs generated by Bernoulli bond percolation[END_REF] introduce an infinite max-linear model to analyze the distribution of extreme opinions in a social network.

Relatively recently the relation between extreme value distributions and conditional independence assumptions has been given theoretical relevance. The earliest is the article of [START_REF] Gissibl | Max-linear models on directed acyclic graphs[END_REF] introducing max-linear models as structural equation models on a DAG, followed by the regularly varying Markov trees in [START_REF] Segers | One-versus multi-component regular variation and extremes of Markov trees[END_REF] and the extremal graphical models in [START_REF] Engelke | Graphical Models for Extremes[END_REF] based on multivariate Pareto distributions. Earlier, [START_REF] Papastathopoulos | Conditional independence among max-stable laws[END_REF] showed that for a max-stable random vector with positive and continuous density, conditional independence implies unconditional independence, thereby concluding that a broad class of max-stable distributions does not exhibit an interesting Markov structure.

A key object of our paper is the multivariate Hüsler-Reiss distribution [START_REF] Hüsler | Maxima of normal random vectors: between independence and complete dependence[END_REF]) with parameter matrix having a particular structure linked to a tree as specified in Eq. ( 5). The structure is motivated by the fact that the max-domain of attraction of the said Hüsler-Reiss distribution contains certain regularly varying Markov trees. The latter property follows from results in [START_REF] Segers | One-versus multi-component regular variation and extremes of Markov trees[END_REF] and sets our work apart from the extremal graphical models in [START_REF] Engelke | Graphical Models for Extremes[END_REF], who impose a non-standard conditional independence relation on the multivariate Pareto distribution associated to a max-stable distribution, but without regard for the latter's max-domain of attraction. Still, it turns out that for trees, the structured Hüsler-Reiss models in [START_REF] Engelke | Graphical Models for Extremes[END_REF] and in our paper are the same, as explained in Appendix A.1. Another structured Hüsler-Reiss distribution based on trees is proposed in [START_REF] Lee | Multivariate extreme value copulas with factor and tree dependence structures[END_REF]. The form they propose is genuinely different from ours, however, as explained in detail in Appendix A.2.

We consider random samples from the distribution of a random vector ξ = (ξ v , v ∈ V ) with continuous margins whose variables are indexed by the node set V = {1, . . . , d} of an undirected tree with edge set E. After marginal standardization to the unit-Pareto distribution, we assume that the random vector is in the max-domain of attraction of the tree-structured Hüsler-Reiss distribution described in the previous paragraph. We emphasize that we do not assume that ξ itself satisfies any conditional independence relations with respect to the tree. The tree only comes into play via the imposed structure on the parameter matrix of the max-stable Hüsler-Reiss distribution containing the distribution of the standardized version of ξ in its max-domain of attraction.

The main result and contribution of our paper is a criterion for identifiability of all d -1 parameters θ e ∈ (0, ∞) for e ∈ E of the tree-structured d-variate Hüsler-Reiss distribution in case some of the d variables are latent (unobservable). To illustrate why the problem of latent variables is relevant, consider again the Seine network on Figure 1. The red dots designate junctions of two river channels (conversely, in a river delta, a channel could split into several ones). No measurement stations being present there, we cannot observe the water levels at those locations. We propose to treat those water levels as latent variables. The question is then whether it is still possible to identify all d -1 parameters. The answer is a surprisingly simple identifiability criterion: it is necessary and sufficient that all nodes indexing latent variables have degree at least three. The important practical implication is that, provided the criterion is met, the latent variables can be included in the model, reflecting the dependence structure more accurately than when they would have been ignored.

Given a random sample from a distribution in the max-domain of attraction of the tree-structured Hüsler-Reiss distribution, we propose three types of estimators of the edge parameters: a first one called method of moments estimator (MME) is based on the estimator proposed in [START_REF] Engelke | Estimation of Hüsler-Reiss distributions and Brown-Resnick processes[END_REF], a second one is based on the composite likelihood function (composite likelihood estimator or CLE) and the third one is essentially the pairwise extremal coefficient estimator (ECE) introduced in [START_REF] Einmahl | A continuous updating weighted least squares estimator of tail dependence in high dimensions[END_REF]. All estimators proposed allow for the fact that some of the d variables are latent, provided the identifiability criterion is met.

We illustrate the method by a detailed analysis of data on high water levels at several locations of the Seine network. The network is represented schematically as a tree with seven nodes indexing five observable variables and two latent ones. As the identifiability criterion is met, we can estimate the six dependence parameters of the tree-structured Hüsler-Reiss distribution, each parameter corresponding to an edge in the tree. For the three proposed estimators we compute parameter estimates and confidence intervals. We assess the goodness-of-fit by comparing the model output with various non-parametric measures of tail dependence. Finally, we compare the fitted tail dependence model incorporating latent variables with a model where the latent variables are ignored.

The outline of the paper is as follows: Section 2 presents some general theory and describes the model to which the identifiability criterion is applied. The latter is the focus of Section 3. Section 4 introduces the three estimators, used for statistical inference and Section 5 is dedicated to the study of high water levels on the Seine network. Concluding remarks and perspectives for further research are discussed in Section 6. The Appendix provides proofs that are not in the text, a numerical comparison between our structured Hüsler-Reiss method and the one of [START_REF] Lee | Multivariate extreme value copulas with factor and tree dependence structures[END_REF], clarification on the relationship between the different objects in our paper and the objects in [START_REF] Engelke | Graphical Models for Extremes[END_REF], some simulation results which aim at comparing the different estimators, and details about some estimation procedures and the data preprocessing.

The model -definition and properties

Preliminaries

Multivariate extremes. Let

V = {1, . . . , d} for some integer d ≥ 2. A d-variate max-stable distribution G is called simple if its margins are unit-Fréchet, that is, a random vector Z with distribution G satisfies P(Z v ≤ x) = exp(-1/x) for x ∈ (0, ∞) and v ∈ V . Let X = (X v , v ∈ V ) be a random vector with unit-Pareto margins, i.e., P(X v ≤ x) = 1 -1/x for x ∈ [1, ∞) and v ∈ V . Let X i = (X v,i , v ∈ V ) for i = 1, .
. . , n be an independent random sample from the distribution of X. We say that X belongs to the max-domain of attraction of the simple max-stable distribution G,

notation X ∈ D(G), if lim n→∞ P max i=1,...,n X v,i ≤ nz v , v ∈ V = G(z), z ∈ (0, ∞) d .
For more background on max-stable distributions and their domains of attractions, we refer to the reader to Resnick (1987, Chapter 5) and de Haan and Ferreira (2007, Chapter 6).

Throughout the paper the stable tail dependence function (stdf) l of G or X ∈ D(G) will appear frequently. It is defined as

l(x) = lim t→∞ t 1 -P(X v ≤ t/x v , v ∈ V ) = -ln G(1/x v , v ∈ V ), x ∈ [0, ∞) d , ( 1
)
with the obvious limit interpretation if x v = 0 for some v ∈ V . The stdf is closely linked to the exponent function of a simple max-stable distribution in Coles and Tawn (1991, Eq. (2.4)). It is introduced and studied in [START_REF] Huang | Statistics of bivariate extremes[END_REF] and [START_REF] Drees | Best attainable rates of convergence for estimators of the stable tail dependence function[END_REF]; see also later literature in de Haan and Ferreira (2007, Chapter 6) and Beirlant et al. (2004, Chapter 8). The stdf evaluated at x J = (1 {j∈J} , j ∈ V ) is known as an extremal coefficient, of which we make use in Sections 4 and 5.

One of the main objects in our paper is the multivariate Hüsler-Reiss distribution. This absolutely continuous max-stable distribution was introduced in [START_REF] Hüsler | Maxima of normal random vectors: between independence and complete dependence[END_REF] and remains a popular parametric model in recent literature [START_REF] Genton | On the likelihood function of a Gaussian max-stable processes[END_REF][START_REF] Huser | Composite likelihood estimation for the Brown-Resnick process[END_REF][START_REF] Asadi | Extremes on river networks[END_REF][START_REF] Engelke | Estimation of Hüsler-Reiss distributions and Brown-Resnick processes[END_REF][START_REF] Einmahl | A continuous updating weighted least squares estimator of tail dependence in high dimensions[END_REF][START_REF] Lee | Multivariate extreme value copulas with factor and tree dependence structures[END_REF]. It arises as the limiting distribution of partial maxima of a triangular array of row-wise independent and identically distributed random vectors from a multivariate normal distribution with correlation matrix ρ(n) depending on the sample size n. In particular, assume that

lim n→∞ 1 -ρ ij (n) ln n = λ 2 ij ∈ (0, ∞)
for every pair of variables i, j ∈ V and let Λ = (λ 2 ij ) i,j∈V denote this limiting matrix. Note that λ 2 ii = 0 for every i ∈ V . For every subset W ⊆ V and any element u ∈ W let Γ W,u (Λ) be the square matrix of size |W | -1 with elements

Γ W,u (Λ) ij = 2(λ 2 iu + λ 2 ju -λ 2 ij ), i, j ∈ W \ u.
(2) [START_REF] Nikoloulopoulos | Extreme value properties of multivariate t copulas[END_REF] and later [START_REF] Genton | On the likelihood function of a Gaussian max-stable processes[END_REF] and [START_REF] Huser | Composite likelihood estimation for the Brown-Resnick process[END_REF] show that the cumulative distribution function (cdf) as deduced by [START_REF] Hüsler | Maxima of normal random vectors: between independence and complete dependence[END_REF] can be written as

H Λ (z) = exp - u∈V 1 z u Φ d-1 ln z v z u + 2λ 2 uv , v ∈ V \ u; Γ V,u (Λ) , z ∈ (0, ∞) d , ( 3
)
where Φ p ( • ; Σ) denotes the p-variate zero mean Gaussian cdf with covariance matrix Σ. The distribution H Λ in (3) is a simple max-stable distribution. In particular, its margins are unit-Fréchet, whereas [START_REF] Hüsler | Maxima of normal random vectors: between independence and complete dependence[END_REF] originally proposed the distribution in terms of Gumbel margins. Multivariate margins of the d-variate Hüsler-Reiss distribution are Hüsler-Reiss distributions too. The corresponding parameter matrix is obtained by selecting the appropriate rows and columns in the original parameter matrix (see, e.g., Engelke and Hitz, 2020, Example 7). In particular, if

X ∈ D(H Λ ) and if U ⊆ V is non-empty, the stdf l U of X U = (X u , u ∈ U ) is l U (x) = u∈U x u Φ |U \u| ln x u x v + 2λ 2 uv , v ∈ U \ u; Γ U,u (Λ) , x ∈ [0, ∞) U . ( 4
)
Here we write U \ u instead of U \ {u}. In case x u = 0 for some u ∈ U , the corresponding term in the sum in (4) vanishes.

Trees. We will need some notions from graph theory. A graph is a pair G = (V, E) where V = {1, . . . , d} is the set of nodes or vertices and E ⊆ {(a, b) ∈ V × V : a = b} is the set of edges. Edges will also be denoted by e = (a, b) ∈ E. The number of vertices in a subset U ⊆ V will be denoted by |U |, while d is reserved for

|V | only. A graph is undirected if (a, b) ∈ E is equivalent to (b, a) ∈ E. A path (u v) from node u to node v is a collection {(u 0 , u 1 ), (u 1 , u 2 ), . . . , (u n-1 , u n )} of distinct, directed
edges such that u 0 = u and u n = v. An undirected tree is an acyclic undirected graph T = (V, E) such that for every pair of distinct nodes a and b there is a unique path (a b).

Model definition

Let T = (V, E) be an undirected tree with node set V = {1, . . . , d} and let ξ = (ξ v , v ∈ V ) be a random vector with joint cdf F and continuous margins

F v (z) = P(ξ v ≤ z) for z ∈ R and v ∈ V . Let the random vector X = (X v , v ∈ V ) be defined as X v = 1/ 1 -F v (ξ v ) for every v ∈ V . Because the functions F v for v ∈ V are continuous, the marginal distributions of X are unit-Pareto.
We assume that X is in the max-domain of attraction of the Hüsler-Reiss distribution H Λ in (3) with Λ = (λ 2 ij ) i,j∈V having the following structure linked to the tree T : there exists a vector θ = (θ e ) e∈E of positive scalars with θ ab = θ ba and such that Λ = Λ(θ) where

Λ(θ) ij = λ 2 ij (θ) = 1 4 e∈(i j)
θ 2 e , i, j ∈ V, i = j.

(5)

The assumption can thus be written compactly as

X ∈ D(H Λ(θ) ) for some θ ∈ (0, ∞) E .
The motivation for the proposed structure is that H Λ(θ) contains in its max-domain of attraction a certain graphical model with respect to T as explained in Section 2.3. Still, it is to be noted that, despite the structure of the parameter matrix, H Λ(θ) itself does not and cannot satisfy any Markov properties with respect to the tree T : by [START_REF] Papastathopoulos | Conditional independence among max-stable laws[END_REF], max-stable distributions with continuous joint densities cannot possess any non-trivial conditional independence properties.

In the parametrization in (5) the extremal dependence in ξ and in X depends on a vector θ = (θ e , e ∈ E) of d -1 free parameters, indexed by the edges of the tree. The main theme in this paper concerns inference on the parameter vector θ in case some of the variables ξ v are latent (unobservable). The first question is whether all edge parameters θ e are still identifiable from (4) when Λ = Λ(θ) and when U V contains the indices of variables that can still be observed. For the Seine network in Figure 1, for instance, there are d = 7 variables in total, of which two are latent. A necessary and sufficient criterion for parameter identifiability is given in Proposition 3.1 below. Provided the criterion is fulfilled, the second question is how to estimate the parameters. Three estimation methods are proposed in Section 4 and illustrated in Section 5.

Note that the random vector ξ itself does not necessarily belong to the max-domain of attraction of some max-stable distribution. The reason is that we do not impose that the marginal distributions of ξ are in the max-domain of attraction of some univariate extreme value distributions. To focus on the tail dependence of ξ, we standardize its margins and formulate the assumption in terms of X.

Motivation of the structured Hüsler-Reiss model

To motivate the structured Hüsler-Reiss parameter matrix Λ(θ) in (5), we construct a graphical model Z * that satisfies the global Markov property with respect to the undirected tree T = (V, E) and such that Z * ∈ D(H Λ(θ) ). Besides serving as a motivation, the auxiliary model Z * plays another important role: in view of Segers (2020, Theorem 2) we are able to project certain asymptotic properties that hold for Z * to X.

For disjoint subsets A, B, C of V , the expression A ⊥ ⊥ T B | C means that C separates A from B in T , also called graphical separation, i.e., all paths from A to B pass through at least one vertex in C. Let Z * be defined on a probability space (Ω, B, P). Conditional independence of Z * A and Z *

B

given Z * C will be denoted by

Z * A ⊥ ⊥ P Z * B | Z * C ; here Z * A = (Z * a , a ∈ A
) and so on. If P = P(Z * ∈ • ) is the law of Z * , we say that the tree T is an independence map (I-map) of P if for any disjoint subsets [START_REF] Koller | Probabilistic Graphical Models: Principles and Techniques[END_REF]. This assumption is equivalent to the assumption that Z * obeys the global Markov property with respect to T [START_REF] Lauritzen | Graphical Models[END_REF]. The law of the random vector Z * = (Z * v , v ∈ V ) is defined by the following two assumptions: (Z1) Z * satisfies the global Markov property (6) with respect to the undirected tree T = (V, E);

A, B, C of V it holds that A ⊥ ⊥ T B | C =⇒ Z * A ⊥ ⊥ P Z * B | Z * C (6)
(Z2) every pair of variables (Z * a , Z * b ) on adjacent nodes (a, b) = e ∈ E has a bivariate Hüsler-Reiss distribution with parameter θ e ∈ (0, ∞) and unit-Fréchet margins, i.e., the special case of (4) with U = {a, b} and λ 2 ab = θ 2 e /4.

The law of Z * is absolutely continuous and its joint density function factorizes in terms of the bivariate Hüsler-Reiss densities along pairs of variables on adjacent nodes through the Hammersley-Clifford theorem; see Appendix A.5 where we describe how to sample from Z * . Moreover, for e = (a, b) ∈ E and if Z has distribution H Λ(θ) , the law of (Z * a , Z * b ) is the same as the one of (Z a , Z b ). However, unless d = 2, the law of Z * is itself not max-stable and thus not equal to the one of Z. One way to see this is to note that by [START_REF] Papastathopoulos | Conditional independence among max-stable laws[END_REF], the law of Z cannot satisfy the global Markov property with respect to T .

Let (M e , e ∈ E) be a random vector of independent lognormal random variables with ln M e ∼ N (-θ 2 e /2, θ 2 e ) for each e ∈ E. In view Theorem 1 and Corollary 1 in Segers (2020), we have the convergence in distribution

(Z * v /Z * u , v ∈ V \ u) | Z * u > x d -→ (Ξ u,v , v ∈ V \ u) = e∈(u v) M e , v ∈ V \ u , x → ∞, (7)
for every u ∈ V . For every u ∈ V the vector (Ξ u,v , v ∈ V \ u) is called a tail tree. The multiplicative structure in (7) goes back to the theory of extremes of Markov chains due to [START_REF] Smith | The extremal index for a Markov chain[END_REF], [START_REF] Perfekt | Extremal behaviour of stationary Markov chains with applications[END_REF], [START_REF] Yun | The extremal index of a higher-order stationary Markov chain[END_REF] and [START_REF] Segers | Multivariate regular variation of heavy-tailed Markov chains[END_REF]. Note that a chain can be seen as a tree with a single branch.

The vector (ln Ξ u,v , v ∈ V \ u) is a linear transformation of a Gaussian random vector and is therefore itself Gaussian. Its mean vector µ V,u (θ) and its covariance matrix Σ V,u (θ) have elements

{µ V,u (θ)} v = - 1 2 e∈(u v) θ 2 e , v ∈ V \ u, ( 8
)
{Σ V,u (θ)} ij = e∈(u i)∩(u j) θ 2 e , i, j ∈ V \ u . ( 9
)
Hence for every u ∈ V and as x → ∞, we have the convergence in distribution

(ln Z * v -ln Z * u , v ∈ V \ u) | Z * u > x d -→ (ln Ξ u,v , v ∈ V \ u) ∼ N |V \u| µ V,u (θ), Σ V,u (θ) , ( 10
)
where N p is the p-variate normal distribution. By construction, Σ V,u (θ) is a covariance matrix and hence positive semi-definite for any θ ∈ (0, ∞) d-1 ; it is actually positive definite since the vector (ln Ξ u,v , v ∈ V \ u) is the result of an invertible linear transformation applied to the vector (ln M e , e ∈ E) of independent and non-degenerate normal random variables. The matrix Σ V,u (θ) is moreover the same as the matrix Γ W,u (Λ) in (2) with W = V and Λ = Λ(θ) in (5):

Σ V,u (θ) ij = e∈(u i)∩(u j) θ 2 e = 1 2   e∈(u i) θ 2 e + e∈(u j) θ 2 e - e∈(i j) θ 2 e   = 2(λ 2 iu + λ 2 ju -λ 2 ij ) = Γ V,u Λ(θ) ij , i, j ∈ V \ u. ( 11
)
In the second equality it is needed to divide by two because the parameters on shared edges are added twice. In addition, the Hüsler-Reiss parameters λ 2 uv are proportional to the means:

2λ 2 uv = 1 2 e∈(u v) θ 2 e = -{µ V,u (θ)} v , v ∈ V \ u. ( 12
) Proposition 2.1. Let T = (V, E) be a tree. If the law of Z * = (Z * v , v ∈ V ) is given by (Z1)-(Z2) above, then Z * ∈ D(H Λ(θ) ) with Λ(θ) in (5).
The proof is given in Appendix A.3 and relies on the properties of Z * mentioned above, in particular on (10). By constructing a graphical model with respect to T in the max-domain of attraction of H Λ(Θ) , we have argued that the latter is a sensible dependence model for extremes of graphical models on trees. Moreover, it follows that any random vector X = (X v , v ∈ V ) with unit-Pareto margins and in the max-domain of attraction of H Λ(θ) shares property (10) with Z * . Corollary 2.2. Let T = (V, E) be a tree and let X = (X v , v ∈ V ) have unit-Pareto margins and belong to D(H Λ(θ) ) with Λ(θ) as in (5) for a vector θ = (θ e , e ∈ E) of positive scalars. Then for every u ∈ V , we have

(ln X v -ln X u , v ∈ V \ u) | X u > t d -→ N |V \u| µ V,u (θ), Σ V,u (θ) , t → ∞. ( 13
)
Proof. The max-domain of attraction condition X ∈ D(H Λ(θ) ) is known to be equivalent to convergence of the measures t P(X/t ∈ • ) as t → ∞ to the exponent measure of H Λ(θ) (Resnick, 1987, Proposition 5.17). Such measure convergence is in turn equivalent to convergence in distribution of X/X u | X u > t as t → ∞ for every u ∈ V to a limit that can be written in terms of the said exponent measure (Segers, 2020, Theorem 2). But for X replaced by Z * , the limiting conditional distribution was found to be a certain multivariate lognormal distribution in (7). The equivalence between ( 7) and ( 10) with Z * replaced by X is clear by the continuous mapping theorem.

The random vector X in Corollary 2.2 does not need to be a graphical model with respect to T . The convergence in (13) appears in Engelke et al. (2015, Theorem 2) for a general random vector with standardized margins and in the max-domain of attraction of a Hüsler-Reiss distribution. With Corollary 2.2 we arrive at the same result but through the properties of the auxiliary model Z * . The convergence in ( 13) is used to build two estimators in the next section.

In [START_REF] Engelke | Graphical Models for Extremes[END_REF], a notion of conditional independence different from the classical one is introduced in the context of multivariate Pareto distributions. When specialized to the Pareto distribution associated to a max-stable Hüsler-Reiss distribution, it yields certain restrictions on the Hüsler-Reiss parameter matrix Λ. In case the conditional independence relations are the ones induced by a tree through graphical separation, the structure of the parameter matrix is the same as the one in (5). We explain the connection in Appendix A.1. Here we just emphasize that in [START_REF] Engelke | Graphical Models for Extremes[END_REF], no graphical model in the classical sense of the term is constructed that belongs to the max-domain of attraction of H Λ(θ) . The way we arrive at the structure of Λ(θ) via the graphical model Z * in (Z1)-( Z2) is thus entirely different from their approach.

Finally, quite another tree-induced structure of the Hüsler-Reiss parameter matrix is proposed in [START_REF] Lee | Multivariate extreme value copulas with factor and tree dependence structures[END_REF]. We provide a comparison in Appendix A.2.

Latent variables and parameter identifiability

A typical application of our model arises in relation to quantities measured on river networks that have a tree-like structure. It is natural to associate a node to an existing measurement station or to locations where two river channels meet (junction) or one channel splits (split) even if there is no measurement station there. Stations are supposed to generate data for the quantity of interest, so for any node associated to a station there is a corresponding variable. In practice, junctions/splits may lack measurements, and this means that there are nodes in the tree with latent variables. Nodes with latent variables are those labelled 2 and 5 in the Seine network in Figure 1.

A naive approach to the presence of latent variables would be to ignore them, that is, to remove the corresponding nodes and all edges incident to them. This will yield a disconnected graph, making it necessary to add edges in some arbitrary way so as to obtain a tree again. In Figure 2 for instance, if node 2 is suppressed, there are three possible ways to reconnect the remaining nodes and form a tree. Each implies a different structured Hüsler-Reiss parameter matrix and thus a different dependence model.

X 1 2 X 3 X 4 θ 12 θ 23 θ 24 X 1 X 4 X 3 β 14 β 13 X 1 X 4 X 3 β 14 β 34 X 1 X 4 X 3 β 13 β 34
Figure 2: The first tree from the left has four nodes where node 2 has a latent variable. If node 2 is suppressed, there are three possible ways to reconnect the three remaining nodes into a tree again.

In this paper we do not modify the original tree but take the latent variables into account. Let T = (V, E) be an undirected tree and consider the Hüsler-Reiss distribution (3) with parameter matrix Λ = Λ(θ) in (5). When there are nodes with latent variables, the question is whether it is still possible to identify the d -1 free edge parameters θ e from the distribution of the subvector of observable variables only. Let U ⊆ V denote the set of indices of the observable variables. On the one hand, Eq. ( 13) implies

(ln X v -ln X u ) v∈U \u | X u > t d -→ N |U \u| µ U,u (θ), Σ U,u (θ) , t → ∞, (14)
with µ U,u (θ) and Σ U,u (θ) as in ( 8) and ( 9) but with V replaced by U . On the other hand, µ U,u (θ) and Σ U,u (θ) together determine the stdf l U of the subvector X U in (4) through the identities ( 11) and ( 12). The question is thus whether the parameter vector θ is still identifiable from the |U \ u|-variate normal distributions on the right-hand side of ( 14), where u ranges over U .

Example. Let X = (X a , X b , X c ) have unit-Pareto margins and suppose that X ∈ D(H Λ(θ) ) where Λ(θ) is as in ( 5) with respect to the chain tree T with nodes V = {a, b, c} and edges between a and b and between b and c. Since a parameter is linked to each (undirected) edge of the graph, the parameter vector is θ = (θ ab , θ bc ). Suppose the variable X b is latent. By ( 14) we have

ln X c -ln X a | X a > t d -→ N -(θ 2 ab + θ 2 bc )/2, (θ 2 ab + θ 2 bc ) , t → ∞.
It is clear that from the limiting normal distribution, we cannot identify θ ab and θ bc .

In this section it is shown that as long as all nodes with missing variables have degree at least three, the parameters associated to the Hüsler-Reiss distribution of the full vector are still identifiable and hence there is no need to change the tree. To this end, note that by ( 8), ( 11) and ( 12) with V replaced by U ⊆ V such that u ∈ U , the mean vectors µ U,u (θ) and covariance matrices Σ U,u (θ) are determined completely by the path sums

p ab = e∈(a b) θ 2 e = 4λ 2 ab , a, b ∈ U, ( 15
)
and that, vice versa, the values of these path sums are determined by the vectors µ U,u (θ) and the matrices Σ U,u (θ). If we know the distribution of X U = (X u , u ∈ U ), we can compute the values of these sums, and if we know these sums, we can compute the stdf l U of X U . The question is thus whether or not the edge parameters θ e are identifiable from the values of the path sums p ab for a, b ∈ U . According to the following proposition, there is a surprisingly simple criterion to decide whether this is the case or not.

Proposition 3.1. Let T = (V, E) be an undirected tree and let X = (X v , v ∈ V ) have unit Pareto margins and be in the max-domain of attraction of the structured Hüsler-Reiss distribution H Λ in (3) with parameter matrix Λ = Λ(θ) in (5). Let U ⊆ V be the set of nodes corresponding to the observable variables. The parameter vector θ is identifiable from X U = (X u , u ∈ U ) if and only if every node u ∈ V \ U has degree at least three.

Proof. Necessity. Assume that the elements of the parameter θ ∈ (0, ∞) d-1 are uniquely identifiable. Let Ū = V \ U = ∅ be the set of nodes with latent variables. We need to show that every v ∈ Ū has degree d(v) at least 3. We will do this by contraposition. As a tree is connected by definition, there cannot be a node of degree zero. First, assume there is v ∈ Ū such that d(v) = 1. The node v must be a leaf node, and in this case there is no path (a b) with a, b ∈ U that passes by v, and thus θ 2 uv , with u the unique neighbor of v, does never appear in the sum (15). Hence θ uv is not identifiable, which is a contradiction to the assumption.

Second, assume there exists v ∈ Ū with d(v) = 2. Then v has exactly two neighbors, i and j, say. Every path sum p ab for a, b ∈ U will contain either the sum of the squared parameters, θ 2 iv + θ 2 jv , or neither of these. Hence, the individual edge parameters θ iv and θ jv are not identifiable, yielding a contradiction. (This generalizes the example given before the statement of the proposition.)

Sufficiency. Assume that all nodes with latent variables are of degree three or more. Let e = (u, v) ∈ E. We will find a linear combination of the path sums (15) equal to θ 2 uv . If u, v ∈ U , then the one-edge path sum p uv = θ 2 uv already meets the condition. Suppose that u ∈ Ū . By assumption, u has at least two other neighbors besides v, say w and x. If v ∈ U , then put v = v. Otherwise, start walking at v away from u until you encounter the first visible node, say v ∈ U . There must always be such a node, since V is finite and since all leaves are observable by assumption. Similarly, let ŵ ∈ U and x ∈ U be the first visible nodes encountered when walking away from u and starting in w and x, respectively. Note that v, ŵ, and x are all different since otherwise the graph would contain a non-trivial cycle, which is not possible in the case of a tree. We can thus observe the sums

p v ŵ = p vu + p u ŵ , p v x = p vu + p ux , p ŵx = p ŵu + p ux .
Since p yz = p zy for every y, z ∈ V , the previous identities constitute three linear equations in three unknowns that can be solved explicitly, producing the values of p uv , p u ŵ, p ux . In particular, summing the first two equations, subtracting the third, and dividing by two, we find

p uv = 1 2 p v ŵ + 1 2 p v x -1 2 p ŵx . If v ∈ U , then v = v
, and (u v) = {e}, so that the above equation shows how to combine path sums in a linear way to extract p uv = θ 2 e . If v ∈ U , then we can repeat the same procedure with u replaced by v. The result is a formula expressing p vv as a linear combination of three visible path sums. Now since

θ 2 e = p uv -p vv ,
we have found a way to extract θ 2 e by a linear combination of at most six visible path sums.

The proof of Proposition 3.1 consists in solving the equations ( 15) with p ab as known and θ 2 e as unknown. Clearly, this is a linear system of equations and the question is thus whether the coefficient matrix defining the system has full column rank. It is an open question how to write down this matrix, which contains only zeroes and ones, in terms of the tree's adjacency matrix in such a way that an algebraic criterion on the latter matrix can be formulated.

The identifiability criterion in Proposition 3.1 allows nodes with latent variables to be adjacent and still counting in the computation of each other's degree. Consider for instance the following tree:

X 1 X 2 X 3 X 4 5 6 θ 15 θ 25 θ 36 θ 46 θ 56
The variables at the adjacent nodes 5 and 6 are latent. Both nodes have degree three and each of the five edge parameters θ e can be solved from the path sums p ab between nodes a, b ∈ {1, . . . , 4}.

The previous example may give the impression that for the identifiability criterion to hold it is actually enough that all variables on leaf nodes are observable. Although the latter property is indeed necessary, it is not sufficient, as illustrated by the example before Proposition 3.1.

Estimation

Let T = (V, E) be an undirected tree with nodes V = {1, . . . , d} and let (ξ v,i , v ∈ V , i = 1, . . . , n) be an independent random sample from the distribution of ξ satisfying the assumptions in Section 2.2. Further, let U ⊆ V be the set of indices of observable variables and assume that every u ∈ V \ U has degree at least three, so that, by Proposition 3.1, the Hüsler-Reiss edge parameters θ = (θ e , e ∈ E) in the definition of Λ(θ) in ( 5) are identifiable from the distribution of the subvector

ξ U = (ξ v , v ∈ U ).
We propose three methods for estimating the parameter vector θ. The first one, called moment estimator (Section 4.1), builds upon the one introduced in [START_REF] Engelke | Estimation of Hüsler-Reiss distributions and Brown-Resnick processes[END_REF]. The second estimator comes from the optimization of a composite likelihood function (Section 4.2). The third estimator, finally, is based on bivariate extremal coefficients (Section 4.3) and on the method in [START_REF] Einmahl | A continuous updating weighted least squares estimator of tail dependence in high dimensions[END_REF]. All estimators are functions of the subvectors (ξ

U,i ) = (ξ v,i , v ∈ U ) for i = 1, . . . , n only.
An important remark for this whole section is related to the fact that X as introduced in Section 2.2 should have unit Pareto margins, obtained after the transformation

X v = 1/(1 -F v (ξ v )) where F v is the marginal distribution function of ξ v for v ∈ V .
It is unrealistic to assume that the functions F v are known, so in practice we use their empirical versions, Fv,n (x) = n i=1 1(ξ v,i ≤ x) /(n + 1). The estimates of the edge parameters will then be based upon the sample X1 , . . . , Xn with coordinates

Xv,i = 1 1 -Fv,n (ξ v,i ) , v ∈ U, i = 1, . . . , n,
considered as a random sample from the distribution of X U = (X u , u ∈ U ).

A variable indexed by the double subscript W, i will denote the i-th observation of variables on nodes belonging to the set W ⊆ U : for instance XW,i = ( Xv,i , v ∈ W ). Such vectors are taken to be column vectors of length |W |. When W = U we just write Xi . [START_REF] Engelke | Estimation of Hüsler-Reiss distributions and Brown-Resnick processes[END_REF] introduce an estimator of the matrix Λ of the Hüsler-Reiss distribution, based on sample counterparts of the matrices Γ W,u (Λ) in (2). Relying on (11) with V replaced by W ⊆ U , we will apply their method to the vector of observable variables and then add a least-squares step to extract the edge parameters θ e .

Method of moments estimator

As a starting point we take the result in ( 14) and as suggested by [START_REF] Engelke | Estimation of Hüsler-Reiss distributions and Brown-Resnick processes[END_REF] for given k ∈ {1, . . . n} we obtain the log-differences

∆ uv,i = ln Xv,i -ln Xu,i , ( 16
)
for u, v ∈ U and for i ∈ I u = {i = 1, . . . , n : Xu,i > n/k}. The proposed estimators of µ U,u and Σ U,u are respectively the sample mean vector

μU,u = 1 |I u | i∈Iu (∆ uv,i , v ∈ U \ u)
and the sample covariance matrix

ΣU,u = 1 |I u | i∈Iu (∆ uv,i -μU,u , v ∈ U \ u)(∆ uv,i -μU,u , v ∈ U \ u) .
To estimate the vector of edge parameters θ = (θ e , e ∈ E), we propose the least squares estimator

θMM n,k = arg min θ∈(0,∞) E u∈U ΣU,u -Σ U,u (θ) 2 F . ( 17
)
where • F is the Frobenius norm. In this way, we take advantage of the empirical covariance matrices ΣU,u for each u ∈ U and thus of each exceedance set I u . In ( 17), for each u ∈ U , we consider the covariance matrix of the log-differences ∆ uv,i for all v ∈ U \ u. However, if v is far away from u in the tree, then the extremal dependence between ξ u and ξ v may be weak and the difference ∆ uv,i may carry little information. Therefore, we propose a modified estimator where, for each u ∈ U , we limit the scope to a subset W u ⊆ U of observable variables indexed by nodes near u, producing the estimator

θMM n,k = arg min θ∈(0,∞) E u∈U ΣWu,u -Σ Wu,u (θ) 2 F . ( 18
)
Besides being simpler to compute, the modified estimator (18) performed better than the one in (17) in Monte Carlo experiments. One possible explanation is that by excluding pairs with weak extremal dependence, the bias of the estimator diminishes.

When choosing the sets W u , care needs to be taken that the parameter vector θ is still identifiable from the collection of covariance matrices Σ Wu,u (θ) for u ∈ U . The set of path sums p ab for a, b ∈ U in Proposition 3.1 is now reduced to the set of the path sums p ab for a, b ∈ W u and u ∈ U . Whether or not these are still sufficient to identify θ needs to be checked on a case-by-case basis. This issue is illustrated in Appendix A.4.

Composite likelihood estimator

The composite likelihood estimator (CLE) is again based on the result in ( 14). This time however we maximize a composite likelihood function with respect to the parameter θ directly. The composite likelihood function consists of multiplication of likelihoods which are defined on subtrees.

As for the method of moments estimator in Section 4.1, we consider for each u ∈ U a set W u ⊆ U of nodes that are close to u in the tree, taking care to include sufficiently many variables so that the edge parameters are still identifiable (Appendix A.4). Recall the log-differences ∆ uv,i in ( 16) and the exceedance set I u right below (16). Let φ p ( • ; Σ) be the density function of the centered p-variate normal distribution with covariance matrix Σ. The composite likelihood estimator θCLE n,k is the maximizer of the composite likelihood

L θ; {∆ uv,i : v ∈ W u \ u, i ∈ I u , u ∈ U } = u∈U i∈Iu φ |Wu\u| (∆ uv,i ) v∈Wu -µ Wu,u (θ); Σ Wu,u (θ) .
We aggregate the likelihoods of the different normal distributions for all u ∈ U treating the samples of log-differences as independent, although they are not. Results from Monte Carlo simulation experiments (Appendix A.5) show that the performance of the CLE is comparable to the one of the moment estimator and the extremal coefficient estimator.

Other estimation methods based on locally defined likelihoods are used by [START_REF] Engelke | Graphical Models for Extremes[END_REF] and [START_REF] Lee | Multivariate extreme value copulas with factor and tree dependence structures[END_REF]. The method of [START_REF] Engelke | Graphical Models for Extremes[END_REF] estimates the parameters associated to each clique separately. For trees this means that there are d -1 one-variate likelihood functions to optimize, a problem which is doable even in trees with many nodes. A problem with this estimator is that it is inapplicable if there are latent variables because there will always be an adjacent pair of variables with one of them being an unobservable, and making it impossible to estimate the corresponding edge parameter. The estimator of Lee and Joe ( 2018) is based on pairwise likelihoods, which can be any pairs, not only adjacent pairs as in the estimator of [START_REF] Engelke | Graphical Models for Extremes[END_REF]. It is obtained by optimizing the composite likelihood which consists of multiplying the pairwise likelihoods. This estimator is applicable when there are latent variables as long as all possible pairs between the observed variables are included in the composite likelihood function. It is close in spirit to the pairwise extremal coefficients estimator considered next.

Pairwise extremal coefficients estimator

The pairwise extremal coefficients estimator (ECE), defined for general tail dependence models in [START_REF] Einmahl | A continuous updating weighted least squares estimator of tail dependence in high dimensions[END_REF], is based on the bivariate stable tail dependence function (stdf) in (4). It minimizes the weighted distance between a non-parametric estimate and the fitted parametric stdf.

Let l be the stdf in (1) and recall that the extremal coefficient associated to a node set J ⊆ V is defined as

l(x J ) = l J (1, . . . , 1) = lim t→∞ t P max j∈J X j > t , ( 19
)
where x J = (1 {j∈J} , j ∈ V ) and where l J is the stdf of the subvector X J . For the Hüsler-Reiss distribution with parameter matrix Λ and for a pair of nodes J = {u, v}, the bivariate extremal coefficient is just l J (1, 1) = 2Φ(λ uv ), with Φ the standard normal cdf. In case Λ = Λ(θ) in ( 5), the pairwise extremal coefficient depends on the path sum

p uv = e∈(u v) θ 2 e via l J (1, 1; θ) = 2Φ( √ p uv /2), J = {u, v}. ( 20
)
The non-parametric estimator of the stdf dates back to [START_REF] Drees | Best attainable rates of convergence for estimators of the stable tail dependence function[END_REF] and yields the following estimator for the extremal coefficient l J (1, . . . , 1) for J ⊆ V :

lJ;n,k (1, . . . , 1) = 1 k n i=1 1 max j∈J n Fj,n (ξ j,i ) > n + 1/2 -k . ( 21
)
Let Q ⊆ {J ⊆ U : |J| = 2} be a collection of pairs of nodes associated to observable variables and put q = |Q|, ensuring that q ≥ |E| = d -1, the number of free edge parameters. The pairwise extremal coefficients estimator (ECE) of θ is

θECE n,k = arg min θ∈(0,∞) E J∈Q lJ;n,k (1, 1) -l J (1, 1; θ) 2 . ( 22
)
If Q is the collection of all possible pairs of nodes in U , then the pairwise extremal coefficients (20) give us access to all path sums p ab for a, b ∈ U , and Proposition 3.1 guarantees we can identify θ. If, however, Q is a smaller set of pairs, then the identifiability of θ from the resulting path sums needs to be checked on the case at hand.

High water levels on the Seine network

We have chosen to present an application that allows us to demonstrate the identifiability criterion outlined in Section 3. Data were collected from http://www.hydro.eaufrance.fr, a web-site of the french Ministry of Ecology, Energy and Sustainable Development, and span the period from January 1987 to April 2019 with gaps for some of the measurement stations. The data represent water levels, in cm, at five locations on the Seine river: Paris, Meaux, Melun, Nemours and Sens. The map on Figure 1 shows part of the actual Seine network. The schematic representation of the graphical model used in the estimation is shown in Figure 3. The tree has d = 7 nodes, two of which are associated to latent variables. Since both these nodes have degree equal to three, Proposition 3.1 guarantees we can still identify all six edge parameters θ 1 , . . . , θ 6 . For more information on the data set, some summary statistics and details on data preprocessing, we refer to Appendix A.6.

Estimates and confidence intervals

We used all three estimators in Section 4 to obtain estimates of the six parameters of extremal dependence. For the pairwise extremal coefficient estimator (ECE) it is possible to calculate standard errors thanks to the asymptotic distribution derived in Einmahl et al. (2018, Theorem 2.2).

Computational details for the standard errors follow in Appendix A.7. The distributions of the MME and CLE are not known so we computed bootstrapped confidence intervals, known as basic bootstrap confidence intervals (Davison and Hinkley, 1997, Chapter 2), by resampling from the data.

The EC estimates and their 95% confidence intervals are displayed in Figure 4 for two of the parameters, namely θ 1 and θ 4 . The confidence intervals using the MME and CLE are narrower as can be seen from Figure 5. The plots for θ 2 , θ 5 , θ 6 are similar to the one for θ 1 : the 95% confidence intervals never include zero, suggesting that the extremal dependence between the corresponding variables is not perfect and hence that the edges cannot be collapsed. In Section 3, we alluded to the possibility of circumventing the issue of latent variables by suppressing nodes and redrawing edges. The fact that the confidence intervals do not include zero indicate that doing so would have produced a misleading picture of extremal dependence.

The plot of θ 3 , similarly to the plot of θ 4 , does contain a segment over k where the lower confidence bound reaches zero: for θ 4 this is approximately k ∈ [260, 360], while for θ 3 it is k ∈ [90, 180]. Although the confidence intervals for θ 3 and θ 4 indicate some instability of the estimated parameters, we believe that collapsing the edges is not advisable, especially in networks with many more unobservable variables. Moreover, the river distance, which is one of the important factors in tail dependence [START_REF] Asadi | Extremes on river networks[END_REF], is rather long between node 2 and Melun and between Melun and node 5, so that there is no physical motivation for collapsing the corresponding edges.

For a point estimate per parameter we need to average out over a range of k. The chosen range per estimator and per parameter need not be the same. As a rule we select a range around the beginning where the estimates start stabilizing around a certain level, omitting the most volatile part for relatively small k. Most of the time we thus consider k ∈ [100,220]. In this way we end up with the point estimates displayed for comparison in Figure 5.

Given the similarities between the MME and CLE, the estimates are pooled in an average of the two for each parameter. Parameters -estimates and 95% confidence intervals index theta q q q q q q q MME CLE ECE

Figure 5: Parameters -estimates and confidence intervals. The confidence intervals of the moment and composite likelihood estimators are bootstrapped, namely θ ∈ [2 θ -q * 0.975 , 2 θ -q * 0.025 ], where q * α is the α-quantile of the bootstrapped distribution of θ.

Considerations on the goodness-of-fit of the model

The Hüsler-Reiss family would not be an appropriate extremal dependence model if some of the variables would exhibit asymptotic independence. As kindly suggested by a Reviewer, we compared non-parametric estimates of multivariate tail dependence coefficients

P(min v∈W X v > t | X u > t) = t P(min v∈W ∪u X v > t)
for W U and u ∈ U \ W at finite thresholds t = n/k with their postulated limits as t → ∞ based on the fitted Hüsler-Reiss stdf. The results (not shown) supported the hypothesis of asymptotic dependence for nearly all subvectors X W ∪u of variables.

To assess how well the model from Section 2.2 fits the data, we compare non-parametric and model-based estimates of quantities describing extremal dependence, such as pairwise and triple-wise extremal coefficients and the Pickands dependence function. For J ⊆ U , recall the extremal coefficient l J (1, . . . , 1) in ( 19) and its non-parametric estimate lJ;n,k (1, . . . , 1) in ( 21), also called empirical extremal coefficient. The extremal coefficient l J (1, . . . , 1) is always between 1 and |J|, corresponding to perfect extremal dependence and to extremal independence, respectively.

Figure 6 compares the model-based extremal coefficients obtained from (4) by plugging in parameter estimates with the empirical counterparts for pairs and triples J ⊆ U . At least visually the fit is quite good for both estimators considered, which are the average of the CLE and MME on the one hand and the ECE on the other hand. Note that the ECE in ( 22) is constructed explicitly to ensure that the model-based pairwise extremal coefficients fit the empirical ones as closely as possible. It is therefore only natural that the extremal coefficients based on the ECE fit the empirical ones 1.0 1.2 1.4 1.6 1.8 2.0 1.0 1.2 1.4 1.6 1.8 2.0

Extremal coefficients -Pairs model based non-parametric q q q q q q q q q q q average of MME and CLE ECE 1.0 1.5 2.0 2.5 3.0 1.0 1.5 2.0 2.5

3.0

Extremal coefficients -Triples model based non-parametric q q q q q q q q q q q average of MME and CLE ECE best. A more comprehensive comparison of the finite-sample performance of MME, CLE and ECE is reported in a simulation study in Appendix A.5.

It should be noted that it is impossible to compute empirical extremal coefficients involving latent variables. Model-based estimates of such extremal coefficients can still be computed however, thanks to the identifiability of the parameter vector θ.

As another visual check of the goodness-of-fit of the assumed model we consider the bivariate Pickands dependence function, usually denoted by A(w) for w ∈ [0, 1]. For the Hüsler-Reiss extreme-value distribution at the pair J = {u, v}, it is equal to

A u,v (w; θ) = l J (1 -w, w; θ) = (1 -w) Φ ln( 1-w w ) + 1 2 p uv √ p uv + w Φ ln( w 1-w ) + 1 2 p uv √ p uv ,
with p uv as in ( 15). Hence the model-based estimator of A u,v (w; θ) is A u,v (w; θn,k ) where θn,k can be the average MME/CLE or the ECE. The non-parametric counterpart of the Pickands dependence function is

Âu,v (w) = 1 k n i=1 1{n Fu,n (ξ u,i ) > n -k(1 -w) + 1/2 or n Fv,n (ξ v,i ) > n -kw + 1/2}.
The model-based Pickands dependence function is compared to the empirical counterpart in Figure 7.

The plot is complemented with non-parametric 95% confidence intervals for A(w) computed by the bootstrap method introduced in Kiriliouk et al. (2018, Section 5). The general idea of the method is to approximate the distribution of

√ k( ln,k -l) by the distribution of √ k( l * n,k -lβ n,k ) where l * n,k
is the empirical stable tail dependence function based on the ranks of a sample of size n from the empirical beta copula and lβ n,k is the stdf based on the empirical beta copula using the ranks of the original sample (ξ v,i , v ∈ U ) for i = 1, . . . , n. A detailed description of the derived bootstrap confidence intervals is provided in Appendix A.8.

Flow-connectedness and tail dependence

In a study by [START_REF] Asadi | Extremes on river networks[END_REF] of data from the Danube, it was found that a key factor for extremal dependence between two locations is whether or not they are flow connected. Two locations are flow connected if one of them is downstream of the other one. Flow connectedness often dominates river distance or Euclidean distance in importance: variables on distant nodes that are flow connected might have stronger tail dependence than variables on nodes that are nearby but not flow connected.

This effect is confirmed in our data too and is illustrated in Figure 8. The cities of Sens and Nemours are not flow connected but the Euclidean and river distance between them is smaller than Sens Nemours q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q qq q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0 Sens Paris q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q Figure 8: Scatterplots of uniform transformed data, Fv,n (ξ v,i ), v ∈ U, i = 1, . . . , n, for two pairs of locations. Sens and Nemours (left) are not flow connected while Sens and Paris (right) are flow connected. It can be seen from the Seine map in Figure 1 that the river and Euclidean distance from Sens to Nemours is much smaller than the one from Sens to Paris. However the tail dependence seems to be stronger for the second pair of locations. the one between the flow connected cities of Sens and Paris. Still, the tail dependence seems to be stronger for the flow connected pair of locations.

Figure 9 illustrates the tail dependence in the Seine network through a heat map of the pairwise extremal coefficients. Pairs which are flow connected are indeed the ones with stronger tail dependence (smaller extremal coefficient). According to both estimators the strongest tail dependence is to be found between Paris and the locations at node 2, node 5 and Melun.

Suppressing latent variables

In Section 3 we alluded to the possibility of suppressing nodes with latent variables. Here we illustrate that method and compare the results with those presented so far. After removing nodes 2 and 5 from the Seine graph in Figure 3, there is no unique way of reconnecting the remaining five nodes into a tree. Two possible structures for the reduced Seine graph are presented in Figure 10. We opt for the right-hand graph and refer to the model associated to that tree as model B. Model A will refer to the one associated to the original graph in Figure 3. In Figure 11, we compare the extremal coefficients induced by model B to the empirical ones and to those induced by model A. The extremal coefficients resulting from both models turn out to be rather close, the little black circles lying almost on the diagonal in all four plots. The reason may be that the original tree is small and that not many nodes have been suppressed, whereas in case of many latent variables, the impact of suppressing them may be large. Furthermore, the results may depend on the particular choice of the reduced tree: out of many possibilities two of which shown in Figure 10 we selected the second one. A final shortcoming of the method of suppressing nodes is that tail dependence cannot be calculated for random vectors involving latent variables. non-parametric/model-based A q q q q q q q q q q q model-based non-parametric/model-based A q q q q q q q q q q q model-based B vs A model-based B vs non-parametric 1.0 1.2 1.4 1.6 1.8 2.0 1.0 1.2 1.4 1.6 1.8 2.0

Melun

Extremal coefficients -Pairs model-based B non-parametric/model-based A q q q q q q q q q q q model-based non-parametric/model-based A q q q q q qq q q q q model-based B vs A model-based B vs non-parametric We conclude in Figure 12 with a depiction of pairwise upper tail dependence in the Seine network. Shown are the complete and reduced trees with edges weighted by the tail dependence coefficients, defined for a pair

J = {u, v} ⊂ V by 2 1 -J (1, 1; θ) (23)
in terms of the pairwise extremal coefficient in (20). In both trees, the strongest tail dependence occurs along the path from Sens to Paris.

Conclusion

We have presented a statistical model suitable for studying extremal dependence within a vector of random variables indexed by the nodes of a tree. The edges between the nodes are meant to indicate links between variables arising from a physical or conceptual network, although we do not impose any conditional independence relations. The main assumption is that, upon marginal standardization, the data-generating distribution is in the max-domain of attraction of a max-stable Hüsler-Reiss distribution whose parameter matrix possesses a certain structure induced by the tree: a free parameter is associated to each edge and each element of the Hüsler-Reiss parameter matrix only depends on the sum of the edge parameters along the path between the two corresponding nodes on the tree. We showed that the max-domain of attraction of this tree-structured Hüsler-Reiss distribution contains a specific distribution that, unlike the max-stable Hüsler-Reiss distribution or the associated multivariate Pareto distribution, satisfies the global Markov property with respect to the tree. This auxiliary model not only motivates the postulated structure, it also allowed us to find extremal dependence properties of any distribution satisfying our main assumption.

The central point and contribution of the paper is related to the identifiability of the edge parameters in case some of the variables are latent (unobservable). This situation occurs for instance in applications on river networks, when measurements on certain locations are missing. We showed that the edge parameters are uniquely identified by the distribution of the observable variables if and only if all nodes indexing latent variables are of degree at least three. Thanks to this result it is possible to quantify tail dependence even between latent variables. The characterization is due to the special structure of the variogram matrix of the Hüsler-Reiss distribution and may not be applicable to other max-stable distributions.

We fitted the model to water level data on the Seine network on a tree with seven variables, two of which were latent. As the corresponding nodes both had degree three, the six edge parameters were still identifiable and could be estimated based on data from the five observable variables. Three different estimators were proposed and implemented, based on the method of moments, on composite likelihood, and on pairwise extremal coefficients. Comparisons of non-parametric and model-based tail dependence quantities confirmed the adequacy of the fitted structured Hüsler-Reiss distribution.

For comparison we estimated a model where the two nodes with latent variables were suppressed and the edges between the affected parts of the network were redrawn in an arbitrary way. Although for the Seine data this reduction did not have a big impact on the fitted tail dependence model of the observable variables, we argued why it is still recommendable to take latent variables into account, provided there is now a sound way to do so.

An open question concerns parameter identifiability criteria in case of latent variables for Hüsler-Reiss distributions with parameter matrices structured in different ways than in this paper. Even the structure itself may be partially unknown. Another interesting direction for further research concerns extensions from the Hüsler-Reiss family to other parametric families of max-stable distributions.

A Appendix

A.1 Relation to extremal graphical models in Engelke and Hitz (2020) [START_REF] Engelke | Graphical Models for Extremes[END_REF] introduce graphical models for extremes in terms of the multivariate Pareto distribution associated to a simple max-stable distribution G. We briefly review their approach and compare it with ours in case G is a Hüsler-Reiss distribution. Let V = {1, . . . , d} and let X = (X v , v ∈ V ) be a random vector with unit-Pareto margins. The condition that

X ∈ D(G) is equivalent to lim t→∞ P(X v ≤ tz v , v ∈ V ) t = G(z), z ∈ (0, ∞) V .
By a direct calculation, it follows that

lim t→∞ P X v /t ≤ z v , v ∈ V max v∈V X v > t = ln G min(z v , 1), v ∈ V -ln G(z) ln G(1, . . . , 1) , (24)
for z ∈ (0, ∞) V , from which

(X v /t, v ∈ V ) | max v∈V X v > t d -→ Y, t → ∞ (25)
where Y = (Y v , v ∈ V ) is a random vector whose distribution function is equal to the right-hand side in (24). The law of Y is a multivariate Pareto distribution, which, upon a change in location, is a special case of the multivariate generalized Pareto distributions arising in [START_REF] Rootzén | Multivariate generalized Pareto distributions[END_REF] and Beirlant et al. (2004, Section 8.3) as limit distributions of multivariate peaks over thresholds. Assuming that Y is absolutely continuous, its support is equal to the L-shaped set {y ∈ (0, ∞) V : max v∈V Y v > 1} or a subset thereof, making conditional independence notions related to density factorizations ill-suited for Y . This is why Engelke and Hitz (2020) study conditional independence relations for the random vector Y u defined in distribution as Y | Y u > 1 for u ∈ V . According to Engelke and Hitz (2020, Definition 2), the law of Y is defined to be an extremal graphical model with respect to some graph G if for all u ∈ V , the law of Y u satisfies the global Markov property with respect to G. Note that Y itself is not required to satisfy the said Markov property.

The multivariate Pareto distribution derived through (24) from the Hüsler-Reiss distribution G = H Λ is referred to in [START_REF] Engelke | Graphical Models for Extremes[END_REF] as the Hüsler-Reiss Pareto distribution. In their article, the term Hüsler-Reiss graphical models is then used for Hüsler-Reiss Pareto distributions that are extremal graphical models.

To show the relation with our approach, note that (25) implies that, for all u ∈ V , we have

(X v /t, v ∈ V ) | X u > t d -→ Y u , t → ∞.
Recall the tail tree (Ξ u,v , v ∈ V \ u) in ( 7) and put Ξ u,u = 0. Equations ( 10) and ( 13) in combination with Theorem 2 in Segers (2020) and the continuous mapping theorem imply that

(X v /t, v ∈ V ) | X u > t d -→ (ζΞ u,v , v ∈ V ), t → ∞,
where ζ is a unit-Pareto random variable, independent of the log-normal random vector (Ξ u,v , v ∈ V ).

Comparing the two previous limit relations, we find that Y u is equal in distribution to (ζΞ u,v , v ∈ V ).

The representation ln Ξ u,v = e∈(u v) ln M e as path sums starting from u over independent Gaussian increments ln M e along the edges implies that the Gaussian vector (ln Ξ u,v , v ∈ V ) satisfies the global Markov property with respect to T . Since ζ is independent of (Ξ u,v , v ∈ U ) this Markov property then also holds for (ζΞ u,v , v ∈ V ) and thus also for Y u . But this means exactly that the multivariate Pareto distribution associated to the Hüsler-Reiss distribution with parameter matrix Λ(θ) in ( 5) is an extremal graphical model with respect to T . By way of comparison, the random vector Z * constructed via properties (Z1)-(Z2) in Section 2.3 is not max-stable nor multivariate Pareto, but it satisfies the global Markov property with respect to T and it belongs to D(H Λ(θ) ), motivating the chosen structure of Λ(θ) in (5). In Section 2.2, our assumption on ξ after transformation to X with unit-Pareto margins is that X ∈ D(H Λ(θ) ). In this sense, we require that the extremal dependence of ξ is like the one of the graphical model Z * . Our approach is thus different from the one in [START_REF] Engelke | Graphical Models for Extremes[END_REF], who postulate a new definition of extremal graphical models for multivariate Pareto vectors, but without regard for the max-domain of attraction of the corresponding max-stable distributions. Still, for graphical models with respect to trees, both methods arrive at the same structure for the Hüsler-Reiss parameter matrix Λ(θ).

A.2 Comparison with the Lee-Joe structured Hüsler-Reiss model [START_REF] Lee | Multivariate extreme value copulas with factor and tree dependence structures[END_REF] already proposed a way to bring structure to the parameter matrix Λ

= (λ 2 ij ) d i,j=1
of a d-variate max-stable Hüsler-Reiss distribution. Recall that [START_REF] Hüsler | Maxima of normal random vectors: between independence and complete dependence[END_REF] studied the asymptotic distribution of the component-wise maxima of a triangular array of row-wise independent and identically distributed Gaussian random vectors, the n-th row having correlation matrix ρ(n).

Assuming (1 -ρ ij (n)) ln(n) → λ 2 ij
as n → ∞, they found the limit to be the distribution bearing their name. Motivated by this property, [START_REF] Lee | Multivariate extreme value copulas with factor and tree dependence structures[END_REF]

propose to set λ 2 ij = (1 -ρ ij )ν where ρ = (ρ ij ) d i,j=1
is a structured correlation matrix and ν > 0 is a free parameter. They then introduce the Hüsler-Reiss distributions that result from imposing on ρ the structure of a factor model or the one of a p-truncated vine. If p = 1, the latter becomes a Markov tree and we can compare their model with ours. In their case, a free correlation parameter α e ∈ (-1, 1) is associated to each edge e ∈ E of the tree on V = {1, . . . , d}. The correlation matrix ρ of the resulting Gaussian graphical model is

ρ ij = e∈(i j) α e , i, j ∈ V.
The Lee-Joe model for the structured Hüsler-Reiss matrix Λ LJ derived from ρ is therefore

λ 2 ij = (1 -ρ ij )ν =   1 - e∈(i j) α e   ν, i, j ∈ V. ( 26
)
The model in ( 26) is to be compared with the one in our Eq. ( 5). The former has (d -1) + 1 = d free parameters, (α e , e ∈ E) and ν, whereas the latter has only d -1 free parameters (θ e , e ∈ E). In Eq. ( 5), the Hüsler-Reiss parameters satisfy

λ 2 ij = e∈(i j) λ 2 e , i, j ∈ V,
where we write λ e = λ ab for e = (a, b) ∈ E. In contrast, the Lee-Joe parameter matrix in Eq. ( 26) only satisfies this additivity relation asymptotically as ν → ∞. For instance, on a tree with d = 3 nodes and edges (1, 2) and (2, 3), i.e., a chain, their and our models satisfy respectively

λ 2 13 = λ 2 12 + λ 2 23 -ν -1 λ 2 12 λ 2 23
for λ 2 ij as in Eq. ( 26), λ 2 13 = λ 2 12 + λ 2 23 for λ 2 ij as in Eq. ( 5).

Since the Lee-Joe parameter ν > 0 takes the role of ln(n) in the Hüsler-Reiss limit relation, we can think of it as being large. In this interpretation, our parametrization becomes a limiting case of the one of [START_REF] Lee | Multivariate extreme value copulas with factor and tree dependence structures[END_REF].

Whereas the Lee-Joe parametrization is motivated from the limit result in [START_REF] Hüsler | Maxima of normal random vectors: between independence and complete dependence[END_REF] for row-wise maxima of Gaussian triangular arrays, ours is motivated as the max-stable attractor of certain regularly varying Markov trees as in [START_REF] Segers | One-versus multi-component regular variation and extremes of Markov trees[END_REF], the vector Z * in Section 2.3 serving as example. A possible advantage of our structure is that the resulting multivariate Pareto vector falls into the framework of conditional independence for such vectors is an extremal graphical model as in Engelke and Hitz (2020, Definition 2), as discussed in Appendix A.1. In general, this is not true for the multivariate Pareto vector induced by the Lee-Joe structure. For the trivariate tree in the preceding paragraph, for instance, the criterion in Proposition 3 in [START_REF] Engelke | Graphical Models for Extremes[END_REF] is easily checked to be verified for our matrix Λ but not for the one of [START_REF] Lee | Multivariate extreme value copulas with factor and tree dependence structures[END_REF].

For the Seine data, we compare the fitted Lee-Joe tail dependence model with ours. In order to avoid possible identifiability issues for the Lee-Joe parameters, we suppress the nodes with latent variables and use the right-hand tree in Figure 10 for the d = 5 observable ones, corresponding to the five locations in the dataset. The estimation method of [START_REF] Lee | Multivariate extreme value copulas with factor and tree dependence structures[END_REF] is based on pairwise copulas and annual maxima via composite likelihood. For year y and for variable j ∈ {1, . . . , d}, let m y,j be the maximum of all observations for that variable and that year, insofar available. These maxima are reported in Table 2 and their availability depends on the variable, i.e., on the location. For Melun there are only 15 such annual maxima in comparison to 33 for Nemours. For each variable j, transform these maxima to uniform margins ûy,j using the empirical cumulative distribution function based on all available maxima for that variable. Note that for this transformation, [START_REF] Lee | Multivariate extreme value copulas with factor and tree dependence structures[END_REF] rely on estimated generalized extreme value distributions instead. For variables i, j ∈ {1, . . . , d}, let Y ij be the set of years y for which annual maxima are available for both variables. For the pair (Paris, Meaux) this is the period 1999-2019 while for the pair (Paris, Nemours) this is 1990-2019. Let c(u, v; λ 2 ) denote the bivariate Hüsler-Reiss copula density with parameter λ 2 . Following [START_REF] Lee | Multivariate extreme value copulas with factor and tree dependence structures[END_REF], we estimate the free parameters in Eq. ( 26) by maximizing a composite likelihood: letting λ 2 ij (α, ν) denote the right-hand side in ( 26), the parameter estimates are

(α, ν) = arg max α∈(-1,1) d-1 ,ν∈(0,∞) d i,j=1 y∈Yij ln c ûi,y , ûj,y ; λ 2 ij (α, ν) .
For the implementation, we relied on the R package CopulaModel [START_REF] Krupskii | CopulaModel: CopulaModel: Dependence Modeling with Copulas[END_REF]. Next, we compute bivariate extremal coefficients and compare them with the non-parametric ones on the one hand and with those obtained using our own model on the other hand. The points in Figure 13 being some distance away from the diagonal, the two methods indeed seem to give somewhat different results. Moreover, there is less concordance between the non-parametric estimates and the ones from the Lee-Joe model than between the non-parametric ones and those resulting from our model: compare the red crosses in Figure 13 with those in the left-hand plots in Figure 11.

A.3 Proof of Proposition 2.1

We show that the stdf l of Z * is equal to l U in (4) with U = V and Λ = Λ(θ). Since the margins of Z * are unit-Fréchet, they are tail equivalent to the unit-Pareto distribution, so that standardization to the unit-Pareto distribution is unnecessary. By the inclusion-exclusion principle,

l(x 1 , . . . , x d ) = lim t→∞ t P(Z * 1 > t/x 1 or . . . or Z * d > t/x d ) = d i=1 (-1) i-1 W ⊆V |W |=i lim t→∞ t P(Z * v > t/x v , v ∈ W ) (27)
for x ∈ (0, ∞) d . For any non-empty W ⊆ V and any u ∈ W , it holds by (7) in combination with Theorem 2 in Segers (2020) that lim

t→∞ t P(Z * v > t/x v , v ∈ W ) = lim t→∞ t 1 t/x u P Z * u t/x u Z * v Z * u > x u x v , v ∈ W \ u Z * u > t x u = x u P(ζΞ uv > x u /x v , v ∈ W \ u),
with ζ a unit-Pareto variable independent of (Ξ uv , v ∈ V \ u). Using the fact that 1/ζ is a uniform variable on [0, 1] and setting Ξ uu = 1 we have

x u P(ζΞ uv > x u /x v , v ∈ W \ u)
1.0 1.2 1.4 1.6 1.8 2.0 1.0 1.2 1.4 1.6 1.8 2.0

Extremal coefficients -Pairs

Lee-Joe model non-parametric/model-based q q q q q q q q q q q Lee-Joe vs non-param Lee-Joe vs ECE Lee-Joe vs MM+CL Figure 13: Bivariate extremal coefficients comparison: using the modelling and estimation method of Lee and Joe (2018, Section 4-5) and those proposed in this paper.

= x u P 1/ζ < min{(x v /x u )Ξ uv , v ∈ W \ u} = x u E[min{1, (x v /x u )Ξ uv , v ∈ W \ u}] = E[min{x v Ξ uv , v ∈ W }] = xu 0 P x v Ξ uv > y, v ∈ W \ u dy = ∞ -ln xu P ln Ξ uv > (-ln x v ) -z, v ∈ W \ u exp(-z) dz
upon a change of variable y = exp(-z). Since (ln Ξ uv , v ∈ V \ u) is multivariate normal with mean vector µ V,u (θ) and covariance matrix Σ V,u (θ), we obtain from (27) that the stdf of Z * is equal toln H Λ(θ) (1/x 1 , . . . , 1/x d ), with H Λ the cumulative distribution function in Eqs. (3.5)-(3.6) in [START_REF] Hüsler | Maxima of normal random vectors: between independence and complete dependence[END_REF], but with unit-Fréchet margins rather than Gumbel ones. By Remark 2.5 in [START_REF] Nikoloulopoulos | Extreme value properties of multivariate t copulas[END_REF], this stdf is equal to the one given in (4), as required.

A.4 Choice of node neighborhoods and parameter identifiability

The MM estimator in (18) and the CL estimator in Section 4.2 involve the choice of subsets W u ⊆ U for u ∈ U . These sets or neighborhoods need to be chosen in such a way that the parameter vector θ is still identifiable from the collection of covariance matrices Σ Wu,u (θ) for u ∈ U and thus from the path sums p a,b for a, b ∈ W u and u ∈ U . Here we illustrate this issue with an example.

Consider the following structure on five nodes where all variables are observable except for the one on node 2:

ξ 1 ξ 2 ξ 3 ξ 4 ξ 5 θ 12 θ 23 θ 34 θ 25
Clearly, the parameter vector θ = (θ 12 , θ 23 , θ 34 , θ 25 ) is identifiable from the distribution of the observable variables because the criterion of Proposition 3.1 is satisfied: the only node whose variable is latent has degree three. First, consider the following subsets W u for u ∈ {1, 3, 4, 5}:

W 1 = {1, 5}, W 3 = {3, 4}, W 4 = {3, 4}, W 5 = {1, 5}.
The four 1 × 1 covariance matrices Σ Wu,u (θ) that correspond to these subsets are We are not able to identify the parameter θ because the set of path sums {p 15 , p 34 } is too small: we have only two equations andfor four unknowns. Second, consider instead the following node sets

Σ W1,1 (θ) = θ 2 12 + θ 2 25 = p 15 , Σ W4
  =   p 13 0 p 23 0 p 34 0 p 23 0 p 35   , Σ W5,5 (θ) = θ 2 12 + θ 2 25 = p 15 .
Clearly, the four edge parameters are identifiable from these covariance matrices.

A.5 Finite-sample performance of the estimators

We assess the performance of the three estimators introduced in Section 4 by numerical experiments involving Monte Carlo simulations. Let ξ = (ξ v , v ∈ V ) be a random vector with continuous joint probability density function and satisfying the global Markov property, (6), with respect to the graph in Figure 14. Let f u (x u ) for any u ∈ V be the marginal density function of the variable ξ u and let x j → f j|v (x j | x v ) be the conditional density function of ξ j given ξ v = x v . For any u ∈ V the joint density function of ξ is

f (x) = f u (x u ) (v,j)∈Eu f j|v (x j | x v ), (28)
with E u ⊆ E the set of edges directed away from u, i.e., (v, j) ∈ E u if and only if v = u or v separates u and j. The joint density f is determined by d -1 bivariate densities f vj . It would seem that the joint density f depends on u, but this is not so, as can be confirmed by writing out the bivariate conditional densities. We make two parametric choices: the univariate margins f u are unit Fréchet densities, f j (x j ) = exp(-1/x j )/x 2 j for x j ∈ (0, ∞), and the bivariate margins for each pair of variables on adjacent vertices j, v are Hüsler-Reiss distributions with parameter θ jv . Hence, ξ corresponds to the vector Z * in Section 2.3.

To generate an observation from the left hand-side of (28) above we use the right hand-side of that equation, proceeding iteratively, walking along paths starting from u using the conditional densities. An observation of ξ j given ξ v = x v is generated via the inverse function of the cdf x j → F j|v (x j | x v ), the conditional cdf of ξ j given ξ v = x v . To do so, the equation F j|v (x j | x v ) -p = 0 is solved numerically as a function in x j for fixed p ∈ (0, 1). The choice of the Hüsler-Reiss bivariate distribution gives the following expression for

F j|v (x j | x v ): Φ θ jv 2 + 1 θ jv ln x j x v • exp - 1 x v Φ θ jv 2 + 1 θ jv ln x j x v -1 - 1 x j Φ θ jv 2 + 1 θ jv ln x v x j .
After generating all the variables (ξ v ) v∈V in this way, independent standard normal noise ε ∼ N d (0, I d ) is added. Although the distribution of ξ = ξ + ε is not necessarily a graphical model with respect to the graph in Figure 14, it is still in the max-domain of attraction of a Hüsler-Reiss distribution with parametric matrix as in (5). Hence the vector ξ is still in the class of models under consideration in Section 2.2. The data on nodes 1 and 3 are discarded and not used in the estimation so as to mimic a model with two latent variables, ξ 1 and ξ 3 ; according to Proposition 3.1, the six dependence parameters are still identifiable. In this way, we generate 200 samples of size n = 1000. The estimators are computed with threshold tuning parameter k ∈ {25, 50, 100, 150, 200, 300}. The bias, standard deviation and root mean squared errors of the three estimators are shown in Figure 15 and Figure 16 for the six parameters. The MME and CLE are computed with the sets W u being W 2 = {2, 4, 5, 6, 7}, W 4 = W 5 = {2, 4, 5}, and W 6 = W 7 = {2, 6, 7}. As is to be expected, the absolute value of the bias is increasing with k, while the standard deviation is decreasing and the mean squared error has a U -shape and eventually increases with k. The MME and CLE have very similar properties. For larger values of the true parameter, e.g. θ 34 = 0.8 and θ 17 = 1.2, all the three estimators perform in a comparable way. The ECE tends to have larger absolute bias and standard deviation for smaller values of the true parameters.

A.6 Seine case study: data preprocessing

The data represent water level in centimeters at the five locations mentioned above and were obtained from Banque Hydro, http://www.hydro.eaufrance.fr, a web-site of the Ministry of Ecology, Energy and Sustainable Development of France providing data on hydrological indicators across the country. The dataset encompasses the period from January 1987 to April 2019 with gaps for some of the stations.

Two major floods in Paris make part of our dataset: the one in June 2016 when the water level was measured at 6.01 m and the one at the end of January 2018 with water levels slightly less than 6 m measured in Paris too. A flood of similar magnitude to the ones in 2016 and 2018 occurred in 1982. By way of comparison, the biggest reported1 flood in Paris is the one in 1910 when the level in Paris reached 8.6 m.

Table 1 shows the average and the maximum water level per station observed in the complete dataset. The maxima of Paris, Meaux, Melun and Nemours occurred either during the floods in June 2016 or the floods in January 2018, which can be seen from From Table 2 it can be observed that for many of the years the dates of maxima occurrence identify a period of several consecutive days during which the extreme event took place. For instance the maxima in 2007 occurred all in the period 4-8 March, which suggests that they make part of one extreme event. Similar examples are the periods 25-31 Dec 2010, 4-12 Feb 2013, 2-4 June 2016, etc. For most of the years this period spans between 3 and 7 days. We will take this into account when forming independent events from the dataset. In particular we choose a window of 7 consecutive calendar days within which we believe the extreme event have propagated through the seven locations. We have experimented with different length of that window, namely 3 and 5 days event period, but we have found that the estimation and analysis results are robust to that choice.

Figure 17 illustrates the water levels attained at the different locations during selected years from Table 2. The maxima of Sens, Nemours and Meaux seem to be relatively homogeneous compared to the maxima in Paris.

For all of the stations water level is recorded several times a day and we take the daily average to form a dataset of daily observations. Accounting for the gaps in the mentioned period (see We highlighted some of the years where there is a clear indication that the dates of the occurrence of the maxima at the different locations form a period of several consecutive days. The maxima attained during this period across stations can thus be considered as one extreme event. The water level in centimeters is rounded to the nearest integer. q q q q q q q q 200 300 400 500 600

Maxima attained during selected events water level (cm)

Dec '99 Jan '03 Jan '04 Mar '07 Dec '10 Feb '13 Jun '16 Jan '18

q Paris Meaux Melun Nemours Sens

Figure 17: Plot of maxima attained at each location during selected events from Table 2.

and Table 2) we end up with a dataset of 3408 daily observations in the period from 1 October 2005 to 8 April 2019. The dataset represents five time series each of length 3408. We consider two sources of non-stationarity: seasonality and serial correlation.

The serial correlation can be due to closeness in time or presence of long term time trend in the observations. We first apply a declustering procedure, similar to the one in [START_REF] Asadi | Extremes on river networks[END_REF] in order to form a collection of supposedly independent events. As a first step each of the series is transformed to ranks and the sum of the ranks is computed for every day in the dataset. The day with the maximal rank is chosen, say d * . A period of 2r + 1 consecutive days, centered around d * is considered and only the observations falling in that period are selected to form the event. Within this period the station-wise maximum is identified and the collection of the station-wise maxima forms one event. Because there is some evidence that the time an extreme event takes to propagate through the seven nodes in our model is about 3-7 days, we choose r = 3, hence we consider that one event lasts 7 days. In this way we obtain 717 observations of supposedly independent events. As it was mentioned the results are robust to the choice of r = {1, 2, 3}.

We test for seasonality and trends each of the series (each having 717 observations). The season factor is significant across all series and the time trend is marginally significant for some of the locations. We used a simple time series model to remove these non-stationarities. The model is based on season indicators and a linear time trend

X t = β 0 + β 1 1 spring t + β 2 1 summert + β 3 1 wintert + αt + t , ( 29
)
where t for t = 1, 2, . . . is a stationary mean zero process. After fitting the model in (29) to each of the five series through ordinary least squares we obtain the residuals and use those in the estimation of the extremal dependence.

A.7 ECE-based confidence interval for the dependence parameters

Let θn,k = θECE n,k denote the pairwise extremal coefficient estimator in ( 22) and let θ 0 denote the true vector of parameters. By Einmahl et al. (2018, Theorem 2) with Ω equal to the identity matrix, the ECE is asymptotically normal,

√ k( θn,k -θ 0 ) d -→ N |E| 0, M (θ 0 ) , n → ∞,
provided k = k n → ∞ such that k/n → 0 fast enough (Einmahl et al., 2012, Theorem 4.6). The asymptotic covariance matrix takes the form

M (θ 0 ) = ( L L) -1 L Σ L L( L L) -1 .
The matrices L and Σ L depend on θ 0 and are described below. For every k and every e ∈ E, an asymptotic 95% confidence interval for the edge parameter θ 0,e is given by θ 0,e ∈ θk,n;e ± 1.96 {M ( θk,n )} ee /k .

First, recall that Q ⊆ {J ⊆ U : |J| = 2} is the set of pairs on which the ECE is based and put q = |Q|. Define the R q -valued map L(θ) = l J (1, 1; θ), J ∈ Q and let L(θ) ∈ R q×|E| be its matrix of partial derivatives. For a pair J = {u, v} and an edge e = (a, b), the partial derivative of l J (1, 1; θ) with respect to θ e is given by

∂l J (1, 1; θ) ∂θ e = φ √ p uv /2 √ p uv θ e 1 {e∈(u v)} ,
where p uv is the path sum as in ( 15) and φ denotes the standard normal density function. The partial derivatives of l J (1, 1; θ) with respect to θ e for every e ∈ E form a row of the matrix L(θ). Second, Σ L (θ 0 ) is the q × q covariance matrix of the asymptotic distribution of the empirical stdf, √ k lj;n,k (1, 1) -l J (1, 1; θ 0 ) m=1,...,q d -→ N q (0, Σ L (θ 0 )), n → ∞.

The elements of the matrix Σ L (θ 0 ) are defined in terms of the stdf evaluated at different coordinates and of the partial derivatives of the stdf l(x; θ) with respect to the elements of x. For details we refer to Einmahl et al. (2018, Section 2.5). Here we note that the partial derivatives just mentioned are

Figure 1 :

 1 Figure 1: Seine network. The data is from the web-site of Copernicus Land Monitoring Service: https://land.copernicus.eu/imagery-in-situ.

Figure 3 :

 3 Figure 3: The Seine network with the tail dependence parameters associated to each edge of the tree.

Figure 4 :

 4 Figure 4: Point estimates and confidence intervals for the pairwise ECE.

Figure 6 :

 6 Figure 6: Non-parametric vs model-based extremal coefficients for pairs (left) and triples (right).

Figure 7 :

 7 Figure 7: The empirical and model-based Pickands dependence function computed using the pooled CL and MM estimates and the EC estimates. The dashed gray lines show the lower limit max(w, 1-w) of any Pickands dependence function A(w).

Figure 9 :

 9 Figure 9: Heat map of the extremal coefficients. The upper diagonal is computed using the pooled MM and CL estimates and the lower diagonal uses the EC estimates. The crosses denote flow connected nodes.

Figure 10 :

 10 Figure 10: Two different versions of the graph of the Seine network in Figure 3 if nodes with latent variables are suppressed and new edges are drawn between the remaining nodes of the affected parts of the tree.

Figure 11 :Figure 12 :

 1112 Figure 11: Comparison of extremal coefficients under model A (latent variables included) and model B (latent variables excluded). Top: combined MM and CL estimates; bottom: EC estimates. Left: pairs; right: triples.

Figure 14 :

 14 Figure 14: Tree used for the graphical model underlying the data-generating process in the simulation study in Appendix A.5. The value of the parameters are θ 12 = 0.1, θ 23 = 0.3, θ 34 = 0.8, θ 35 = 0.5, θ 16 = 0.2 and θ 17 = 1.2. Variables ξ 1 and ξ 3 are latent.

Figure 15 :

 15 Figure 15: Bias (left), standard deviation (middle) and root mean squared error (right) of the method of moment estimator (MME), composite likelihood estimator (CLE) and pairwise extremal coefficient estimator (ECE) of the parameters θ 12 (top), θ 23 (middle), and θ 34 (bottom) as a function of the threshold parameter k. Model and settings as described in Appendix A.5.

Figure 16 :

 16 Figure 16: Bias (left), standard deviation (middle) and root mean squared error (right) of the method of moment estimator (MME), composite likelihood estimator (CLE) and pairwise extremal coefficient estimator (ECE) of the parameters θ 35 (top), θ 16 (middle), and θ 17 (bottom) as a function of the threshold parameter k. Model and settings as described in Appendix A.5.

  ∂l J (x u , x v ; confidence interval for the Pickands dependence function For assessing the goodness-of-fit of the proposed model (Section 5.2), we construct non-parametric 95% confidence intervals for A(w) = l(1 -w, w) for w ∈ [0, 1]. As shown inKiriliouk et al. (2018, Section 5) this can be achieved by resampling from the empirical beta copula. For every fixed w ∈ [0, 1] we seek with a(w) and b(w) such thatP a(w) ≤ ln,k (1 -w, w) -l(1 -w, w) ≤ b(w) = 0.95 ,where ln,k is the non-parametric estimator of the stdf. For a(w) and b(w) satisfying the above expression, a point-wise confidence interval is given byA(w) ∈ ln,k (1 -w, w) -b(w), ln,k (1 -w, w) -a(w) . (30)Let (Y * v,i ) v∈U , for i = 1, . . . , n, be a random sample from the empirical beta copula drawn according to steps A1-A4 ofKiriliouk et al. (2018, Section 5). Let the function lβ n,k be the empirical beta stdf based on the original data and let the function l * n,k be the non-parametric estimate of the stdf using the bootstrap sample.We use the distribution of l * n,k -lβ n,k conditionally on the data as an estimate of the distribution of ln,k -l. Hence, we estimate a(w) and b(w) by a * (w) and b * (w) respectively defined implicitly by0.95 = P * a * (w) ≤ l * n,k (1 -w, w) -lβ n,k (1 -w, w) ≤ b * (w) = P * a + lβ n,k (1 -w, w) ≤ l * n,k (1 -w, w) ≤ b + lβ n,k (1 -w, w) .We further estimate the bootstrap distribution of l * n,k by a Monte Carlo approximation obtained by N = 1000 samples of size n from the empirical beta copula. As a consequence, the lower and upper bounds for l * n,k (1 -w, w) above are equated to the empirical 0.025-and 0.975-quantiles, respectively, yielding l *0.025 (w, 1 -w) = a * (w) + lβ n,k (w, 1 -w), l * 0.975 (w, 1 -w) = b * (w) + lβ n,k (w, 1 -w). (31)Replacing a(w) and b(w) in (30) by a * (w) and b * (w) respectively as solved from (31) yields the bootstrapped confidence interval for A(w) shown in Figure7.

Table 1 :

 1 Table 2 which displays the annual maxima at the five locations and the date of occurrence. Average and maximum water level per station in the whole dataset.

	Station	Paris	Meaux	Melun	Nemours	Sens
	Period 1 Jan 1990 -1 Nov 1999 -1 Oct 2005 -16 Jan 1987 -1 Jan 1990 -
		9 Apr 2019	9 Apr 2019	9 Apr 2019	9 Apr 2019	9 Apr 2019
	(#obs)	(10,621)	(6,287)	(4,443)	(10,154)	(9,159)
	Mean (cm)	139.11	275.85	296.61	210.07	133.46
	Max (cm)	601.95	468.70	545.48	439.03	333.80

Table 1

 1

Table 2 :

 2 Annual maxima for all stations.

	Year	Paris	Meaux	Melun	Nemours	Sens
		date	cm	date	cm	date	cm	date	cm	date	cm
	1987	n/a	n/a	n/a	n/a	n/a	n/a 15/11 221	n/a	n/a
	1988	n/a	n/a	n/a	n/a	n/a	n/a 13/02 247	n/a	n/a
	1989	n/a	n/a	n/a	n/a	n/a	n/a 04/03 213	n/a	n/a
	1990 17/02 254	n/a	n/a	n/a	n/a 03/07 217 18/02 183
	1991 10/01 339	n/a	n/a	n/a	n/a 23/04 212 04/01 175
	1992 06/12 293	n/a	n/a	n/a	n/a 15/01 218 06/12 170
	1993 28/12 377	n/a	n/a	n/a	n/a 26/09 217 26/12 184
	1994 11/01 478	n/a	n/a	n/a	n/a 19/10 253 09/01 260
	1995 30/01 500	n/a	n/a	n/a	n/a 21/03 277 28/01 259
	1996 04/12 324	n/a	n/a	n/a	n/a 03/12 219 04/12 194
	1997 28/02 313	n/a	n/a	n/a	n/a 03/07 214	n/a	n/a
	1998 02/05 358	n/a	n/a	n/a	n/a 21/12 216	n/a	n/a

According to the report of the Organisation for Economic Co-operation and Development (OECD) Preventing the flooding of the Seine in the Paris -Ile de France region -p.4.

Acknowledgements

We wish to thank the two Reviewers and the Associate Editor for a wealth of valuable remarks and suggestions, which helped us to enhance the understanding of the place of our models in the existing literature and of our contribution to it. Stefka Asenova is also grateful to David Lee for the clarifications and indications on the model in Lee and Joe (2018).