# Partial Tail-Correlation Coefficient Applied to Extremal-Network Learning

Yan Gong, Peng Zhong, Thomas Opitz, Raphaël Huser

<figure><img src="image_1.png" /></figure>

## To cite this version:

Yan Gong, Peng Zhong, Thomas Opitz, Raphaël Huser. Partial Tail-Correlation Coefficient Applied to Extremal-Network Learning. Technometrics, 2024, pp.1-25. 10.1080/00401706.2024.2304334. hal-04643301

HAL Id: hal-04643301

https://hal.inrae.fr/hal-04643301v1

Submitted on 10 Jul 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Partial Tail-Correlation Coefficient Applied to Extremal-Network Learning

Yan Gong, Peng Zhong, Thomas Opitz & Raphaël Huser

**To cite this article:** Yan Gong, Peng Zhong, Thomas Opitz & Raphaël Huser (28 Feb 2024): Partial Tail-Correlation Coefficient Applied to Extremal-Network Learning, Technometrics, DOI: 10.1080/00401706.2024.2304334

**To link to this article:** https://doi.org/10.1080/00401706.2024.2304334

<figure><img src="image_3.png" /></figure>

© 2024 The Author(s). Published with license by Taylor & Francis Group, LLC.

<figure><img src="image_4.png" /></figure>

View supplementary material 🌐

<figure><img src="image_5.png" /></figure>

Published online: 28 Feb 2024.

<figure><img src="image_6.png" /></figure>

Submit your article to this journal 🌐

<figure><img src="image_7.png" /></figure>

Article views: 495

<figure><img src="image_8.png" /></figure>

View related articles 🌐

<figure><img src="image_9.png" /></figure>

View Crossmark data 🌐

CrossMark
<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /></figure>

OPEN ACCESS

<figure><img src="image_12.png" /></figure>

Check for updates

# Partial Tail-Correlation Coefficient Applied to Extremal-Network Learning

Yan Gong<sup>a,d</sup>, Peng Zhong<sup>b</sup>, Thomas Opitz<sup>c</sup>, and Raphaël Huser<sup>a,◎</sup>

<sup>a</sup>Statistics Program, Computer, Electrical and Mathematical Sciences and Engineering (CEMSE) Division, King Abdullah University of Science and Technology (KAUST), Thuwal, Saudi Arabia; <sup>b</sup>Data Science Hub, School of Mathematics and Statistics, University of New South Wales, Sydney, Australia; <sup>c</sup>Biostatistics and Spatial Processes, INRAE, Avignon, France; <sup>d</sup>Harvard T.H. Chan School of Public Health, Boston, MA

## ABSTRACT

We propose a novel extremal dependence measure called the partial tail-correlation coefficient (PTCC), in analogy to the partial correlation coefficient in classical multivariate analysis. The construction of our new coefficient is based on the framework of multivariate regular variation and transformed-linear algebra operations. We show how this coefficient allows identifying pairs of variables that have partially uncorrelated tails given some other variables in a random vector. Unlike other recently introduced conditional independence frameworks for extremes, our approach requires minimal modeling assumptions and can thus be used in exploratory analyses to learn the structure of extremal graphical models. Similarly to traditional Gaussian graphical models where edges correspond to the nonzero entries of the precision matrix, we can exploit classical inference methods for high-dimensional data, such as the graphical Lasso with Laplacian spectral constraints, to efficiently learn the extremal network structure via the PTCC. We apply our new method to study extreme risk networks in two different datasets (extreme river discharges and historical global currency exchange data) and show that we can extract interesting extremal structures with meaningful domain-specific interpretations.

## ARTICLE HISTORY

Received November 2022

Accepted January 2024

## KEYWORDS

Extreme event; Graphical Lasso; Multivariate regular variation; Network structure learning; Partial tail dependence

## 1. Introduction

Characterizing the extremal dependence of complex stochastic processes (e.g., in spatial, temporal, and spatio-temporal settings) is fundamental for both statistical modeling and applications, such as risk assessment in environmental and financial contexts. Important applications include the modeling of precipitation extremes (Huser and Davison 2014; Opitz et al. 2018; Bacro et al. 2020; Saunders, Stephenson, and Karoly 2021; Richards, Tawn, and Brown 2022), heatwaves (Winter and Tawn 2016; Zhong, Huser, and Opitz 2022), and air pollution (Vettori, Huser, and Genton 2019; Vettori et al. 2020), as well as financial risk assessment (Bassamboo, Juneja, and Zeevi 2008; Ferro and Stephenson 2011; Marcon, Padoan, and Antoniano-Villalobos 2016; Bekiros and Uddin 2017; Yan, Wu, and Zhang 2019; Gong and Huser 2022). Often, applications illustrate the benefits of methodological innovations, such as the application of the extremal dependence measure in Larsson and Resnick (2012), which is a key ingredient of the present work, in the analysis of financial data.

Models for extremal dependence traditionally rely on asymptotic frameworks, such as max-stable processes for block maxima or r-Pareto processes for threshold exceedances of a summary functional of the process over a high threshold (Davison and Huser 2015; Davison, Huser, and Thibaud 2019). Recently, more advanced models have been proposed

to further improve flexibility, especially toward modeling of asymptotically independent data with dependence vanishing at the most extreme levels, such as inverted max-stable processes (Wadsworth and Tawn 2012), max-mixture models (Ahmed et al. 2020), random scale mixtures of Gaussian processes (Opitz 2016; Huser, Opitz, and Thibaud 2017) or of more general processes (Wadsworth et al. 2017; Engelke, Opitz, and Wadsworth 2019; Huser and Wadsworth 2019), max-infinitely divisible processes (Bopp, Shaby, and Huser 2021), and conditional spatial extremes models (Wadsworth and Tawn 2022); for a comprehensive review, see Huser and Wadsworth (2022). Specifics of serial extremal dependence have been studied by Davis, Mikosch, and Zhao (2013), among others.

In the study of stochastic dependence structures, networks and graphs are natural tools to represent dependence relationships in multivariate data. Conditional independence, sparsity, and parsimonious representations are key concepts in graph-based approaches for random vectors. Recently, graph-based tools have also been developed for extremal dependence, where variants of conditional independence apply to variables not directly connected by the edges of the graph. For example, Huang et al. (2019) provide an exploratory tool, called the χ-network, for modeling extremal dependence, and they use it to analyze maximum precipitation during the hurricane season

<figure><img src="image_13.png" /></figure>
<figure><img src="image_14.png" /></figure>

in the United States (US) Gulf Coast and in surrounding areas. In their approach, however, the χ-network does not remove the effect of confounding variables, so it does not naturally lead to sparse extremal dependence representations. More recently, Engelke and Hitz (2020) introduce a notion of conditional independence adapted to multivariate Pareto distributions arising for limiting multivariate threshold exceedances, and they use it to develop parametric graphical models for extremes based on the Hüsler-Reiss distribution. Similarly, Gissibl (2018) and Klüppelberg and Krali (2021) propose max-linear constructions for modeling maxima on tree-like supports, and Tran, Buck, and Klüppelberg (2021) propose QTree, a simple and efficient algorithm to solve the “Latent River Problem” for the important case of extremes on trees. In the same vein, Engelke and Volgushev (2022) develop a data-driven methodology for learning the graphical structure in the setting of Engelke and Hitz (2020), whereas Röttger, Engelke, and Zwiernik (2023) further propose Hüsler-Reiss graphical models under the assumption of multivariate total positivity of order two (MTP₂), which allows estimating sparse graphical structures. Finally, Engelke and Ivanovs (2021) review the recent developments in sparse representations, dimension reduction approaches, and graphical models for extremes. Overall, existing graphical representations for extremes from the literature often rely on rather stringent asymptotically justified models, sometimes leading to issues when dealing with relatively high-dimensional problems or when specific graph structure assumptions (e.g., trees) are required. A recent exception is Engelke, Lalancette, and Volgushev (2021), who develop theory and methods for learning extremal graphical structures in high dimensions based on $L₁$ regularized optimization, though their methodology still assumes a parametric extremal dependence structure of Hüsler-Reiss type.

By contrast, rather than restricting ourselves to a strict parametric modeling framework, we adopt a more pragmatic and empirical approach. Specifically, our goal is to extend and enrich existing approaches by defining the new concept of *partial tail-correlation* as an extreme-value analog of the notion of partial correlation widely used in classical multivariate analysis, and by introducing a new coefficient that enables estimation of general extremal networks under minimal modeling assumptions. In the same way that correlation does not imply independence in general, our concept of *partial tail-uncorrelatedness* is a weaker assumption than *conditional tail independence*. However, we shall show that it still provides relevant insights into various forms of extremal dependence structures and helps in guiding modeling choices at an exploratory stage.

As a novel extremal dependence measure, we propose the *partial tail-correlation coefficient (PTCC)* as an equivalent of the partial correlation coefficient in the non-extreme setting. In the classical setting, the Pearson correlation coefficient between two random variables can give misleading interpretations when there are confounding variables that influence both variables, whereas the partial correlation coefficient measures the residual degree of association between two variables after the linear effects of a set of other variables have been removed. To compute the partial correlation between two variables of interest, we regress each of these variables onto the set of covariates given

by all the other variables in the multivariate random vector and then compute the correlation between the residuals from the two fitted linear regressions. In the Gaussian setting, a partial correlation of zero is equivalent to conditional independence between two variables (Lawrance 1976; Baba, Shibata, and Sibuya 2004), and the elements of the inverse of the covariance matrix (i.e., the *precision matrix*) of the full vector are known to characterize this conditional (in)dependence structure. In this article, we adopt a similar strategy to define the PTCC, namely by computing a suitable measure of tail dependence between residuals obtained by regressing variables using transformed-linear operations that do not alter tail properties. While classical linear regression only makes sense for Gaussian-like data, such transformed-linear operations can be used for tail regression of multivariate regularly-varying random vectors, which is a fundamental assumption characterizing asymptotically dependent extremes.

