# Graphical models for multivariate extremes

Sebastian Engelke  Université de Genève<EMAIL_ADDRESS>Manuel
Hentschel  Université de Genève<EMAIL_ADDRESS>Michaël Lalancette
Université du Québec à Montréal<EMAIL_ADDRESS>Frank Röttger  TU
Eindhoven<EMAIL_ADDRESS>
###### Abstract

Graphical models in extremes have emerged as a diverse and quickly expanding
research area in extremal dependence modeling. They allow for parsimonious
statistical methodology and are particularly suited for enforcing sparsity in
high-dimensional problems. In this work, we provide the fundamental concepts
of extremal graphical models and discuss recent advances in the field.
Different existing perspectives on graphical extremes are presented in a
unified way through graphical models for exponent measures. We discuss the
important cases of nonparametric extremal graphical models on simple graph
structures, and the parametric class of Hüsler–Reiss models on arbitrary
undirected graphs. In both cases, we describe model properties, methods for
statistical inference on known graph structures, and structure learning
algorithms when the graph is unknown. We illustrate different methods in an
application to flight delay data at US airports.

##  1 Introduction

Conditional independence for a random vector
$\mathbf{X}=(X_{i}:i\in\mathcal{V})$ with index set
$\mathcal{V}=\\{1,\dots,d\\}$ forms the basis for many tools of modern
statistical inference. A probabilistic graphical model represents the
components of the vector $\mathbf{X}$ and their conditional independencies
through a graph structure $\mathcal{G}=(\mathcal{V},\mathcal{E})$, where
variables are indexed by nodes $\mathcal{V}$ and connected by edges
$\mathcal{E}$. This notion naturally decomposes the joint distribution into
smaller models and thus simplifies statistical inference. Both undirected and
directed graphical models offer principled and interpretable ways to define
sparse statistical models, and they have lead to many successful applications.
For instance, they provide sparse representations of the joint dependence
structure of Gaussian models (Wainwright and Jordan, 2008) and several
extensions, such as skew-normal and Gaussian copula models (Capitanio et al.,
2003; Liu et al., 2009). Multiple graphical models have been developed for
continuous exponential family distributions (Inouye et al., 2016; Yang et al.,
2015), including dependence models for multivariate angular data (Klein et
al., 2020). The Ising model, an extensively studied model in particle physics,
is a graphical model for multivariate binary data (Ravikumar et al., 2010).
For more general categorical data, graphical models on log-linear interaction
models are commonly used for contingency tables (Darroch et al., 1980;
Lauritzen, 1996). Other graphical models for discrete data include the
discrete square root and Poisson graphical models (Inouye et al., 2016; Yang
et al., 2013). Bayesian networks encode conditional independence and density
factorization via separation statements for directed acyclic graphs (Pearl,
2009). They link to causal inference through structural causal models, where
the graphical structure allows for a causal interpretation (Maathuis et al.,
2019; Peters et al., 2017).

Multivariate extreme value theory is concerned with the theoretical
understanding and the statistical modeling of the distributional tail of the
random vector $\mathbf{X}$. Mathematically, this tail can be described by
three different, equivalent approaches: the point process approach, the block
maxima method leading to max-stable distributions, and the threshold
exceedances approach with limiting multivariate Pareto distributions. For all
of these methods, for a sample of size $n$, the effective sample size $k\ll n$
is much smaller since only the largest $k$ observations of $\mathbf{X}$ carry
relevant information on its tail and are used for inference. Therefore, sparse
extreme value models that make efficient use of the data are of utmost
importance (Engelke and Ivanovs, 2021). A natural question is thus how
powerful tools from graphical modeling can be employed in the framework of
statistics of extremes.

Compared to the classical applications of graphical models described above,
the focus here is different. Indeed, while probabilistic conditional
independence mainly concerns the bulk of the distribution, an appropriate
notion of an extremal graphical model should describe the dependence structure
in the tail of $\mathbf{X}$. Conceptually, there are at least two different
ways of approaching this problem: (a) one can assume that $\mathbf{X}$ is a
classical graphical model and study its extremal limit; or (b) one may
directly define graphical models for the limiting objects in the three
approaches to extremes described above, relying on a tailor-made notion of
extremal conditional independence. Eventually, both of the resulting theories,
where they overlap, yield the same graphical models for extremes.

Perspective (a) goes back to the analysis of extremes of time series, where
the extremal limit is a multiplicative random walk on a chain graph (Smith,
1992; Perfekt, 1994; Bortot and Coles, 2003; Segers, 2007). In the same
spirit, this random walk theory may be extended to trees (Segers, 2020), block
graphs (Asenova and Segers, 2023) and certain directed graphs (Segers and
Asenova, 2022); see Figure 1 for examples of graphs. A limitation of this
perspective is that technical assumptions are required on the distribution of
$\mathbf{X}$ to guarantee multivariate regular variation (2). Moreover, no
results exist that go beyond simple graph structures where the limit would no
longer have a random walk structure.

(a) Tree

(b) Block

(c) Cycle

(d) DAG

Figure 1:  Examples of graph structures on the set of nodes
$\mathcal{V}=\\{1,\dots,4\\}$.

The main difficulty of perspective (b) is that classical conditional
independence is not useful for any of the three approaches to extreme value
modeling. Indeed, considering the point process approach, a random set of
points is a non-standard object for graphical modeling. For block maxima,
given a max-stable vector $\mathbf{Z}$ with positive continuous density and
disjoint subsets $A,B,C\subset\mathcal{V}$, Papastathopoulos and Strokorb
(2016) show that classical conditional independence already implies
unconditional independence,

$\displaystyle{\mathbf{Z}_{A}\perp\\!\\!\\!\perp\mathbf{Z}_{B}\mid\mathbf{Z}_{C}}\quad\Rightarrow\quad\mathbf{Z}_{A}\perp\\!\\!\\!\perp\mathbf{Z}_{B}.$
(1)

Consequently, no nontrivial conditional independencies are possible in this
model class and densities of max-stable distributions therefore cannot
factorize on a graph into lower dimensional functions. For threshold
exceedances, the multivariate Pareto distribution $\mathbf{Y}$ is defined on a
non-product space and classical conditional independence is not applicable to
this random vector since even the support depends on the conditioning event.
Consequently, in all three cases, a different notion of conditional
independence is required. Engelke and Hitz (2020) introduce the concept of
extremal conditional independence for the class of multivariate Pareto
distributions, and Engelke et al. (2022) extend it to a unifying notion
involving the so-called exponent measure that covers all of the three
approaches.

Since then, a lot of progress has been made at the interface of graphical
models and extreme value theory, providing a better understanding of the
theory and a tool box of statistical methods. In this work we discuss the
fundamental ideas of extremal graphical models and give an overview over the
statistical methodology. While we concentrate on probabilistic graphical
models where graphs encode conditional independence properties (Lauritzen,
1996), there exist also other ways of using graphs to build parsimonious
statistical models (Lee and Joe, 2018; Vettori et al., 2020; Kiriliouk et al.,
2023). Our focus is the limiting model perspective (b) under the classical
asymptotic dependence framework (Resnick, 2008); graphical models for
asymptotically independent models are still fairly unexplored (Kulik and
Soulier, 2015; Papastathopoulos et al., 2017; Casey and Papastathopoulos,
2023).

Section 2 provides background on multivariate extremes and classical graphical
models. We introduce extremal notions of conditional independence and
graphical models in Section 3 and discuss implications for the limiting models
arising in the three approaches to extremes mentioned above. The remaining
sections present fundamental properties and statistical methodology for
extremal graphical models, with a focus on threshold exceedances. In
particular, we describe estimation methods if the underlying graph structure
$\mathcal{G}$ is known, and algorithms for structure learning if the graph is
unknown. Section 4 focuses on non-parametric methods on simple graph
structures, while Section 5 allows for general graphs but for the parametric
class of Hüsler–Reiss models, an analogue to Gaussian distributions in
extremes. Section 6 showcases many of the methods on a flight delay dataset
from the R package graphicalExtremes (Engelke et al., 2022).

##  2 Background

###  2.1 Multivariate extreme value theory

Multivariate extreme value theory studies the dependence between the largest
observations of the components of a random vector
$\mathbf{X}=(X_{i}:i\in\mathcal{V})$ with index set
$\mathcal{V}=\\{1,\dots,d\\}$. In order to focus only on the dependence
structure and abstract away the marginal distribution, we assume throughout
that $\mathbf{X}$ is normalized to have standard Pareto margins. In order to
mathematically describe this extremal dependence structure, a common
assumption is multivariate regular variation, which requires the existence of
a Borel measure $\Lambda$ on $E^{+}=[0,\infty]^{d}\setminus\\{\mathbf{0}\\}$,
such that for all Borel subsets $A\subset E^{+}$ bounded away from the origin
with $\Lambda(\partial A)=0$ the limit

$\displaystyle\lim_{n\to\infty}n\mathbb{P}(\mathbf{X}/n\in A)=\Lambda(A),$ (2)

exists; see Resnick (2008) for details. We say that $\mathbf{X}$ is in the
domain of attraction of the exponent measure $\Lambda$, which then has the
homogeneity property $\Lambda(tB)=t^{-1}\Lambda(B)$ for any $t>0$ and Borel
set $B\subset E^{+}$. If $\Lambda$ possesses a Lebesgue density $\lambda$,
then it satisfies the two properties

  1. (P1)

normalization: for any $m\in\\{1,\dots,d\\}$ it holds that
$\int_{\\{y_{m}>1\\}}\lambda(\mathbf{y}){d}\mathbf{y}=1$;

  2. (P2)

homogeneity: $\lambda(t\mathbf{y})=t^{-d-1}\lambda(\mathbf{y})$ for any $t>0$
and $\mathbf{y}\in E^{+}$.

Conversely, any positive function on $E^{+}$ satisfying these two properties
uniquely defines an exponent measure. This can be verified by equating the
normalization property (P1) to the marginal constraint that is known to
characterize the class of all possible spectral measures (Beirlant et al.,
2006, Section 8.2.3). For any subset $I\subset\mathcal{V}$, the marginal
exponent measure $\Lambda_{I}$ and its density $\lambda_{I}$ are defined by
integrating out the components $\mathcal{V}\setminus I$.

While (2) is a purely mathematical condition, it is in fact equivalent to
either of the three main approaches to multivariate extremes, which we
describe in the sequel. Let $\mathbf{X}_{1},\dots,\mathbf{X}_{n}$ be
independent copies of a random vector $\mathbf{X}$ with standard Pareto
margins.

  * •

(Point process) As $n\to\infty$, the random point set
$\sum_{i=1}^{n}\delta_{\mathbf{X}_{i}/n}$ converges in distribution in the
space of point process to the Poisson point process $\Pi_{\Lambda}$ on $E^{+}$
with intensity measure $\Lambda$; see (Resnick, 2008, Proposition 3.21) for
details.

  * •

(Block maximum) As $n\to\infty$, the componentwise maximum
$n^{-1}\max_{i=1,\dots,n}\mathbf{X}_{i}$ converges in distribution to a max-
stable random vector $\mathbf{Z}$ with distribution function

$\mathbb{P}(\mathbf{Z}\leq\mathbf{z})=\exp\left\\{-V(\mathbf{z})\right\\},\quad\mathbf{z}\in[0,\infty)^{d},$

where $V(\mathbf{y})=\Lambda(E^{+}\setminus[\mathbf{0},\mathbf{y}])$ for any
$\mathbf{y}\in E^{+}$.

  * •

(Peaks-over-threshold) As the threshold $u\to\infty$, the exceedances
$u^{-1}\mathbf{X}\mid\max_{j=1,\dots,{d}}X_{j}>u$ converge in distribution to
a multivariate Pareto vector $\mathbf{Y}$ with distribution function

$\mathbb{P}(\mathbf{Y}\leq\mathbf{y})=\frac{V(\min\\{\mathbf{y},\mathbf{1}\\})-V(\mathbf{y})}{V(\mathbf{1})},\quad\mathbf{y}\in\mathcal{L},$

where $\min\\{\mathbf{y},\mathbf{1}\\}$ denotes the componentwise minimum.

Figure 2:  All panels show $n$ samples of the random vector $\mathbf{X}$ and
illustrate the three main approaches to multivariate extremes: the blue points
in the left panel are the random point set that is approximated by the Poisson
point process $\Pi_{\Lambda}$; the solid blue point in the center panel is the
componentwise maximum with approximate max-stable distribution; the blue
points in the right panel are the threshold exceedances that are approximate
samples of a multivariate Pareto distribution.

######  Remark 2.1.

We emphasize that while the three approaches above appear fairly different,
they capture the same information on the tail of the random vector
$\mathbf{X}$. Indeed, either of the three limiting objects is in one-to-one
correspondence to the exponent measure $\Lambda$. When $\mathbf{X}$ is
attracted to $\Lambda$, we may equivalently say that it is in the domain of
attraction of the associated point process $\Pi_{\Lambda}$, max-stable
distribution $\mathbf{Z}$, or multivariate Pareto distribution $\mathbf{Y}$.
Another equivalent way of representing the extremal dependence in $\Lambda$ is
via the extreme value copula that attracts the copula of $\mathbf{X}$
(Gudendorf and Segers, 2010).

A widely used summary of extremal dependence is the extremal correlation

$\chi_{ij}=\lim_{u\to\infty}\mathbb{P}(X_{j}>u\mid X_{i}>u)\in[0,1],$

where the limit always exists if $\mathbf{X}$ is in the domain of attraction
of $\Lambda$. It can be reexpressed in terms of any of the limiting objects in
the above approaches:

$\chi_{ij}=2-V_{ij}(1,1)=2-\log\mathbb{P}(Z_{i}\leq 1,Z_{j}\leq
1)=\mathbb{P}(Y_{j}>1\mid Y_{i}>1).$

The set of $d\times d$ extremal correlation matrices $\chi$ with entries
$\chi_{ij}$, $i,j=1,\dots,d$, is a strict subset of the set of correlation
matrices with positive entries (Fiebig et al., 2017).

Another aspect of the extremal dependence of $\mathbf{X}$ is given by
considering the sub-faces of $E^{+}$ defined as

$E^{+}_{I}=\\{\mathbf{y}\in E^{+}\,:\,y_{j}\neq 0\,\,\forall\,j\in I\,\text{
and }\,y_{j}=0\,\forall\,j\notin I\\},$

for a subset $I\subset V$. If $\Lambda(E^{+}_{I})>0$, then the variables
$\mathbf{X}_{I}$ can be concomitantly extreme while the other variables
$\mathbf{X}_{\mathcal{V}\setminus I}$ are small. There are several statistical
methods in the literature that aim to detect such sub-faces with mass
(Chiapino et al., 2019; Simpson et al., 2020; Meyer and Wintenberger, 2023);
for a review see Engelke and Ivanovs (2021).

There are numerous parametric and non-parametric models for multivariate
extremes. Any model defined for one of the three approaches above directly
induces the associated models for the other two approaches. We discuss some
classical models and refer to Gudendorf and Segers (2010) for a detailed
overview.

######  Example 2.1 (Max-linear model).

Let ${\varepsilon}_{j}$, $j=1,\dots,p$ be independent standard Fréchet
variables and $A=(a_{ij})$ be a ${d}\times p$ matrix of non-negative
coefficients. Define a ${d}$-dimensional random vector $\mathbf{Z}$ with
entries

$\displaystyle Z_{i}=\max_{j=1,\dots,p}a_{ij}{\varepsilon}_{j},\quad
i=1,\dots,{d}.$ (3)

Assume that the rows of $A$ sum to $1$, ensuring that all $Z_{i}$ have
standard Fréchet distributions. The max-linear model is max-stable with
exponent measure $\Lambda$ supported on $p$ rays specified by the columns of
$A$, that is, the angles of the rays take $p$ possible values $a_{\cdot
j}/\|a_{\cdot j}\|$ with mass proportional to $\|a_{\cdot j}\|$.

The max-linear model is somewhat degenerate since neither the max-stable
distribution nor the corresponding exponent measure $\Lambda$ possess Lebesgue
densities. On the other hand, many continuous models that admit densities have
been proposed in the literature, such as the extremal logistic (Gumbel, 1960),
the extremal Dirichlet (Coles and Tawn, 1991) or extremal $t$-distribution
(Schlather, 2002; Opitz, 2013).

A highly flexible and very popular parametric model is the Hüsler–Reiss model
(Hüsler and Reiss, 1989). Similarly to a Gaussian distribution, it is
parameterized by a symmetric $(d\times d)$-matrix of bivariate dependence
parameters. Hüsler–Reiss distributions also arise as finite-dimensional
marginals of Brown–Resnick processes (Brown and Resnick, 1977; Kabluchko et
al., 2009), which are ubiquitous models in spatial extreme value modeling.
This model will be the focus of Section 5.

######  Example 2.2 (Hüsler–Reiss model).

The Hüsler–Reiss model is parameterized by a symmetric, strictly conditionally
negative definite matrix $\Gamma$, also called variogram, in the cone

$\mathcal{C}^{d}=\left\\{\Gamma\in\mathbb{[}0,\infty)^{d\times
d}:\Gamma=\Gamma^{\top},\;\;{\rm diag}{\Gamma}=\mathbf{0},\;\;v^{\top}\Gamma
v<0\;\forall\,\mathbf{0}\neq v\perp\mathbf{1}\right\\}.$

It is defined through its exponent measure density

$\lambda(\mathbf{y})\propto\prod_{i=1}^{d}y_{i}^{-1-1/d}\exp\Big{\\{}-\frac{1}{2}(\log
y-\mu_{\Gamma})^{\top}\Theta(\log y-\mu_{\Gamma})\Big{\\}},\quad\mathbf{y}\in
E^{+},$ (4)

where $\Theta=(P(-\Gamma/2)P)^{+}$, $P=I_{d}-\mathbf{1}\mathbf{1}^{\top}/d$ is
the projection matrix onto the orthogonal complement of $\mathbf{1}$, $A^{+}$
is the Moore–Penrose pseudoinverse of a matrix $A$, and
$\mu_{\Gamma}=P(-\Gamma/2)\mathbf{1}$; see Hentschel et al. (2022) for
details.

###  2.2 Conditional independence and graphical models

Before discussing graph separation and its role in probabilistic graphical
models, we summarize the necessary graph theory. We give complete details on
undirected graphs, whereas some more technical notions on directed graphs are
deferred to Section A.

A directed graph is a pair $\mathcal{G}=(\mathcal{V},\mathcal{E})$ consisting
of a vertex set $\mathcal{V}$ and an edge set
$\mathcal{E}\subseteq\mathcal{V}\times\mathcal{V}$, containing ordered pairs
of distinct vertices; the terms vertex and node are generally used
interchangeably. Without loss of generality we shall consider that
$\mathcal{V}=\\{1,\dots,d\\}$. An undirected graph is a graph in which the
presence of an edge $(i,j)\in\mathcal{E}$ also implies $(j,i)\in\mathcal{E}$;
the two edges are then considered equivalent and graphically represented by a
single undirected connection. Figure 1 shows examples of undirected graphs in
the first three panels and a directed graph in the fourth panel.

A path is a sequence of distinct vertices $(i_{0},i_{1},\ldots,i_{L})$ such
that $(i_{k-1},i_{k})\in\mathcal{E}$ or $(i_{k},i_{k-1})\in\mathcal{E}$ for
all $k=1,\ldots,L$. A cycle is a sequence of vertices
$(i=i_{0},i_{1},\ldots,i_{L}=i)$ such that $(i_{k-1},i_{k})\in\mathcal{E}$ or
$(i_{k},i_{k-1})\in\mathcal{E}$ for all $k=1,\ldots,L$, and where no vertex is
visited more than once except the start and end point $i$. Paths and cycles
are called directed if for each $k$, $(i_{k-1},i_{k})\in\mathcal{E}$. A
directed graph that contains no directed cycles is called a directed acyclic
graph (DAG). Two vertices $i,j$ are said to be connected if there exists a
(not necessarily directed) path $(i,\ldots,j)$ between them. If all pairs of
vertices in a graph are connected, the graph is referred to as connected,
otherwise it can be partitioned into connected components. If there is an edge
$(i,j)$ between vertices $i$ and $j$, they are called adjacent. The set of all
vertices adjacent to $i$ is called its neighborhood and denoted
$\text{ne}(i)$. A maximal subset of nodes where every pair of nodes is
adjacent is called a clique.

In order to connect graphical structures to probability distributions and make
them amenable for statistical modeling, the notion of separation in a graph is
crucial.

######  Definition 2.1 (Undirected separation).

Let $\mathcal{G}=(\mathcal{V},\mathcal{E})$ be an undirected graph and let
$A,B,C\subseteq\mathcal{V}$ be disjoint sets. The set $A$ is separated from
$B$ by $C$ if any path between $A$ and $B$ in $\mathcal{G}$ has a non-empty
intersection with $C$.

We call the smallest set of nodes that separates $A$ from $B$ a minimal
separator. An undirected graph where each minimal separator between any pair
of connected nodes is a clique is called decomposable. For DAGs, the
ubiquitous notion of separation is termed $d$-separation. It is considerably
more involved, so its formal definition is deferred to Appendix A. In short we
write $A\perp_{\mathcal{G}}B\mid C$ whenever $A$ is graph separated from $B$
by $C$ in an undirected graph or DAG $\mathcal{G}$. When some of these sets
are singletons, we typically write $a\perp_{\mathcal{G}}b\mid c$ for
$\\{a\\}\perp_{\mathcal{G}}\\{b\\}\mid\\{c\\}$.

######  Example 2.3.

The graph separation statements implied by the undirected graph in Figure 1(c)
are $1\perp_{\mathcal{G}}3\mid\\{2,4\\}$ and
$2\perp_{\mathcal{G}}4\mid\\{1,3\\}$. In contrast, the DAG in Figure 1(d) does
not imply that $2\perp_{\mathcal{G}}4\mid\\{1,3\\}$, but instead that
$2\perp_{\mathcal{G}}4\mid 1$.

Probabilistic graphical models on both undirected graphs and DAGs are
stochastic models that are constrained to satisfy certain conditional
independence relations. In this section we discuss conditional independence of
a random vector $\mathbf{X}=(X_{i}:i\in\mathcal{V})$ in the classical sense
and write $\mathbf{X}_{A}\perp\\!\\!\\!\perp\mathbf{X}_{B}\mid\mathbf{X}_{C}$
when $\mathbf{X}_{A}$ is conditionally independent of $\mathbf{X}_{B}$ with
respect to $\mathbf{X}_{C}$, for subsets $A,B,C\subset\mathcal{V}$; a
different notion of conditional independence for extremes is introduced in
Section 3. The conditional independencies are determined by through a Markov
property that links them to graph separation statements in $\mathcal{G}$.

######  Definition 2.2 (Global Markov property and graphical models).

We say that a random vector $\mathbf{X}$ indexed by $\mathcal{V}$ satisfies
the global Markov property with respect to the undirected graph or DAG
$\mathcal{G}$ when

$A\perp_{\mathcal{G}}B\mid
C\quad\Rightarrow\quad\mathbf{X}_{A}\perp\\!\\!\\!\perp\mathbf{X}_{B}\mid\mathbf{X}_{C}$

for all disjoint $A,B,C\subseteq\mathcal{V}$. In this case, we call
$\mathbf{X}$ a (probabilistic) graphical model with respect to graph
$\mathcal{G}$.

While the global Markov property allows a unified definition of graphical
models, it is customary and often more convenient to utilize other Markov
properties that are equivalent under mild assumptions. For example, for an
undirected graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$ the pairwise Markov
property requires that

$X_{i}\perp\\!\\!\\!\perp
X_{j}\mid\mathbf{X}_{\mathcal{V}\setminus\\{i,j\\}},\quad\forall(i,j)\not\in\mathcal{E}.$

It is equivalent to the global Markov property on $\mathcal{G}$ under the
assumption of a positive density with respect to a product measure, for
instance. We refer to Lauritzen (1996) for a discussion of this and other
Markov properties.

Both conditional independence and graph separation are ternary relations for
disjoint index sets $A,B,C\subseteq\mathcal{V}$. Such ternary relations have
been generalized to a more abstract notions of independence, which can be
studied through axiomatization. It has been shown that conditional
independence of probability distributions, undirected graph separation and
$d$-separation for DAGs all conform with a popular collection of axioms called
a semi-graphoid (Pearl, 1988; Lauritzen, 1996; Geiger et al., 1990). For more
details we further refer to Lauritzen (1996); Lauritzen and Sadeghi (2018).

##  3 Graphical modeling in extremes

###  3.1 Extremal graphical models

Extremal conditional independence was originally defined in Engelke and Hitz
(2020) as an extension of classical conditional independence to the framework
of threshold exceedances. It was used to define extremal graphical models for
multivariate Pareto distributions that admit a density. In Section 2.1 we
discussed that the extreme observations of some random vector $\mathbf{X}$ can
be described by the corresponding exponent measure $\Lambda$, regardless of
whether the perspective of point processes, block maxima or threshold
exceedances is taken.

We therefore follow the more general notion of extremal conditional
independence on the level of the exponent measure $\Lambda$ in Engelke et al.
(2022). This approach does not require the existence of densities and
naturally encompasses all three perspectives on extremes. The challenge is
that $\Lambda$ is not a probability measure since it explodes at the origin
and therefore has infinite mass on the space $E^{+}$. To circumvent this
issue, we define the set of all rectangles with positive $\Lambda$-mass that
are bounded away from the origin as

$\displaystyle\mathcal{R}(\Lambda)=\bigg{\\{}R=[\mathbf{a},\mathbf{b}]\,:\mathbf{a}\leq\mathbf{b},\,\Lambda(R)>0,\,\mathbf{0}\notin
R\bigg{\\}}.$

The following definition of conditional independence for $\Lambda$ essentially
requires the corresponding independence statement to hold for any restriction
of $\Lambda$ to rectangles with finite $\Lambda$-mass.

######  Definition 3.1.

Let $A,B,C\subseteq\mathcal{V}$ be a partition of $\mathcal{V}$. The exponent
measure $\Lambda$ is said to admit extremal conditional independence of $A$
and $B$ given $C$,

$\displaystyle{A\perp_{\Lambda}B\mid C},$ (5)

if we have classical conditional independence on any restriction of $\Lambda$
on rectangles bounded away from the origin, that is,

$\displaystyle{\mathbf{U}_{A}\perp\\!\\!\\!\perp\mathbf{U}_{B}\mid\mathbf{U}_{C}}\text{
for }\mathbf{U}\sim\Lambda(\cdot)/\Lambda(R)\quad\text{ for all
}R\in\mathcal{R}(\Lambda).$ (6)

This is trivially true for $A$ or $B$ being empty, and for $C=\emptyset$ we
say that $\Lambda$ admits extremal independence of $A$ and $B$, and write

${A\perp_{\Lambda}B}.$

