Graphical Models for Extremes

Introduction

Evaluation of the risk related to heat waves, extreme flooding, financial crises, or other rare events requires the quantification of their small occurrence probabilities. Empirical estimates are unreliable since the regions of interest are in the tail of the distribution and typically contain few or no data points. Extreme value theory provides the theoretical foundation for extrapolations to the distributional tail of a $d$-dimensional random vector $\boldsymbol{X}$. The univariate case $d=1$ is well-studied and the generalized extreme value and Pareto distributions are widely applied in areas such as hydrology (Katz et al., 2002), climate science (Min et al., 2011) and finance (McNeil et al., 2015); see also Embrechts et al. (1997) and Beirlant et al. (2004).

In the multivariate setting, $d\geq 2$, the result of the extrapolation strongly depends on the strength of extremal dependence between the components of $\boldsymbol{X}$. Most current statistical models assume multivariate regular variation for $\boldsymbol{X}$ (Resnick, 2008) since this entails mathematically elegant descriptions of the asymptotic tail distribution. Similar to the univariate setting, two different but closely related approaches exist. Max-stable distributions arise as limits of normalised maxima of independent copies of $\boldsymbol{X}$ and have been extensively studied and applied in multivariate and spatial risk problems (cf., de Haan, 1984; Gudendorf and Segers, 2010; Davison et al., 2012). On the other hand, multivariate Pareto distributions describe the random vector $\boldsymbol{X}$ conditioned on the event that at least one component exceeds a high threshold; see Rootzén and Tajvidi (2006), Rootzén et al. (2018) and Kiriliouk et al. (2018) for their construction, stability properties and statistical inference.

A drawback of the current multivariate models is their limitation to rather moderate dimensions $d$, and the construction of tractable parametric models in higher dimensions is challenging, both for max-stable and multivariate Pareto distributions. Sparse multivariate models require the notion of conditional independence (Dawid, 1979), which is not easy to define for tail distributions. In fact, Papastathopoulos and Strokorb (2016) show that if $(Z_{1},Z_{2},Z_{3})$ is a max-stable random vector with positive continuous density, then the conditional independence of $Z_{1}\perp\!\!\!\perp Z_{3}\mid Z_{2}$ already implies the independence $Z_{1}\perp\!\!\!\perp Z_{3}$; see also Dombry and Éyi-Minko (2014). Meaningful conditional independence structures can thus only be obtained for max-stable distributions with discrete spectral measure (Gissibl and Klüppelberg, 2018). Since these models do not admit densities, this excludes most of the currently used parametric families.

In this work we take the perspective of threshold exceedances and introduce a new notion of conditional independence for a multivariate Pareto distribution $\boldsymbol{Y}=(Y_{1},\dots,Y_{d})$, which we denote by $\perp_{e}$ to stress that it is designed for extremes. It is different from classical conditional independence since the support of $\boldsymbol{Y}$ is not a product space, but the homogeneity property of $\boldsymbol{Y}$ can be used to show that it is well-defined. Conditional independence is tightly linked to graphical models. For an undirected graph $\mathcal{G}=(V,E)$ with nodes $V=\{1,\dots,d\}$ and edge set $E$, we say that $\boldsymbol{Y}$ is an extremal graphical model if it satisfies the pairwise Markov property

The main advantage of conditional independence and graphical models is that they imply a simple probabilistic structure and possibly sparse patterns in multivariate random vectors (Lauritzen, 1996; Wainwright and Jordan, 2008). For extremal graphical models on decomposable graphs, we prove a Hammersley–Clifford type theorem stating that (1) is equivalent to the factorization of the density $f_{\boldsymbol{Y}}$ of $\boldsymbol{Y}$ into lower-dimensional marginals. This underlines that our notion of conditional independence is in fact natural for multivariate Pareto distributions.

Applications of this result are manifold. From a probabilistic perspective, we analyse models in the literature regarding their graphical properties in the sense of our definition (1). Extremal graphical models whose underlying graph is a tree have a particularly simple multiplicative stochastic representation in terms of extremal functions, a notion that is known from the simulation of max-stable processes (Dombry et al., 2016). In multivariate extremes, one may argue that the family of Hüsler–Reiss distributions (Hüsler and Reiss, 1989) takes a similar role as Gaussian distributions in the non-extreme world. Instead of covariance matrices, they are parameterized by a variogram matrix $\Gamma$. We show that the extremal graphical structure of a Hüsler–Reiss distribution can be identified by zero patterns on matrices derived from $\Gamma$.

Extremal graphical models enable the construction of parsimonious models for multivariate Pareto distributions $\boldsymbol{Y}$, which further enjoy the advantage of interpretability in terms of the underlying graph. Thanks to the factorization of the densities, statistical inference can be efficiently carried out on lower-dimensional marginals. For decomposable graphs with singleton separator sets, so-called block graphs, this allows the use of multivariate Pareto models in fairly high dimensions. In many cases the underlying graphical structure is unknown and has to be learned from data. We discuss how a maximum likelihood tree can be obtained using standard algorithms by Kruskal (1956) or Prim (1957), and how the best model can be selected among different extremal graphical models.

There is previous work on the construction of parsimonious extreme value models. A large body of literature studies spatial max-stable random fields (Schlather, 2002; Kabluchko et al., 2009; Opitz, 2013). Such models have small parameter dimension but they rely on strong assumptions on stationarity and cannot be applied to multivariate, non-spatial data without information on an underlying space. Other approaches include constructions through factor copulas (Lee and Joe, 2018), ensembles of trees combining bivariate copulas (Yu et al., 2017), graphical models for large censored observations (Hitz and Evans, 2016) and eigendecompositions (Cooley and Thibaud, 2018). Closely related to our concept of conditional independence is the work of Coles and Tawn (1991) and Smith et al. (1997) who propose a Markov chain model where all bivariate marginals are extreme value distributions. This can be seen as a special case of our approach when the graph has the simple structure of a chain. Similar limiting objects also arise as the tail chains in the theory of extremes of stationary Markov chains with regularly varying marginals (Smith, 1992; Basrak and Segers, 2009; Janssen and Segers, 2014). This theory has recently been extended to regularly varying Markov trees (Segers, 2019). Gissibl and Klüppelberg (2018) and Gissibl et al. (2018) study the causal structure of directed acyclic graphs for max-linear models, and they develop methods for model identification based on tail dependence coefficients. Their work is in some sense complementary to ours, since their models do not possess densities whereas we will explicitly assume the existence of densities.

To the best of our knowledge, our work is the first principled attempt to define conditional independence for general multivariate extreme value models that naturally extends to the factorization of densities, sparsity and graphical models. Section 2 introduces background on extreme value theory and graphical models needed throughout the paper. The new notion of conditional independence is defined in Section 3 and equivalent properties are derived. Section 4 contains the main probabilistic results on extremal graphical models, the representation of trees and the characterization for Hüsler–Reiss distributions. Statistical models on block graphs and their estimation, simulation and model selection are discussed in Section 5. In these graphical models the dependence is modeled directly between lower-dimensional subsets of variables, whereas the global dependence is implicitly implied by the conditional independence structure of the graph. There are many potential applications of extremal graphical models. In Section 6, we illustrate the advantages of this structured approach compared to classical extreme value models on a data set related to flooding on a river network in the upper Danube basin (cf., Asadi et al., 2015). The interpretation of the graphical structures obtained in this application is particularly interesting since there is a seemingly natural underlying tree associated to the flow-connections. Our conditional independence is formulated for multivariate Pareto distributions, but the results in this paper have implications for max-stable distributions. This point and further research directions will be addressed in the discussion in Section 7. The Appendix contains proofs and some additional results.

An implementation for R (R Core Team, 2019) is available in the package graphicalExtremes (Engelke et al., 2019). The code for the simulation study and application can be found in the supplementary material.

Background

2 Multivariate extreme value theory

The tail behavior of the random vector $\boldsymbol{X}=(X_{1},\dots,X_{d})$ can be described through two different approaches, one based on componentwise maxima and the other one on threshold exceedances. We briefly discuss both approaches and the close link between them.

The standardized vector $\boldsymbol{X}$ is said to be in the max-domain of attraction of the random vector $\boldsymbol{Z}=(Z_{1},\dots,Z_{d})$ if for any $\boldsymbol{z}=(z_{1},\dots,z_{d})$

where the exponent measure $\Lambda$ is a Radon measure on the cone $\mathcal{E}=[0,\infty)^{d}\setminus\{\boldsymbol{0}\}$, and $\Lambda\left(\boldsymbol{z}\right)$ is shorthand for $\Lambda\left(\mathcal{E}\setminus[\boldsymbol{0},\boldsymbol{z}]\right)$. If $\Lambda$ is absolutely continuous with respect to Lebesgue measure on $\mathcal{E}$, its Radon–Nikodym derivative, denoted by $\lambda$, has the following properties:

homogeneity of order $-(d+1)$, i.e., $\lambda(t\boldsymbol{y})=t^{-(d+1)}\lambda(\boldsymbol{y})$ for any $t>0$ and $\boldsymbol{y}\in\mathcal{E}$;

normalised marginals, i.e., for any $i=1,\dots,d$,