To be more precise, we here define the PTCC by building upon the framework of Cooley and Thibaud (2019), who developed a customized transformed-linear algebra on the positive orthant that preserves multivariate regular variation and is therefore well adapted to develop “linear” methods for joint extremes. Cooley and Thibaud (2019) used this framework for principal component analysis of extremes based on decompositions of the so-called tail pairwise dependence matrix (TPDM), which conveniently summarizes information about extremal dependence in random vectors and possesses appealing properties for such decompositions; Rohrbeck and Cooley (2023) applied this method of principal components for simulating extreme flood event sets. The TPDM can be thought of as an analogy of the classical covariance matrix but tailored for multivariate extremes. In follow-up work, Mhatre and Cooley (2020) then developed nonnegative regularly varying time series models with autoregressive moving average (ARMA) structure using the transformed-linear operations for time series extremes. For spatial extremes, Fix, Cooley, and Thibaud (2021) extended the simultaneous autoregressive (SAR) model under the transformed-linear framework and developed an estimation method to minimize the discrepancy between the TPDM of the fitted model and an empirically estimated TPDM. Furthermore, Lee and Cooley (2021) recently introduced transformed-linear prediction methods for extremes. In the aforementioned papers, the TPDM always plays a central role. Similar to covariances, the entries of the TPDM are tail dependence measures giving insights into the bivariate extremal dependence structure without removing the influence of other confounding variables. However, just as the covariance matrix does not reflect partial correlations, the TPDM does not directly inform us about partial associations among extremes. In this work, we fill this gap with our new proposed PTCC. Thanks to its definition in terms of best prediction using transformed-linear operations, we show that the PTCC inherits several appealing features of the classical partial correlation coefficient. In particular, the PTCC between two components $X_i$ and $X_k$ from a random vector $X$ is such that the $(i, k)$th entry of the inverse TPDM matrix of $X$ equals zero if and only if the corresponding PTCC for these two variables is also equal to zero. In other words, *partial tail-uncorrelatedness* can be conveniently read off from the zero elements of the inverse
<figure><img src="image_15.png" /></figure>

TPDM, similar to classical Gaussian graphical models. We exploit this property to define a new class of extremal graphical models based on the PTCC and then propose efficient inference methods to learn the extremal network structure from high-dimensional data based on state-of-the-art techniques from graph theory (e.g., the graphical Lasso with or without Laplacian spectral constraints). Our focus is on studying undirected graph structures, which is different from causal inference, where causal relationships can be encoded using directed graph edges.

The remainder of this article is organized as follows. In Section 2, we first review the necessary background on multivariate regular variation and transformed-linear algebra, as introduced in Cooley and Thibaud (2019). Then, we define the new PTCC and the related notion of partial tail-uncorrelatedness. In Section 3, we present methods for learning general extremal network structures from the PTCC in a high-dimensional data setting, and we discuss two particularly appealing approaches, namely the graphical Lasso and Laplacian spectral constraint-based methods. Section 4 presents a simulation study to illustrate the effectiveness of the above two graph inference methods, for various structured undirected graphs and sample sizes. In Section 5, we apply these new tools to explore the risk networks formed by river discharges observed at a collection of monitoring stations in the upper Danube basin, and by historical global currency exchange rate data from different historical periods, covering different economic cycles, the COVID-19 pandemic, and the 2022 military conflict in Ukraine.

## 2. Transformed-Linear Algebra for Multivariate Extremes

Before introducing the partial tail-correlation coefficient (PTCC) and the related notion of partial tail-uncorrelatedness, we first briefly review the multivariate regular variation framework, which is our main assumption for defining the PTCC. We also summarize the foundations of transformed-linear algebra, which allows defining an inner product on the finite-dimensional vector space spanned by all transformed-linear combinations of certain independent univariate regularly varying variables. Thanks to the inner product, we can define “linear” prediction in such vector spaces as the projection of a given variable onto the space spanned by other variables; and this in turn allows us to define the PTCC by analogy with the classical partial correlation coefficient. Importantly, results from the literature highlight that any TPDM can arise from such a transformed-linear combination in a finite-dimensional vector space with an appropriately chosen dimension. Consequently, as we shall see, the TPDM of a given random vector encapsulates all information we need to determine which variables from this random vector are partially tail-uncorrelated.

### 2.1. Regular Variation Framework and Transformed Linear Algebra

A random vector is multivariate regularly varying (Resnick 2007) (i.e., jointly heavy-tailed) if its joint tail decays like a power function. Precisely, we say that a p-dimensional random vector

$X \in \mathbb{R}_+^p = [0, \infty)^p$ with $p \in \mathbb{N}$ is regularly varying if there exists a sequence $b_n \to \infty$ such that

$$n \Pr(b_n^{-1}X \in \cdot) \xrightarrow{\nu} v_X(\cdot), \quad n \to \infty, \qquad (1)$$

where $\nu$ denotes vague convergence to the non-null limit measure $v_X$, a Radon measure defined on the space $[0, \infty]^p \setminus \{0\}$. This measure has the scaling property $r^\alpha v_X(rB) = v_X(B)$ for $r > 0$ and Borel sets $B \subset [0, \infty]^p \setminus \{0\}$, where $\alpha > 0$ controls the tail decay (with $1/\alpha$ commonly called the tail index). For this reason, the measure can be further decomposed into a radial measure and an angular measure $H_X$ on the unit sphere $\mathbb{S}_{p-1}^+ = \{w \in \mathbb{R}_+^p : \|w\|_2 = 1\}$, where $\| \cdot \|_2$ is the $L_2$ norm, such that $v_X(\{x \in [0, \infty]^p \setminus \{0\} : \|x\|_2 \ge r, x/\|x\|_2 \in B_H\}) = r^{-\alpha} \times H_X(B_H)$ for $r > 0$ and Borel subsets $B_H$ of $\mathbb{S}_{p-1}^+$. The normalizing sequence $b_n$ is not uniquely determined but must satisfy $b_n = L(n)n^{1/\alpha}$, where $L(n)$ is a slowly varying function (at infinity), that is, $L(n) > 0$ and $L(rn)/L(n) \to 1$ for any $r > 0$, as $n \to \infty$. We use the short-hand notation $X \in \mathbb{RV}_+^p(\alpha)$ for a regularly varying vector $X$ with tail index $1/\alpha$.

Cooley and Thibaud (2019) introduced the transformed-linear algebra framework to construct an inner-product space on an open set (the so-called target space) via a suitable transformation, where the distribution of the random vector $X$ has support within this set. Our use will mainly concern transformation toward the target space $\mathbb{R}_+^p$ from the space $\mathbb{R}^p$, but we first present the general approach. Let $t$ be a bijective transformation from $\mathbb{R}$ onto some open set $X \subset \mathbb{R}$, and let $t^{-1}$ be its inverse. For a p-dimensional vector $y \in \mathbb{R}^p$, we define $x = t(y) \in \mathbb{R}^p$ componentwise. Then, arithmetic operations among elements of the target space are carried out in the space $\mathbb{R}^p$ before transforming back to the target space. We define vector addition in $\mathbb{R}^p$ as $x_1 \oplus x_2 = t(t^{-1}(x_1) + t^{-1}(x_2))$, and scalar multiplication with a factor $a \in \mathbb{R}$ as $a \circ x = t(at^{-1}(x))$. The additive identity in $\mathbb{R}^p$ is set to $0_{\mathbb{R}^p} = t(0)$, and the additive inverse of $x \in \mathbb{R}^p$ is given as $\ominus x = t(-t^{-1}(x))$. A valid inner product between two elements $x_1 = (x_{1,1}, \dots, x_{1,p})^T$, $x_2 = (x_{2,1}, \dots, x_{2,p})^T \in \mathbb{R}^p$ from the target space is then obtained by applying the usual scalar product in $\mathbb{R}^p$, that is, we set $\langle x_1, x_2 \rangle = \sum_{j=1}^p t^{-1}(x_{1,j})t^{-1}(x_{2,j})$. To obtain an inner-product space on the positive orthant for which arithmetic operations preserve multivariate regular variation, thus, having an asymptotically negligible effect on large values, we follow Cooley and Thibaud (2019) and define the specific transformation $t: \mathbb{R} \mapsto (0, \infty)$ given by

$$t(y) = \log\{1 + \exp(y)\},$$

although there are other possibilities. We have $y/t(y) \to 1$ as $y \to \infty$, such that the upper tail behavior of a random vector $Y = t(X)$ is preserved through $t$. To prevent the lower tail from interfering with the upper tail (Cooley and Thibaud 2019), we also need $\exp(y)/t(y) \to 1$ as $y \to -\infty$, which is clearly satisfied here. The inverse transformation is $t^{-1}(x) = \log(\exp(x) - 1)$, $x > 0$. Algebraic operations done in the vector space induced by the above transformation $t$ are commonly called transformed-linear operations, and we can exploit this framework to extend classical linear algebra methods (e.g., principal component analysis) to the multivariate extremes setting,
<figure><img src="image_16.png" /></figure>

where heavy-tailed vectors and models are often conveniently expressed on the positive orthant, $\mathbb{R}_+^p$. Our main assumption of multivariate regular variation is useful for data that can potentially exhibit asymptotic tail dependence, and it implies homogeneity of the limit measure in (1), but it does not impose further parametric structural assumptions such as with Pareto models of Hüsler-Reiss type (Hüsler and Reiss 1989).

## 2.2. Inner-Product Space of Regularly Varying Random Variables

With transformed-linear operations, we can use a vector of iid random variables to construct new regularly varying random vectors on the positive orthant that possess tail dependence. Suppose that $Z = (Z_1, \dots, Z_q)^T > 0$ is a vector of $q \in \mathbb{N}$ iid regularly varying random variables with tail index $1/\alpha$, such that there exists a sequence $\{b_n\}$ that yields

$$n \Pr(Z_j > b_n z) \to z^{-\alpha}, \quad n \Pr\{Z_j \le \exp(-kb_n)\} \to 0,$$

for $k > 0$, and $j = 1, \dots, q$, where the first condition is equivalent to regular variation (1) in dimension $p = 1$. The random vector $Z$ with independent components has a limit measure of multivariate regular variation characterized by $v_Z\{[0, z]^C\} = \sum_{j=1}^q z_j^{-\alpha}$ for $z = (z_1, \dots, z_q)^T > 0$. Then, we can construct new regularly varying $p$-dimensional random vectors $X = (X_1, \dots, X_p)^T$ by exploiting transformed-linear operations, via a matrix product with a deterministic matrix $A = (a_1, \dots, a_q) \in \mathbb{R}_+^{p \times q}$, with columns $a_j = (a_{1j}, \dots, a_{pj})^T \in \mathbb{R}_+^p$; precisely, the $i$th component of the random vector $X$ may be expressed as follows:

$$X_i = \bigoplus_{j=1}^{q} a_{ij} \circ Z = (a_{i1} \circ Z_1) \oplus \dots \oplus (a_{iq} \circ Z_q), \quad (2)$$

for $i = 1, \dots, p$, and we write $X = A \circ Z \in \text{RV}_+^p(\alpha)$. This construction ensures that the multivariate regular variation property is preserved with the same index $\alpha$ (Cooley and Thibaud 2019, Corollary 1). Furthermore, we require $A$ to have a full row-rank. Based on the construction (2), it is possible to define a (different) inner-product space spanned by the random variables obtained by transformed-linear operations on $Z$, where some but not all of the components of $a_j$ are further allowed to be nonpositive. Following Lee and Cooley (2021), an inner product of $\langle X_i, X_k \rangle$ on the space spanned by all possible transformed-linear combinations of the elements of the random vector $X$ constructed as in (2), may be defined as follows:

$$\langle X_i, X_k \rangle = \sum_{j=1}^{q} a_{ij} a_{kj},$$

where $a_{ij}$ refers to the entry in row $i$ of the column $j$ of the matrix $A$ for $i \in \{1, \dots, p\}$ and $j \in \{1, \dots, q\}$, and the corresponding norm becomes $\|X\| = \sqrt{\langle X, X \rangle}$. The metric induced by the inner product is $d(X_i, X_k) = \|X_i \ominus X_k\| = [\sum_{j=1}^q (a_{ij} - a_{kj})^2]^{1/2}$, for $i, k = 1, \dots, p$.

## 2.3. Generality of the Framework

In practice, given a random vector $X$ for which we assume $X \in \text{RV}_+^p(\alpha)$, we will further assume that it allows for a stochastic representation as in (2). Since the constructions of type (2) form a dense subclass of the class of multivariate regularly varying vectors (if $q$ is not fixed but allowed to tend to infinity, that is, $q \to \infty$), this assumption is not restrictive; see Fougères, Mercadier, and Nolan (2013) and Cooley and Thibaud (2019), Proposition 4.

Thanks to their flexibility, representations akin to the transformed-linear random vectors in (2) have recently found widespread interest in statistical learning for extremes. The fundamental model structure used in the causal discovery framework for extremes developed by Gnecco et al. (2021) is essentially based on a variant of (2). In the setting of max-linear models, in particular the graphical models of Gissibl (2018), we can use (2) to construct random vectors $X$ possessing the same limit measure $v_X$ as the max-linear vectors. Finally, low-dimensional representations of extremal dependence in random vectors obtained through variants of the $k$-means algorithm can be shown to be equivalent to the extremal dependence induced by construction (2); see Janßen and Wan (2020).

## 2.4. Tail Pairwise Dependence Matrix

The tail pairwise dependence matrix (TPDM; Cooley and Thibaud 2019, sec. 4) is defined to summarize the pairwise extremal dependence of a regularly varying random vector using the second-order properties of its angular measure. Let $\alpha = 2$, which ensures desirable properties; in practice, this condition can be ensured through appropriate marginal pre-standardization of data (Cooley and Thibaud 2019). Then, the TPDM $\Sigma$ of $X \in \text{RV}_+^p(2)$ is defined as follows:

$$\Sigma = (\sigma_{ik})_{i,k=1,\dots,p}, \quad \text{with} \quad \sigma_{ik} := \int_{\mathbb{S}_{p-1}^+} w_i w_k dH_X(w),$$

where $H_X$ is the angular measure on $\mathbb{S}_{p-1}^+ = \{w \in \mathbb{R}_+^p : \|w\|_2 = 1\}$ as introduced in Section 2.1. The matrix $\Sigma$ is an extreme-value analog of the covariance matrix, and it has similar useful properties. It is positive semidefinite and completely positive, that is, there always exists a finite $p \times q$ matrix $A$ with nonnegative entries such that the TPDM can be factorized as $\Sigma = AA^T$ (Cooley and Thibaud 2019, Proposition 5). The matrix $A$ is not unique, in particular, if we do not impose nonnegative entries. Specifically, for random vectors $X$ obtained by the transformed-linear construction (2), the entries of $\Sigma$ correspond to the values of the inner product $\sigma_{ik} = \langle X_i, X_k \rangle$. In the following, we further assume that $\Sigma$ is positive definite, which guarantees the existence of the inverse matrix of the TPDM.

We emphasize that the special case where $\sigma_{ik} = 0$ is equivalent to asymptotic tail independence of the components $X_i$ and $X_k$ (see Cooley and Thibaud 2019), meaning that the conditional exceedance probability $\Pr(F_{X_i}(X_i) > u | F_{X_k}(X_k) > u)$ tends to zero as $u \to 1$ (Sibuya 1960; Ledford and Tawn 1996).

By exploiting the property that the TPDM is completely positive, we can construct new transformed-linear random vectors that have the same TPDM as a given random vector $X$. Since we
<figure><img src="image_17.png" /></figure>

can always factorize the TPDM as $\Sigma = AA^T$ for some matrix $A$ of dimension $p \times q$, we can then use the construction (2) by multiplying $A$ with a (new) random vector $Z \in RV_+^p(2)$ of independent regularly varying random variables, and the resulting vector will still have TPDM $\Sigma$. It is also worth noting that as $q \to \infty$, the angular measure of the new random vector can be arbitrarily close to that of $X$ thanks to the denseness property of discrete angular measures. While $A$ is not unique, we note that it is very important that the inner product depends only on the entries of the matrix $\Sigma = AA^T$, such that the specific choice of $A$ does not matter.

An estimator of the TPDM was proposed by Cooley and Thibaud (2019). For an iid sequence of vectors $x_t$, with $t = 1, \dots, n_{\text{samp}}$, that is, independent realizations from a random vector $X \in RV_+^p(2)$, define

$$ \hat{\sigma}_{ik} = \hat{m} \int_{S_{p-1}^+} w_t w_k d\hat{N}_X(w) = \hat{m} n_{\text{ext}}^{-1} \sum_{t=1}^{n_{\text{samp}}} w_{ti} w_{tk} \mathbb{1}(r_t > r_0), \quad (3) $$

where $r_t = \|x_t\|_2$ and $w_t = x_t/r_t$ are the observed radial and angular components, respectively, $r_0$ is a high threshold for the radial component, $n_{\text{ext}} = \sum_{t=1}^{n_{\text{samp}}} \mathbb{1}(r_t > r_0)$ refers to the number of threshold exceedances, and the probability measure $N_X(\cdot) = m^{-1}H_X(\cdot)$ (with $m = H_X(S_{p-1}^+)$) is obtained by normalizing $H_X$. Moreover, $\hat{m}$ denotes an estimate of $H_X(S_{p-1}^+)$, and $\hat{N}_X$ is the empirical counterpart of $N_X$. We note that when the data are preprocessed to have a common scale (for instance, Fréchet(2)), we can set $m = p$ and there is no need to estimate the normalizing factor Cooley and Thibaud (2019, sec. 7.1). The estimator (3) was discussed by Larsson and Resnick (2012) in the bivariate case.

## 2.5. Partial Tail-Correlation Coefficient

In this section, we introduce our new measure, the partial tail-correlation coefficient (PTCC), which is analogous to the partial correlation coefficient but tailored to heavy-tailed random vectors.

Let $X = A \circ Z \in RV_+^p(2)$ be a $p$-dimensional vector constructed as in (2), with TPDM $\Sigma$. We write $X_{ik} = (X_i, X_k)^T$; $X_{\text{rest}}$ for the $(p-2)$-dimensional random vector obtained by removing the two components $X_i$ and $X_k$ from $X$; $A_{ik}$ for the matrix comprising the $i$th and $k$th rows of $A$; and $A_{\text{rest}}$ for the matrix $A$ without its $i$th and $k$th rows. Moreover, we define the $p$-dimensional random vector $X'$ by reordering the components of $X$ as $X' = (X_{ik}^T, X_{\text{rest}}^T)^T = (A_{ik}^T, A_{\text{rest}}^T)^T \circ Z$. It is straightforward to show that the best transformed-linear predictor of $X_{ik}$ given $X_{\text{rest}}$ can be obtained as $\hat{X}_{ik} = B \circ X_{\text{rest}} = (\hat{X}_{ik}^i, \hat{X}_{ik}^k)^T$, where $B = (b_1, b_2)^T$ is a $2 \times (p-2)$ matrix with $b_1, b_2 \in \mathbb{R}_+^{p-2}$ chosen such that $d(X_i, \hat{X}_{ik}^i)$ and $d(X_k, \hat{X}_{ik}^k)$ attain their minimum, respectively. Suppose that the TPDM of $X'$ is of the block-matrix form