If the sets $A$, $B$ and $C$ are not a partition of $\mathcal{V}$, then the
above definition remains the same with the test class in (6) replaced by
$\mathcal{R}(\Lambda_{A\cup B\cup C})$.

This definition of conditional independence is tailor-made for extremes and
turns out to be very natural. First, it satisfies the so-called semi-graphoid
axioms (Lauritzen, 1996), a set of desirable criteria for conditional
independence notions; see Engelke et al. (2022) for details. Second, if
$\Lambda$ admits a positive Lebesgue density $\lambda$, then (5) is equivalent
to the factorization of this density into $\lambda(\mathbf{y})=\lambda_{A\cup
B}(\mathbf{y}_{A\cup B})\lambda_{B\cup C}(\mathbf{y}_{B\cup
C})/\lambda_{C}(\mathbf{y}_{C})$ for all $\mathbf{y}\in E^{+}$.

It is worthwhile to note that extremal conditional independence has different
properties than those known from the classical, probabilistic notion for
random vectors. For instance, the case of extremal independence
${A\perp_{\Lambda}B}$ between two sets $A$ and $B$ that form a partition of
$\mathcal{V}$ is not equivalent to the factorization of the exponent measure
density into the marginal densities; in fact, in this case, a Lebesgue density
of $\Lambda$ can not exist. Instead, extremal independence is equivalent to
the fact that $\Lambda$ has only mass on the sub-faces $E^{+}_{A}$ and
$E^{+}_{B}$, that is,

$\displaystyle{A\perp_{\Lambda}B}\quad\Leftrightarrow\quad\Lambda(\mathbf{y}_{A}\neq\mathbf{0}_{A}\text{
and }\mathbf{y}_{B}\neq\mathbf{0}_{B})=0.$ (7)

This fact is highly desirable since it means that this notion is compatible
with the traditional notion of asymptotic independence; indeed, the latter is
said to hold for sub-vectors $\mathbf{X}_{A}$ and $\mathbf{X}_{B}$ of a random
vector $\mathbf{X}$ in the domain of attraction of $\Lambda$ whenever the
right-hand side of (7) is satisfied (Strokorb, 2020; Engelke et al., 2022). As
a consequence of (7), the existence of a density $\lambda$ implies the graph
$\mathcal{G}$ to be connected, a surprising property of extremal graphical
models.

We now define extremal graphical models through an extremal global Markov
property, which links graph separation for undirected graphs and DAGs
(Definitions 2.1 and A.1) with the notion of extremal conditional
independence.

######  Definition 3.2.

The exponent measure $\Lambda$ satisfies the global Markov property with
respect to an undirected graph or DAG $\mathcal{G}=(\mathcal{V},\mathcal{E})$
if

$\displaystyle A\perp_{\mathcal{G}}B\mid
C\quad\Rightarrow\quad{A\perp_{\Lambda}B\mid C}$

for any disjoint $A,B,C\subset\mathcal{V},$ where a (possibly empty) set $C$
separates $A$ from $B$. In this case, we say that $\Lambda$ is an extremal
graphical model with respect to $\mathcal{G}$.

This definition includes both undirected graphical models and directed
graphical models. For those two classes, we next discuss some fundamental
properties and give examples.

####  3.1.1 Undirected extremal graphical models

If $\mathcal{G}$ is an undirected graph, the global Markov property implies
the pairwise Markov property

$(i,j)\notin\mathcal{E}\quad\Rightarrow\quad{i\perp_{\Lambda}j\mid
V\setminus\\{i,j\\}},$

which allows us to read off extremal conditional independencies directly from
missing edges in $\mathcal{G}$. Extremal graphical models also contain
information on which variables can be extreme at the same time. More
precisely, for any $D\subset\mathcal{V}$ such that the sub-graph of
$\mathcal{G}$ restricted to $D$ is disconnected, it holds that
$\Lambda{}(E^{+}_{D})=0$ (Engelke et al., 2022, Corollary 6.3).

If the exponent measure $\Lambda$ admits a positive Lebesgue density on
$E^{+}$ then pairwise and global Markov properties are equivalent. If the
graph $\mathcal{G}$ is decomposable then, by a Hammersley–Clifford theorem,
this is further equivalent to the factorization of the exponent measure
density into marginal densities on the cliques $\mathcal{C}$ of the graph

$\displaystyle\lambda(\mathbf{y})=\prod_{C\in\mathcal{C}}\lambda_{C}(\mathbf{y}_{C})\Big{/}\prod_{D\in\mathcal{D}}\lambda_{D}(\mathbf{y}_{D}),\quad\mathbf{y}\in
E^{+},$ (8)

where $\mathcal{D}$ is a multiset containing intersections between the cliques
called separator sets (Engelke and Hitz, 2020, Theorem 1). This property is
the basis for efficient statistical modeling and inference, by specifying and
estimating the lower-dimensional exponent measure densities associated to
cliques and separator sets (Engelke and Hitz, 2020).

######  Example 3.1.

Let $\Lambda$ be an undirected extremal graphical model with respect to the
decomposable graph in Figure 1(b) that admits a positive Lebesgue density
$\lambda$. Then for all $\mathbf{y}\in E^{+}$,

$\lambda(\mathbf{y})=\lambda_{123}(\mathbf{y}_{123})\lambda_{34}(\mathbf{y}_{34})/\lambda_{3}(y_{3}).$

####  3.1.2 Directed extremal graphical models

If $\Lambda$ is an extremal graphical model with respect to a DAG
$\mathcal{G}$, the global Markov property again implies an extremal directed
factorization property. In particular, if an extremal graphical model
$\Lambda$ admits a positive Lebesgue density $\lambda$, then this density
factorizes on $\mathcal{G}$ as

$\displaystyle\lambda(\mathbf{y})$
$\displaystyle=\prod_{v\in\mathcal{V}}\lambda(y_{v}\mid\mathbf{y}_{\text{pa}(v)}),\quad\mathbf{y}\in
E^{+},$ (9)

where we define the conditional exponent measure density as
$\lambda(y_{v}\mid\mathbf{y}_{\text{pa}(v)})=\lambda(y_{v},\mathbf{y}_{\text{pa}(v)})/\lambda(\mathbf{y}_{\text{pa}(v)})$
and $\text{pa}(v)$ denotes the parent set of node $v$; see Appendix A. We note
that by the normalization property (P1) in Section 2.1 the conditional density
$\lambda(\cdot\mid\mathbf{y}_{\text{pa}(v)})$ is a probability density for
almost all $\mathbf{y}_{\text{pa}(v)}\in[0,\infty)^{|\text{pa}(v)|}$ and that,
by the homogeneity property (P2) it satisfies

$\displaystyle\lambda(ty_{v}\mid
t\mathbf{y}_{\text{pa}(v)})=t^{-1}\lambda(y_{v}\mid\mathbf{y}_{\text{pa}(v)}),$
(10)

for any $\mathbf{y}\in E^{+}$ and $t>0$.

######  Example 3.2.

Let $\Lambda$ be an extremal graphical model with respect to the DAG in Figure
1(d) that admits a positive Lebesgue density $\lambda$. Then,

$\lambda(\mathbf{y})=\lambda(y_{1})\lambda(y_{2}|y_{1})\lambda(y_{4}|y_{1})\lambda(y_{3}|y_{2},y_{4}),\quad\mathbf{y}\in
E^{+}.$

Extremal graphical models imply factorizations of the exponent measure density
$\lambda$ as in (8) and (10). This directly links to the point process
approach to extremes since the density $\lambda$ is the intensity of the
limiting Poisson point process $\Pi_{\Lambda}$ in this framework. In the next
two sections we discuss the implications of extremal conditional independence
and graphical models for the two other approaches: the peak-over-threshold and
the block maxima method.

###  3.2 Multivariate Pareto distributions

Of the three approaches to extreme values described in Section 2.1, the
framework of threshold exceedances and their limiting multivariate Pareto
distributions $\mathbf{Y}$ is the most naturally suited to graphical modeling.
This is due to the fact that their probability measure is proportional to the
corresponding exponent measure, in the sense that

$\displaystyle\mathbb{P}(\mathbf{Y}\in A)=\Lambda(A)/V(\mathbf{1}),\quad
A\in\mathcal{B}(\mathcal{L}),$

and, if it exists, so is the probability density function
$f(\mathbf{y})=\lambda(\mathbf{y})/V(\mathbf{1})$. In view of (8), if
$\Lambda$ is an extremal graphical model on the undirected graph
$\mathcal{G}$, the density $f$ factorizes on $\mathcal{G}$.

Since the support of $\mathbf{Y}$ is not a product space, we need to restrict
the random vector to suitable rectangles to obtain classical conditional
independence statements. For $m\in\mathcal{V}$, choose the rectangle
$\mathcal{L}^{(m)}=\\{x\in\mathcal{L}\mid x_{m}>1\\}$ in Definition 3.1 and
note that the corresponding random vector with law
$\Lambda(\cdot)/\Lambda(\mathcal{L}^{(m)})$ has the same distribution as

$\displaystyle\mathbf{Y}^{(m)}=\mathbf{Y}\mid\\{\mathbf{Y}_{m}>1\\}.$ (11)

This random vector has the interpretation as the limit when the $m$th
component is extreme, and it exhibits classical conditional independence
according to the underlying graph $\mathcal{G}$; in fact, this was the
original definition of extremal conditional independence in Engelke and Hitz
(2020). Working with $\mathbf{Y}^{(m)}$ often allows us to apply statistical
methods tailored to classical conditional independence models. Together their
supports cover the entire set $\mathcal{L}$ and a joint estimator based on
observations from each $\mathbf{Y}^{(m)}$, $m\in\mathcal{V}$, uses all
observations from $\mathbf{Y}$ itself. An important summary statistic based on
these random vectors is the extremal variogram rooted at node $m\in V$ defined
as the matrix $\Gamma^{(m)}$ with entries

$\displaystyle\Gamma^{(m)}_{ij}=\text{Var}\left\\{\log Y_{i}^{(m)}-\log
Y_{j}^{(m)}\right\\},\quad i,j\in V,$ (12)

whenever the right-hand side is finite (Engelke and Volgushev, 2022, Section
3). For many parametric models, this matrix has a one-to-one correspondence to
the model parameters, as for instance for the logistic and the Hüsler–Reiss
models.

######  Example 3.3.

If $\mathbf{Y}$ follows a Hüsler–Reiss multivariate Pareto distribution with
parameter matrix $\Gamma\in\mathcal{C}^{d}$, for $m\in\mathcal{V}$, then the
random vector in (11) has the representation

$\displaystyle\left(\log Y^{(m)}_{i}-\log Y^{(m)}_{m}\right)_{i\neq m}\sim
N(-\text{diag}(\Sigma^{(m)})/2,\Sigma^{(m)}),$ (13)

where $\Sigma^{(m)}$ is a $(d-1)\times(d-1)$-dimensional positive definite
covariance matrix obtained from $\Gamma$ via
$\Sigma_{ij}^{(m)}=\frac{1}{2}\left(\Gamma_{im}+\Gamma_{jm}-\Gamma_{ij}\right)$
for $i,j\neq m$. This implies that for Hüsler–Reiss distributions, all
extremal variograms rooted at nodes $m\in\mathcal{V}$ coincide and are equal
to the parameter matrix, that is,

$\Gamma=\Gamma^{(1)}=\dots=\Gamma^{(d)}.$

The extremal conditional independence
${\mathbf{Y}_{A}\perp_{\Lambda}\mathbf{Y}_{B}\mid\mathbf{Y}_{C}}$ has another
intuitive interpretation in terms of classical conditional independence
(Hentschel et al., 2022). Suppose that we observe an extreme event of any
variable in $C$, that is, $\max_{j\in C}Y_{j}>1$. Then, conditionally on
$\mathbf{Y}_{C}$, the sub-vectors $\mathbf{Y}_{A}$ and $\mathbf{Y}_{B}$ are
independent in the usual sense. For an extremal graphical structure
$\mathcal{G}$, this translates into an extremal local prediction property of
the form

$\left(Y_{i}\mid\mathbf{Y}_{\setminus\\{i\\}}=\mathbf{y}_{\setminus\\{i\\}}\right)\stackrel{{\scriptstyle
d}}{{=}}\left(Y_{i}\mid\mathbf{Y}_{\text{ne}(i)}=\mathbf{y}_{\text{ne}(i)}\right),\qquad\mathbf{y}\in
E^{+},\;\max_{j\in\text{ne}(i)}y_{j}>1,$

where $\text{ne}(i)$ denotes the neighborhood of vertex $i$ in $\mathcal{G}$;
see Section 2. Thus, knowing that an extreme event has occurred in the
neighborhood of $Y_{i}$, it suffices to know the neighbors of $Y_{i}$ for
predicting its value.

###  3.3 Max-stable distributions

Following (1) in the introduction, max-stable distributions seem not to be
suited for graphical modeling since their densities can only factorize
trivially. On the other hand, suppose that the exponent measure satisfies the
extremal conditional independence ${A\perp_{\Lambda}B\mid C}$. What does this
mean for the corresponding max-stable distribution? While there is no density
factorization, this statement still introduces sparsity in the max-stable
distribution function that is characterized by the exponent measure $\Lambda$.
This allows the natural construction of parsimonious max-stable statistical
models and therefore simplifies statistical inference. Moreover, in virtually
all exact simulation methods for a max-stable vector $\mathbf{Z}$, many
samples from densities proportional to the exponent measure density $\lambda$
have to be drawn (Dieker and Mikosch, 2015; Dombry et al., 2016; Liu et al.,
2019). Sparsity in form of extremal conditional independence statements of
$\Lambda$ therefore speed up simulation significantly (Engelke and Hitz, 2020,
Section 5.4).