The two properties follow from the max-stability and the standard Fréchet marginals of $\boldsymbol{Z}$, respectively. For a non-empty subset $I\subset\{1,\dots,d\}$, we define the marginal of $\lambda$ by

and note that it is homogeneous of order $-(|I|+1)$. In particular, if $I=\{i\}$ for some $i=1,\dots,d$, then $\lambda_{\{i\}}(y_{i})=1/y_{i}^{2}$ as a consequence of (L1) and (L2). Conversely, any positive and continuous function $\lambda$ satisfying (L1) and (L2) defines a valid density of an exponent measure $\Lambda(\boldsymbol{z})$ by integration over $\mathcal{E}\setminus[\boldsymbol{0},\boldsymbol{z}]$, $\boldsymbol{z}\in\mathcal{E}$, that satisfies similar homogeneity and normalization properties as $\lambda$. By (4) this also defines a max-stable distribution.

Another perspective on multivariate extremes is through threshold exceedances. By Proposition 5.17 in Resnick (2008), the convergence in (3) is equivalent to

Consequently, the multivariate distribution of the threshold exceedances of $\boldsymbol{X}$ satisfies

The distribution of the limiting random vector $\boldsymbol{Y}$ is called a multivariate Pareto distribution (cf., Rootzén and Tajvidi, 2006). It is defined through the exponent measure $\Lambda$ of the max-stable distribution $\boldsymbol{Z}$, with support on the $L$-shaped space $\mathcal{L}=\{\boldsymbol{x}\in\mathcal{E}:\|\boldsymbol{x}\|_{\infty}>1\}$. We say that $\boldsymbol{Z}$ and $\boldsymbol{Y}$ are associated, since their distributions mutually determine each other.

Multivariate Pareto distributions are the only possible limits in (6) and, owing to the homogeneity of the exponent measure, they enjoy certain stability properties (cf., Rootzén et al., 2018). The exponent measure $\Lambda$, and hence the distribution of $\boldsymbol{Y}$, may place mass on some lower-dimensional faces of the space $\mathcal{E}$. For the remainder of this paper we exclude this case to avoid technical difficulties. We further assume that the distribution of $\boldsymbol{Y}$ admits a positive and continuous density $f_{\boldsymbol{Y}}$ on $\mathcal{L}$, which is

since $\Lambda(\boldsymbol{y}\wedge\boldsymbol{1})$ is always constant along at least one coordinate for $\boldsymbol{y}\in\mathcal{L}$. The density $f_{\boldsymbol{Y}}$ is thus proportional to the density $\lambda$ of the exponent measure $\Lambda$. By the homogeneity of $\lambda$, $f_{\boldsymbol{Y}}$ is also homogeneous of order $-(d+1)$. The normalization constant $\Lambda(\boldsymbol{1})\in[1,d]$ is known as the $d$-variate extremal coefficient (cf., Schlather and Tawn, 2003). The assumption of a positive and continuous density $f_{\boldsymbol{Y}}$ implies that the multivariate Pareto distributions we study are models for asymptotic extremal dependence, and all $p$-variate extremal coefficients, $1\leq p\leq d$, are strictly smaller than their upper limit $p$.

where $\Lambda_{I}$ is the exponent measure of $\boldsymbol{Z}_{I}$, and $\lambda_{I}$ is the density of $\Lambda_{I}$.

The extremal logistic distribution with parameter $\theta\in(0,1)$ induces a multivariate Pareto distribution with density

The matrix $\Sigma^{(k)}$ is strictly positive definite; see Appendix B for details. The representation of the density in (9) seems to depend on the choice of $k$, but, in fact, the value of the right-hand side of this equation is independent of $k$. The Hüsler–Reiss multivariate Pareto distribution has density $f_{\boldsymbol{Y}}(\boldsymbol{y})=\lambda(\boldsymbol{y})/\Lambda(\mathbf{1})$ and the strength of dependence between the $i$th and $j$th component is parameterized by $\Gamma_{ij}$, ranging from complete dependence for $\Gamma_{ij}=0$ and independence for $\Gamma_{ij}=+\infty$. In the bivariate case $d=2$ we have

and $\Lambda(1,1)=2\Phi\left(\sqrt{\Gamma_{12}}/2\right)$, where $\Phi$ is the standard normal distribution function. The extension of Hüsler–Reiss distributions to random fields are Brown–Resnick processes (Brown and Resnick, 1977; Kabluchko et al., 2009), which are widely used models for spatial extremes.

The above is a general construction principle, since every valid exponent measure density can be obtained in this way. The bivariate Hüsler–Reiss distribution in (11) corresponds to the case of log-normal $U^{1}_{2}$ and $U^{2}_{1}$, but many other parametric and non-parametric examples are available (e.g., Boldi and Davison, 2007; Cooley et al., 2010; Ballani and Schlather, 2011; de Carvalho and Davison, 2014).

3 Graphical models

A graph $\mathcal{G}=(V,E)$ is defined as a set of nodes $V=\{1,\dots,d\}$, also called vertices, together with a set of edges $E\subset V\times V$ of pairs of distinct nodes. The graph is called undirected if for two nodes $i,j\in V$, $(i,j)\in E$ if and only if $(j,i)\in E$. For notational convenience, for undirected graphs we sometimes represent edges as unordered pairs $\{i,j\}\in E$. When counting the number of edges, we count $\{i,j\}\in E$ such that each edge is considered only once. A subset $C\subset V$ of nodes is called complete if it is fully connected in the sense that $(i,j)\in E$ for all $i,j\in C$. We denote by $\mathcal{C}$ the set of all cliques, that is, the complete subsets that are not properly contained within any other complete subset.

and we write $\boldsymbol{X}_{A}\perp\!\!\!\perp\boldsymbol{X}_{C}\mid\boldsymbol{X}_{B}$. If $B=\emptyset$, then (13) amounts to independence of $\boldsymbol{X}_{A}$ and $\boldsymbol{X}_{C}$.

The random vector $\boldsymbol{X}$ is said to be a probabilistic graphical model on the graph $\mathcal{G}=(V,E)$ if its distribution satisfies the pairwise Markov property relative to $\mathcal{G}$, that is, $X_{i}\perp\!\!\!\perp X_{j}\mid\boldsymbol{X}_{\setminus\{i,j\}}$ for all $(i,j)\notin E$. If in addition, for any disjoint subsets $A,B,C\subset V$ such that $B$ separates $A$ from $C$ in $\mathcal{G}$, $\boldsymbol{X}_{A}\perp\!\!\!\perp\boldsymbol{X}_{C}\mid\boldsymbol{X}_{B}$, then $\boldsymbol{X}$ is said to obey the global Markov property relative to $\mathcal{G}$. Since $f_{\boldsymbol{X}}$ is positive and continuous, it follows from the Hammersley–Clifford theorem (cf., Lauritzen, 1996, Theorem 3.9) that the two Markov properties are equivalent, and they are further equivalent to the factorization of the density

for suitable functions $\psi_{C}$ on $\times_{i\in C}\mathcal{X}_{i}$. If the graph $\mathcal{G}$ is decomposable, then this factorization can be rewritten in terms of marginal densities

where $\mathcal{D}$ is a multiset containing intersections between the cliques called separator sets; see Lauritzen (1996) and Appendix A for the definition of decomposability and separator sets.

We recall that for a normal distribution $\boldsymbol{W}=(W_{i})_{i\in V}$ with invertible covariance matrix $\Sigma$, the precision matrix $\Sigma^{-1}$ contains the conditional independencies, or equivalently the graph structure, since for $i,j\in V$,

Conditional independence for threshold exceedances

The notion of conditional independence has not been exploited in extreme value theory. In fact, for max-stable distributions it only leads to trivial probabilistic structures (Papastathopoulos and Strokorb, 2016). An exception are directed acyclic graphs for max-linear models studied in Gissibl and Klüppelberg (2018) and Gissibl et al. (2018), which do however not admit densities.

We therefore approach the problem from the perspective of threshold exceedances. Since the notion of independence is only defined on product spaces, the meaning of conditional independence is not straightforward for a multivariate Pareto distribution $\boldsymbol{Y}=(Y_{i})_{i\in V}$, $V=\{1,\dots,d\}$, with support on the $L$-shaped space $\mathcal{L}=\{\boldsymbol{x}\in\mathcal{E}:\|\boldsymbol{x}\|_{\infty}>1\}$. In this section we show that there is nevertheless a natural definition of conditional independence for $\boldsymbol{Y}$. To this end, we restrict $\boldsymbol{Y}$ to product spaces. For any $k\in V$, we consider the random vector $\boldsymbol{Y}^{k}$ defined as $\boldsymbol{Y}$ conditioned on the event that $\{Y_{k}>1\}$. Clearly, $\boldsymbol{Y}^{k}$ has support on the product space $\mathcal{L}^{k}=\{\boldsymbol{x}\in\mathcal{L}:x_{k}>1\}$ (cf., Figure 1) and it admits the density

For any set $I\subset V$ with $k\in I$, the marginal $\boldsymbol{Y}_{I}^{k}$ has density