$$ \Sigma_{X'} = \begin{bmatrix} \Sigma_{ik,ik} & \Sigma_{ik,\text{rest}} \\ \Sigma_{\text{rest},ik} & \Sigma_{\text{rest},\text{rest}} \end{bmatrix} $$

where $\Sigma_{ik,ik} \in \mathbb{R}_+^{2 \times 2}$ is the TPDM restricted to $X_{ik}$, that is, $\Sigma_{ik,ik} = \text{TPDM}(X_{ik}) = \begin{bmatrix} \Sigma_{ii} & \Sigma_{ik} \\ \Sigma_{ki} & \Sigma_{kk} \end{bmatrix}$; $\Sigma_{\text{rest},\text{rest}} \in \mathbb{R}_+^{(p-2) \times (p-2)}$

is the TPDM restricted to $X_{\text{rest}}$, and the off-diagonal block-matrix $\Sigma_{ik,\text{rest}} = \Sigma_{\text{rest},ik}^T = (\Sigma_{i,\text{rest}}^T, \Sigma_{k,\text{rest}}^T)^T \in \mathbb{R}_+^{2 \times (p-2)}$ is the cross-TPDM between $X_{ik}$ and $X_{\text{rest}}$. Then, based on the projection theorem for vector spaces with inner products, we have that

$$ B = \Sigma_{ik,\text{rest}} \Sigma_{\text{rest},\text{rest}}^{-1} $$

and the inner product matrix of the projection error $e = X_{ik} \ominus \hat{X}_{ik}$ can be expressed as

$$ \begin{bmatrix} \Sigma_{ii} - \Sigma_{i,\text{rest}} \Sigma_{\text{rest},\text{rest}}^{-1} \Sigma_{\text{rest},i} & \Sigma_{ik} - \Sigma_{i,\text{rest}} \Sigma_{\text{rest},\text{rest}}^{-1} \Sigma_{\text{rest},k} \\ \Sigma_{ki} - \Sigma_{k,\text{rest}} \Sigma_{\text{rest},\text{rest}}^{-1} \Sigma_{\text{rest},i} & \Sigma_{kk} - \Sigma_{k,\text{rest}} \Sigma_{\text{rest},\text{rest}}^{-1} \Sigma_{\text{rest},k} \end{bmatrix} = \Sigma_{ik,ik} - \Sigma_{ik,\text{rest}} \Sigma_{\text{rest},\text{rest}}^{-1} \Sigma_{\text{rest},ik}. \quad (4) $$

The diagonal elements of this 2-by-2 matrix are always nonnegative but its off-diagonal elements can be negative. Therefore, this matrix is in general not completely positive and therefore not a TPDM in general. The possibility of negative entries is a consequence of the structure of the inner product that does not allow excluding negative entries of the coefficient matrix $A$. In the following, we further assume that the diagonal elements in (4) are positive, such that the (degenerate) case where a component $X_i$ is an exact transformed-linear combination of the components different from $X_i$ and $X_k$ is excluded.

We can now define the partial tail-covariance and the partial tail-correlation coefficients.

**Definition 1.** The partial tail-covariance coefficient of two random variables $X_i$ and $X_k$ is defined as the off-diagonal coefficient in the inner product matrix of the bivariate residual vector $e$ shown in (4), that is, $\Sigma_{ik} - \Sigma_{i,\text{rest}} \Sigma_{\text{rest},\text{rest}}^{-1} \Sigma_{\text{rest},k}$, such that transformed-linear dependence with respect to all other random variables is removed.

**Definition 2.** The partial tail-correlation coefficient (PTCC) is defined as the standardized partial tail-covariance coefficient, that is

$$ \text{PTCC}_{ik} = \frac{\Sigma_{ik} - \Sigma_{i,\text{rest}} \Sigma_{\text{rest},\text{rest}}^{-1} \Sigma_{\text{rest},k}}{\sqrt{\Sigma_{ii} - \Sigma_{i,\text{rest}} \Sigma_{\text{rest},\text{rest}}^{-1} \Sigma_{\text{rest},i} \Sigma_{kk} - \Sigma_{k,\text{rest}} \Sigma_{\text{rest},\text{rest}}^{-1} \Sigma_{\text{rest},k}}} \in [-1, 1]. $$

From the way the PTCC is constructed, it is obvious that this new measure is only informative about the joint tail structure but is invariant to the properties of the bulk of the distribution (as it can be obtained from the angular measure), and that it is a scale-free copula-based metric, thus, invariant to monotone marginal transformations (this follows from the pre-standardization to a common scale with fixed $\alpha = 2$).

**Definition 3.** Let $X_{ik} = (X_i, X_k)^T$ and $X_{\text{rest}}$ be defined as above. Given $X_{\text{rest}}$, we say that $X_i$ and $X_k$ are partially tail-uncorrelated if the PTCC of $X_i$ and $X_k$ (given $X_{\text{rest}}$) is equal to zero, which is equivalent to $\Sigma_{ik} - \Sigma_{i,\text{rest}} \Sigma_{\text{rest},\text{rest}}^{-1} \Sigma_{\text{rest},k} = 0$.

The following proposition links tail-uncorrelatedness to the entries of the inverse of the TPDM of $X$.
<figure><img src="image_18.png" /></figure>

**Proposition 1.** Given the representation of the TPDM of $X'$ as a 3 × 3 block matrix as follows,

$$ \Sigma_{X'} = \text{TPDM}\{(X_i, X_k, \mathbf{X}_{\text{rest}}^T)^T\} \quad (5) \\ = \begin{bmatrix} \Sigma_{ii} & \Sigma_{ik} & \Sigma_{i,\text{rest}} \\ \Sigma_{ki} & \Sigma_{kk} & \Sigma_{k,\text{rest}} \\ \Sigma_{\text{rest},i} & \Sigma_{\text{rest},k} & \Sigma_{\text{rest},\text{rest}} \end{bmatrix}, $$

where the dimensions of submatrices are as above, the following two statements are equivalent:

$$ (1) \Sigma_{ik} - \Sigma_{i,\text{rest}} \Sigma_{\text{rest},\text{rest}}^{-1} \Sigma_{\text{rest},k} = 0, $$

$$ (2) (\Sigma^{-1})_{ik} = 0 $$

where $\Sigma$ is the TPDM of the original vector $X$.

Since $\Sigma_{X'}$ is a positive definite and therefore an invertible covariance matrix, this result is a direct consequence of the equivalence of statements (i') and (ii) of Proposition 1 of Speed and Kiiveri (1986). In the notation of Speed and Kiiveri (1986), we have the following index sets: $a = \{k, \text{rest}\}$, $b = \{i, \text{rest}\}$, and $ab = a \cap b = \text{rest}$. A direct consequence of Proposition 1 and Definitions 2 and 3 is that two variables $X_i$ and $X_k$ are partially tail-uncorrelated given $X_{\text{rest}}$ (i.e., $\text{PTCC}_{ik} = 0$) if and only if the $(i, k)$th entry of the inverse TPDM is zero, that is, $Q_{ik} = 0$ with $Q = \Sigma^{-1}$ where $\Sigma$ is the TPDM of the vector $X$.

# 3. Learning Extremal Networks for High-Dimensional Extremes

Using the PTCC, we can now explore the partial tail-correlation structure of multivariate random variables under the framework of multivariate regular variation. In this section, we define new graphical models to represent extremal dependence for extremes. Thanks to the transformed-linear framework exposed in the previous section, we can proceed as for classical graphical models by replacing the classical covariance matrix with the TPDM.

## 3.1. Graphical Models for Extremes

Let $G = (V, E)$ be a graph, where $V = \{1, \dots, p\}$ represents the node set and $E \subset V \times V$ the edge set. We call $G$ an undirected graph if for two nodes $i, k \in V$, the edge $(i, k)$ is in $E$ if and only if the edge $(k, i)$ is also in $E$. We show how to estimate graphical structures for extremes for any type of undirected graph in which we have no edge $(i, k)$ if and only if the variables $X_i$ and $X_k$ are partially tail-uncorrelated given all the other variables in the graph, which we write as

$$ X_i \perp_p X_k | \mathbf{X}_{\text{rest}}. $$

Our methods work for general undirected graph structures, including trees, decomposable graphs, and nondecomposable graphs, see the example illustrations in Figure 1. Note, however, that our general method cannot restrict the estimated graph to be of a specific type, such as a tree.

## 3.2. Sparse Representation of High-Dimensional Extreme Networks

For high-dimensional extremes with a relatively large number of components $p$ (e.g., up to tens and hundreds of variables), a graphical representation of the extremal dependence structure is desirable for reasons of parsimony and interpretability. We now introduce two efficient inference methods to learn extremal networks from high-dimensional data via the PTCC based on two state-of-the-art graphical methods: the extremal graphical Lasso, and the structured graph learning method via Laplacian spectral constraints (Friedman, Hastie, and Tibshirani 2008; Kumar et al. 2019), which were initially designed for learning Gaussian graphical models. These two methods both work efficiently in high-dimensional settings, and Ravikumar et al. (2011) and Pavez (2022) established their consistency under broad conditions. In our context, they return an estimate of the underlying extremal dependence with sparse structure, that is, with the cardinality of $E$ being of the same order as that of $V$.

### 3.2.1. Extremal Graphical Lasso

Given an empirical estimator $\hat{\Sigma}$ of the TPDM and a penalty parameter $\lambda \ge 0$, the optimization carried out by the extremal graphical Lasso method is expressed as follows (Friedman, Hastie, and Tibshirani 2008):

$$ \hat{Q}_{\lambda} = \arg \max_{\Theta \ge 0} \left[ \log \det \Theta - \text{tr}(\hat{\Sigma} \Theta) - \lambda \sum_{i \neq k} |\Theta_{ik}| \right] $$

where $\ge$ indicates positive-semidefiniteness and $\hat{Q}_{\lambda}$ is a $L_1$-regularized estimate of the precision matrix $Q = \Sigma^{-1}$. Note that thanks to the $L_1$ regularization, the estimate $\hat{Q}_{\lambda}$ will tend to be sparse (with exact zero entries) and thus contains information on the extremal graph structure. A larger $\lambda$ enforces a larger proportion of zeros in $\hat{Q}_{\lambda}$ and hence fewer edges in the graph. Choosing an appropriate value for $\lambda$ is thus critical. On the one hand, we want to enforce sparsity in the graph (for interpretation purposes), so only significant connections are maintained in the network. On the other hand, $\hat{Q}_{\lambda}$ should not be overly sparse so that it remains nonsingular (for numerical stability purposes) and relatively unbiased in order to appropriately reflect the true dependence structure in the observed data. In our river discharge application in Section 5.1 we use a voting procedure to select the best value for $\lambda$, while in our global currency application in Section 5.2 we set the sparsity level to a pre-defined level for interpretation purposes.

### 3.2.2. Structured Graph Learning via Laplacian Spectral Constraints

As an alternative to the graphical Lasso approach, we can seek to include more structural information into the graph by using the structured graph Laplacian (SGL) method of Kumar et al. (2019), which assumes that the signal residing on the graph changes “smoothly” between connected nodes. Essentially, the method enforces the precision matrix $Q$ to be a Laplacian matrix, that is, $Q \in S_{\Omega} = \{\Omega \mid \Omega_{ij} = \Omega_{ji} \le 0, i \neq j; \Omega_{ii} = -\sum_{j \neq i} \Omega_{ij}\}$. Imposing such “Laplacian constraints” allows us to better balance the sparsity and connectedness of the estimated precision matrix, as the eigenvalues of the Laplacian matrix encode information about the graph structure. For instance,
<figure><img src="image_19.png" /></figure>

if exactly one eigenvalue is zero and all other eigenvalues are positive, then the graph is connected. Here, we set $Q = \mathcal{L}w$ with $\mathcal{L}$ the linear operator that maps a nonnegative set of edge weights $w \in \mathbb{R}_+^{p(p-1)/2}$ to the matrix $Q$ with Laplacian constraints. The Laplacian matrix $Q$ can then be factorized as $Q = U\text{diag}(\eta)U^T$ (with an orthogonal matrix $U$ and eigenvalues $\eta$) to enforce appropriate constraints on $\eta$. Therefore, the optimization problem can be formulated as follows:

$$ (\hat{\eta}, \hat{U}) = \arg \max_{\lambda, U} \max_{w} \left( \log \text{gdet}(\text{Diag}(\eta)) - \text{tr}(\hat{\Sigma} \mathcal{L} w) - \beta \|\mathcal{L}w\|_1 - \frac{\gamma}{2} \|\mathcal{L}w - U \text{Diag}(\eta) U^T\|_F^2 \right), \\ \text{subject to } w \ge 0, \eta \in \mathcal{H}, \text{ and } U^T U = I, $$

where $\mathcal{H}$ denotes the set of constrained eigenvalues, $\| \cdot \|_F$ is the Frobenius norm, and gdet is the generalized determinant defined as the product of all positive values in $\eta$. We then obtain the estimate as $\hat{Q} = \hat{U}\text{diag}(\hat{\eta})\hat{U}^T$. The optimization problem can be viewed as a penalized likelihood if data have been generated from a Gaussian Markov random field; in more general cases such as ours, it still provides meaningful graphical structures since it can be viewed as a so-called penalized log-determinant Bregman divergence problem. Therefore, this method can be seen as an extension of the graphical Lasso with additional spectral constraints used to control the graph structure. A larger value of $\beta$ increases the sparsity level of the graph. The hyperparameter $\gamma \ge 0$ additionally controls the level of connectedness, and a larger value of $\gamma$ enforces a higher level of connectedness of the estimated graph structure.

## 4. Simulation Study

This section presents two simulation studies aimed at illustrating the efficacy of the proposed graph learning methods. In Section 4.1, we demonstrate the versatility of our methods across various graph types in small dimensions. Subsequently, in Section 4.2, we simulate data from a higher-dimensional tree structure that emulates a real-world example considered in one of our data applications.

### 4.1. Small-Dimensional Setting with General Extremal Graph Structure

We start with a simulation study, in which we consider three small-dimensional extremal dependence graph structures of different types; the cases we choose are displayed in Figure 1, that is, a tree (Case 1), a decomposable graph (Case 2), and a nondecomposable graph (Case 3), each with $p = 4$ variables.

The simulated models are constructed as follows. We simulate datasets of iid random variables $R_1, R_2, R_3, R_4 \stackrel{\text{iid}}{\sim} \text{Fréchet}(2)$ and construct the random vector $X = (X_1, X_2, X_3, X_4)^T$ according to the following three cases, for which we know the exact true TPDM $\Sigma$ and its inverse $Q = \Sigma^{-1}$:

**Case 1:**

$$ \left\{ \begin{aligned} X_1 &= R_1 \\ X_2 &= R_1 \oplus R_2 \\ X_3 &= R_1 \oplus R_3, \\ X_4 &= R_1 \oplus R_4 \end{aligned} \right. , \quad \Sigma = \begin{bmatrix} 1 & 1 & 1 & 1 \\ 1 & 2 & 1 & 1 \\ 1 & 1 & 2 & 1 \\ 1 & 1 & 1 & 2 \end{bmatrix}, \quad Q = \begin{bmatrix} 4 & -1 & -1 & -1 \\ -1 & 1 & 0 & 0 \\ -1 & 0 & 1 & 0 \\ -1 & 0 & 0 & 1 \end{bmatrix}. $$

**Case 2:**

$$ \left\{ \begin{aligned} X_1 &= R_1 \\ X_2 &= R_1 \oplus R_2 \\ X_3 &= R_1 \oplus R_3 \\ X_4 &= R_1 \oplus (2 \circ R_3) \oplus R_4 \end{aligned} \right. , \quad \Sigma = \begin{bmatrix} 1 & 1 & 1 & 1 \\ 1 & 2 & 1 & 1 \\ 1 & 1 & 2 & 3 \\ 1 & 1 & 3 & 6 \end{bmatrix}, \quad Q = \begin{bmatrix} -4 & -1 & -3 & 1 \\ -1 & 1 & 0 & 0 \\ -3 & 0 & 5 & -2 \\ 1 & 0 & -2 & 1 \end{bmatrix}. $$

**Case 3:**

$$ \left\{ \begin{aligned} X_1 &= R_1 \\ X_2 &= R_1 \oplus (3/\sqrt{6} \circ R_2) \\ X_3 &= R_1 \oplus (1/\sqrt{6} \circ R_2) \oplus (2/\sqrt{3} \circ R_3) \\ X_4 &= R_1 \oplus (\sqrt{6}/3 \circ R_2) \oplus (1/\sqrt{3} \circ R_3) \oplus R_4 \end{aligned} \right. , \quad \Sigma = \begin{bmatrix} 1 & 1 & 1 & 1 \\ 1 & 2.5 & 1.5 & 2 \\ 1 & 1.5 & 2.5 & 2 \\ 1 & 2 & 2 & 3 \end{bmatrix}, \quad Q = \begin{bmatrix} 2 & -0.5 & -0.5 & 0 \\ -0.5 & 1 & 0 & -0.5 \\ -0.5 & 0 & 1 & -0.5 \\ 0 & -0.5 & -0.5 & 1 \end{bmatrix}. $$

In each case, we generate datasets of size $n = 1000, 4000$, and 50,000 following the above construction and proceed as follows to infer the extremal graph structure. First, we estimate the TPDM $\Sigma$ of $X$ based on the estimator $\hat{\Sigma}$ specified through (3) with different choices of threshold $r_0$, which determines the number of the threshold exceedances to estimate the TPDM. To assess whether our proposed method is consistent, we choose $r_0$ to be a high empirical $u$-quantile of the observed vector of $L_2$ norms with $u = 0.9, 0.95, 0.99$ when $n = 1000, 4000, 50,000$, leading to $n_e = 100, 200, 500$ exceedances, respectively. Then we

<figure><img src="image_20.png" /><figcaption>Figure 1. Examples of undirected graph structures: a tree (left), a decomposable graph (middle), and a nondecomposable graph (right).</figcaption></figure>
<figure><img src="image_21.png" /></figure>

apply the extremal graphical Lasso and SGL methods to estimate the unknown graph structure. In each setting, the range of $\lambda$ and $\{\beta, \gamma\}$ penalty values is chosen to span a wide range of graphical structures, from fully connected to fully sparse (no connection). The simulation is repeated $N = 300$ times to estimate the graph recovery rate.

Table 1 (extremal graphical Lasso) and Table 2 (SGL) report the percentage of times (out of $N = 300$ simulations) that the method correctly recovers the true graph structure in each setting, when the sparsity level is correctly specified (i.e., when the number of estimated edges matches the number of edges of the true underlying graph). We find that, in these simple settings, the two methods perform well at recovering the true underlying graph structure, and that the SGL method is uniformly better than the extremal graphical Lasso. Interestingly, we also see that the results improve as the sample size increases (together with the threshold), which is in line with consistency.

To better understand how the estimated graph varies as a function of the chosen penalty parameter, Figure 2 displays the estimated extremal structure for one dataset following Case 3 (general nondecomposable “square graph”) with $n = 50,000$, $u = 0.99$, when the extremal graphical Lasso method is used and the penalty parameter $\lambda$ ranges in a suitable interval. The heading of each display shows the range of $\lambda$ values that leads to the estimated graph shown below. The penalty

parameter $\lambda$ can be tuned to control the number of edges in the graph, that is, the sparsity level: when $\lambda$ decreases, the number of edges increases, and vice versa. Interestingly, the proposed method always retrieves true connections in this case (i.e., it never yields wrong connections) whenever the estimated graph is as sparse, or sparser than the true graph. This experiment shows that in simple settings (e.g., small-dimensional graph, and large sample size), our method is generally able to retrieve the true extremal dependence graph structure, provided the penalty parameter is well specified.

## 4.2. Higher-Dimensional Setting with Extremal Tree Structure

To assess the performance of our methods in a higher-dimensional and more realistic setting, we further conduct a simulation study with a graph of dimension $p = 31$. The true graph structure is set to be the tree structure from our river flow data application in Section 5.1 (see the river network in Figure 5). Data are simulated according to a transformed-linear causal structural model that follows the river network, where a downstream variable depends on its parent (i.e., upstream) variables and a heavy-tailed noise term. Specifically, for a single replicate, we first simulate iid random innovations $R_i \sim \text{Fréchet}(2)$, $i = 1, \dots, 31$; then, for each source node $X_i$ in the graph (i.e., with no parent variable), we set $X_i = R_i$, whereas for each non-source variable $X_i$ (i.e., with $K_i \neq 0$ parent variables), we set

$$X_i = \{(\theta/K_i)^{1/2} \circ R_i\} \oplus \left[ \{1 - (\theta/K_i)\}/K_i \right]^{1/2} \circ X_k,$$

where $Pa(i)$ is the set of parent nodes of node $i$, $K_i = |Pa(i)| \in \{1, 2, \dots\}$ is their number, and $\theta \in [0, 1]$ is a weight that balances the contributions of the innovation and the parent variables. This construction ensures that each variable $X_i$ has the tail of a Fréchet(2) variable with unit scale parameter; moreover, as is expected with true river networks, the influence of the random innovation noise term diminishes when there are multiple incoming flows. The overall strength of extremal dependence along the graph is determined by $\theta$ and we set it here to $\theta = 0.5$. As in the small-dimensional simulation example, we generate $N = 300$ datasets of size $n = 1000, 4000, 50,000$ following the above construction, and estimate the TPDM in each case using the empirical $u$-quantile of the observed $L_2$ norms with $u = 0.9, 0.95, 0.99$ when $n = 1000, 4000, 50,000$, leading to $n_e = 100, 200, 500$ exceedances, respectively. The extremal graphical Lasso method is then applied with a suitable

<table><caption>**Table 1.** Simulation results for the three types of graphs with $p = 4$ (rows), for three sample sizes and thresholds (columns).</caption><thead><tr><th></th><th>n = 1000, u = 0.9<br>n<sub>e</sub> = 100</th><th>n = 4000, u = 0.95<br>n<sub>e</sub> = 200</th><th>n = 50,000, u = 0.99<br>n<sub>e</sub> = 500</th></tr></thead><tbody><tr><td>Case 1</td><td>99.67%</td><td>100%</td><td>100%</td></tr><tr><td>Case 2</td><td>99.33%</td><td>100%</td><td>100%</td></tr><tr><td>Case 3</td><td>65.67%</td><td>81%</td><td>97.33%</td></tr></tbody></table>

NOTE: Numbers in each cell are the percentage of the times (out of $N = 300$ simulations) that the extremal graphical Lasso recovers the true graph extremal structure, when the sparsity level is correctly specified.

<table><caption>**Table 2.** Simulation results for the three types of graphs with $p = 4$ (rows), for three sample sizes and thresholds (columns).</caption><thead><tr><th></th><th>n = 1000, u = 0.9<br>n<sub>e</sub> = 100</th><th>n = 4000, u = 0.95<br>n<sub>e</sub> = 200</th><th>n = 50,000, u = 0.99<br>n<sub>e</sub> = 500</th></tr></thead><tbody><tr><td>Case 1</td><td>100%</td><td>100%</td><td>100%</td></tr><tr><td>Case 2</td><td>100%</td><td>100%</td><td>100%</td></tr><tr><td>Case 3</td><td>99%</td><td>100%</td><td>100%</td></tr></tbody></table>

NOTE: Numbers in each cell are the percentage of the times (out of $N = 300$ simulations) that the SGL method recovers the true graph extremal structure, when the sparsity level is correctly specified.

<figure><img src="image_22.png" /><figcaption>**Figure 2.** Estimated extremal graph structures using the extremal graphical Lasso method based on the PTCC for a single dataset following Case 3, with $n = 50,000$, $u = 0.99$, plotted as a function of the penalty parameter $\lambda$ (shown at the top of each display).</figcaption></figure>
<figure><img src="image_23.png" /></figure>

<figure><img src="image_24.png" /><figcaption>**Figure 3.** Boxplots (based on *N* = 300 simulations) of the percentage of correctly identified edges in the high-dimensional simulation study mimicking the river network example with *p* = 31, for three sample sizes and threshold levels (left to right), using the extremal graphical Lasso with a correctly specified sparsity level (setting *p* − 1 = 30 edges).</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>**Figure 4.** Edgewise precision (i.e., percentage of the *N* = 300 simulations that a specific edge is correctly identified) for each edge of the true graph, assumed to be the river network of Section 5 with *p* = 31, for two sample sizes and threshold levels (Left: *n* = 1000, *u* = 0.9; right: *n* = 4000, *u* = 0.95), using the extremal graphical Lasso with a correctly specified sparsity level (setting *p* − 1 = 30 edges).</figcaption></figure>

penalty value λ set to span a wide range of sparsity levels. We do not consider the SGL method in this simulation, because it was found to outperform the extremal graphical Lasso in the small-dimensional setting, so we expect it to be better here, as well. Figure 3 displays, for each sample size, boxplots containing 300 values (one value for each of the *N* simulations) showing the proportion of edges among those in the estimated graph that are indeed true edges. Overall, the extremal graphical Lasso works quite well at recovering the true extremal graph structure and improves as the sample size increases, again supporting consistency in this higher-dimensional case.

To have a closer look at the edgewise precision, Figure 4 displays the percentage of times (out of the *N* = 300 simulations) that each true connection is correctly identified by our extremal graphical Lasso method when we set the number of estimated edges to be exactly 30, for the first two settings of Figure 3 (with sample size *n* = 1000 and 4000, respectively). The results for the last setting can be found in the supplementary material.

As before, we find that the edgewise precision increases for each true edge as *n*, *u*, and *n*<sub>e</sub> increase simultaneously.

We finally note that with our distribution-free approach, it is not obvious how to estimate the true (unknown) sparsity level and optimally select the penalty parameters. However, we can use expert knowledge or problem-specific criteria to facilitate

interpretation of the estimated graph (and thus set the penalty parameters); see the data applications in Section 5.

**5. Real Data Applications**

Risk networks are useful in quantitative risk management to elucidate complex extremal dependence structures in collections of random variables. We show two examples of both environmental and financial risk analysis. First, we study river discharge data of the upper Danube basin (Asadi, Davison, and Engelke 2015), which has become a benchmark dataset for learning extremal networks in the recent literature. The true underlying physical river flow network is available, which can be used as a benchmark to compare the performance of our method with other existing approaches. Second, we apply our method to historical global currency exchange rate data from different historical periods, including different economic cycles, the COVID-19 period, and the period of the 2022 military conflict opposing Russia and Ukraine (from its start until September 26, 2022).

**5.1. Extremal Network Estimation for a River Network**

We apply our method to study the dependence structure of extreme discharges on the river network of the upper Danube basin (see the left panel of Figure 5 for the topographic map).
<figure><img src="image_26.png" /></figure>

<figure><img src="image_27.png" /><figcaption>**Figure 5.** Left: Topographic map of the upper Danube basin (from Asadi, Davison, and Engelke 2015), showing 31 sites of gauging stations (red circles) and the altitudes of the region. Right: The true physical river flow connections; the arrows show the flow directions.</figcaption></figure>

<figure><img src="image_28.png" /><figcaption>**Figure 6.** Left: Estimated TPDM of the river discharge data from the upper Danube basin. Right: Votes (%) of the edges selected based on the extremal graphical Lasso method. Darker red cells indicate that the corresponding edge has been selected more often by the extremal graphical Lasso.</figcaption></figure>

This region has been regularly affected by severe flooding events in its history, which have caused losses of human lives and damage to material goods. The original daily discharge data from 1960 to 2009 were provided by the Bavarian Environmental Agency (http://www.gkdl.bayern.de), and Asadi, Davison, and Engelke (2015) preprocessed the data, which now include $n = 428$ approximately independent events $X_1, \dots, X_n \in \mathbb{R}^p$ recorded at $p = 31$ gauging stations located on the river network from three summer months (June, July, and August), obtained using declustering methods. The data were later also studied by Engelke and Hitz (2020) and Hu, Peng, and Segers (2022) among others, using graphical models for extremes based on a conditional independence notion adapted to multivariate Pareto distributions. The true physical river flow connections and directions are represented by a directed graph shown in the right panel of Figure 5, where the arrows indicate the flow directions. This can serve as an accurate benchmark of the “true” conditional independence structure, against which we can compare the results from our proposed

extremal graphical structure learning methods based on the PTCC.

### 5.1.1. Graph Structure Learning Using the Extremal Graphical Lasso

To learn the extremal dependence structure of the river network, we first perform a nonparametric empirical transformation of the data to satisfy Fréchet($\alpha = 2$) margins, for each station separately. Next, we estimate the TPDM, $\Sigma$, using the proposed estimator, $\hat{\Sigma}$, defined through (3). In particular, we choose $m = p = 31$ because the margins are preprocessed to have a common Fréchet scale and $r_0 = 11.4$, which corresponds to the empirical 90% quantile. Therefore, $n_{exc} = 43$ is the number of extreme observations (i.e., threshold exceedances) that are used to estimate $\Sigma$. The left panel of Figure 6 displays the estimated TPDM of the river discharge data from the upper Danube basin, while the right panel displays the votes (in percentage) of the edges selected based on the extremal graphical Lasso method, obtained from multiple fits with a range of $\lambda$ values producing
<figure><img src="image_29.png" /></figure>

<figure><img src="image_30.png" /><figcaption>**Figure 7.** Left: Number of edges selected by the SGL method as a function of the penalty parameters $\beta$ and $\gamma$. Lighter blue cells correspond to parameter combinations producing sparser graphs. Right: Votes (%) of the edges selected based on the SGL method. Darker red cells indicate that the corresponding edge has been selected more often by the SGL method.</figcaption></figure>

different dependence structures, from fully connected to fully sparse graphs.

### 5.1.2. Graph Structure Learning Using the SGL Method

To enhance connectedness, we further explore the SGL method, which learns sparse graph structures under additional spectral constraints. In particular, as described in Section 3.2.2, we can control both the sparsity level and the graph connectedness by tuning the two penalty parameters $\beta$ and $\gamma$, respectively. As shown in the left panel of Figure 7, the overall graph sparsity varies for different combinations of $\beta$ and $\gamma$. The right panel of Figure 7 displays the votes (in percentage) of the edges selected based on the SGL method, obtained by fitting a large number of models for each of the $\{\beta, \gamma\}$ combinations shown in the left panel.

### 5.1.3. Results: Estimated Extremal River Discharge Networks

For both the extremal graphical Lasso and the SGL method, it is important to carefully select the penalty parameters, $\lambda$ and $\{\beta, \gamma\}$, respectively. These tuning parameters impact both the sparsity level and the connectedness of the resulting graph structure. Ideally, we would like to obtain a sparse graph while keeping it connected. However, there is a tradeoff between these two requirements. Our approach is to control the overall sparsity level while imposing a soft connectedness condition: we start from the fully sparse graph (no edges) and sequentially add edges between nodes according to the ranking of the votes shown in Figures 6 and 7, until no node is left alone (i.e., each node has at least one connection with another node). The estimated graph structure thus prioritizes edges that are most often selected and is obtained by blending the results from several model fits, which also makes it less sensitive to specific values of the penalty parameters.

Figure 8 displays the estimated extremal river discharge network using our approach for both the extremal graphical Lasso and the SGL method, respectively. The edge thickness is proportional to the votes shown in Figures 6 and 7.

Recall that when two nodes are not connected, it means that they are partially tail-uncorrelated in our framework. The estimation based on the extremal graphical Lasso method has a few more edges than the true physical river flow network. The extra links could be interpreted as being due to extremal dependence induced by regional weather events, though this is not very clear. By contrast, the estimation result based on the SGL method matches most of the true river flow connections, while the votes (shown in terms of the edge thickness) represent the strength of the extremal dependence connections. Interestingly, this dependence strength seems well-aligned with the actual physical strength of the river flow. Compared with existing methods (Engelke and Hitz 2020; Klüppelberg and Krali 2021), our estimated extremal network based on the SGL method looks more realistic (thus, more easily interpretable), as it is closer to the true flow structure of the river network, though the recent results from Engelke, Lalancette, and Volgushev (2021) are quite similar to ours (with a few extra connections).

## 5.2. Extremal Network Estimation for Global Currency Exchange Rate Network

We now apply our method to explore historical global currency exchange rate data for 20 currencies from different historical periods, including two different global economic cycles (2009–2014 and 2015–2019, where the segmentation is determined from the world GDP cycles illustrated in the supplementary material, COVID-19 (January 01, 2020–February 23, 2022), and the period from the beginning of the 2022 military conflict between Russia and Ukraine until a most recent date when we downloaded the data (February 24, 2022–September 26, 2022). Historical data were downloaded from Yahoo Finance. We chose the currencies from all G20 countries and also added currencies from Ukraine and Kazakhstan. The list of selected currencies and their corresponding symbols can be found in the supplementary material. As the unit of all currencies is the US Dollar
<figure><img src="image_31.png" /></figure>

<figure><img src="image_32.png" /><figcaption>**Figure 8.** Extremal river discharge network estimated using the extremal graphical Lasso (left) and the SGL method (right). The edge thickness is proportional to the votes shown in Figures 6 and 7.</figcaption></figure>

<figure><img src="image_33.png" /><figcaption>**Figure 9.** Extremal currency exchange rate risk network estimated for 2009–2014 (left) and 2015–2019 (right). Different edge types (respectively thickest, thick, and thin) with different colors (respectively red, blue, and gray) indicate a frequency of > 90%, 70%–90%, and 50%–70%, respectively, to be selected among the 300 bootstrap fitted models, while absent edges indicate that this frequency is less than 50%.</figcaption></figure>

(USD), USD is not considered in our list of currencies under study.

First, we preprocess the historical daily closing prices of the currencies. An ARMA(1,1)-GARCH(1,1) time series model is fitted to the negative log return time series of each currency, and then standardized residuals are extracted and transformed marginally to Fréchet margins with shape parameter $\alpha = 2$. The extremal dependence graph structure governing negative log returns therefore represents partial associations among extreme losses, shedding light on the integration and/or vulnerability of major economies in periods of stress. To estimate this risk network, we follow the same procedure as before, estimating first the TPDM using the estimator in (3) with $r_0$ as the empirical 90% quantile, and then applying the SGL method.

In this analysis, we use only the SGL method, because it has shown better performance than the extremal graphical Lasso method in our simulation study and the Danube River application. Furthermore, for comparability among different historical periods, we here fix the sparsity level to 80% (i.e., with only

about 38 edges), rather than using the method based on votes. This approach yields a specific combination of {$\beta, \gamma$} penalty parameters (when there are multiple combinations available, we choose the one with the minimal $\beta$ value), which then yields the final estimated graph structure. To assess the estimation uncertainty of the graph structure, we have further conducted 300 bootstrap simulations, whereby standardized residuals are resampled with replacement and the TPDM is then re-estimated, as well as the graph structure fixing the same sparsity level (i.e., potentially with different selected {$\beta, \gamma$} values for each bootstrap simulation). The bootstrap results are shown in Figures 9 and 10 for the first two, and last two periods, respectively, where different edge types represent the “significance” of the displayed connections: the thickest edges (in red) indicate a frequency of at least 90% to be included among the 300 bootstrap fitted models; thick edges (in blue) indicate a frequency between 70% and 90% to be selected; thin edges (in gray) indicate a frequency between 50% and 70% to be selected; absent edges have been selected less than 50% of the time. Moreover, the size of nodes in the
<figure><img src="image_34.png" /></figure>

<figure><img src="image_35.png" /><figcaption>**Figure 10.** Extremal currency exchange rate risk network estimated for the COVID-19 period (left) and the period since the 2022 military conflict between Russia and Ukraine (right). Different edge types (respectively thickest, thick, and thin) with different colors (respectively red, blue, and gray) indicate a frequency of > 90%, 70%–90%, and 50%–70%, respectively, to be selected among the 300 bootstrap fitted models, while absent edges indicate that this frequency is less than 50%.</figcaption></figure>

displayed graphs is proportional to their degree, that is, to their number of connections. The bigger a node, the more connected it is, giving an idea of the centrality of a currency in the risk network.

To show that our methods can be helpful tools at an exploratory stage, we now provide our own interpretation of the estimated risk networks, assuming that more strongly connected nodes tend to be more vulnerable/exposed to network risks, or are important currencies that strongly determine the behavior of other currencies in the monetary system. While we here attempt to back up some of our findings using historical events and economic development regimes for different countries, we also note that further research (through more advanced models and insightful economic expertise) would be required to fully confirm the conclusions drawn below. From the left panel in Figure 9, representing the economic cycle 2009–2014, the strongest-connected currencies within the estimated risk network are GBP (United Kingdom) and ARS (Argentina), whereas KZT (Kazakstan) and UAH (Ukraine) are less connected from the major clusters. By contrast, in the right panel of Figure 9, representing the economic cycle 2015–2019, the strongest-connected currencies from the estimated risk network are KRW (Republic of Korea), AUD (Australia) and EUR (EU), whereas ARS, KZT, and MXN (Mexico) are relatively isolated. ARS and GBP are two examples of currencies that were much weaker connected during the second period. Historical economic data show that the Argentinian economy went through a major crisis around 2000, establishing strong economic links to other countries and international financial markets in the following years. This seems to clearly transpire from the estimated graph for the 2009–2014 period following the global financial crises of 2008, but the monetary interconnections of ARS to other currencies have drastically been reduced in the later 2015–2019 period. As to GBP, the exit of Great Britain from the European Union (Brexit), and the

political and economic decisions accompanying it, could be the origin of this change.

The left panel of Figure 10 shows that during the COVID-19 pandemic, the ARS currency is again the strongest-connected one, by far, in the risk network, whereas TRY (Turkey) and CAD (Canada) are isolated and most of the other currencies are also less interconnected. As for ARS, our result might be due to the fact that Argentina was struck by COVID-19 during an already fragile economic situation, and its aforementioned strong international monetary ties might have led to strong extremal connectedness to other countries. Finally, the right panel of Figure 10 shows the risk network for currency exchange rates from the beginning of the 2022 military conflict between Russia and Ukraine until September 26, 2022. Since the length of the time period is shorter, there is higher estimation uncertainty, and no edges have been selected more than 90% of the time among the 300 bootstrap model fits. Nevertheless, we can still make the interesting observation that RUB (Russia) is the only isolated currency from the rest of the network, which might be due to the antagonism and strong economic sanctions imposed by Western countries, as well as the dramatic changes in monetary policies.

## 6. Conclusions

We have proposed partial tail-correlation as a novel notion of extremal dependence measure that removes the effect of confounding variables through transformed-linear operations. Unlike other approaches from the recent literature, our new partial tail-correlation coefficient (PTCC) assumes multivariate regular variation but does not rely on any further strict parametric assumptions. Furthermore, the PTCC has appealing theoretical properties and it can be used to define a new class of extremal graphical models, where the absence of edges indicates partial tail-uncorrelatedness between variables (i.e., when
<figure><img src="image_36.png" /></figure>

the PTCC of the corresponding edges equals zero). We have shown that the zero PTCC values between variable pairs can be retrieved by identifying the zero entries in the inverse tail pairwise dependence matrix (TPDM). This convenient property, which is akin to classical Gaussian graphical models, allows us to efficiently learn high-dimensional extremal networks defined in terms of the PTCC by exploiting state-of-the-art methods from graph theory, such as the graphical Lasso and structured graph learning via Laplacian spectral constraints.

Our graph-inference approach is flexible, can be applied to general undirected graphs, and easily scales to high dimensions. We demonstrate the effectiveness of our method as an exploratory tool for interpretable extremal network estimation. In our first application to river discharge data from the upper Danube basin, we show that the proposed method outperforms other existing methods by realistically capturing most physical flow connections, together with the strength of the connections, while largely avoiding spurious connections. In our second application to historical currency exchange rate data, we obtain interesting findings based on the estimated risk network for four recent periods. While further insights could be obtained from domain experts (e.g., economists), our exploratory results seem to make sense in view of certain recent historical events.

We have identified some theoretical and methodological challenges for future research. First, we have not obtained formal theoretical guarantees that our graph learning method can consistently estimate the unknown graph structure. However, we note that our proposed method relies on an estimator of the TPDM, whose consistency was established by Davis and Mikosch (2009) and Larsson and Resnick (2012), and on widely-used graph-learning techniques, whose consistency was established by Ravikumar et al. (2011) and Pavez (2022) under broad conditions. While these separate results give us some reasonable confidence that our proposed extremal graphical learning methods are also consistent under mild conditions, it would be interesting in future research to formally verify it. Nevertheless, the results of our simulation study show that we can recover the true extremal graph quite well when the sparsity level is correctly specified, and they are in line with the proposed procedure being consistent. Moreover, our real data applications have provided convincing evidence that the proposed methods can extract insightful structural information from extreme observations. While the true sparsity level is unknown in practice and often decided by expert knowledge or set to a fixed percentage to facilitate interpretation (as we do in our currency application), it would also be interesting in future research to develop a more principled data-driven way of optimally specifying the number of edges. Second, we have focused in this article on graph learning through the extremal graphical Lasso, and the SGL method that imposes further Laplacian constraints. It would be interesting to extend these methods, in order to let the graphical structure depend on well-chosen covariates (e.g., temperature in our river network application, or the time period index in our currency network application) so that the estimated (nonstationary) risk network can then be interpreted through the glasses of these covariates. Moreover, unlike our currency application, where separate model fits were obtained for each time period under consid-

eration, introducing covariates in the graph learning procedure would allow estimating the networks simultaneously from the combined dataset, thus, gaining efficiency for higher accuracy of the estimated graph structure. A possibility could be to extend the Gaussian graphical regression approach proposed by Zhang and Li (2023) to our PTCC-based extremes framework. Third, by analogy with the widely used Gaussian Markov random field framework, one could imagine constructing Markov random fields for extremes where the distribution of connected variables may be fitted with an asymptotically justified model for extremes (e.g., of multivariate Pareto type), whereas the factorization of the likelihood could be pre-specified by the extremal network learned from the data in a preliminary step using our proposed method. Finally, another direction to investigate concerns the geometric representation of multivariate extremes (Simpson, Wadsworth, and Tawn 2021; Nolde and Wadsworth 2022). It would be interesting to see if the new PTCC can be defined through this geometric approach and to explore its links with other popular measures of extremal dependence.

## Supplementary Materials

**PDF Supplement:** This supplement contains the list of currencies used in our application and their symbols, a figure that shows the world GDP and how we identified economic cycles, as well as further simulation results for the high-dimensional tree example in the large sample size case (PDF file)

**R Code:** This supplement contains a working example implemented in R to reproduce the results of our simulation study in one setting (R file)

## Acknowledgments

We point out that there has been independent and parallel work by Lee and Cooley (2021), who also investigate partial tail-correlation for extremes. Our understanding is that inference in their work focuses rather on hypothesis testing (with the goal of checking if the partial tail-correlations for given pairs of variables are significantly different from zero), and not on learning extremal networks with structural constraints. The authors would like to thank the Editor, Associate Editor, and two anonymous reviewers for their constructive comments and suggestions that helped improve the quality of the manuscript.

## Disclosure Statement

The authors report there are no competing interests to declare.

## Funding

This publication is based upon work supported by the King Abdullah University of Science and Technology (KAUST) Office of Sponsored Research (OSR) under Award No. OSR-CRG2020-4394.

## ORCID

Raphaël Huser http://orcid.org/0000-0002-1228-2071

## References

Ahmed, M., Maume-Deschamps, V., Ribereau, P., and Vial, C. (2020), “Spatial Risk Measures for Max-Stable and Max-Mixture Processes,” *Stochastics*, 92, 1005–1020. [1]
<figure><img src="image_37.png" /></figure>

Asadi, P., Davison, A. C., and Engelke, S. (2015), “Extremes on River Networks,” *The Annals of Applied Statistics*, 9, 2023–2050. [9,10]

Baba, K., Shibata, R., and Sibuya, M. (2004), “Partial Correlation and Conditional Correlation as Measures of Conditional Independence,” *Australian & New Zealand Journal of Statistics*, 46, 657–664. [2]

Bacro, J.-N., Gaetan, C., Opitz, T., and Toulemonde, G. (2020), “Hierarchical Space-Time Modeling of Asymptotically Independent Exceedances with an Application to Precipitation Data,” *Journal of the American Statistical Association*, 115, 555–569. [1]

Bassamboo, A., Juneja, S., and Zeevi, A. (2008) Portfolio credit risk with extremal dependence: Asymptotic analysis and efficient simulation. *Operations Research*, 56, 593–606. [1]

Bekiros, S., and Uddin, G. S. (2017), “Extreme Dependence Under Uncertainty: An Application to Stock, Currency and Oil Markets,” *International Review of Finance*, 17, 155–162. [1]

Bopp, G. P., Shaby, B. A., and Huser, R. (2021), “A Hierarchical Max-Infinitely Divisible Spatial Model for Extreme Precipitation,” *Journal of the American Statistical Association*, 116, 93–106. [1]

Cooley, D., and Thibaud, E. (2019), “Decompositions of Dependence for High-Dimensional Extremes,” *Biometrika*, 106, 587–604. [2,3,4,5]

Davis, R. A., and Mikosch, T. (2009), “The Extremogram: A Correlogram for Extreme Events,” *Bernoulli*, 15, 977–1009. [14]

Davis, R. A., Mikosch, T., and Zhao, Y. (2013), “Measures of Serial Extremal Dependence and their Estimation,” *Stochastic Processes and their Applications*, 123, 2575–2602. [1]

Davison, A. C., and Huser, R. (2015), “Statistics of Extremes,” *Annual Review of Statistics and its Application*, 2, 203–235. [1]

Davison, A. C., Huser, R., and Thibaud, E. (2019), “Spatial Extremes,” in *Handbook of Environmental and Ecological Statistics*, eds. A. E. Gelfand, M. Fuentes, J. A. Hoeting, and R. L. Smith, pp. 711–744, Boca Raton, FL: CRC Press. [1]

Engelke, S., and Hitz, A. S. (2020), “Graphical Models for Extremes,” *Journal of the Royal Statistical Society*, Series B, 82, 871–932. [2,10,11]

Engelke, S., and Ivanovs, J. (2021), “Sparse Structures for Multivariate Extremes,” *Annual Review of Statistics and Its Application*, 8, 241–270. [2]

Engelke, S., Lalancette, M., and Volgushev, S. (2021), “Learning Extremal Graphical Structures in High Dimensions,” arXiv preprint arXiv:2111.00840. [2,11]

Engelke, S., Opitz, T., and Wadsworth, J. (2019), “Extremal Dependence of Random Scale Constructions,” *Extremes*, 22, 623–666. [1]

Engelke, S., and Volgushev, S. (2022), “Structure Learning for Extremal Tree Models,” *Journal of the Royal Statistical Society*, Series B, 84, 2055–2087. [2]

Ferro, C. A., and Stephenson, D. B. (2011), “Extremal Dependence Indices: Improved Verification Measures for Deterministic Forecasts of Rare Binary Events,” *Weather and Forecasting*, 26, 699–713. [1]

Fix, M. J., Cooley, D. S., and Thibaud, E. (2021), “Simultaneous Autoregressive Models for Spatial Extremes,” *Environmetrics*, 32, e2656. [2]

Fougères, A.-L., Mercadier, C., and Nolan, J. P. (2013), “Dense Classes of Multivariate Extreme Value Distributions,” *Journal of Multivariate Analysis*, 116, 109–129. [4]

Friedman, J., Hastie, T., and Tibshirani, R. (2008), “Sparse Inverse Covariance Estimation with the Graphical Lasso,” *Biostatistics*, 9, 432–441. [6]

Gissibl, N. (2018), “Graphical Modeling of Extremes: Max-Linear Models on Directed Acyclic Graphs,” Ph.D. thesis, Universitätsbibliothek der TU München. [2,4]

Gnecco, N., Meinshausen, N., Peters, J., and Engelke, S. (2021), “Causal Discovery in Heavy-Tailed Models,” *The Annals of Statistics*, 49, 1755–1778. [4]

Gong, Y., and Huser, R. (2022), “Asymmetric Tail Dependence Modeling, with Application to Cryptocurrency Market Data,” *The Annals of Applied Statistics*, 16, 1822–1847. [1]

Hu, S., Peng, Z., and Segers, J. (2022), “Modelling Multivariate Extreme Value Distributions via Markov Trees,” arXiv preprint arXiv:2208.02627. [10]

Huang, W. K., Cooley, D. S., Ebert-Uphoff, I., Chen, C., and Chatterjee, S. (2019), “New Exploratory Tools for Extremal Dependence: χ Networks and Annual Extremal Networks,” *Journal of Agricultural, Biological and Environmental Statistics*, 24, 484–501. [1]

Huser, R., and Davison, A. C. (2014), “Space-Time Modelling of Extreme Events,” *Journal of the Royal Statistical Society*, Series B, 76, 439–461. [1]

Huser, R., Opitz, T., and Thibaud, E. (2017), “Bridging Asymptotic Independence and Dependence in Spatial Extremes Using Gaussian Scale Mixtures,” *Spatial Statistics*, 21, 166–186. [1]

Huser, R., and Wadsworth, J. L. (2019), “Modeling Spatial Processes with Unknown Extremal Dependence Class,” *Journal of the American Statistical Association*, 114, 434–444. [1]

——— (2022), “Advances in Statistical Modeling of Spatial Extremes,” *Wiley Interdisciplinary Reviews: Computational Statistics*, 14, e1537. [1]

Hüsler, J., and Reiss, R.-D. (1989), “Maxima of Normal Random Vectors: Between Independence and Complete Dependence,” *Statistics & Probability Letters*, 7, 283–286. [4]

Janßen, A., and Wan, P. (2020), “k-means Clustering of Extremes,” *Electronic Journal of Statistics*, 14, 1211–1233. [4]

Klöppelberg, C., and Krali, M. (2021), “Estimating an Extreme Bayesian Network via Scalings,” *Journal of Multivariate Analysis*, 181, 104672. [2,11]

Kumar, S., Ying, J., de Miranda Cardoso, J. V., and Palomar, D. (2019), “Structured Graph Learning via Laplacian Spectral Constraints,” in *Advances in Neural Information Processing Systems* (Vol. 32). [6]

Larsson, M., and Resnick, S. I. (2012), “Extremal Dependence Measure and Extremogram: The Regularly Varying Case,” *Extremes*, 15, 231–256. [1,5,14]

Lawrance, A. (1976), “On Conditional and Partial Correlation,” *The American Statistician*, 30, 146–149. [2]

Ledford, A. W. and Tawn, J. A. (1996), “Statistics for Near Independence in Multivariate Extreme Values,” *Biometrika*, 83, 169–187. [4]

Lee, J., and Cooley, D. (2021), “Transformed-Linear Prediction for Extremes,” arXiv preprint arXiv:2111.03754. [2,4,14]

Marcon, G., Padoan, S. A., and Antoniano-Villalobos, I. (2016), “Bayesian Inference for the Extremal Dependence,” *Electronic Journal of Statistics*, 10, 3310–3337. [1]

Mhatre, N., and Cooley, D. (2020), “Transformed-Linear Models for Time Series Extremes,” arXiv preprint arXiv:2012.06705. [2]

Nolde, N., and Wadsworth, J. L. (2022), “Linking Representations for Multivariate Extremes via a Limit Set,” *Advances in Applied Probability*, 54, 688–717. [14]

Opitz, T. (2016), “Modeling Asymptotically Independent Spatial Extremes based on Laplace Random Fields,” *Spatial Statistics*, 16, 1–18. [1]

Opitz, T., Huser, R., Bakka, H., and Rue, H. (2018), “INLA Goes Extreme: Bayesian Tail Regression for the Estimation of High Spatio-Temporal Quantiles,” *Extremes*, 21, 441–462. [1]

Pavez, E. (2022), “Laplacian Constrained Precision Matrix Estimation: Existence and High Dimensional Consistency,” in *International Conference on Artificial Intelligence and Statistics*, pp. 9711–9722. [6,14]

Ravikumar, P., Wainwright, M. J., Raskutti, G., and Yu, B. (2011), “High-Dimensional Covariance Estimation by Minimizing $\ell_1$-Penalized Log-Determinant Divergence,” *Electronic Journal of Statistics*, 5, 935–980. [6,14]

Resnick, S. I. (2007), *Heavy-Tail Phenomena: Probabilistic and Statistical Modeling*, New York: Springer. [3]

Richards, J., Tawn, J. A., and Brown, S. (2022), “Modelling Extremes of Spatial Aggregates of Precipitation Using Conditional Methods,” *The Annals of Applied Statistics*, 16, 2693–2713. [1]

Rohrbeck, C., and Cooley, D. (2023), “Simulating Flood Event Sets Using Extremal Principal Components,” *The Annals of Applied Statistics*, 17, 1333–1352. [2]

Röttger, F., Engelke, S., and Zwiernik, P. (2023), “Total Positivity in Multivariate Extremes,” *The Annals of Statistics*, 51, 962–1004. [2]

Saunders, K., Stephenson, A., and Karoly, D. (2021), “A Regionalisation Approach for Rainfall based on Extremal Dependence,” *Extremes*, 24, 215–240. [1]

Sibuya, M. (1960), “Bivariate Extreme Statistics, I,” *Annals of the Institute of Statistical Mathematics*, 11, 195–210. [4]

Simpson, E. S., Wadsworth, J. L., and Tawn, J. A. (2021), “A Geometric Investigation into the Tail Dependence of Vine Copulas,” *Journal of Multivariate Analysis*, 184, 104736. [14]
<figure><img src="image_38.png" /></figure>

Speed, T. P., and Kiiveri, H. T. (1986), “Gaussian Markov Distributions over Finite Graphs,” *The Annals of Statistics*, 14, 138–150. [6]

Tran, N. M., Buck, J., and Klüppelberg, C. (2021), “Estimating a Latent Tree for Extremes,” arXiv preprint arXiv:2102.06197. [2]

Vettori, S., Huser, R., and Genton, M. G. (2019), “Bayesian Modeling of Air Pollution Extremes Using Nested Multivariate Max-Stable Processes,” *Biometrics*, 75, 831–841. [1]

Vettori, S., Huser, R., Segers, J., and Genton, M. G. (2020), “Bayesian Model Averaging over Tree-based Dependence Structures for Multivariate Extremes,” *Journal of Computational and Graphical Statistics*, 29, 174–190. [1]

Wadsworth, J., Tawn, J. A., Davison, A., and Elton, D. (2017), “Modelling Across Extremal Dependence Classes,” *Journal of the Royal Statistical Society*, Series B, 79, 149–175. [1]

Wadsworth, J. L., and Tawn, J. (2022), “Higher-Dimensional Spatial Extremes via Single-Site Conditioning,” *Spatial Statistics*, 51, 100677. [1]

Wadsworth, J. L., and Tawn, J. A. (2012), “Dependence Modelling for Spatial Extremes,” *Biometrika*, 99, 253–272. [1]

Winter, H. C., and Tawn, J. A. (2016), “Modelling Heatwaves in Central France: A Case-Study in Extremal Dependence,” *Journal of the Royal Statistical Society, Series C*, 65, 345–365. [1]

Yan, X., Wu, Q., and Zhang, W. (2019), “Cross-sectional Learning of Extremal Dependence among Financial Assets,” in *Advances in Neural Information Processing Systems* (Vol. 32). [1]

Zhang, J., and Li, Y. (2023), “High-Dimensional Gaussian Graphical Regression Models with Covariates,” *Journal of the American Statistical Association*, 118, 2088–2100. [14]

Zhong, P., Huser, R., and Opitz, T. (2022), “Modeling Nonstationary Temperature Maxima based on Extremal Dependence Changing with Event Magnitude,” *The Annals of Applied Statistics*, 16, 272–299. [1]