Importantly, implication (1) only holds for continuous distributions
$\mathbf{Z}$, and consequently, a max-stable distribution that does not admit
a density may allow for nontrivial conditional independence structures. The
authors in Gissibl and Klüppelberg (2018) therefore define recursive max-
linear models on a DAG $\mathcal{G}$ by

$Z_{i}=\max_{k\in\text{pa}(i)}\max\\{c_{ki}Z_{k},c_{ii}{\varepsilon}_{i}\\},\quad
i=1,\dots,d,$ (14)

with independent, positive noise terms ${\varepsilon}_{i}$ and edge weights
$c_{ki}>0$, $i\in V$, $k\in\text{pa}(i)\cup\\{i\\}$. The maximum over an empty
set is defined as $0$ and thus, for a source node $i$ for which $\text{pa}(i)$
is empty, we have $Z_{i}=c_{ii}{\varepsilon}_{i}$. While Gissibl and
Klüppelberg (2018) do not fix the noise distribution, we follow later papers
and assume that ${\varepsilon}_{i}$ are standard Fréchet distributed. The
random vector $\mathbf{Z}$ in (14) is then a special case of the max-linear
model in Example 2.1 with $p=d$ noise terms and the coefficient matrix $A$ is
defined by $a_{ii}=c_{ii}$ and

$a_{ij}=c_{jj}\max_{\omega\in P_{ji}}\prod_{(u,v)\in\omega}c_{uv}\quad i,j\in
V,i\neq j,$

(Gissibl and Klüppelberg, 2018, Theorem 2.2), where $P_{ji}$ is the set of all
directed paths from $j$ to $i$ in the DAG $\mathcal{G}$; here a path $\omega$
is to be understood as a set of edges and the product ranges over these edges.
While the coefficient matrix $A$ is fully identifiable from the joint
distribution of $\mathbf{X}$, the DAG and edge weights in the formulation of
(14) are not: different graph structures and sets of edge weights can lead to
the same coefficient matrix $A$, hence the same observational distribution
(Gissibl et al., 2021). A notable exception is the case where $\mathcal{G}$ is
a directed tree. Its structure is then fully identifiable and can be learned
via Chu–Liu/Edmonds’ algorithm (Tran et al., 2021).

Nevertheless, for recursive max-linear models on an arbitrary DAG, Klüppelberg
and Krali (2021) describe algorithms to determine a causal order implied by
the unknown graph $\mathcal{G}$, and to learn this order by estimating the
max-linear coefficients $a_{ij}$ using scalings. Building on this, Krali et
al. (2023) determine sufficient conditions for the edge weights $c_{ki}$ to be
identifiable even in the presence of unobserved variables in the graph, a
problem reminiscent of that considered in Gnecco et al. (2021) for heavy-
tailed linear models.

Like all structural equation models on DAGs the model (14) satisfies the
global Markov property with respect to $\mathcal{G}$, so
${\mathbf{Z}_{A}\perp\\!\\!\\!\perp\mathbf{Z}_{B}\mid\mathbf{Z}_{C}}$ whenever
$A$ and $B$ are $d$-separated by $C$ in $\mathcal{G}$ (Pearl, 2009, Theorem
1.4.1). However, due to the special structure of recursive max-linear models,
they satisfy additional conditional independence relations that are not
enforced by $d$-separation. In fact, Améndola et al. (2022) show that
conditional independence in such models is characterized by $*$-separation, a
graphical criterion strictly weaker than $d$-separation. This is illustrated
by the Cassiopeia graph $\mathcal{G}$ in Figure 3; while nodes $1$ and $3$ are
$d$-connected relative to $\\{4,5\\}$, they are $*$-separated by $\\{4,5\\}$,
and any recursive max-linear model on $\mathcal{G}$ satisfies
${Z_{1}\perp\\!\\!\\!\perp Z_{3}\mid\mathbf{Z}_{\\{4,5\\}}}$. Interestingly,
$*$-separation and $d$-separation share the same Markov equivalence classes
(Améndola et al., 2021), meaning that two DAGs imply the same conditional
independence relations on arbitrary distributions if and only if they do so in
recursive max-linear models.

Figure 3: Cassiopeia graph

Since recursive max-linear models are max-stable distributions, it is
interesting to understand the relationship between extremal conditional
independence statements on the corresponding exponent measure $\Lambda$ in the
sense of Definition 3.1, and conditional independence in the classical sense
discussed above. As a partial result in this direction, Engelke et al. (2022)
show the exponent measure of any recursive max-linear model is also a directed
extremal graphical model on the same DAG in the $\Lambda$ sense. On the other
hand, it is still unknown whether $d$-separation fully characterizes extremal
conditional independence or if it is strictly stronger. It is also an open
question whether in these models, extremal conditional independence is
stronger than, weaker than, equivalent to (or none of the above)
$*$-separation.

##  4 Simple graph structures

###  4.1 Properties

The complexity of the graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$ determines
the difficulty of statistical inference for the corresponding extremal
graphical models. Here, complexity does not only refer to the number of edges
$|\mathcal{E}|$, but also to the structural properties of the graph
$\mathcal{G}$. In this section we discuss simple graph structures, which often
allow statistical methods without any parametric assumption on the exponent
measure $\Lambda$.

A tree $\mathcal{T}=(\mathcal{V},\mathcal{E})$ with nodes $\mathcal{V}$ and
edge set $\mathcal{E}$ is a connected undirected graph without cycles, such as
in Figure 1(a). We say that $\Lambda$ is an extremal tree model if it is an
extremal graphical model with respect to a tree. Trees are the sparsest
connected graphs; they contain $|\mathcal{E}|=d-1$ edges, which is much
smaller than the number $d(d-1)/2$ of all possible edges in an undirected
graph. A convenient property of trees is the fact that between any two nodes
$h,l\in\mathcal{V}$ there is a unique path on $\mathcal{T}$ denoted by
$\text{ph}(hl;\mathcal{T})$, again to be understood as a set of edges. To get
a first intuition on how the graph structure implies certain properties of the
extremal dependence structure, we can consider the extremal correlation
coefficient. For an extremal tree model $\Lambda$ on a tree $\mathcal{T}$ it
can be shown (Engelke and Volgushev, 2022, Proposition 5) that

$\displaystyle\chi_{hl}\leq\chi_{ij}\quad\forall(i,j)\in\text{ph}(hl;\mathcal{T}).$
(15)

This result means that extremal dependence decays with the distance on the
tree: the closer two nodes on a path on the tree, the larger the extremal
correlation and the stronger the dependence. These inequalities between
different extremal correlations can be used to learn the tree structure in a
data driven way. It turns out however, that other coefficients, namely the
extremal variograms $\Gamma^{(m)}$ as defined in (12), are much more natural
for this purpose. The reason is that they contain information on the
underlying tree more directly through the tree metric property

$\displaystyle\Gamma_{hl}^{(m)}=\sum_{(i,j)\in\text{ph}(hl;\mathcal{T})}\Gamma_{ij}^{(m)},$
(16)

that is, the entry $\Gamma_{hl}^{(m)}$ of two non-adjacent nodes can be
computed by summing up the entries of $\Gamma^{(m)}$ on the path between these
nodes (Engelke and Hitz, 2020; Asenova et al., 2021; Engelke and Volgushev,
2022). Consequently, the whole matrix $\Gamma^{(m)}$ is determined by the
$d-1$ entries corresponding to the edges $\mathcal{E}$ of the tree; this
principle is illustrated in Figure 4. As long as the extremal variograms
exist, this property is independent of any parametric assumptions on the
exponent measure.

The above results suggest that for trees, the bivariate properties completely
imply the multivariate dependence structure. This is formalized in terms of
the density of the exponent measure which factorizes by (8) into the bivariate
marginal densities (Engelke and Hitz, 2020, Theorem 1) as

$\displaystyle\lambda(\mathbf{y})=\prod_{(ij)\in\mathcal{E}}\frac{\lambda_{ij}(y_{i},y_{j})}{y_{i}^{-2}y_{j}^{-2}}\prod_{i\in
V}y_{i}^{-2}\quad\mathbf{y}\in E^{+},$ (17)

where the terms $y_{i}^{-2}$, $i\in\mathcal{V}$, correspond to the univariate
marginal distributions of $\lambda$. A valid extremal tree model on
$\mathcal{T}$ can therefore be constructed by choosing, for each edge
$(i,j)\in\mathcal{E}$, an arbitrary bivariate exponent measure density (e.g.,
logistic, Hüsler–Reiss) and combining them as described above. For chain
graphs, a subset of tree structures corresponding to Markov chains, this
modeling approach was proposed by Coles and Tawn (1991); Smith et al. (1997).

A slightly more general class than trees are block graphs, which are defined
as connected, decomposable graphs where the separator sets are singletons; see
Figure 1(b) for an example. In particular, there is a unique shortest path
between any two nodes of a block graph. Therefore, some of the properties
above generalize naturally to this class of graphs, using this unique shortest
path in place of the unique path. For instance, a factorization similar to
(17) allows the construction of parametric extremal block graph models by
specifying exponent measure densities on the blocks (Engelke and Hitz, 2020,
Section 5.1). Moreover, the additivity property of extremal variograms has
been shown to continue to hold if each sub-model on the blocks is Hüsler–Reiss
(Engelke and Hitz, 2020; Asenova and Segers, 2023).

Trees and block graphs can be seen as simple ways of combining lower-
dimensional distributions into a sparse, higher-dimensional model. For this
reason, they are a natural starting point to go beyond multivariate regular
variation and build sparse asymptotic independence models. For instance,
extremes of graphical models that allow for asymptotic independence have been
studied in Papastathopoulos et al. (2017, 2023) for Markov chains and for
trees and block graphs in Casey and Papastathopoulos (2023). Engelke et al.
(2022) use the general $\Lambda$ conditional independence in Definition 3.1 to
define a Hüsler–Reiss tree model which allows for mass on sub-faces of
$E^{+}$, and they propose a construction principle for asymptotic independence
trees.

$\displaystyle\Gamma=\left(\begin{matrix}0&6&6&\underline{12}&\underline{10}&\underline{22}\\\
6&0&4&\underline{10}&\underline{8}&\underline{20}\\\
6&4&0&6&4&\underline{16}\\\ \underline{12}&\underline{10}&6&0&10&10\\\
\underline{10}&\underline{8}&4&10&0&\underline{20}\\\
\underline{22}&\underline{20}&\underline{16}&10&\underline{20}&0\end{matrix}\right)$

Figure 4:  Example of a block graph and an associated extremal variogram of a
Hüsler–Reiss distribution. Edges are labelled with their corresponding
variogram entry. Underlined entries, corresponding to non-edges, can be
computed using (16), for example $\Gamma_{16}=6+6+10=22$.

###  4.2 Parameter estimation

Parameter estimation for trees and block graphs is particularly easy. Instead
of estimating a full $|\mathcal{V}|$-dimensional model, we can estimate the
lower-dimensional models on the cliques of the graph and then combine them
through the Hammersley–Clifford theorem; recall that the cliques in trees are
all two-dimensional. Parameter estimation for extremal block graph models
therefore boils down to the estimation of low-dimensional, parametric extreme
value models. The latter is a classical problem in extreme value statistics
and many statistical methods have been proposed.

We concentrate here on the non-parametric, empirical estimator of the extremal
variogram $\Gamma^{(m)}$ in (12), since it will serve as input for most graph
structure learning algorithms. Let $\mathbf{X}_{1},\dots,\mathbf{X}_{n}$ be
independent copies of the $d$-dimensional random vector $\mathbf{X}$ in the
domain of attraction of a multivariate Pareto distribution $\mathbf{Y}$. Based
on the idea of extremal increments (Heffernan and Tawn, 2004; Engelke et al.,
2014, 2015), the empirical variogram $\hat{\Gamma}^{(m)}$ is a method of
moments estimator with entries

$\hat{\Gamma}_{ij}^{(m)}=\widehat{\text{Var}}\Big{(}\log(1-\hat{F}_{i}(X_{ti}))-\log(1-\hat{F}_{j}(X_{tj})):\hat{F}_{m}(X_{tm})\geq
1-k/n\Big{)},$ (18)

where $\widehat{\text{Var}}$ denotes the sample variance, $k=k_{n}$ is an
intermediate sequence and $\hat{F}_{i}$ is the empirical distribution function
of $X_{1i},\dots,X_{ni}$. A close look at this estimator reveals that, first,
all marginals are empirically normalized to standard Pareto distributions
through the transformations $1/(1-\hat{F}_{i})$. Then, $k$ exceedances in the
$m$th component are selected to obtain approximate samples of
$\mathbf{Y}^{(m)}$, and finally these samples are used to empirically estimate
the variance in (12). Under the assumption $k/n\to 0$ as $n\to\infty$ and mild
conditions on the underlying data generation, this estimator is consistent for
$\Gamma_{ij}^{(m)}$ (Engelke and Volgushev, 2022), and enjoys subexponential
concentration properties (Engelke et al., 2022). In many applications it makes
sense to consider the joint estimator

$\displaystyle\hat{\Gamma}=\frac{1}{d}\sum_{m=1}^{d}\hat{\Gamma}^{(m)},$ (19)

in order to combine information from exceedances in all variables.