which is homogeneous of order $-(|I|+1)$ on $\mathcal{L}_{I}^{k}=\{\boldsymbol{x}_{I}\in\mathcal{L}_{I}:x_{k}>1\}$; see (5). This is however not the case if $k\notin I$, since integration over $\boldsymbol{y}_{\setminus I}$ then includes $y_{k}$ whose domain is $(1,\infty)$ in $\mathcal{L}^{k}$, and thus in general $f^{k}_{I}(\boldsymbol{y}_{I})\neq\lambda_{I}(\boldsymbol{y}_{I})$, $\boldsymbol{y}_{I}\in[0,\infty)^{|I|}.$

Suppose that $\boldsymbol{Y}$ is multivariate Pareto and admits a positive and continuous density $f_{\boldsymbol{Y}}$ on $\mathcal{L}$, and let $A,B,C\subset V$ be non-empty disjoint subsets whose union is $V=\{1,\ldots,d\}$. We say that $\boldsymbol{Y}_{A}$ is conditionally independent of $\boldsymbol{Y}_{C}$ given $\boldsymbol{Y}_{B}$ if

In this case we write $\boldsymbol{Y}_{A}\perp_{e}\boldsymbol{Y}_{C}\mid\boldsymbol{Y}_{B}$.

In fact, this definition can be shown to be equivalent to a slightly weaker condition, and to a factorization of the exponent measure density $\lambda$.

Let $f_{\boldsymbol{Y}}$ and the sets $A,B,C$ be as in the above definition, then $\boldsymbol{Y}_{A}\perp_{e}\boldsymbol{Y}_{C}\mid\boldsymbol{Y}_{B}$ is equivalent to any of the following two conditions.

The density of the exponent measure factorizes as

A natural question is whether one can extend the definition of $\boldsymbol{Y}_{A}\perp_{e}\boldsymbol{Y}_{C}\mid\boldsymbol{Y}_{B}$ to the case where $B=\emptyset$, meaning that $\boldsymbol{Y}_{A}$ and $\boldsymbol{Y}_{C}$ are independent on $\mathcal{L}$. In terms of the original definition, that would mean that for any $k\in V$, $f^{k}(\boldsymbol{y})=f^{k}_{A}(\boldsymbol{y}_{A})f^{k}_{C}(\boldsymbol{y}_{C})$ for all $\boldsymbol{y}\in\mathcal{L}^{k}$. Without losing generality, suppose $k\in A$, then $f^{k}_{C}(\boldsymbol{y}_{C})=\lambda(\boldsymbol{y}_{A},\boldsymbol{y}_{C})/\lambda_{A}(\boldsymbol{y}_{A})$ for any $\boldsymbol{y}_{A}\in\mathcal{L}^{k}_{A}$ and $\boldsymbol{y}_{C}\in[0,\infty)^{|C|}$. Therefore $f^{k}_{C}$ would be homogeneous of order $-|C|$ and thus not integrable on $[0,\infty)^{|C|}$, a contradiction. In the next section we show that this property implies that all graphical models defined in terms of the conditional independence $\perp_{e}$ must be connected.

Graphical models for threshold exceedances

The notion of conditional independence allows us to define graphical models for threshold exceedances. As before, let $f_{\boldsymbol{Y}}$ be the positive and continuous density on $\mathcal{L}$ of a multivariate Pareto distribution $\boldsymbol{Y}$, proportional to the density $\lambda$ of the exponent measure $\Lambda$, and homogeneous of order $-(d+1)$. Let $\mathcal{G}=(V,E)$ be an undirected graph with nodes $V=\{1,\dots,d\}$ and edge set $E$. Similarly to the classical probabilistic graphical models described in Section 2.3, we say that $\boldsymbol{Y}$ satisfies the pairwise Markov property on $\mathcal{L}$ relative to $\mathcal{G}$ if

that is, $Y_{i}$ and $Y_{j}$ are conditionally independent in the sense of Definition 1 given all other nodes whenever there is no edge between $i$ and $j$ in $\mathcal{G}$. By definition, this is equivalent to saying that $\boldsymbol{Y}^{k}$ satisfies the usual pairwise Markov property on $\mathcal{L}^{k}$ relative to $\mathcal{G}$ for all $k\in V$. The global Markov property for $\boldsymbol{Y}$ is defined similarly.

Let $\mathcal{G}=(V,E)$ be an undirected graph. If the multivariate Pareto distribution $\boldsymbol{Y}$ with positive and continuous density $f_{\boldsymbol{Y}}$ satisfies the pairwise Markov property (20) relative to $\mathcal{G}$, we call the distribution of $\boldsymbol{Y}$ an extremal graphical model with respect to $\mathcal{G}$.

For a decomposable graph $\mathcal{G}$ we obtain a factorization of the density $f_{\boldsymbol{Y}}$ similar to the classical Hammersley–Clifford theorem, showing that the Definition 1 of conditional independence is natural for multivariate Pareto distributions. Let $\mathcal{C}$ and $\mathcal{D}$ be the sequences of cliques and separators of $\mathcal{G}$, respectively, satisfying the running intersection property (44) in Appendix A.

Let $\mathcal{G}=(V,E)$ be a decomposable graph and suppose that $\boldsymbol{Y}$ has a multivariate Pareto distribution with positive and continuous density $f_{\boldsymbol{Y}}$ on $\mathcal{L}$. Denote the corresponding exponent measure and its density by $\Lambda$ and $\lambda$, respectively. Then the following are equivalent.

The distribution of ${\boldsymbol{Y}}$ satisfies the pairwise Markov property relative to $\mathcal{G}.$

The distribution of ${\boldsymbol{Y}}$ satisfies the global Markov property relative to $\mathcal{G}.$

The density $f_{\boldsymbol{Y}}$ factorizes according to $\mathcal{G}$, that is,

where the marginals $\lambda_{I}$ are positive, continuous and homogeneous of order $-(|I|+1)$ for any $I\subset V$.

In all cases, the graph $\mathcal{G}$ is necessarily connected.

The above theorem shows that only connected extremal graphical models can arise. This is related to the assumption of multivariate regular variation in (3) that implies asymptotic dependence between all components. Loosely speaking, unconnected components would correspond to asymptotically independent components.

If the graph $\mathcal{G}$ in the above theorem is non-decomposable, it is expected that the density $f_{\boldsymbol{Y}}$ still factorizes into factors on the cliques of the graph. These factors can however no longer be identified with marginal densities, and it is an open problem whether they can be chosen to be homogeneous.

Since $\mathcal{L}$ is not a product space, unlike in the classical Hammersley–Clifford theorem for decomposable graphs in (15), the factors in the factorization of the density $f_{\boldsymbol{Y}}$ in (21) are not the marginals $f_{I}$ but the marginals of the exponent measure density $\lambda_{I}$. It holds however that $f_{I}(\boldsymbol{y}_{I})=\lambda_{I}(\boldsymbol{y}_{I})/\Lambda_{I}(\boldsymbol{1})$ for all $\boldsymbol{y}_{I}\in\mathcal{L}_{I}\subset\{\boldsymbol{x}_{I}:\boldsymbol{x}\in\mathcal{L}\}$.

As a first application, the above theorem allows us to formally analyse the conditional independencies and graphical structures of models in the multivariate extreme value literature.

One of the simplest examples of a graph is a chain, that is,

Coles and Tawn (1991) proposed a model that factorizes with respect to this chain where all bivariate marginals are logistic (cf., Example 1). This was extended to general bivariate marginals in Smith et al. (1997). More generally, in the study of extremes of stationary Markov chains the limiting objects are so-called tail chains. The latter induce multivariate Pareto distributions that can readily be seen to factorize with respect to a chain; see Smith (1992) Basrak and Segers (2009) and Janssen and Segers (2014).

It turns out that many of the multivariate models in the literature do not have any conditional independencies, that is, their underlying graph is fully connected. For instance, this holds for the logistic multivariate Pareto distribution in Example 1, the Dirichlet mixture model in Boldi and Davison (2007), and the pairwise beta distribution in Cooley et al. (2010).

Similar to Gaussian distributions, an appealing property of the Hüsler–Reiss model is its stability under taking marginals. Indeed, for any $I\subset V$ and $k\in I$ the marginal density of the exponent measure is

with the notation of Example 2, where $\Sigma^{(k)}_{I}$ is the matrix in (10) induced by the submatrix $\Gamma_{I}$. Thus, $f_{I}(\boldsymbol{y}_{I})=\lambda_{I}(\boldsymbol{y}_{I})/\Lambda_{I}(\mathbf{1})$ is the density of the $|I|$-dimensional Hüsler–Reiss Pareto distribution with parameter matrix $\Gamma_{I}$.

By Theorem 1, the density of a Hüsler–Reiss distribution that satisfies the pairwise Markov property relative to some decomposable graph $\mathcal{G}$ factorizes into lower-dimensional Hüsler–Reiss distributions. The explicit formula is given in Corollary 2 in Appendix C.

Theorem 1 can also be seen as a construction principle to build new classes of extreme value distributions in a modular way by combining lower-dimensional marginals. The following corollary shows how a multivariate Pareto distributions can be defined that factorizes according to a desired underlying graphical structure. This is particularly useful in high-dimensional problems to ensure model sparsity.

Let $\mathcal{G}$ be a decomposable and connected graph and suppose that $\{\lambda_{I}:I\in\mathcal{C}\cup\mathcal{D}\}$ is a set of valid, positive and continuous exponent measure densities in the sense of (L1) and (L2) in Section 2.2. For $D\subset C$, $D\in\mathcal{D},$ $C\in\mathcal{C}$, assume that they satisfy the consistency constraint

The density of a valid $d$-dimensional exponent measure $\Lambda$ is then given by

and the function $f_{\boldsymbol{Y}}(\boldsymbol{y})=\lambda(\boldsymbol{y})/\Lambda(\mathbf{1})$, $\boldsymbol{y}\in\mathcal{L}$, is the density of a multivariate Pareto distribution satisfying the pairwise Markov property relative to $\mathcal{G}$.

A tree is a special case of a decomposable graphical model that is connected and has no cycles. All cliques are then of size two and the separators contain only one node. Let $\mathcal{T}=(V,E)$ be an undirected tree with nodes $V=\{1,\dots,d\}$ and edge set $E\subset V\times V$. Suppose that $\boldsymbol{Y}=(Y_{i})_{i\in V}$ follows a multivariate Pareto distribution on $\mathcal{L}$ with density $f_{\boldsymbol{Y}}$ that is an extremal graphical model with respect to the tree $\mathcal{T}$. Theorem 1 yields the factorization

where $\lambda_{ij}=\lambda_{\{i,j\}}$ are the bivariate marginals of the exponent measure density $\lambda$ corresponding to $\boldsymbol{Y}$. The formula (23) allows the extension of the modelling approach by Smith et al. (1997) described in Example 5 from time series to general tree structures. Such tree models are able to represent more complex dependencies and, moreover, are suitable beyond temporal data for multivariate or spatial applications.

Thanks to the relatively simple structure of a tree, more explicit results can be derived than for general graphical models. To this end, we define a new, directed tree $\mathcal{T}^{k}=(V,E^{k})$ rooted at an arbitrary but fixed node $k\in V$. The edge set $E^{k}$ consist of all edges $e\in E$ of the tree $\mathcal{T}$ pointing away from node $k$; see Figure 2 for an example with $k=2$. For the resulting directed tree we define a set $(U_{e})_{e\in E^{k}}$ of independent random variables, where for $e=(i,j)$, the distribution of $U_{e}=U^{i}_{j}$ is the extremal function of $\lambda_{ij}$ at coordinate $j$, relative to coordinate $i$; see (12) in Example 3 for the definition of extremal functions. The following stochastic representation of the random vectors $\boldsymbol{Y}^{k}$, $k\in V$, provides a better understanding of the stochastic structure of multivariate Pareto distributions factorizing on a tree, and it is the main ingredient for efficient simulation in Section 5.4.

Let $\boldsymbol{Y}$ be a multivariate Pareto distribution with positive and continuous density on $\mathcal{L}$ that factorizes with respect to the tree $\mathcal{T}$. With the notation above, and for a standard Pareto distribution $P$, we have the joint stochastic representation for $\boldsymbol{Y}^{k}$ on $\mathcal{L}^{k}$

where $\operatorname{ph}(ki)$ denotes the set of edges on the unique path from node $k$ to node $i$ on the tree $\mathcal{T}^{k}$.

The same object as in (24) has been obtained in Segers (2019) as the limit of regularly varying random vectors that satisfy a Markov condition on a tree. In analogy to the tail chains in Example 5, they term it a tail tree.

Suppose all bivariate marginals $\lambda_{ij}$ for $\{i,j\}\in E$ of a tree Pareto model are of logistic type with parameter $\theta_{ij}\in(0,1)$ as defined in Example 1. This tree logistic model is a generalization of the chain logistic model considered in Coles and Tawn (1991). In this symmetric case, the extremal functions $U^{i}_{j}$ and $U^{j}_{i}$ have the same distribution with stochastic representation $F/G$, where $F$ follows a Fréchet$(1/\theta,c_{\theta})$ distribution with scale parameter $c_{\theta}=\Gamma(1-\theta)^{-1}$ and $(G/c_{\theta})^{-1/\theta}$ follows a Gamma$(1-\theta,1)$ distribution, where we abbreviated $\theta=\theta_{ij}$ and $\Gamma$ is the Gamma function.

Similarly we can define a Hüsler–Reiss tree model, or use asymmetric models for $\lambda_{ij}$ such as the Dirichlet model in Boldi and Davison (2007) for some or all of the edges $\{i,j\}\in E$. In asymmetric models, the extremal functions $U^{i}_{j}$ and $U^{j}_{i}$ have in general different distributions. We refer to Section 4 in Dombry et al. (2016) for explicit formulas for extremal function distributions of commonly used model classes.

2 Hüsler–Reiss graphical models

In many ways, the class of Hüsler–Reiss distributions introduced in Example 2 can be seen as the natural analog of Gaussian distributions in the world of asymptotically dependent extremes. They are parameterized by the variogram of Gaussian distributions, and their statistical inference (Wadsworth and Tawn, 2014; Engelke et al., 2015) and exact simulation (Dombry et al., 2016) involves tools that are closely related to the corresponding methods for normal models.

Despite the similarities to Gaussian distributions, there are subtle but important differences that render analysis and statistical inference of Hüsler–Reiss distributions more difficult. In order to characterise conditional independence and graphical structures in these models, we first recall some results related to the original construction. The max-stable Hüsler–Reiss distribution has a stochastic representation as componentwise maxima

Importantly, this implies that the representation in (25) is not unique since any centred, possibly degenerate normal distribution $\boldsymbol{W}$ with variogram matrix $\Gamma$ leads to the same max-stable Hüsler–Reiss distribution. Let

be the set of all covariance matrices that correspond to the same variogram matrix $\Gamma$; see Appendix B. The Hüsler–Reiss Pareto distribution $\boldsymbol{Y}$ associated with $\boldsymbol{Z}$ is defined by its density in Example 2, which is also parameterized by $\Gamma$. We recall that for a normal distribution $\boldsymbol{W}$ with invertible covariance matrix $\Sigma$, the precision matrix $\Sigma^{-1}$ contains the conditional independencies; see Example 4. A first, naive guess would be that the graph structure of $\boldsymbol{W}$ used in the construction of $\boldsymbol{Z}$ directly translates into the extremal graph structure of the Hüsler–Reiss Pareto distribution $\boldsymbol{Y}$. This is however not the case.

We consider three examples for $\boldsymbol{W}$ in the representation (25) with $d=4$.

Let $W_{i}$, $i=1,\dots,4$, be independent standard normal distributions, then $\Sigma^{-1}=\operatorname{diag}(1,\dots,1)$ and $\Gamma_{ij}=2$ if $i\neq j$ and zero otherwise. The graph underlying the distribution of $\boldsymbol{W}$ is the graph with four unconnected nodes. The graph of the corresponding Hüsler–Reiss Pareto distribution $\boldsymbol{Y}$ turns out to be the fully connected graph on the left-hand side of Figure 3.

Consider the centred normal distribution $\boldsymbol{W}$ with precision matrix and variogram matrix

respectively. The Gaussian graphical model is the graph in the centre of Figure 3 with an additional edge between the nodes $2$ and $3$. On the contrary, the corresponding Hüsler–Reiss model factorizes according to the graph in the centre of Figure 3.

Consider the centred normal distribution $\boldsymbol{W}$ with precision matrix and variogram matrix

respectively. It can be checked that both the Gaussian and the Hüsler–Reiss graphical model are as in the right-hand side of Figure 3. Also note that this graph is not decomposable.

We denote the precision matrix of $\Sigma^{(k)}$ by $\Theta^{(k)}=(\Sigma^{(k)})^{-1}$. For notational convenience, the indices of the matrices $\Sigma^{(k)}$ and $\Theta^{(k)}$ range in $\{1,\dots,d\}\setminus\{k\}$ instead of $\{1,\dots,d-1\}$.

For $k,k^{\prime}\in V$, $k\neq k^{\prime}$, the precision matrices $\Theta^{(k)}$ and $\Theta^{(k^{\prime})}$ satisfy for $i,j\in V\setminus\{k^{\prime}\}$

The above lemma is of independent interest since it explains the link between the precision matrices $\Theta^{(k)}$ for different different $k\in V$. The proof uses blockwise inversion of the precision matrices. This result is the crucial ingredient to characterise conditional independence in Hüsler–Reiss models.

For a Hüsler–Reiss Pareto distribution $\boldsymbol{Y}$ with parameter matrix $\Gamma$, it holds for $i,j\in V$ with $i\neq j$, and for any $k\in V$, that

For any $k\in V$, the single matrix $\Theta^{(k)}$ contains all information on conditional independence of $\boldsymbol{Y}$. Conditional independence concerning the $k$th component is encoded in the row and column sums of $\Theta^{(k)}$, and it might sometimes be easier to switch to another representation $\Theta^{(k^{\prime})}$, $k^{\prime}\neq k$, where it simply figures as a zero entry. In Example 9 we can now easily determine the graphical model $\mathcal{G}=(V,E)$ for each of the three Hüsler–Reiss Pareto distributions. For a given $\Sigma$ we first compute the matrix $\Gamma$ as in (26), then transform it by (10) to obtain $\Sigma^{(k)}$ for any $k\in V$ and use Proposition 3 to decide whether $(i,j)\in E$ for all $i,j\in V$. These transformations are implemented in our R-package graphicalExtremes (Engelke et al., 2019).