For block graphs, for each clique $C\subset\mathcal{V}$, the empirical
extremal variogram $\hat{\Gamma}_{CC}$ of the sub-model can be estimated using
only the observations of the components $\mathbf{X}_{C}=(X_{i}:i\in C)$, and
in many parametric models this directly implies estimates of the parameters.
Alternatively, Engelke and Hitz (2020) propose maximum likelihood estimation
for the parameters on each clique, and Asenova et al. (2021) use the pairwise
extremal coefficients estimator of Einmahl et al. (2018). We stress again that
any method from extreme value statistics can be used here, the graphical
structure providing us with a way to combine the estimated models on the
cliques to a valid, $d$-dimensional model on the block graph.

###  4.3 Structure learning

The particular properties of trees allow for more efficient statistical
structure learning methods than in the case of general graphs. Most of these
methods rely on the notion of the minimum spanning tree. For a set of
symmetric weights $\rho_{ij}>0$ associated with each pair of nodes
$i,j\in\mathcal{V}$, $i\neq j$, the latter is defined as the tree structure
that minimizes the sum of distances on that tree, that is,

$\displaystyle\mathcal{T}_{\text{mst}}=\mathop{\operatorname{arg\,min}}_{\mathcal{T}=(\mathcal{V},\mathcal{E})}\sum_{(i,j)\in\mathcal{E}}\rho_{ij},$
(20)

where $\mathcal{T}$ ranges over all spanning trees on $\mathcal{V}$. Given the
set of weights, there exist the well-known greedy algorithms by Prim (Prim,
1957) and Kruskal (Kruskal, 1956) that constructively and efficiently solve
this minimization problem. The crucial ingredient for minimum spanning trees
are the weights $\rho_{ij}$, which ideally should be chosen in such a way that
$\mathcal{T}_{\text{mst}}$ recovers the true underlying tree structure that
represents the conditional independence relations. In the classical case of
Gaussian graphical models, such tree recovery is possible by using weights
$\rho_{ij}=\log(1-r_{ij}^{2})/2$, where $r_{ij}$ are the correlation
coefficients (Drton and Maathuis, 2017). Interestingly, the assumption of
Gaussianity is crucial and the result no longer holds outside this specific
parametric class.

For extremal tree models, consistent tree recovery turns out to hold more
generally without the need for any parametric assumption on the model class.
Indeed, there are two natural summary statistics for the strength of extremal
dependence as candidates for weights in the minimum spanning tree: the
extremal correlation in (15) and the extremal variogram in (16). Suppose that
$\mathbf{Y}$ is an extremal tree model on the tree $\mathcal{T}$ and let
$\hat{\chi}$ and $\hat{\Gamma}$ be the empirical estimators based on $n$
observations of $\mathbf{X}$ in the domain of attraction of $\mathbf{Y}$. We
further denote by $\hat{\mathcal{T}}_{\chi}$ and $\hat{\mathcal{T}}_{\Gamma}$
the minimum spanning trees with weights $\rho_{ij}=-\log\hat{\chi}_{ij}$ and
$\rho_{ij}=\hat{\Gamma}_{ij}$, respectively. Then, under the conditions for
consistency of the empirical estimators, we have

$\displaystyle\mathbb{P}(\hat{\mathcal{T}}_{\chi}=\mathcal{T})\to
1,\quad\mathbb{P}(\hat{\mathcal{T}}_{\Gamma}=\mathcal{T})\to 1,\quad
n\to\infty;$ (21)

for details see (Engelke and Volgushev, 2022, Theorem 2). This result is
surprising since it guarantees non-parametric tree recovery, which is stronger
than in the non-extreme setting. Let us give some intuition why this is
possible. A tree graphical model is characterized by the bivariate
distributions on the edges $(i,j)\in\mathcal{E}$. In the case of extremal tree
models, these bivariate distributions are represented by the bivariate
exponent measure densities $\lambda_{ij}(y_{i},y_{j})$ in the factorization
(17). The crucial property of these densities is the homogeneity described in
Section 2.1, which essentially reduces this bivariate distribution to the
univariate distribution of the angular part. This implies that any extremal
tree model has a random walk structure on the underlying tree $\mathcal{T}$ as
shown in (Segers, 2020, Theorem 1) and (Engelke and Volgushev, 2022,
Proposition 1). This explains why tree learning is easier in the extremal than
in the non-extreme setting.

The tree recovery result in (21) can be further extended to the setting where
the number $d(n)$ of nodes in the tree grows with the sample size. In
particular, based on concentration bounds for the empirical variogram
$\hat{\Gamma}$ in Engelke et al. (2022), (Engelke and Volgushev, 2022, Theorem
4) show that extremal trees can still be consistently learned in high
dimensions where $d(n)$ grows exponentially faster than $k$.

##  5 Hüsler–Reiss graphical models

###  5.1 Properties

As described in the previous section, for simple graph structures tailor-made
and often non-parametric methods yield efficient statistical inference. For
general graphs, a suitable parametric model class is typically required to
specify arbitrary conditional independence structures and to perform parameter
estimation and structure learning. For extremal graphical models, the
Hüsler–Reiss distribution (Hüsler and Reiss, 1989) constitutes the only
parametric family with the desired properties. Most importantly, in this
class, extremal conditional independence is encoded in a parametric way,
namely as the zero pattern of a certain Hüsler–Reiss precision matrix. This is
one of the reasons why this distribution can be seen as an analogue of the
Gaussian distribution in extremes.

Recall from Example 2.2 that an extreme value model with Hüsler–Reiss
distribution is parameterized by a variogram matrix $\Gamma$ and defined
through the density of its exponent measure. The matrix $\Theta$ that appears
in this density is called the Hüsler–Reiss precision matrix and plays a
crucial role for graphical modeling. It is in one-to-one correspondence with
the variogram matrix through the mapping $\Theta=(P(-\Gamma/2)P)^{+}$, which
is homeomorphic between $\mathcal{C}^{d}$ and the set
$\mathcal{S}^{d}_{\mathbf{1}}$ of symmetric, positive semidefinite matrices
with zero row sums and rank ${d}-1$ (Hentschel et al., 2022).

The Hüsler–Reiss distribution is widely used in applications due to its
flexibility for statistical modeling (Engelke et al., 2019; de Fondeville and
Davison, 2018; Thibaud et al., 2013). It turns out to be the only class of
exponent measure densities that exhibits the structure of a pairwise
interaction model (Lalancette, 2023), an interpretable and computationally
advantageous exponential family that is ubiquitous in high dimensional
dependence and graphical modeling (Yang et al., 2015; Lin et al., 2016; Klein
et al., 2020). Besides these properties, the main importance of the
Hüsler–Reiss family stems from the fact that extremal conditional independence
and the extremal graphical structure can directly be read off from the
precision matrix as

${i\perp_{\Lambda}j\mid\mathcal{V}\setminus\\{i,j\\}\quad\Leftrightarrow\quad\Theta_{ij}=0};$

(Engelke and Hitz, 2020; Hentschel et al., 2022). This fact is the basis for
virtually all estimation and structure learning algorithms that we present in
the sequel.

For Hüsler–Reiss tree models, it is enough to specify the parameter matrix
$\Gamma$ on the edges of the graph and the remaining entries are implied by
the graphical structure through property (16). This result can in fact be
generalized to arbitrary graphs $\mathcal{G}=(\mathcal{V},\mathcal{E})$. Let
$\mathring{\Gamma}$ be a partially specified matrix with given entries for
$(i,j)\in\mathcal{E}$ such that all (fully specified) sub-matrices
corresponding to cliques are valid variograms. The matrix completion problem
of finding a valid variogram matrix $\Gamma\in\mathcal{C}^{d}$ with
corresponding Hüsler–Reiss precision matrix $\Theta$ such that

$\displaystyle\begin{cases}\Gamma_{ij}=\mathring{\Gamma}_{ij}&(i,j)\in\mathcal{E},\\\
{\Theta}=0&(i,j)\not\in\mathcal{E},\\\ \end{cases}$ (22)

then has a unique solution (Hentschel et al., 2022), where, if the graph is
non-decomposable, we need the additional assumption that any valid, possibly
non-graph structured solution exists. In practice, the solution $\Gamma$ is
obtained through a recursive algorithm, which for decomposable graphs finishes
after finitely many steps. Figure 5 illustrate the matrix completion problem
and its solution on an example of a decomposable graph.

It is worthwhile to highlight some similarities and difference of Hüsler–Reiss
and Gaussian distributions. At first glance, (4) appears to be the density of
a multivariate lognormal distribution. A crucial distinction is in that the
precision matrix $\Theta$ is rank deficient, and in the “direction” of the
vector $\mathbf{1}$, the density decays at a Pareto rather than lognormal
rate. This ensures homogeneity of the exponent measure density $\lambda$, and
also explains why it is not integrable on $E^{+}$. To make the connection to
Gaussian distributions more explicit, we consider the multivariate Pareto
distribution $\mathbf{Y}$ that corresponds to a Hüsler–Reiss model. On the
log-scale, it admits the stochastic representation

$\displaystyle\log\mathbf{Y}\mid\\{\mathbf{1}^{\top}\log\mathbf{Y}>0\\}=R\mathbf{1}+\mathbf{W},$
(23)

where $R$ is a standard exponential random variable and, independently,
$\mathbf{W}$ has a degenerate Gaussian distribution with covariance matrix
$\Sigma=\Theta^{+}$. Compared to classical Gaussian graphical models, the main
difficulty of working with Hüsler–Reiss graphical models, both in terms of
theory and statistical methodology, is the zero row sums property on the
precision matrix $\Theta$, as well as its rank which is constrained to be
$d-1$.

$\displaystyle\left(\begin{matrix}0&10&4&3&\mathrm{?}\\\
10&0&18&\mathrm{?}&\mathrm{?}\\\ 4&18&0&3&\mathrm{?}\\\ 3&\mathrm{?}&3&0&6\\\
\mathrm{?}&\mathrm{?}&\mathrm{?}&6&0\end{matrix}\right)\mapsto\left(\begin{matrix}0&10&4&3&\underline{9}\\\
10&0&18&\underline{15}&\underline{21}\\\ 4&18&0&3&\underline{9}\\\
3&\underline{15}&3&0&6\\\
\underline{9}&\underline{21}&\underline{9}&6&0\end{matrix}\right)$

Figure 5:  Example of a decomposable graph and an associated partial and
completed extremal variogram. Underlined entries, corresponding to non-edges,
are uniquely recovered from the known variogram entries on edges. Note that
unlike in the case of a block graph (Figure 4), the solution is not simply
obtained by summing the entries along paths. For example, $\Gamma_{24}=15$ is
neither equal to $\Gamma_{21}+\Gamma_{14}=13$, nor to
$\Gamma_{23}+\Gamma_{34}=21$.

###  5.2 Parameter estimation

####  5.2.1 Maximum likelihood on known graph

Similarly to Section 4.2, suppose that we have $n$ independent copies of a
random vector $\mathbf{X}$ in the domain of attraction of a Hüsler–Reiss
distribution with parameter matrix $\Gamma$. In order to use properties of
Gaussian maximum likelihood theory we use the transformation
$\mathbf{Y}^{(m)}$ of the corresponding multivariate Pareto distribution
$\mathbf{Y}$, whose increments have a lognormal distribution with mean vector
$-{\rm diag}(\Sigma^{(m)})/2$ and covariance matrix $\Sigma^{(m)}$; see
Example 3.3. Similarly to the empirical variogram estimator in (18), out of
the $n$ samples, we choose the $k$ largest observations in the $m$th component
and apply the log increments transformation in (13). The resulting dataset
then has approximately the surrogate positive semidefinite Gaussian log-
likelihood

$\displaystyle\ell(\Theta;\overline{\Gamma})=\log\operatorname{Det}(\Theta)+\frac{1}{2}{\rm
tr}(\Theta\overline{\Gamma}),$ (24)

which we parameterize in terms of Hüsler–Reiss precision matrices
$\Theta\in\mathcal{S}^{d}_{\mathbf{1}}$, and $\overline{\Gamma}$ is an
estimator of $\Gamma$; see Röttger et al. (2023); Hentschel et al. (2022) for
a detailed derivation. It is a surrogate likelihood since the mean is
empirically normalized to zero and $\overline{\Gamma}$ can be any estimator of
$\Gamma$, and it is approximate since the exceedances of $\mathbf{X}$ only
follow a multivariate Pareto distribution in the limit. Nevertheless, this
likelihood turns out to be computationally tractable and accurate enough to
provide estimators with good statistical properties. In order to use all data,
the most common choice is to combine data from conditioning on all components
$m\in\mathcal{V}$ and use the joint empirical variogram
$\overline{\Gamma}=\hat{\Gamma}$ from (19).

Suppose now that the Hüsler–Reiss distribution is an extremal graphical model
on the known graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$. In order to obtain
a first estimator that is consistent for the entries of $\Theta$ and has the
zeros at the correct entries, we can maximize the log-likelihood (24) under
the constraint that

$\displaystyle\Theta_{ij}=0\quad(i,j)\not\in\mathcal{E}.$ (25)

Surprisingly, the maximizer $\hat{\Theta}_{\mathcal{G}}$ of this optimization
problem is identical to the solution of the completion problem (22) with
partial variogram matrix chosen as $\mathring{\Gamma}_{ij}=\overline{\Gamma}$
for $(i,j)\in\mathcal{E}$. In practice, we can therefore proceed as follows:
choose an estimator of the matrix $\Gamma$ such as the joint empirical
variogram $\hat{\Gamma}$; use the entries of $\hat{\Gamma}$ on the edges of
$\mathcal{G}$ in the partial matrix $\mathring{\Gamma}$ as input for the
matrix completion problem together with (25); apply the algorithm of Hentschel
et al. (2022) to obtain the completion $\hat{\Gamma}_{\mathcal{G}}$ and note
that it exists and is unique since $\mathring{\Gamma}$ has the valid
completion $\hat{\Gamma}$. The resulting estimator
$\hat{\Gamma}_{\mathcal{G}}$ is consistent and its precision
$\hat{\Theta}_{\mathcal{G}}$ has extremal graph structure $\mathcal{G}$
(Hentschel et al., 2022, Section 5.2).