In this section, we have so far not required that the underlying graph $\mathcal{G}$ is decomposable. If this is the case then, as shown in Example 7, Theorem 1 implies that the density of the Hüsler–Reiss graphical model factorizes into lower-dimensional Hüsler–Reiss densities; see Corollary 2 in Appendix C.

Statistical inference for block graphs

The notion of conditional independence and graphical models for multivariate Pareto distributions allows the construction of new statistical models with two major advantages. First, sparsity can be imposed on the model, which is a crucial ingredient for tractable and parsimonious models in higher dimensions. Second, under certain graphical structures, the model parameters can be estimated separately on lower-dimensional subsets of the data.

We consider here, and throughout the rest of the paper, decomposable and connected graphs $\mathcal{G}=(V,E)$ with clique set $\mathcal{C}$ and separator set $\mathcal{D}$, where all separators in $\mathcal{D}$ are single nodes. Such graph structures with singleton separator sets are known as block graphs (cf., Harary, 1963) and have already been seen to have appealing properties for discrete data (Loh and Wainwright, 2013). In our case, they are a convenient way of restricting the model complexity in order to obtain a tractable class of extremal graphical models. In fact, Corollary 1 provides a simple construction principle for multivariate Pareto distributions that factorize with respect to the block graph $\mathcal{G}$.

For each clique $C\in\mathcal{C}$, choose possibly different parametric families of valid exponent measure densities $\{\lambda_{C}(\cdot;\theta_{C}):\theta_{C}\in\Omega_{C}\}$ for suitable parameter spaces $\Omega_{C}$. If $\mathcal{G}$ is a tree $\mathcal{T}$, then this reduces to choosing $d-1$ bivariate exponent measure densities $\lambda_{ij}$, for each $\{i,j\}\in E$; see Example 3 for a general representation of such densities.

Since all separator sets consist of a single node, the consistency constraint (22) is trivially fulfilled as a consequence of (L1) and (L2) in Section 2.2 and the fact that $\lambda_{D}(y_{D})=y_{D}^{-2}$ for all $D\in\mathcal{D}$.

For any fixed combination of parameters $\theta=(\theta_{C})_{C\in\mathcal{C}}\in\Omega=\times_{C\in\mathcal{C}}\Omega_{C}$, the product of the normalised lower-dimensional exponent measure densities,

defines a valid $d$-variate Pareto distribution factorizing according to the graph $\mathcal{G}$, which is a member of the parametric family parameterized by $\theta\in\Omega$. For a tree $\mathcal{T}$, this reduces to the density in (23).

Concrete examples for this construction are tree logistic or tree Hüsler–Reiss models as described in Example 8, where all cliques have the same type of distributions. The above construction is much more flexible, as it allows us to use different distribution families for the different cliques. Moreover, some, or even all of the cliques may be modeled by non-parametric methods; see Lafferty et al. (2012) for non-parametric tree models in the non-extreme case. In this direction, there is a line of research on kernel-based estimation of exponent measure densities (cf., de Carvalho and Davison, 2014; Marcon et al., 2017; Kiriliouk et al., 2018) that could be used as clique models. We will not follow this approach here.

In the graphical models above, the dependence inside each clique is modeled directly, whereas dependence between components from different cliques is implicitly implied by the conditional independence structure of the graph. Even if all cliques are modeled with the same type of parametric family, the joint distribution (30) is typically not of this distribution type. For a tree logistic distribution, for instance, this can easily be seen by comparing its density (23) with that of $d$-variate logistic distribution in Example 1. The latter only has one parameter governing the whole $d$-dimensional dependence structure, whereas the tree has $d-1$ logistic parameters $\{\theta_{ij};\{i,j\}\in E\}$ and thus much higher flexibility.

An important exception is the family of Hüsler–Reiss distributions, which is stable under taking marginal distributions; see Example 7. The following proposition shows that for a given graphical structure as above, if all cliques have Hüsler–Reiss distributions, then so has the full $d$-dimensional model. This is the converse of Corollary 2 in Appendix C.

Let $\mathcal{G}=(V,E)$ be a block graph as above, and suppose that on each clique $C\in\mathcal{C}$, $\boldsymbol{Y}$ has a $|C|$-variate Hüsler–Reiss distribution with exponent measure density $\lambda_{C}(\cdot;\Gamma^{(C)})$ parameterized by a $|C|\times|C|$-dimensional variogram matrix $\Gamma^{(C)}$. Then there exists a unique solution to the problem:

with the notation from Proposition 3. The corresponding $d$-variate Hüsler–Reiss distribution factorizes according to the graph $\mathcal{G}$ into the lower-dimensional Hüsler–Reiss densities on the cliques.

This is a matrix completion problem for variograms similar to what Dempster (1972) introduced for covariance matrices. In our case, the graph is decomposable and the above result relates to the marginal problem studied in Kellerer (1964) and Dawid and Lauritzen (1993). For Hüsler–Reiss marginals on block graphs we even see that the implied $d$-dimensional distribution is again Hüsler–Reiss. We give a direct, constructive proof in Appendix F. This provides a method to construct high-dimensional Hüsler–Reiss distributions out of many low-dimensional ones. The full $d$-variate Hüsler–Reiss model without any conditional independencies has $d(d-1)/2$ parameters. A Hüsler–Reiss distribution as in Proposition 4 that factorizes on a block graph with click set $\mathcal{C}$ has only

parameters, which can be much smaller than $d(d-1)/2$.

2 Estimation

Extremal graphical models can be used to build parsimonious statistical models for the tail of a multivariate random vector. In this section we discuss how the model parameters can be estimated efficiently by considering each clique distribution separately.

Let $\boldsymbol{X}=(X_{j})_{j\in V}$, $V=\{1,\dots,d\}$, be a random vector in the max-domain of attraction of the max-stable random vector $\boldsymbol{Z}$ as in (3), with marginal distribution $X_{j}$ in the max-domain of attraction of a generalized extreme value distribution with shape parameter $\xi_{j}$, $j\in V$. Equivalently, there exist a sequence of high thresholds $\boldsymbol{t}_{u}=(t_{u1},\dots,t_{ud})$ with $t_{uj}$ tending to the upper endpoint of $X_{j}$ as $u\to\infty$, and positive normalizing functions $\sigma_{u}=(\sigma_{u1},\dots,\sigma_{ud})$, such that the distribution of exceedances converges weakly

where $\boldsymbol{Y}$ is the multivariate Pareto distribution associated with $\boldsymbol{Z}$. We assume $\boldsymbol{Y}$ to be in the model class of the previous section with density (30), and for now we suppose that the underlying graph $\mathcal{G}=(V,E)$ is known and fixed. The conditional density of $\boldsymbol{X}-\boldsymbol{t}_{u}$ given that $\|\boldsymbol{X}/\boldsymbol{t}_{u}\|_{\infty}>1$ is then approximated by

This density can be used to estimate jointly the marginal parameters $(\sigma_{uj},\xi_{j})$, $j\in V$, and the dependence parameter vector $\theta=(\theta_{C})_{C\in\mathcal{C}}$ of $f_{\boldsymbol{Y}}$.

In the sequel we concentrate on estimation of the dependence, and we therefore assume that the marginal parameters are known or have been estimated separately. As described in Section 2.2, we can then normalise $\boldsymbol{X}$ to standard Pareto marginals, in which case $\xi_{j}=1$, $t_{uj}=u$ and $\sigma_{uj}=u$ for all $j\in V$. We recover the standardized setting of (6) considered throughout the paper, where $\boldsymbol{X}/u$ given that $\|\boldsymbol{X}\|_{\infty}>u$ converges to $\boldsymbol{Y}$, whose likelihood is proportional as a function of $\theta$ to

Direct maximization of the likelihood with contributions (34) for each data point is tedious since the normalizing constant $Z_{\theta}$ contains all parameters and does not factorize. Fortunately the class of block graphs has the property that we can estimate the parameters $\theta_{C}$ of each $\lambda_{C}$ separately, without having to enforce the consistency constraints at the separator sets. In fact, we use the following observation. If $\boldsymbol{X}$ is in the domain of attraction of the family of multivariate Pareto distributions $\{f_{\boldsymbol{Y}}(\cdot;\theta):\theta\in\Omega\}$, then for a fixed clique $C\in\mathcal{C}$, the subvector $\boldsymbol{X}_{C}$ is in the domain of attraction of $\{f_{C}(\cdot;\theta_{C}):\theta_{C}\in\Omega_{C}\}$, and the distribution of the normalised exceedance $\boldsymbol{X}_{C}/u\mid\|\boldsymbol{X}_{C}\|_{\infty}>u$ is approximated for large $u$ by $\boldsymbol{Y}_{C}$ with density

see (7) in Section 2.2. We can therefore obtain an estimate of $\theta_{C}$ based only on data of the components in $C$, whose dimension is typically much smaller than the dimension $d$ of the full graph. Estimating the cliques separately might in principle result in a loss of estimation efficiency compared to using the joint likelihood (34). The normalizing constant $Z_{\theta}$ does however not contain much information on the parameter $\theta$ and the maximum likelihood estimate using $f_{\boldsymbol{Y}}(\boldsymbol{y};\theta)$ is generally very close to the estimate obtained by maximizing separate likelihoods based on (35). We discuss this point in the simulation study in Section 5.5.

In practice, some components of $\boldsymbol{X}$ might not have converged to the limiting distribution $\boldsymbol{Y}$. In order to avoid biased estimates of the dependence parameters $\theta_{C}$, it has become a standard approach to apply censoring to the data; see Ledford and Tawn (1997), Smith et al. (1997). For a data point $\boldsymbol{X}_{C}$ with $\|\boldsymbol{X}_{C}\|_{\infty}>u$ for a high threshold $u>0$, define $J$ to be the set of indices $j\in C$ such that $Y_{j}<1$, i.e., $X_{j}<u$. For this data point we use the censored likelihood contribution

which uses for all $j\in J$ only the information that this component of $\boldsymbol{Y}_{C}$ is smaller than $1$, but not its exact value. For explicit forms of the censored likelihoods for many parametric models see Dombry et al. (2017) and Kiriliouk et al. (2018).

For $n$ independent data $\boldsymbol{y}^{(h)}\in\mathcal{L}$, $h=1,\dots,n$, of $\boldsymbol{X}/u\mid\|\boldsymbol{X}\|_{\infty}>u$, for each clique $C$ we define $\widehat{\theta}_{C}$ as the maximizer of the censored log-likelihood

where $\mathcal{L}_{C}=\{\boldsymbol{y}\in\mathcal{L}:\exists j\in C\text{ s.t. }y_{j}>1\}$, and each $\boldsymbol{y}^{(h)}_{C}$ has its own censoring set $J^{(h)}\subset C$.

Maximum likelihood estimation is only one possibility to infer the parameters $\theta_{C}$ based on exceedances of $\boldsymbol{X}_{C}$ and the limiting distribution (35). Alternative methods use $M$-estimators (Einmahl et al., 2012; Einmahl et al., 2016) or proper scoring rules (de Fondeville and Davison, 2018).

3 Model selection

Up to now we have assumed that a graphical structure $\mathcal{G}$ was a priori given and we analysed models that factorize with respect to this structure. In many applications the underlying graph structure is unknown and should be learned in a data-driven way. Theorem 1 implies that all extremal graphical structures are connected, and a simple and flexible class of connected graphs are trees; see Section 4.1. It is thus natural to first build a suitable tree as a baseline model, and then extend the tree by adding additional edges in order to obtain more complex graphs.

Since trees are a special case of general graphical models, there are specific methods to learn these simpler structures. The notion of a minimum spanning tree is crucial (Kruskal, 1956). Let $\mathcal{G}_{0}=(V,E_{0})$ be the fully connected graph on $V=\{1,\dots,d\}$ with edge set $E_{0}=\{(i,j):i,j\in V\}$. Suppose that a positive weight $w_{ij}>0$ is attached to each edge $(i,j)\in E_{0}$ of $\mathcal{G}_{0}$. This number can be seen as the length of the edge $(i,j)$ or the distance between nodes $i$ and $j$, and it is assumed that $w_{ij}=w_{ji}$ and $w_{ii}=0$, $i,j\in V$. The minimum spanning tree is the tree $\mathcal{T}_{\operatorname{mst}}=(V,E_{\operatorname{mst}})$ with $E_{\operatorname{mst}}\subset E_{0}$, that minimizes the sum of weights on that tree, i.e.,

If all edges of $\mathcal{G}_{0}$ have distinct lengths, then $\mathcal{T}_{\operatorname{mst}}$ is unique. This minimization problem can be solved efficiently by the greedy algorithms proposed in Kruskal (1956) or Prim (1957).

The weights $w_{ij}$ determine the tree structure and should be chosen carefully. A common approach in graphical modelling is to search the conditional independence structure that maximizes the likelihood, (cf., Cowell et al., 2006, Chapter 11). Such a tree is also called a Chow–Liu tree (Chow and Liu, 1968). We fix a parametric family of bivariate Pareto distributions that is used for all pairs of nodes $\{f(\cdot;\theta_{ij}):\theta_{ij}\in\Omega\}$. For $n$ independent data $\boldsymbol{y}^{(h)}$, $h=1,\dots,n$, the maximal log-likelihood of a fixed tree within this parametric class is essentially the sum over the maximized clique log-likelihoods in (37) over all edges of this tree. In order to find the tree that maximizes the log-likelihood over all trees and all distributions in this parametric family, we therefore find the minimum spanning tree in (38) with weights

where we include the censored marginal densities $y_{i}^{-2}$ and $y_{j}^{-2}$ in (30) for the clique $\{i,j\}$, since now the edges are no longer fixed but parameters of the optimization. The resulting tree ${\mathcal{T}}_{\operatorname{mst}}$ is the baseline model for the data. If the model fit is not satisfactory, it is possible to extend this tree to graphs with more complex structures by adding additional edges. The family of Hüsler–Reiss distributions is particularly appealing since the bivariate marginals remain in the same class. We illustrate this model extension through a greedy forward selection in Section 5.5.

The different multivariate Pareto models can then be compared by the Akaike information criterion (Kiriliouk et al., 2018),

where $p$ is the number of parameters in the respective model, and the second term is twice the negative log-likelihood based on the censored version of (34), evaluated at the optimized parameters of each clique.

4 Exact simulation

Exact simulation of a max-stable random vector $\boldsymbol{Z}$ relies on the notion of extremal functions (Dombry and Éyi-Minko, 2013). The extremal function of $\boldsymbol{Z}$, or of its associated multivariate Pareto distribution $\boldsymbol{Y}$, relative to coordinate $k\in V$ is the $d$-dimensional random vector $\boldsymbol{U}^{k}$ with $U^{k}_{k}=1$ such that the exponent measure density of $\boldsymbol{Z}$ can be written as

The distributions of the extremal functions $\boldsymbol{U}^{k}$, $k\in V$, for most commonly used models have explicit forms and are derived in Section 4 of Dombry et al. (2016). Theorem 2 in the same paper relates the distribution of the so-called spectral measure to these extremal functions. Together with the following representation of $\boldsymbol{Y}$, this enables simulation of multivariate Pareto distributions by rejection sampling. Recall that for any $k\in V$, the random vector $\boldsymbol{Y}^{k}$ is defined as $\boldsymbol{Y}$ conditioned on the event that $\{Y_{k}>1\}$.

The distribution of the extremal function $\boldsymbol{U}^{k}$ of $\boldsymbol{Y}$ relative to coordinate $k\in V$ is given by the distribution of $\boldsymbol{Y}^{k}/Y^{k}_{k}$. Independently, let $P$ be a standard Pareto random variable and $T$ uniformly distributed on $\{1,\dots,d\}$. We then have for any Borel set $A\subset\mathcal{L}$

The above representation yields a simple algorithm for exact simulation of $\boldsymbol{Y}$; see also de Fondeville and Davison (2018).

1. Simulate a standard Pareto random variable $P$. 2. Simulate $T$ uniformly on $\{1,\dots,d\}$ and sample a realization of the extremal function $\boldsymbol{U}^{T}$ relative to coordinate $T$. 3. If $\max\{P\|\boldsymbol{U}^{T}\|_{\infty}/\|\boldsymbol{U}^{T}\|_{1}\}>1$, 3. return $\boldsymbol{Y}=P\boldsymbol{U}^{T}/\|\boldsymbol{U}^{T}\|_{1}$ as realization of the multivariate Pareto distribution. 4. Else, 4. reject the simulation and go back to step 1.

The complexity of this simulation algorithm as a function of the dimension $d$ of the vector $\boldsymbol{Y}$ is driven by the number of times one has to sample from one of the extremal functions $\boldsymbol{U}^{1},\dots,\boldsymbol{U}^{d}$, since simulation of the variables $P$ and $T$ requires much less computational effort. Let $C_{\boldsymbol{Y}}(d)$ denote the number of extremal functions that have to be simulated in the above algorithm. The random variable $C_{\boldsymbol{Y}}(d)$ follows a geometric distribution and from (50) in the proof of Lemma 2 its expectation is

The complexity measures $C_{\boldsymbol{Y}}(d)$ and $C_{\boldsymbol{Z}}(d)$ only consider the number of extremal functions required for one exact simulation of $\boldsymbol{Y}$ and $\boldsymbol{Z}$, respectively. The computational effort of sampling $\boldsymbol{U}^{k}$ can however be significantly lower if $\boldsymbol{Y}$ has a sparse structure. If $\boldsymbol{Y}$ factorizes according to a graph, then, by the Definition 1 of conditional independence, the $\boldsymbol{Y}^{1},\dots,\boldsymbol{Y}^{d}$ inherit the sparsity of this graph structure. This is particularly important in the case of trees and for Hüsler–Reiss distributions, as shown in the examples below. It is important to note that more efficient simulation of the extremal functions speeds up exact simulation of the multivariate Pareto distribution $\boldsymbol{Y}$, but also of the max-stable distribution $\boldsymbol{Z}$.