####  5.2.2 Colored graphical models

In addition to the graphical constraint (25), it might be desirable to further
reduce the dimensionality of the parameter matrix. A potential approach are
parameter symmetries, which can be visualized through colored graphs (Röttger
et al., 2023). For an undirected graph $\mathcal{G}=(\mathcal{V},\mathcal{E})$
let $\lambda:\mathcal{E}\to\\{1,\dots,r\\}$ be an edge coloring function that
maps every edge to a natural number $1,\ldots,r$. In inference, graph
colorings can for example be determined using clustering algorithms. A
Hüsler–Reiss restricted concentration (RCON) model is a linear submodel of
Hüsler–Reiss graphical models on $\mathcal{G}$ with additional symmetries

$\Theta_{st}=\Theta_{uv}\text{ whenever }\lambda(st)=\lambda(uv)\text{ for
}(s,t),(u,v)\in\mathcal{E}.$

Such models allow for a simple surrogate maximum likelihood estimator, as all
constraints are affine in the natural parameter of the underlying exponential
family. Alternatively, one may impose the symmetries on the level of the
variogram matrix $\Gamma$, resulting in the Hüsler–Reiss restricted variogram
(RVAR) model with

$\Gamma_{st}=\Gamma_{uv}\text{ whenever }\lambda(st)=\lambda(uv)\text{ for
}(s,t),(u,v)\in\mathcal{E}.$

The RVAR model usually does not allow an affine description in the precision
matrix $\Theta$, such that simple surrogate maximum likelihood estimation may
be difficult to apply. As an alternative, Röttger et al. (2023) employ the
two-step estimation procedure of Lauritzen and Zwiernik (2022), which allows
to efficiently impose linear constraints in a mixed parametrization of an
exponential family. Depending on the number of edge color classes $r$, both
RCON and RVAR models result in a drastic parameter reduction while showing
good performance in applications (Röttger et al., 2023).

####  5.2.3 Positive dependence

Multivariate total positivity ($\text{MTP}_{2}$) and positive association are
properties of a multivariate random vector that describe more or less strong
notions of positive dependence. Many datasets naturally exhibit an intrinsic
positivity and statistical methodology based on this notion can therefore
improve accuracy and interpretability; see Lauritzen and Zwiernik (2022) and
references therein.

For multivariate Pareto distributions, Röttger et al. (2023) introduce
extremal $\text{MTP}_{2}$ ($\text{EMTP}_{2}$) and extremal association by
requiring that for all $m\in\mathcal{V}$ the vectors $\mathbf{Y}^{(m)}$ in
(11) be $\text{MTP}_{2}$ and associated in the usual sense, respectively. For
Hüsler–Reiss distributions, $\text{EMTP}_{2}$ is equivalent to the parametric
constraint that the parameter matrix $\Theta$ is a graph Laplacian, that is,

$\displaystyle\Theta_{ij}\leq 0\quad\forall i\neq j.$ (26)

In particular, this implies that all Hüsler–Reiss tree models are
$\text{EMTP}_{2}$. To estimate a Hüsler–Reiss model that is $\text{EMTP}_{2}$,
Röttger et al. (2023) propose to maximize likelihood (24) under the linear
constraint (26). The resulting estimator $\hat{\Theta}_{+}$ typically contains
many zeros, without enforcing this explicitly in the optimization problem. If
the underlying Hüsler–Reiss distribution is $\text{EMTP}_{2}$ and an extremal
graphical model on $\mathcal{G}=(\mathcal{V},\mathcal{E})$, then the estimated
graph $\mathcal{G}_{+}=(\mathcal{V},\hat{\mathcal{E}}_{+})$ corresponding to
$\hat{\Theta}_{+}$ is asymptotically a supergraph of $\mathcal{G}$, that is,
$\mathbb{P}(\mathcal{E}\subseteq\hat{\mathcal{E}}_{+})\to 1$ as the sample
size $n\to\infty$.

A sufficient condition for extremal association for Hüsler–Reiss distributions
is the metric property, which requires the triangle inequalities

$\Gamma_{ij}\leq\Gamma_{ik}+\Gamma_{jk}$

to hold for any triple $i,j,k\in\mathcal{V}$. An $\text{EMTP}_{2}$
Hüsler–Reiss distribution always satisfies this property, while the inverse is
only true for ${d}=3$. While $\text{EMTP}_{2}$ is by construction a global
property, it can be a too strong assumption for applications with localized
positive dependence structure. For a given Hüsler–Reiss graphical model on
$\mathcal{G}$, Röttger and Schmitz (2023) therefore propose a model with local
metric property where the triangle inequalities are only required for all
triples $i,j,k\in C$, for all cliques $C$ of the graph $\mathcal{G}$. They
estimate this model based on mixed convex exponential families (Lauritzen and
Zwiernik, 2022) in a two-step procedure.

####  5.2.4 Score matching

Score matching is an alternative to likelihood inference that proposes to
minimize the Fisher information distance (Hyvärinen, 2005; Drton and Maathuis,
2017). The approach was originally proposed because its loss function does not
depend on the normalizing constant of the underlying distribution, but also
became popular for multivariate Gaussian models for its attractive
computational properties.

In graphical extremes, Lederer and Oesting (2023) propose a computationally
efficient score matching estimator for a family of distributions that
generalizes the Hüsler–Reiss model. Specifically, they introduce a class of
pseudo-densities that is parameterized by a superset of all Hüsler–Reiss
precision matrices and a nuisance location parameter. This function class
contains all Hüsler–Reiss multivariate Pareto densities, but also improper and
proper (though not multivariate Pareto) densities. Since the score matching
objective is well-defined even for non-integrable densities, they are able to
define an $L^{1}$ penalized score matching estimator in this class. A main
contribution in Lederer and Oesting (2023) is in providing an estimator of the
Hüsler–Reiss precision matrix that, despite not always being a valid
Hüsler–Reiss precision matrix itself, enjoys sharp concentration guarantees in
high dimensions. Nevertheless, the regularization scheme inherently produces
sparsity in the estimated precision matrix.

###  5.3 Structure learning

####  5.3.1 An extremal graphical lasso

The representation of conditional independence as zeros in the precision
matrix and the surrogate log-likelihood in (24) suggest the use of $L^{1}$
penalization methods to enforce graph sparsity. The Gaussian graphical lasso
(Yuan and Lin, 2007; Banerjee et al., 2008; Friedman et al., 2008) could then
be adapted to estimate sparse Hüsler–Reiss graphical models, giving rise to an
extremal graphical lasso problem of the form

$\mathop{\rm
minimize}_{\Theta\in\mathcal{S}^{d}_{\mathbf{1},+}}-\ell(\Theta;\overline{\Gamma})+\lambda\|\Theta\|_{1,\rm{off}},$
(27)

where $\|\Theta\|_{1,\rm{off}}:=\mathop{\sum\sum}_{i\neq j}|\Theta_{ij}|$ and
$\lambda\geq 0$ is a penalty parameter that controls the level of sparsity in
the solution. The matrix $\overline{\Gamma}$ can be a general estimator of
$\Gamma$ but will typically be chosen as the joint empirical variogram
$\hat{\Gamma}$. The Gaussian graphical lasso has implementations that are fast
and stable in high dimensions, and simultaneously learns a graph structure and
yields a model estimate on the learned graph. An estimator obtained by solving
(27) would share this property. However, an issue with this optimization
problem is the positive semi-definiteness of the Hüsler–Reiss precision
matrix. Indeed, the interaction between the $L^{1}$ penalization and
constraint of zero row sums in $\Theta$ has been shown to lead to dense
estimated models where every parameter is estimated to be nonzero (Ying et
al., 2020, 2021). We describe two approaches that circumvent this issue, and
which can be regarded as approximations of the ideal extremal graphical lasso
problem in (27).

####  5.3.2 Graph learning through majority voting

In Engelke et al. (2022), the authors consider the matrices
$\Theta^{(1)},\dots,\Theta^{(d)}\in\mathbb{R}^{(d-1)\times(d-1)}$, indexed by
$\mathcal{V}\setminus\\{m\\}$, where $\Theta^{(m)}$ is obtained by removing
the $m$th row and column of $\Theta$; these are the inverses of the covariance
matrices in (13). The authors leverage the fact that for each $m$,
$\Theta^{(m)}$ parameterizes the Hüsler–Reiss model and encodes most of the
extremal graphical structure in its sparsity pattern: by the properties of
$\Theta$, for any $i,j\neq m$,

${i\perp_{\Lambda}j\mid\mathcal{V}\setminus\\{i,j\\}}\quad\Leftrightarrow\quad\Theta_{ij}^{(m)}=0;$

see (Engelke and Hitz, 2020, Section 4.3). Conditional independencies
involving node $m$ are instead encoded as zero row sums in $\Theta^{(m)}$.
These matrices are usual positive definite precision matrices, and thus, their
support can be estimated by the Gaussian graphical lasso, neighborhood
selection (Meinshausen and Bühlmann, 2006) or similar algorithms. The
estimated support of $\Theta^{(m)}$ offers information on the presence of any
possible edge not connected to node $m$ in the graph $\mathcal{G}$. By
aggregating the estimated supports of $\Theta^{(1)},\dots,\Theta^{(d)}$
through the majority voting meta-algorithm EGlearn, Engelke et al. (2022)
obtain a full estimated extremal graph $\hat{\mathcal{G}}$. The method is
moreover shown to be model selection consistent (sparsistent) in exponentially
growing dimension under usual conditions. EGlearn is a pure structure learning
method. In combination with the matrix completion of the empirical variogram
$\hat{\Gamma}$ on the estimated graph $\hat{\mathcal{G}}$ in a second step,
the method provides a consistent and sparsistent estimator of the Hüsler–Reiss
parameter matrix (Hentschel et al., 2022).

####  5.3.3 Graph learning through parameter shift

In Wan and Zhou (2023), it is proposed to estimate the shifted parameter
matrix $\Theta^{*}:=\Theta+c\mathbf{1}\mathbf{1}^{\top}$, for an explicit
positive constant $c$. As opposed to $\Theta$, $\Theta^{*}$ is a precision
matrix of full rank and its inverse $\Sigma^{*}$ can be estimated from data.
The authors then solve the graphical lasso type optimization problem

$\mathop{\rm minimize}_{\Theta^{*}}-\log\det(\Theta^{*})+{\rm
tr}(\Theta^{*}\hat{\Sigma}^{*})+\lambda\|\Theta^{*}-c\mathbf{1}\mathbf{1}^{\top}\|_{1,\rm{off}},$

where $\Theta^{*}$ ranges over the symmetric positive definite matrices and
where $\hat{\Sigma}^{*}$ is an estimator of $\Sigma^{*}$. The idea is that the
first two terms in this objective function can be seen as a negative Gaussian
pseudolikelihood, which admits a minimum close to $\Theta^{*}$ if the
estimator $\hat{\Sigma}^{*}$ is sufficiently precise. The nonstandard penalty
term enforces sparsity in $\Theta$ rather than in $\Theta^{*}$. The method is
fast and has the advantage of simultaneously learning a graph and producing an
estimator of $\Theta$. However, the estimator is not guaranteed to be a valid
Hüsler–Reiss precision matrix. Nevertheless, partially relying on
concentration results proved in Engelke et al. (2022), Wan and Zhou (2023)
establish consistency in terms of both model selection and parameter
estimation.

Note that both this method and EGlearn estimate the sparsity pattern of
surrogate, full rank precision matrices. By doing so, they avoid the ill
behavior of positive semi-definite $L^{1}$ penalized maximum likelihood
estimation. They sacrifice, however, the special structure of the zero row
sums of Hüsler–Reiss precision matrices. For this reason, they share an
important shortcoming: neither is guaranteed to produce a connected estimated
graph, despite the fact Hüsler–Reiss graphical models are by definition
connected. Solving the ideal extremal graphical lasso problem in (27) would
remedy this issue, in addition to producing a valid parameter estimate
$\hat{\Theta}\in\mathcal{S}^{d}_{\mathbf{1},+}$.

##  6 Application

###  6.1 Dataset

Applications of extremal graphical models range from modeling river flows for
flood risk assessment (Asadi et al., 2015; Engelke and Hitz, 2020; Asenova et
al., 2021; Wan and Zhou, 2023), over financial applications (Klüppelberg and
Krali, 2021; Engelke and Volgushev, 2022) to food nutrition (Buck and
Klüppelberg, 2021; Krali et al., 2023).

In this section we illustrate the methods presented in Sections 5 and 4 by
applying them to the flight delay dataset introduced in Hentschel et al.
(2022), which is available from the R package graphicalExtremes (Engelke et
al., 2022). The goal is to qualitatively compare the methods and to illustrate
their specific properties. The dataset contains daily accumulated (positive)
flight delays in minutes at major US airports from 2005 to 2020. We consider
the airports in the “Texas Cluster” obtained in Hentschel et al. (2022) as the
result of a $k$-medoids clustering algorithm. The resulting dataset consists
of $n=$ 3603 observations at $d=$ 29 airports. It can be obtained as follows:

library(graphicalExtremes)delays <\- getFlightDelayData(airportFilter =
’tcCluster’, dateFilter = ’tcAll’)The airports and the graph that contains an
edge whenever there is at least one monthly flight between two airports is
shown in Figure 6(a). We split the data into a training set (2005-01-01 to
2010-12-31) and a test set (2011-01-01 to 2020-12-31). The former is used to
perform structure learning and parameter estimation with the different methods
and possibly different hyperparameters. The latter is used to compare the
fitted models in terms of the test likelihoods and other summary statistics.

###  6.2 Methods

Throughout, we assume a Hüsler–Reiss distribution to model the extremal
dependence structure of this dataset. We choose the empirical extremal
variogram $\hat{\Gamma}$ defined in (18) as input for all methods that we
apply. Recall that the marginal distributions are empirically normalized to
standard Pareto distributions; see Section 4.2. As a threshold we use the
empirical marginal $p=1-k/n=95\%$ quantiles so that the effective sample size
in each direction is $k=88$ in the training set. We use the Hüsler–Reiss model
with parameter matrix given by the empirical variogram as a first benchmark
and note that it corresponds to a fully connected extremal graph. To obtain
two other benchmarks with non-trivial graph structures, we choose the “Flight
Graph” obtained by considering all connections with at least one monthly
flight, and a randomly generated graph structure. For both graphs, we fit a
Hüsler–Reiss model by completing the empirical variogram using the matrix
completion in Section 5.2.

The “EGlearn” method in Section 5.3.2 and the extremal minimum spanning tree
“EMST” in Section 4.3 only provide an estimator for the underlying extremal
graphical structure $\mathcal{G}$, but not for the parameter matrix $\Gamma$.
In this application, the estimated parameter matrix of the “Parameter Shift”
method in Section 5.3.3 does not return a valid $\Gamma$ matrix, and we
therefore only use the estimated graph. In the case of “EGlearn”, the
algorithm used to estimate the support of the precision matrices
$\Theta^{(m)}$ is neighborhood selection with common penalty parameter $\rho$
(Engelke et al., 2022). For these three methods, we obtain a full model by
completing the empirical variogram on the estimated graphs as explained in
Section 5.2.

The “$\text{EMTP}_{2}$” method in Section 5.2.3 and the “Score Matching”
method in Section 5.2.4 directly yield estimates of the parameter matrix
$\Gamma$ and a graph structure $\mathcal{G}$. The method “Colored Graph” fits
an RVAR model where we choose the “$\text{EMTP}_{2}$” estimate as the
underlying graph and use a clustering algorithm to assign edge colorings.

Some of above methods are available in the R package graphicalExtremes:

p <\- 0.95G_hat <\- emp_vario(delays, p = p)models_eglearn <\- eglearn(delays,
p)model_emst <\- emst(delays, p)G_empt2 <\- emtp2(G_hat)

###  6.3 Results

Figure 6 shows the graph structures estimated by the different methods and
Table 1 contains the corresponding edge counts and log-likelihoods on the test
data. Figure 8 compares the variogram entries implied by the fitted models to
the corresponding entries of the empirical variograms on the test data.

By construction, the sparsest estimated graph is the “EMST” estimate. On this
data, the test likelihood and the consistent underestimation of the non-
adjacent empirical variograms indicate that a denser graph is necessary to
capture the full dependence structure. On the other side of the spectrum is
the full empirical variogram corresponding to a dense, fully connected graph.
The test likelihood suggest that this model is not sparse enough and overfits
to the training data. The random graph with an ad-hoc number of 100 edges is
in between the two methods in terms of test likelihood. The additional edges
compared to the tree allow for more flexible modeling, but since the graph is
not learned in a data-driven way, it is not competitive with the other
methods. The “Flight Graph” uses a similarly sparse graph, but the edges are
chosen in a more reasonable way according to the existence of flight
connections. This results in a better test likelihood than the previous
methods.

The remaining methods, which all enforce sparsity in the extremal graph in an
automatic, data-driven way, perform significantly better on this data set. The
methods “Score Matching”, “EGlearn”, and “Parameter Shift” each have a penalty
parameter, and “Colored Graph” has the number of colors as hyperparameter.
These methods were fitted for a range of different tuning parameters and the
best-performing model was chosen by considering the log-likelihood on the test
data, as shown in Figure 7. Furthermore, Figure 9 shows the convergence of
entries in $\Theta$ to zero as the penalization parameter is increased in the
three structure learning methods, illustrating their similarity to the
classical graphical lasso. Overall, the “Colored Graph” seems to have a
slightly better test likelihood than the other methods on this data set.
Interestingly, this is not reflected in the bivariate summary statistics in
Figure 8, which would hint at a slight overestimation of extremal dependence.

| Edges | Log-likelihood
---|---|---
Flight Graph | 129 | -9.23
Full Variogram | 406 | -293.68
Random Graph | 100 | -713.79
EMST | 28 | -6435.20
EGlearn | 101 | 1123.45
$\text{EMTP}_{2}$ | 126 | 1212.97
Parameter Shift | 142 | 1192.02
Colored Graph | 126 | 1341.17
Score Matching | 173 | 1144.97
Table 1:  Edge counts and test log-likelihoods for the different methods.

(a) Flight Graph

(b) Full Variogram

(c) Random Graph

(d) EMST

(e) EGlearn

(f) $\text{EMTP}_{2}$

(g) Parameter Shift

(h) Colored Graph

(i) Score Matching

Figure 6:  Graph structures estimated by the different methods. In (6(h)),
edge colors correspond to variogram entries with the same value.

(a)

(b)

(c)

(d)

Figure 7: Test log-likelihood for different tuning parameters for (7(a))
“EGlearn”, (7(b)) “Parameter Shift”, (7(c)) “Colored Graph”, and (7(d)) “Score
Matching”. Methods without tuning parameters are indicated by horizontal lines
(“Random Graph” dotted, “Full Variogram” short-dashed, “Flow Graph” solid,
“$\text{EMTP}_{2}$” long-dashed), “EMST” is below the axis and not shown.

##  Appendix A $d$-separation

Let $\mathcal{G}=(\mathcal{V},\mathcal{E})$ be a DAG. The set of parents
$\text{pa}(i)$ of node $i$ is defined as all nodes $j$ such that
$(j,i)\in\mathcal{E}$. The set of descendants $\text{des}(i)$ of some node
$i\in\mathcal{V}$ contains all nodes $j\in\mathcal{V}$ for which the graph
contains a directed path $(i,\ldots,j)\in\mathcal{E}$.

######  Definition A.1 (Pearl’s $d$-separation).

Let $\mathcal{G}=(\mathcal{V},\mathcal{E})$ be a directed acyclic graph. A
path $(i_{1},\ldots,i_{m})$ is blocked by a set $C$ when it contains a node
$i_{k}$ such that one of the following holds:

  * •

$i_{k}\in C$ and $i_{k-1}\rightarrow i_{k}\rightarrow i_{k+1}$ or
$i_{k-1}\leftarrow i_{k}\leftarrow i_{k+1}$ or $i_{k-1}\leftarrow
i_{k}\rightarrow i_{k+1}$,

  * •

$i_{k}\not\in C$ and $\text{des}(i_{k})\cap S=\emptyset$ and
$i_{k-1}\rightarrow i_{k}\leftarrow i_{k+1}$.

Disjoint sets $A,B\subseteq\mathcal{V}$ are called $d$-separated by a disjoint
set $C$ when every path between nodes in $A$ and $B$ is blocked by $C$.

Note that there are equivalent ways to define $d$-separation that rely on the
notion of collider. A collider is a vertex $i\in\mathcal{V}$ that has at least
two parents. The moralized graph of $\mathcal{G}$ is defined by adding an
edge, if not already present, between any two common parents of every collider
(the orientation of that new edge is irrelevant). The undirected skeleton of
this moralized graph is then obtained by remove the orientation of every edge.
It can be shown that $d$-separation between $A$ and $B$ by $C$ is equivalent
to graph separation (in the undirected sense) in the skeleton of the moralized
graph of $\mathcal{G}$.

##  Appendix B Additional Figures

(a) Flight Graph

(b) Full Variogram

(c) Random Graph

(d) EMST

(e) EGlearn

(f) $\text{EMTP}_{2}$

(g) Parameter Shift

(h) Colored Graph

(i) Score Matching

Figure 8:  Extremal Correlation based on the empirical and fitted extremal
variogram. In (8(h)), colored points correspond to edges of the same color.

(a) EGlearn

(b) Parameter Shift

(c) Score Matching

Figure 9:  Convergence of Entries in $\Theta$ to zero for penalization based
methods. Only those entries that eventually converge to zero are shown.