Suppose that $\boldsymbol{Y}$ factorizes according to a tree $\mathcal{T}=(V,E)$. It follows from Proposition 2 and Lemma 2 that the extremal function $\boldsymbol{U}^{k}$ relative to coordinate $k\in V$ is

For exact simulation of $\boldsymbol{Y}$ it therefore suffices to simulate the univariate random variables $U_{e}$. This is feasible even in very large dimensions.

If $\boldsymbol{Y}$ has a Hüsler–Reiss distribution that factorizes on the graph $\mathcal{G}=(V,E)$, then it follows from (28) that the extremal function $\boldsymbol{U}^{k}$ relative to coordinate $k\in V$ is

The exact simulation algorithms for both multivariate Pareto and max-stable distributions are implemented in our R-package graphicalExtremes (Engelke et al., 2019).

5 Simulation study

We assess the efficiency of parameter estimation and model selection in the framework of graphical models for extremes described in the previous sections. We fix a dimension $d$ of variables or nodes $V=\{1,\dots,d\}$ and a block graph $\mathcal{G}=(V,E)$ as in Section 5.1. In this study we simulate samples directly from the limiting distribution $\boldsymbol{Y}$ using the exact Algorithm 1, but we use the censored estimation since this is common practice in applications.

We first choose $d=5$ and let $\mathcal{G}$ be the undirected version of the tree in Figure 2. We simulate $n\in\{100,200\}$ samples $\boldsymbol{y}^{(1)},\dots,\boldsymbol{y}^{(n)}$ of a Hüsler–Reiss distribution with parameter matrix $\Gamma$ that factorizes according to $\mathcal{G}$. The entries of $\Gamma$ need to be specified only on the submatrices $\Gamma^{(C)}$ for all cliques $C\in\mathcal{C}$ of $\mathcal{G}$, since the solution to the matrix completion problem in Proposition 4 then yields the unique variogram matrix $\Gamma$. In this simulation we set

where we only specified the four parameters $\Gamma_{ij}$ for $(i,j)\in E$, $i<j$, to the values in bold, and the rest of the matrix is implied by the graph structure.

In this dimension we can still maximize the censored version of the joint likelihood (34) to obtain an estimate $\widehat{\Gamma}_{ij}^{\text{joint}}$, $\{i,j\}\in E$, of the parameters corresponding to the four edges of the tree. We also obtain estimates $\widehat{\Gamma}_{ij}$, $\{i,j\}\in E$, of the parameters of each clique separately by maximizing the censored clique likelihood (37). In both cases, the four estimated parameters yield estimates $\widehat{\Gamma}^{\text{joint}}$ and $\widehat{\Gamma}$ of the whole variogram matrix $\Gamma$ through the graph structure. We repeat the simulation and estimation $200$ times and compare the efficiency of both approaches in Figure 4, displaying only the four free parameters that have actually been estimated.

The difference in efficiency between the joint and clique likelihoods seems to be small or even negligible. This is due to two reasons. For non-censored points the two likelihoods only differ by the normalizing constant $Z_{\theta}$. Since this constant only measures the global strength of dependence and does not depend on the data, it seems not very sensitive to changes in the parameter $\theta$. The second difference between the two approaches is that they use slightly different data. Consider a clique $C\in\mathcal{C}$ and the corresponding model parameter $\theta_{C}$. The joint likelihood uses all data $\boldsymbol{Y}$ in the space $\mathcal{L}=\{\boldsymbol{y}\in\mathcal{E}:\exists j\in V\text{ s.t. }y_{j}>1\}$, but censors all components with $y_{j}\leq 1$. On the other hand, the clique likelihood uses the marginals $\boldsymbol{Y}_{C}$ of all data $\boldsymbol{Y}$ in $\mathcal{L}_{C}=\{\boldsymbol{y}\in\mathcal{L}:\exists j\in C\text{ s.t. }y_{j}>1\}$. Consequently, the additional data used in the joint likelihood is in $\mathcal{L}\setminus\mathcal{L}_{C}=\{\boldsymbol{y}\in\mathcal{L}:y_{j}\leq 1\text{ for all }j\in C\}$. But the contribution to the joint likelihood of data in this set with regard to the parameter $\theta_{C}$ is completely censored and does therefore not add significant additional information. These two considerations underline that estimating the parameters for each clique separately does not result in significant efficiency losses. This is one of the main advantages of graphical models, namely that the distribution is defined locally by the cliques and extends globally by the conditional independence structure. In terms of computational aspects, the joint likelihood becomes infeasible even in moderate dimensions, whereas the clique likelihood is applicable in high dimensions as long as the cliques have small enough sizes. Moreover, the computations for different cliques can be easily parallelized.

For the second experiment we take $d=16$ and let $\mathcal{G}$ be the graph on the left-hand side of Figure 5, which is not a tree. We simulate $n=100$ samples of a Hüsler–Reiss distribution with parameter matrix $\Gamma$ that factorizes according to $\mathcal{G}$. The parameters of the $p=18$ edges are independently sampled from a uniform distribution on $(0.5,1)$, under the constraint that $\Gamma$ is conditionally negative definite on cliques with three nodes. We illustrate how we can choose the best graphical model, where we restrict to block graphs as in Section 5.1 with cliques of sizes two and three. We first construct the minimum spanning tree as described in Section 5.3 within the class of Hüsler–Reiss distributions. The estimated edge set of this tree is denoted by $E_{1}$. The $15$ parameter estimates $\widehat{\Gamma}_{ij}$, $\{i,j\}\in E_{1}$ obtained by fitting the clique likelihoods of each clique of the tree yield a unique estimate $\widehat{\Gamma}$ of the $d\times d$-dimensional variogram matrix; see Proposition 4. This tree model does not contain all edges of the true underlying graph. We therefore perform a greedy forward selection in order to add additional edges and improve the model. In each step, we define an enlarged edge set $E_{m+1}=E_{m}\cup\{i,j\}$, $m=1,2,\dots$, restricting to those edges $\{i,j\}$, $i,j\in V$, that still yield a block graph with cliques of maximal size three. We continue this process until no more edge can be added in this way. For the same parameter matrix $\Gamma$, we repeat the simulation and model selection 100 times. The right-hand side of Figure 5 shows the graph with the selected edges, where the line width of each edge indicates the number of times it has been selected among the first $18$ edges. It can be seen that the graph structure is generally very well identified. For each model and each repetition we also compute the resulting $\operatorname{AIC}$ according to (40). The proportion of times that the model with $\{15,\dots,20\}$ edges has the smallest $\operatorname{AIC}$ are $\{0.01,0.11,0.23,0.39,0.23,0.03\}$. Even though the $\operatorname{AIC}$ is a criterion built for model estimation and not for identification (cf., Arlot and Celisse, 2010), it seems to be well suited to select the correct degree of sparsity for this extremal graphical model.

Application

We illustrate the applicability of extremal graphical models at the example of river discharges in the upper Danube basin, a region that is prone to serious flooding. The data are provided by the Bavarian Environmental Agency (http://www.gkd.bayern.de) and we use $d=31$ gauging stations with $50$ years of common daily data from 1960–2009. The tree induced by the physical flow-connections at these stations is shown on the left-hand side of Figure 6, where the path $10\to 9\to\dots\to 1$ is on the Danube and the other branches are tributaries. The spatial extremal dependence structure of this data set has been studied in Asadi et al. (2015) and we follow their preprocessing steps to make the results comparable. Out of all daily data only the three months June, July and August are considered since the most severe floods occur in this period and are caused by heavy summer rain (Böhm and Wetzel, 2006). The $50\times 92=4600$ observations in these months are declustered in time in order to remove temporal dependence and to match slightly shifted peak flows at different locations. We refer to Asadi et al. (2015) for more details on the data, the declustering method and exploratory analysis concerning stationarity and asymptotic dependence; see also Keef et al. (2009, 2013) for other approaches to flood risk assessment.

The max-stable Brown–Resnick model in Asadi et al. (2015) corresponds to a parametric family of Hüsler–Reiss Pareto distributions $\{f_{\boldsymbol{Y}}(\cdot;\theta):\theta\in\Omega\}$ at the $31$ gauging stations. The dependence model is tailor-made for this particular application to river extremes and uses several covariates such as distance on the river network, catchment sizes and altitudes. In terms of our new notion of extremal graphical models it is readily checked using the results of Proposition 3 that for any parameter value $\theta\in\Omega$ their model does not exhibit conditional independencies.

We propose a different Hüsler–Reiss model that factorizes according to a sparse graph and does not require any domain knowledge or additional covariates. In fact, we propose a sequence of models

where $\theta^{(l)}=(\theta^{(l)}_{C})_{C\in\mathcal{C}^{(l)}}$, and $\mathcal{C}^{(l)}$ is the set of all cliques of the $l$th extremal graphical model $\mathcal{G}^{(l)}$ according to which the model family $M^{(l)}$ factorizes. As simplest model we take $\mathcal{G}^{(1)}$ to be the minimum spanning tree within the family of Hüsler–Reiss models as described in Section 5.3. Similarly as in the simulation study in Section 5.5, we obtain $\mathcal{G}^{(2)},\dots,\mathcal{G}^{(L)}$ by successively adding edges to the tree $\mathcal{G}^{(1)}$ in a greedy way while restricting the model class to block graphs with cliques of size at most three. The estimated tree $\mathcal{G}^{(1)}$ is shown on the left-hand side of Figure 9 in Appendix D. It is very similar to the tree in Figure 6 that corresponds to the tree induced by the flow-connections of the river network. There are however differences, and it is important to note that the flow-connection tree is not necessarily the optimal tree structure in terms of extreme river discharges. Appendix D also contains a sensitivity analysis of the tree structure for different thresholds $u$, and a comparison to a Gaussian tree model fitted to non-extremal data.

Figure 7 shows the $\operatorname{AIC}$ values for the different models $M^{(1)},\dots,M^{(L)}$. The forward selection is a greedy approach and it does not guarantee to find the optimal graph. We therefore also initialize the forward selection with the simplest model $\mathcal{G}^{(1)}$ being the flow-connection tree on the left-hand side of Figure 6. This tree must have a larger $\operatorname{AIC}$ than the minimum spanning tree, but interestingly, the left panel of Figure 7 shows that by adding additional edges the optimal $\operatorname{AIC}$ is better than the previous optimal $\operatorname{AIC}$. In this particular case, we thus choose the graph initiated with the flow-connection tree with $9$ additional edges. In general, a tree structure appears to be too simple for this application. The reason is that only part of the extremal dependence of discharges at different locations can be explained by flow-connections. Additional dependence may arise even between flow-unconnected locations due to proximity of their catchments that are affected by the same spatial precipitation events. Asadi et al. (2015) model this explicitly through a variogram with two parts, one for the dependence on the river network and one for the spatial, meteorological dependence. The $9$ additional edges of the graphical model on the right-hand side of Figure 6, which minimizes the $\operatorname{AIC}$, partly improve the model in terms of this spatial dependence between flow-unconnected stations, but also strengthen it between some flow-connected locations. This best graphical model has $39$ edges and an $\operatorname{AIC}$ of $5269.43$. It significantly outperforms the simpler tree models with $30$ edges and the spatial model of Asadi et al. (2015), which has only six parameters but an $\operatorname{AIC}$ of $5291.34$, which is indicated by the dashed orange line in the left panel of Figure 7.

A popular summary statistic for extremal dependence between $Y_{i}$ and $Y_{j}$, $i,j\in V$, is the tail correlation (cf., Coles et al., 1999), which can be expressed as $\chi_{ij}=2-\Lambda_{ij}(1,1)$. The centre and right panels of Figure 7 compare empirical estimates of these statistics for all pairs of stations with those implied by the fitted models. In terms of this bivariate summary, both models seem to fit the data well, even though the graphical model seems to be slightly less biased than the spatial model. There are also versions of $\chi$ that assess how a model captures the higher-order extremal dependence structure. In Figure 11 in Appendix E we compare the trivariate empirical $\chi$ coefficients with those implied from the fitted spatial and graphical model. Both models fit well the trivariate dependence, again with a slightly lower bias of the graphical model.

In this application we have only considered block graphs, which are particularly convenient in terms of statistical inference as seen in the previous sections. In general it should be assessed whether this sparse model class is justified for the data. In our case, the bivariate and trivariate $\chi$ coefficients indicate that block graphs are flexible enough to capture the extremal dependence structure of the river data. This is further supported by the fact that the AIC curve in Figure 7 attains its minimum even before the maximal number of edges is added in this model class. It is an important question for future research how extremal graphical models with more complicated structures can be estimated.

Discussion

The conditional independence relation $\perp_{e}$ introduced in this paper is natural for a multivariate Pareto distribution $\boldsymbol{Y}$ as it explains the factorization of its density $\boldsymbol{f}_{\boldsymbol{Y}}$ into lower-dimensional marginals (cf., Theorem 1). This establishes a link of extreme value statistics to the broad field of graphical models, and it opens the door to define sparsity and to perform structure learning for tail distributions. In this work we have studied the probabilistic structure and statistical inference for some important models, with the main purpose of modelling the extremal dependence structure. Many subsequent research directions are possible. Directed acyclic graphs as in Gissibl and Klüppelberg (2018) for max-linear models may be formulated in our setting and would yield different factorizations than for undirected graphs, and this would form the basis to extend work on causal inference for extremes (Naveau et al., 2018; Mhalla et al., 2019; Gnecco et al., 2019) to continuous extreme value distributions. The models in this paper are well-suited for asymptotic dependence. Another line of research focuses on multivariate tail models under asymptotic independence (Ledford and Tawn, 1997; Heffernan and Tawn, 2004; Wadsworth et al., 2017). Conditional independence and graphical models have not been studied in this framework, except for the special case of Markov chains (Kulik and Soulier, 2015; Papastathopoulos et al., 2017).

Conditional independence for $\boldsymbol{Y}$ does not carry over to factorization of the density of the associated max-stable distribution $\boldsymbol{Z}$. By Proposition 1, the conditional independence relation $\perp_{e}$ does however imply the factorization of the exponent measure density $\lambda$ of $\boldsymbol{Z}$, which is the key object in simulation (Dombry et al., 2016) and full likelihood estimation (Thibaud et al., 2016; Dombry et al., 2017; Huser et al., 2019) of max-stable processes. Thus, sparsity in our notion for multivariate Pareto distributions also facilitates inferential tasks for max-stable distributions, a fact that has been briefly discussed for simulation in Section 5.4 but deserves further investigation.

Acknowledgments

We thank Robin J. Evans and Nicola Gnecco for helpful discussions. We are grateful to the editorial team and the referees for knowledgeable comments that improved the paper. Financial support by the Swiss National Science Foundation (S. Engelke) and by the Berrow Foundation (A. S. Hitz) is gratefully acknowledged. The paper was completed while S. Engelke was a visitor at the Department of Statistical Sciences, University of Toronto.

Appendix

Let $\mathcal{G}=(V,E)$ be an undirected graph with node set $V=\{1,\dots,d\}$ and edge set $E\subset V\times V$; see Section 2.3. We define the notion decompositions and decomposability for the graph $\mathcal{G}$ (cf., Lauritzen, 1996, Definition 2.1).

A triplet $(A,B,C)$ of disjoints subsets of $V$ is said to form a decomposition of $\mathcal{G}$ into the components $\mathcal{G}_{A\cup B}$ and $\mathcal{G}_{B\cup C}$ if $V=A\cup B\cup C$ and

$B$ separates $A$ from $C$ (i.e., every path from $A$ to $C$ intersects $B$);

The decomposition is called proper if $A$ and $C$ are both non-empty. A graph $\mathcal{G}$ is decomposable if it is complete or if there exists a proper decomposition $(A,B,C)$ into decomposable subgraphs $\mathcal{G}_{A\cup B}$ and $\mathcal{G}_{B\cup C}.$ Decomposable graphs are also known as triangulated or chordal graphs.

For instance, $(\{1,2,3,4,5\},\{4,5\},\{4,5,6\})$ is a proper decomposition of the decomposable graph in Figure 8.

For a connected, decomposable graph $\mathcal{G}$, we can order the set of the cliques $\mathcal{C}=\{C_{1},\dots,C_{m}\}$ such that for all $i=2,\dots,m$,

a condition called the running intersection property; cf., Lauritzen (1996, Chapter 2) and Green and Thomas (2013). The sets $D_{i}$, $i=2,\dots,m$, are called separators of the graph, and both $\mathcal{C}$ and the collection of separators $\mathcal{D}=\{D_{2},\dots,D_{m}\}$ are uniquely determined up to different orderings. The separators may not all be distinct, and we say that $\mathcal{D}$ is a multiset. A possible enumeration of cliques and separators for the graph in Figure 8 that satisfies the running intersection property is

From (44) we note that the clique $C_{m}$ intersects the other cliques only in $D_{m}$. Consider the connected, decomposable subgraph $\mathcal{G}_{m-1}$ of $\mathcal{G}$ with node set $V_{m-1}=V\setminus(C_{m}\setminus D_{m})$ and corresponding induced edge set. The property (44) then holds for $\mathcal{G}_{m-1}$, which has one clique less. Continuing this process, we note that each $C_{j}$ intersects the subgraph $\mathcal{G}_{j}$ only in $D_{j}$, $j=2,\dots,m$, and $\mathcal{G}_{1}$ with nodes $V_{1}=C_{1}$ is complete.

B Link between variogram and covariance matrices

For any $k\in V$, there is a bijection $\varphi_{k}:\mathcal{D}_{d}\to\mathcal{P}_{d-1}^{k}$ given by

using the fact that $\Gamma$ is symmetric and $\Gamma_{ii}=0$ for all $i\in V$. The assertion then follows; see also the proof of Lemma 3.2.1 in Berg et al. (1984). ∎

C Hüsler–Reiss densities on decomposable graphs

Let $\mathcal{G}=(V,E)$ be a decomposable and connected graph, and suppose that $\boldsymbol{Y}$ is a Hüsler–Reiss Pareto distribution that satisfies the pairwise Markov property