## References

  * Améndola et al. (2021) Améndola, C., B. Hollering, S. Sullivant, and N. Tran (2021).  Markov equivalence of max-linear Bayesian networks.  In Uncertainty in Artificial Intelligence, pp. 1746–1755. PMLR.
  * Améndola et al. (2022) Améndola, C., C. Klüppelberg, S. Lauritzen, and N. M. Tran (2022).  Conditional independence in max-linear Bayesian networks.  Ann. Appl. Probab. 32(1), 1–45.
  * Asadi et al. (2015) Asadi, P., A. C. Davison, and S. Engelke (2015).  Extremes on river networks.  Ann. Appl. Statist. 9(4), 2023–2050.
  * Asenova et al. (2021) Asenova, S., G. Mazo, and J. Segers (2021).  Inference on extremal dependence in the domain of attraction of a structured Hüsler–Reiss distribution motivated by a Markov tree with latent variables.  Extremes 24, 461–500.
  * Asenova and Segers (2023) Asenova, S. and J. Segers (2023).  Extremes of Markov random fields on block graphs: max-stable limits and structured Hüsler-Reiss distributions.  Extremes 26(3), 433–468.
  * Banerjee et al. (2008) Banerjee, O., L. El Ghaoui, and A. d’Aspremont (2008).  Model selection through sparse maximum likelihood estimation for multivariate Gaussian or binary data.  J. Mach. Learn. Res. 9, 485–516.
  * Beirlant et al. (2006) Beirlant, J., Y. Goegebeur, J. Segers, and J. L. Teugels (2006).  Statistics of extremes: theory and applications.  John Wiley & Sons.
  * Bortot and Coles (2003) Bortot, P. and S. Coles (2003).  Extremes of Markov chains with tail switching potential.  J. R. Stat. Soc. Ser. B Stat. Methodol. 65(4), 851–867.
  * Brown and Resnick (1977) Brown, B. M. and S. I. Resnick (1977).  Extreme values of independent stochastic processes.  J. Appl. Probab. 14, 732–739.
  * Buck and Klüppelberg (2021) Buck, J. and C. Klüppelberg (2021).  Recursive max-linear models with propagating noise.  Electron. J. Stat. 15(2), 4770–4822.
  * Capitanio et al. (2003) Capitanio, A., A. Azzalini, and E. Stanghellini (2003).  Graphical models for skew-normal variates.  Scand. J. Statist. 30(1), 129–144.
  * Casey and Papastathopoulos (2023) Casey, A. and I. Papastathopoulos (2023).  Decomposable tail graphical models.  Available from https://arxiv.org/abs/2302.05182.
  * Chiapino et al. (2019) Chiapino, M., A. Sabourin, and J. Segers (2019).  Identifying groups of variables with the potential of being large simultaneously.  Extremes 22, 193–222.
  * Coles and Tawn (1991) Coles, S. G. and J. A. Tawn (1991).  Modelling extreme multivariate events.  J. Roy. Statist. Soc. Ser. B 53(2), 377–392.
  * Darroch et al. (1980) Darroch, J. N., S. L. Lauritzen, and T. P. Speed (1980).  Markov fields and log-linear interaction models for contingency tables.  Ann. Statist. 8(3), 522–539.
  * de Fondeville and Davison (2018) de Fondeville, R. and A. C. Davison (2018).  High-dimensional peaks-over-threshold inference.  Biometrika 105, 575–592.
  * Dieker and Mikosch (2015) Dieker, A. B. and T. Mikosch (2015).  Exact simulation of Brown-Resnick random fields at a finite number of locations.  Extremes 18(2), 301–314.
  * Dombry et al. (2016) Dombry, C., S. Engelke, and M. Oesting (2016).  Exact simulation of max-stable processes.  Biometrika 103, 303–317.
  * Drton and Maathuis (2017) Drton, M. and M. H. Maathuis (2017).  Structure learning in graphical modeling.  Annu. Rev. Stat. Appl. 4, 365–393.
  * Einmahl et al. (2018) Einmahl, J. H. J., A. Kiriliouk, and J. Segers (2018).  A continuous updating weighted least squares estimator of tail dependence in high dimensions.  Extremes 21(2), 205–233.
  * Engelke et al. (2019) Engelke, S., R. de Fondeville, and M. Oesting (2019).  Extremal behaviour of aggregated data with an application to downscaling.  Biometrika 106(1), 127–144.
  * Engelke and Hitz (2020) Engelke, S. and A. S. Hitz (2020).  Graphical models for extremes (with discussion).  J. R. Stat. Soc. Ser. B Stat. Methodol 82(4), 871–932.
  * Engelke et al. (2022) Engelke, S., A. S. Hitz, N. Gnecco, and M. Hentschel (2022).  graphicalExtremes: Statistical Methodology for Graphical Extreme Value Models.  Available from https://github.com/sebastian-engelke/graphicalExtremes, R package version 0.1.0.9000.
  * Engelke and Ivanovs (2021) Engelke, S. and J. Ivanovs (2021).  Sparse structures for multivariate extremes.  Annu. Rev. Stat. Appl. 8, 241–270.
  * Engelke et al. (2022) Engelke, S., J. Ivanovs, and K. Strokorb (2022).  Graphical models for infinite measures with applications to extremes and Lévy processes.  Available from https://arxiv.org/abs/2211.15769.
  * Engelke et al. (2022) Engelke, S., M. Lalancette, and S. Volgushev (2022).  Learning extremal graphical structures in high dimensions.  Available from https://arxiv.org/abs/2111.00840.
  * Engelke et al. (2015) Engelke, S., A. Malinowski, Z. Kabluchko, and M. Schlather (2015, may).  Estimation of Hüsler-Reiss distributions and Brown-Resnick processes.  J. R. Stat. Soc. Ser. B Stat. Methodol 77(1), 239–265.
  * Engelke et al. (2014) Engelke, S., A. Malinowski, M. Oesting, and M. Schlather (2014).  Statistical inference for max-stable processes by conditioning on extreme events.  Adv. in Appl. Probab. 46(2), 478–495.
  * Engelke and Volgushev (2022) Engelke, S. and S. Volgushev (2022).  Structure learning for extremal tree models.  J. R. Stat. Soc. Ser. B. Stat. Methodol. 84(5), 2055–2087.
  * Fiebig et al. (2017) Fiebig, U.-R., K. Strokorb, and M. Schlather (2017).  The realization problem for tail correlation functions.  Extremes 20(1), 121–168.
  * Friedman et al. (2008) Friedman, J., T. Hastie, and R. Tibshirani (2008, 12).  Sparse inverse covariance estimation with the graphical lasso.  Biostatistics 9(3), 432–441.
  * Geiger et al. (1990) Geiger, D., T. Verma, and J. Pearl (1990).  Identifying independence in Bayesian networks.  Networks 20(5), 507–534.
  * Gissibl and Klüppelberg (2018) Gissibl, N. and C. Klüppelberg (2018).  Max-linear models on directed acyclic graphs.  Bernoulli 24(4A), 2693–2720.
  * Gissibl et al. (2021) Gissibl, N., C. Klüppelberg, and S. Lauritzen (2021).  Identifiability and estimation of recursive max-linear models.  Scand. J. Stat. 48(1), 188–211.
  * Gnecco et al. (2021) Gnecco, N., N. Meinshausen, J. Peters, and S. Engelke (2021).  Causal discovery in heavy-tailed models.  Ann. Statist. 49(3), 1755 – 1778.
  * Gudendorf and Segers (2010) Gudendorf, G. and J. Segers (2010).  Extreme-value copulas.  In Copula Theory and Its Applications: Proceedings of the Workshop Held in Warsaw, 25-26 September 2009, pp. 127–145. Springer.
  * Gumbel (1960) Gumbel, E. J. (1960).  Distributions des valeurs extrêmes en plusieurs dimensions.  Publ. Inst. Statist. Univ. Paris 9, 171–173.
  * Heffernan and Tawn (2004) Heffernan, J. E. and J. A. Tawn (2004).  A conditional approach for multivariate extreme values.  J. R. Stat. Soc. Ser. B Stat. Methodol 66(3), 497–546.
  * Hentschel et al. (2022) Hentschel, M., S. Engelke, and J. Segers (2022).  Statistical inference for hüsler-reiss graphical models through matrix completions.
  * Hüsler and Reiss (1989) Hüsler, J. and R.-D. Reiss (1989, February).  Maxima of normal random vectors: Between independence and complete dependence.  Statist. Prob. Letters 7(4), 283–286.
  * Hyvärinen (2005) Hyvärinen, A. (2005).  Estimation of non-normalized statistical models by score matching.  J. Mach. Learn. Res. 6, 695–709.
  * Inouye et al. (2016) Inouye, D., P. Ravikumar, and I. Dhillon (2016).  Square root graphical models: Multivariate generalizations of univariate exponential families that permit positive dependencies.  In International conference on machine learning, pp. 2445–2453. PMLR.
  * Kabluchko et al. (2009) Kabluchko, Z., M. Schlather, and L. de Haan (2009).  Stationary max-stable fields associated to negative definite functions.  Ann. Probab. 37(5), 2042 – 2065.
  * Kiriliouk et al. (2023) Kiriliouk, A., J. Lee, and J. Segers (2023).  X-vine models for multivariate extremes.
  * Klein et al. (2020) Klein, N., J. Orellana, S. L. Brincat, E. K. Miller, and R. E. Kass (2020).  Torus graphs for multivariate phase coupling analysis.  Ann. Appl. Stat. 14(2), 635–660.
  * Klüppelberg and Krali (2021) Klüppelberg, C. and M. Krali (2021).  Estimating an extreme Bayesian network via scalings.  J. Multivariate Anal. 181, Paper No. 104672, 23.
  * Krali et al. (2023) Krali, M., A. C. Davison, and C. Klüppelberg (2023).  Heavy-tailed max-linear structural equation models in networks with hidden nodes.  Available from https://arxiv.org/abs/2306.15356.
  * Kruskal (1956) Kruskal, Jr., J. B. (1956).  On the shortest spanning subtree of a graph and the traveling salesman problem.  Proc. Amer. Math. Soc. 7, 48–50.
  * Kulik and Soulier (2015) Kulik, R. and P. Soulier (2015).  Heavy tailed time series with extremal independence.  Extremes 18(2), 273–299.
  * Lalancette (2023) Lalancette, M. (2023).  On pairwise interaction multivariate Pareto models.  Stat 12, Paper No. e613, 10.
  * Lauritzen and Sadeghi (2018) Lauritzen, S. and K. Sadeghi (2018).  Unifying Markov properties for graphical models.  Ann. Statist. 46(5), 2251–2278.
  * Lauritzen and Zwiernik (2022) Lauritzen, S. and P. Zwiernik (2022).  Locally associated graphical models and mixed convex exponential families.  Ann. Statist. 50(5), 3009–3038.
  * Lauritzen (1996) Lauritzen, S. L. (1996).  Graphical models, Volume 17 of Oxford statistical science series.  Oxford: Clarendon Press.
  * Lederer and Oesting (2023) Lederer, J. and M. Oesting (2023).  Extremes in high dimensions: Methods and scalable algorithms.  Available from https://arxiv.org/abs/2303.04258.
  * Lee and Joe (2018) Lee, D. and H. Joe (2018).  Multivariate extreme value copulas with factor and tree dependence structures.  Extremes 21(1), 147–176.
  * Lin et al. (2016) Lin, L., M. Drton, and A. Shojaie (2016).  Estimation of high-dimensional graphical models using regularized score matching.  Electron. J. Stat. 10(1), 806–854.
  * Liu et al. (2009) Liu, H., J. Lafferty, and L. Wasserman (2009).  The nonparanormal: semiparametric estimation of high dimensional undirected graphs.  J. Mach. Learn. Res. 10, 2295–2328.
  * Liu et al. (2019) Liu, Z., J. H. Blanchet, A. B. Dieker, and T. Mikosch (2019).  On logarithmically optimal exact simulation of max-stable and related random fields on a compact set.  Bernoulli 25(4A), 2949–2981.
  * Maathuis et al. (2019) Maathuis, M., M. Drton, S. Lauritzen, and M. Wainwright (Eds.) (2019).  Handbook of graphical models.  Chapman & Hall/CRC Handbooks of Modern Statistical Methods. CRC Press, Boca Raton, FL.
  * Meinshausen and Bühlmann (2006) Meinshausen, N. and P. Bühlmann (2006).  High-dimensional graphs and variable selection with the lasso.  Ann. Stat. 34(3), 1436–1462.
  * Meyer and Wintenberger (2023) Meyer, N. and O. Wintenberger (2023).  Multivariate sparse clustering for extremes.  J. Am. Stat. Assoc. (forthcoming), 1–23.
  * Opitz (2013) Opitz, T. (2013).  Extremal $t$ processes: Elliptical domain of attraction and a spectral representation.  J. Multivariate Anal. 122, 409–413.
  * Papastathopoulos et al. (2023) Papastathopoulos, I., A. Casey, and J. A. Tawn (2023).  Hidden tail chains and recurrence equations for dependence parameters associated with extremes of higher-order Markov chains.  Available from https://arxiv.org/abs/1903.04059.
  * Papastathopoulos and Strokorb (2016) Papastathopoulos, I. and K. Strokorb (2016).  Conditional independence among max-stable laws.  Statistics & Probability Letters 108, 9–15.
  * Papastathopoulos et al. (2017) Papastathopoulos, I., K. Strokorb, J. A. Tawn, and A. Butler (2017).  Extreme events of Markov chains.  Adv. in Appl. Probab. 49(1), 134–161.
  * Pearl (1988) Pearl, J. (1988).  Probabilistic reasoning in intelligent systems: networks of plausible inference.  Morgan Kaufmann.
  * Pearl (2009) Pearl, J. (2009).  Causality (Second ed.).  Cambridge University Press, Cambridge.  Models, reasoning, and inference.
  * Perfekt (1994) Perfekt, R. (1994).  Extremal behaviour of stationary Markov chains with applications.  Ann. Appl. Probab. 4(2), 529–548.
  * Peters et al. (2017) Peters, J., D. Janzing, and B. Schölkopf (2017).  Elements of causal inference.  Adaptive Computation and Machine Learning. MIT Press, Cambridge, MA.  Foundations and learning algorithms.
  * Prim (1957) Prim, R. C. (1957).  Shortest connection networks and some generalizations.  Bell System Technical Journal 36, 1389–1401.
  * Ravikumar et al. (2010) Ravikumar, P., M. J. Wainwright, and J. D. Lafferty (2010).  High-dimensional Ising model selection using $\ell_{1}$-regularized logistic regression.  Ann. Statist. 38(3), 1287–1319.
  * Resnick (2008) Resnick, S. I. (2008).  Extreme Values, Regular Variation and Point Processes.  New York: Springer.
  * Röttger et al. (2023) Röttger, F., S. Engelke, and P. Zwiernik (2023).  Total positivity in multivariate extremes.  Ann. Statist. 51(3), 962–1004.
  * Röttger et al. (2023) Röttger, F., J. I. Coons, and A. Grosdos (2023).  Parametric and nonparametric symmetries in graphical models for extremes.  Available from https://arxiv.org/abs/2306.00703.
  * Röttger and Schmitz (2023) Röttger, F. and Q. Schmitz (2023).  On the local metric property in multivariate extremes.  Available from https://arxiv.org/abs/2212.10350.
  * Schlather (2002) Schlather, M. (2002).  Models for stationary max-stable random fields.  Extremes 5(1), 33–44.
  * Segers (2007) Segers, J. (2007).  Multivariate regular variation of heavy-tailed Markov chains.  Available from https://arxiv.org/abs/math/0701411.
  * Segers (2020) Segers, J. (2020).  One- versus multi-component regular variation and extremes of Markov trees.  Adv. in Appl. Probab. 52, 855–878.
  * Segers and Asenova (2022) Segers, J. and S. Asenova (2022).  Max-linear graphical models with heavy-tailed factors on trees of transitive tournaments.  Available from https://arxiv.org/abs/2209.14938.
  * Simpson et al. (2020) Simpson, E. S., J. L. Wadsworth, and J. A. Tawn (2020, 05).  Determining the dependence structure of multivariate extremes.  Biometrika 107(3), 513–532.
  * Smith (1992) Smith, R. L. (1992).  The extremal index for a Markov chain.  J. Appl. Probab. 29(1), 37–45.
  * Smith et al. (1997) Smith, R. L., J. A. Tawn, and S. G. Coles (1997).  Markov chain models for threshold exceedances.  Biometrika 84(2), 249–268.
  * Strokorb (2020) Strokorb, K. (2020).  Extremal independence old and new.  Available from https://arxiv.org/abs/2002.07808.
  * Thibaud et al. (2013) Thibaud, E., R. Mutzner, and A. C. Davison (2013).  Threshold modeling of extreme spatial rainfall.  Water Resources Research 49(8), 4633–4644.
  * Tran et al. (2021) Tran, N. M., J. Buck, and C. Klüppelberg (2021).  Estimating a directed tree for extremes.  Available from https://arxiv.org/abs/2102.06197.
  * Vettori et al. (2020) Vettori, S., R. Huser, J. Segers, and M. G. Genton (2020).  Bayesian model averaging over tree-based dependence structures for multivariate extremes.  J. Comput. Graph. Statist. 29(1), 174–190.
  * Wainwright and Jordan (2008) Wainwright, M. and M. Jordan (2008).  Graphical models, exponential families, and variational inference.  Foundations and Trends in Machine Learning 1(1–2), 1–305.
  * Wan and Zhou (2023) Wan, P. and C. Zhou (2023).  Graphical lasso for extremes.  Available from https://arxiv.org/abs/2307.15004.
  * Yang et al. (2015) Yang, E., P. Ravikumar, G. I. Allen, and Z. Liu (2015).  Graphical models via univariate exponential family distributions.  J. Mach. Learn. Res. 16, 3813–3847.
  * Yang et al. (2013) Yang, E., P. K. Ravikumar, G. I. Allen, and Z. Liu (2013).  On poisson graphical models.  Advances in neural information processing systems 26.
  * Ying et al. (2021) Ying, J., J. M. Cardoso, and D. Palomar (2021).  Minimax estimation of Laplacian constrained precision matrices.  In Int. Conf. Artif. Intell. Stat., pp. 3736–3744. PMLR.
  * Ying et al. (2020) Ying, J., J. V. d. M. Cardoso, and D. P. Palomar (2020).  Does the $\ell_{1}$-norm learn a sparse graph under Laplacian constrained graphical models?  Available from https://arxiv.org/abs/2006.14925.
  * Yuan and Lin (2007) Yuan, M. and Y. Lin (2007).  Model selection and estimation in the Gaussian graphical model.  Biometrika 94(1), 19–35.