# Estimating probabilities of multivariate failure sets based on pairwise tail
dependence coefficients

Anna<EMAIL_ADDRESS>[ Chen
<EMAIL_ADDRESS>[ Namur Institute for Complex Systems (Naxys),
Université de Namur, Belgiumpresep=, ]e1  Econometric Institute, Erasmus
University Rotterdam presep=, ]e2

###### Abstract

An important problem in extreme-value theory is the estimation of the
probability that a high-dimensional random vector falls into a given extreme
failure set. This paper provides a parametric approach to this problem, based
on a generalization of the tail pairwise dependence matrix (TPDM). The TPDM
gives a partial summary of tail dependence for all pairs of components of the
random vector. We propose an algorithm to obtain an approximate completely
positive decomposition of the TPDM. The decomposition is easy to compute and
applicable to moderate to high dimensions. Based on the decomposition, we
obtain parameters estimates of a max-linear model whose TPDM is equal to that
of original random vector. We apply the proposed decomposition algorithm to
industry portfolio returns and maximal wind speeds to illustrate its
applicability.

Extreme-value theory,

Tail dependence,

Max-linear model,

###### keywords:

and

##  1 Introduction

Extreme value statistics deal with the characterization of extreme events,
such as stock market crashes, that occur with low frequency but with a
potentially severe impact. A typical extreme value question is to estimate the
(low) probability of such an extreme event. In the context of multivariate
extremes, we are interested in estimating the probability that a
$d$-dimensional random vector $\bm{X}=(X_{1},\ldots,X_{d})$ falls into some
“extreme region” $C$ (also called failure region). For example, in the example
of portfolio risk, the components of $\bm{X}$ represent negative log-returns
of $d$ stocks in a portfolio. Another example is in catastrophic climate
events, where $\bm{X}$ may represent maximum daily rainfall amounts at $d$
monitoring stations. To estimate such a probability $\mathbb{P}[\bm{X}\in C]$,
it is important to take the tail dependence across the margins of $\bm{X}$
into account, i.e., the tendency of extremes to occur in several components
simultaneously.

In practice, the extreme region of interest might contain few or even no
observed data points, especially when the dimension $d$ is large. One
potential solution is to employ appropriate parametric models for multivariate
extremes and to estimate the model parameters using intermediate level
observations. However, as the dimension of the vector $\bm{X}$ increases, such
a parametric approach becomes increasingly complex, which leads to the
following dilemma. On the one hand, models with a small number of parameters
tend to oversimplify the tail dependence structure. On the other hand, models
with many parameters are difficult to estimate and may even lead to
overfitting.

This paper provides a parametric model for the case where the underlying
random vector is high dimensional. The model can capture the so-called tail
dependence coefficients for all pairs of components of the $d-$dimensional
random vector. Even though it has a large number of parameters (up to
$d(d+1)/2$), they are easily obtained based on a fast algorithm, allowing one
to approximate the failure probability using the estimated model. Other
estimators of (multivariate) failure probabilities are proposed in Drees and
de Haan (2015) (for the bivariate case) and Valk (2016), among others.

In the existing literature, tail dependence of $\bm{X}$ is characterized in
various ways, such as by the exponent measure, the angular measure, and the
stable tail dependence function, among others (Beirlant et al., 2004; de Haan
and Ferreira, 2006). While such objects fully characterize the extremal
behavior of a random vector, their estimation in high dimensions suffers from
the “curse of dimensionality”. Instead, one may consider a partial summary of
tail dependence by focusing on the tail dependence coefficients $\chi$ or the
extremal coefficients $\theta$ for all pairs, see e.g. Coles, Heffernan and
Tawn (1999). This is similar to summarizing the dependence structure of a
random vector by the covariance or correlation matrix.

In this paper, we focus on a generalization of the so-called _tail pairwise
dependence matrix_ (TPDM), which collects all pairwise tail dependence
coefficients. This dependence matrix was originally proposed by Larsson and
Resnick (2012) and has gained in popularity since the work of Cooley and
Thibaud (2019); for example, Klüppelberg and Krali (2021) use it to estimate
the parameters of an extreme Bayesian network, Fomichov and Ivanovs (2020)
propose a clustering algorithm inspired by spherical $k$-means, where the
means are replaced by the Perron–Frobenius eigenvectors of the TPDM, and Kim
and Kokoszka (2022) study the extremal dependence between vectors of scores of
functional data.

The TPDM is particularly appropriate when $\bm{X}$ is multivariate regularly
varying with tail index $\alpha=2$. By contrast, our approach allows for
multivariate regularly variation with a general tail index $\alpha$, without
standardizing the data to $\alpha=2$. Since marginal standardization does not
affect dependence modeling, it is common in practice to separate marginal and
dependence modeling by first estimating the marginal distributions
(parametrically or non-parametrically), standardizing each marginal to a
common distribution, and then modeling the dependence structure only. However,
such an approach can be suboptimal in some applications. For example, in the
portfolio risk example, if one intends to estimate the Value-at-Risk or the
Expected Shortfall of a portfolio loss $L=\sum_{j=1}^{d}v_{j}X_{j}$, where
$v_{1},\ldots,v_{d}$ represent weights of the stocks (see, e.g., Mainik and
Embrechts, 2013), the failure region is of the form
$\\{\sum_{j}v_{j}X_{j}>x\\}$ with $x$ large. In this case, we need to directly
model the original vector $\bm{X}$: a nonparametric transformation for the
marginals will lead to untractable failure regions. Our approach without
standardizing a priori is therefore necessary. Such failure regions are also
of interest in climate applications, where $\sum_{j=1}^{d}v_{j}X_{j}$
represents the spatial average over a given region of some environmental
variable (see, e.g. Kiriliouk and Naveau, 2020).

The tail pairwise dependence matrix is closely related to the max-linear
model, a simple but rather flexible multivariate extreme-value model. In such
a model, each component of a $d$-dimensional vector can be interpreted as the
maximum shock among a set of $q$ independent heavy-tailed factors. The max-
linear model has $d\times q$ parameters (or $d\times(q-1)$ when margins are
standardized), usually stocked in a parameter matrix $A\in\mathbb{R}^{d\times
q}$. Max-linear models are, among other applications, used to model extremes
on directed acyclic graphs (Gissibl and Klüppelberg, 2018) or latent factor
structures for financial returns (Cui and Zhang, 2018). If the number of
parameters is limited, one can use the minimum distance estimation (Einmahl,
Kiriliouk and Segers, 2018) or spherical $k$-means (Janßen and Wan, 2020) to
estimate a max-linear model. Kiriliouk (2020) presents hypothesis tests aiding
in deciding the size of $q$.

Theoretically, Fougeres, Mercadier and Nolan (2013) showed that the max-linear
model is dense in the class of $d$-dimensional multivariate extreme-value
distributions. Moreover, Cooley and Thibaud (2019) showed that there exists a
finite $q\in\mathbb{N}$ such that the tail pairwise dependence matrix
$\Sigma_{\bm{X}}$ of any multivariate regularly varying random vector $\bm{X}$
with $\alpha=2$ is equal to that of a max-linear model with $q$ factors via
the relation $\Sigma_{\bm{X}}=AA^{T}$. In practice, the parameter matrix $A$
of this max-linear model can be estimated through a completely positive
decomposition of $\Sigma_{\bm{X}}$. A matrix $S$ is called _completely
positive_ if there exists a non-negative matrix $A$ such that $S=AA^{T}$;
finding this matrix $A$ amounts to finding a _completely positive
decomposition_ of $S$. Positive decompositions are well studied in the linear
algebra literature (Groetzner and Dür, 2020; Boţ and Nguyen, 2021), where it
is shown that the maximum number of columns of $A$ is $d(d+1)/2-4$. The
positive decomposition can be obtained via a decomposition algorithm,
nevertheless, the algorithm is computationally intensive.

In this paper, we propose an algorithm for an _approximate_ completely
positive decomposition $\Sigma_{\bm{X}}\approx AA^{T}$ with
$A\in\mathbb{R}^{d\times d}$. The decomposition is easy to compute and
applicable to moderate to high dimensions. In addition, it is fast enough to
obtain various decompositions corresponding to a single TPDM and hence to
characterize the model uncertainty. We show that the decomposition is exact
when $\bm{X}$ follows a max-linear model with a triangular coefficient matrix
$A$. Moreover, in all other cases, the decomposition is exact for all off-
diagonal entries of $\Sigma_{\bm{X}}$, but may overestimate its diagonal
entries. We apply the proposed decomposition algorithm to industry portfolio
returns and maximal wind speeds to illustrate its applicability. A third
application, concerning rainfall in Switzerland, has been deferred to Appendix
B.

##  2 Background

###  2.1 Multivariate regular variation and the tail pairwise dependence
matrix

Let $\bm{X}\in[0,\infty)^{d}$ be a random vector with joint cumulative
distribution function $F$ and continuous marginal distributions
$F_{1},\ldots,F_{d}$. Write
$\mathbb{E}_{0}=[0,\infty]^{d}\setminus\\{\bm{0}\\}$. We assume that $\bm{X}$
is _multivariate regularly varying_ , i.e., there exists a sequence
$b_{n}\rightarrow\infty$ and a limit measure $\nu_{\bm{X}}$ such that

$n\,\mathbb{P}\left[b_{n}^{-1}\bm{X}\in\,\cdot\,\right]\overset{v}{\rightarrow}\nu_{\bm{X}}(\,\cdot\,),\qquad\text{as
}n\rightarrow\infty,$ (1)

where $\overset{v}{\rightarrow}$ denotes vague convergence in the space of
non-negative Radon measures on $\mathbb{E}_{0}$ and $\nu_{\bm{X}}$ is called
the _exponent measure_ (Resnick, 2007, Chapter 6). Then there exists an
$\alpha>0$, called the _tail index_ of $\bm{X}$, such that the limit measure
$\nu_{\bm{X}}$ satisfies the homogeneity property of order $\alpha$, i.e.,
$\nu_{\bm{X}}(tC)=t^{-\alpha}\nu_{\bm{X}}(C)$ for any constant $t>0$ and any
($\nu_{\bm{X}}$-measurable) Borel set $C\subset\mathbb{E}_{0}$. The sequence
$b_{n}$ can be shown to satisfy $b_{n}\sim L(n)n^{1/\alpha}$, where $L$
denotes a slowly varying function.

The homogeneity property suggests that it is possible to express $\bm{X}$ and
correspondingly the limit relation in (1) in polar coordinates. Let
$\lVert{\,\cdot\,}\rVert$ denote any norm on $\mathbb{R}^{d}$ and consider the
transformation $T:\mathbb{E}_{0}\mapsto(0,\infty)\times\mathbb{S}_{d-1}$,
where $\mathbb{S}_{d-1}=\\{\bm{w}\in\mathbb{E}_{0}:\lVert{\bm{w}}\rVert=1\\}$,
and

$T(\bm{X})=\left(\lVert{\bm{X}}\rVert,\frac{\bm{X}}{\lVert{\bm{X}}\rVert}\right)=:(R,\bm{W}).$

Then (1) is equivalent to the existence of a finite measure $H_{\bm{X}}$
(called the _angular measure_) on $\mathbb{S}_{d-1}$ and a sequence
$b_{n}\rightarrow\infty$ such that

$n\,\mathbb{P}\left[\left(b_{n}^{-1}R,\bm{W}\right)\in\,\cdot\,\right]\overset{v}{\rightarrow}\mu_{\alpha}\times
H_{\bm{X}}(\,\cdot\,),\qquad\text{as }n\rightarrow\infty,$ (2)

in the space of non-negative Radon measures on
$(0,\infty]\times\mathbb{S}_{d-1}$, where $\mu_{\alpha}$ is a measure on
$[0,\infty)$ given by $\mu_{\alpha}\left((x,\infty]\right)=x^{-\alpha}$ for
$x>0$ (Resnick, 2007, Theorem 6.1). Let $B\subset\mathbb{S}_{d-1}$ be
$H_{\bm{X}}$-measurable and consider the set

$C_{s,B}:=\\{\bm{x}\in\mathbb{E}_{0}:\lVert{\bm{x}}\rVert>s,\bm{x}/\lVert{\bm{x}}\rVert\in
B\\}.$

Equation (2) implies $\nu_{\bm{X}}(C_{s,B})=s^{-\alpha}H_{\bm{X}}(B)$ and
$\nu_{\bm{X}}(\,\textnormal{d}r\times\,\textnormal{d}\bm{w})=\alpha
r^{-\alpha-1}\,\textnormal{d}r\,\textnormal{d}H_{\bm{X}}(\bm{w})$. The angular
measure $H_{\bm{X}}$ is related to the probability measure $S_{\bm{X}}$ in
Resnick (2007, Chapter 6) through $H_{\bm{X}}=mS_{\bm{X}}$ for
$m:=\nu(C_{1,\mathbb{S}_{d-1}})=H_{\bm{X}}(\mathbb{S}_{d-1})>0$. The sequence
$b_{n}$ in (2) is determined by

$n\,\mathbb{P}[R>b_{n}x]\rightarrow mx^{-\alpha},\qquad\text{ as
}n\rightarrow\infty,$

and the marginals of $\bm{X}$ satisfy

$\lim_{n\rightarrow\infty}n\,\mathbb{P}\left[X_{j}>b_{n}x\right]=x^{-\alpha}\int_{\mathbb{S}_{d-1}}w_{j}^{\alpha}\,\textnormal{d}H_{\bm{X}}(\bm{w}).$
(3)

To summarize the dependence of a multivariate regularly-varying vector
$\bm{X}$, we define the _tail pairwise dependence matrix (TPDM)_
$\Sigma_{\bm{X}}$ as

$\Sigma_{\bm{X}}=(\sigma_{\bm{X}_{jk}})_{j,k=1,\ldots,d},\qquad\text{with
}\,\,\,\sigma_{\bm{X}_{jk}}=\int_{\mathbb{S}_{d-1}}w_{j}^{\alpha/2}w_{k}^{\alpha/2}\,\textnormal{d}H_{\bm{X}}(\bm{w}).$
(4)

Note that this constitutes a generalization of the matrix introduced in
Larsson and Resnick (2012) whose integrand does not depend on $\alpha$; see
also Samorodnitsky and Taqqu (1994). We will write $\bm{X}\sim\Sigma_{\bm{X}}$
if $\Sigma_{\bm{X}}$ is the TPDM of $\bm{X}$. When no confusion can arise, we
will write the elements of $\Sigma_{\bm{X}}$ as $\sigma_{jk}$. Two variables
$X_{j},X_{k}$ are _tail dependent_ if and only if $\sigma_{jk}>0$. From (3) we
observe that

$\lim_{n\rightarrow\infty}n\,\mathbb{P}\left[X_{j}>b_{n}x\right]=x^{-\alpha}\sigma_{{jj}}.$

Hence, the (limiting) scale of $X_{j}$ is
$\left(\sigma_{{jj}}\right)^{1/\alpha}$, and we say that the margins of
$\bm{X}$ are _standardized_ with index $\alpha$ whenever $\sigma_{jj}=1$ for
all $j=1,\ldots,d$.

Finally, let $\lVert{\,\cdot\,}\rVert_{\alpha}$ denote the $L_{\alpha}$ norm,
i.e., for $\bm{x}\in[0,\infty)^{d}$,
$\lVert{\bm{x}}\rVert_{\alpha}=\left(\sum_{j=1}^{d}\lvert
x_{j}\rvert^{\alpha}\right)^{1/\alpha}$ for $\alpha\geq 1$. In the following,
we will always use the $L_{\alpha}$ norm in the definition of the angular
measure $H_{\bm{X}}$. The total mass of the angular measure then equals

$m=H_{\bm{X}}(\mathbb{S}_{d-1})=\int_{\mathbb{S}_{d-1}}\lVert{\bm{w}}\rVert_{\alpha}^{\alpha}\,\textnormal{d}H_{\bm{X}}(\bm{w})=\sum_{j=1}^{d}\sigma_{{jj}}.$

######  Proposition 2.1 (Generalization of Cooley and Thibaud (2019)).

The tail pairwise dependence matrix $\Sigma_{\bm{X}}$, as defined in (4), is
positive semi-definite.

###  2.2 Max-linear model

We are interested in calculating the probability that $\bm{X}$ falls into some
extreme failure region $C\subset\mathbb{E}_{0}$. In practice, to estimate such
probabilites $\mathbb{P}[\bm{X}\in C]$, we need to assume a parametric model
for the exponent measure $\nu$, or, equivalently, for the angular measure $H$.
A simple yet effective parametric model for the goal of estimating failure
probabilities is the _max-linear model_. Let $A$ denote a $d\times q$ matrix
with non-negative entries, with columns $\bm{a}_{1},\ldots,\bm{a}_{q}$, and
suppose that $\max(\bm{a}_{l})=\max(a_{1l},\ldots,a_{dl})>0$ for all
$l\in\\{1,\ldots,q\\}$. Let $Z_{1},\ldots,Z_{q}$ be independent random
variables whose distribution is Fréchet with scale 1 and shape $\alpha$,
$\mathbb{P}[Z_{l}\leq z]=\exp\left\\{-z^{-\alpha}\right\\}$, and define

$\displaystyle
A\times_{\max}\bm{Z}:=\left(\max_{l=1,\ldots,q}a_{1l}Z_{l},\ldots,\max_{l=1,\ldots,q}a_{dl}Z_{l}\right)^{T},$

for $\bm{Z}=(Z_{1},\ldots,Z_{q})$. Denoting $\bm{Y}=A\times_{\max}\bm{Z}$ and
choosing $b_{n}=n^{1/\alpha}$, the angular measure of $\bm{Y}$ is

$H_{\bm{Y}}(\,\cdot\,)=\sum_{l=1}^{q}\lVert{\bm{a}_{l}}\rVert_{\alpha}^{\alpha}\,\delta_{\bm{a}_{l}/\lVert{\bm{a}_{l}}\rVert_{\alpha}}(\,\cdot\,).$
(5)

Note that the order of the columns of $A$ does not matter. The marginal
distribution of $Y_{j}$ is Fréchet with scale
$\left(\sum_{l=1}^{q}a_{jl}^{\alpha}\right)^{1/\alpha}$ and shape $\alpha$ for
$j=1,\ldots,d$. The TPDM of $\bm{Y}$ has elements

$\sigma_{\bm{Y}_{jk}}=\sum_{l=1}^{q}\lVert{\bm{a}_{l}}\rVert_{\alpha}^{\alpha}\left(\frac{a_{jl}}{\lVert{\bm{a}_{l}}\rVert_{\alpha}}\right)^{\alpha/2}\left(\frac{a_{kl}}{\lVert{\bm{a}_{l}}\rVert_{\alpha}}\right)^{\alpha/2}=\sum_{l=1}^{q}a_{jl}^{\alpha/2}a_{kl}^{\alpha/2},$

and hence

$\displaystyle\Sigma_{\bm{Y}}=A_{*}A^{T}_{*},\qquad\text{where
}\,A_{*}:=\left(a_{jk}^{\alpha/2}\right)_{j,k=1,\ldots,d}$

Cooley and Thibaud (2019) show that the class of max-linear angular measures,
as defined in (5) is dense in the class of possible angular measures (for any
norm), i.e., any multivariate regularly varying vector $\bm{X}$ can be
approximated by a max-linear model with a $(d\times q)$ parameter matrix as
$q\rightarrow\infty$. Moreover, if attention is restricted to the TPDM, a max-
linear model with _finite_ $q$ is sufficient to exactly match
$\Sigma_{\bm{X}}$; this result has been proven in Cooley and Thibaud (2019)
for $\alpha=2$, but can be easily generalized to any $\alpha>0$.

######  Proposition 2.2 (Cooley and Thibaud (2019), Propositions 4 & 5).

  1. 1.

Let $Z_{1},\ldots,Z_{q}$ be Fréchet$(1,\alpha)$ random variables. Given any
angular measure $H$, there exists a sequence of non-negative $(d\times q)$
matrices $\\{A_{q}\\}$, for $q=1,2,\ldots$ such that
$H_{A_{q}\times_{\max}\bm{Z}}\rightarrow H$ weakly as $q\rightarrow\infty$.

  2. 2.

If $\bm{X}\in\mathbb{R}^{d}$ has TPDM $\Sigma_{\bm{X}}$, there exists a non-
negative $(d\times\widetilde{q})$ matrix $A_{*}$ with $\widetilde{q}<\infty$,
such that $\Sigma_{\bm{X}}=A_{*}A_{*}^{T}=\Sigma_{A\times_{\max}\bm{Z}}$.

Whenever a matrix $\Sigma$ can be decomposed as $\Sigma=AA^{T}$, and the
matrix $A$ has non-negative entries, we say that $\Sigma$ is _completely
positive_. The second part of Proposition 2.2 not only shows that the TPDM is
completely positive, but also that for any vector $\bm{X}$, there exists a
max-linear vector $\bm{Y}$ whose TPDM is equal to that of $\bm{X}$. This
result is key for the estimation of failure probabilites.

###  2.3 The probability of a failure region under the max-linear model

We intend to estimate the failure probability $\mathbb{P}[\bm{Y}\in C]$ for a
max-linear vector $\bm{Y}$. Indeed, whenever the set $C\subset\mathbb{E}_{0}$
is “extreme enough”, we have

$\mathbb{P}\left[\bm{Y}\in
C\right]=\frac{1}{n}\left\\{n\,\mathbb{P}\left[n^{-1/\alpha}\bm{Y}\in
n^{-1/\alpha}C\right]\right\\}\approx\frac{1}{n}\nu_{\bm{Y}}\left(n^{-1/\alpha}C\right).$
(6)

The choice of the region $C$ in (6) depends on the application of interest.
Nevertheless, for some practically relevant failure regions, the probability
can be calcuated directly from the parameters of the max-linear model. For
environmental data, e.g., if $\bm{Y}$ represents a vector of maximum daily
wind speeds, rainfall amounts, or temperatures in $d$ locations, the region

$\displaystyle C_{(\max)}(\bm{x})$
$\displaystyle=\\{\bm{y}\in\mathbb{E}_{0}:y_{1}>x_{1}\text{ or }\ldots\text{
or }y_{d}>x_{d}\\},$

is a common choice if an extreme event at at least one location can lead to a
climatological catastrophe. The exponent measure for this region equals

$\displaystyle\nu(C_{(\max)}(\bm{x}))=\sum_{l=1}^{q}\max_{j=1,\ldots,d}\left(\frac{a_{jl}}{x_{j}}\right)^{\alpha}.$

On the other hand, if $\bm{Y}$ represents for example the negative log-returns
of $d$ stocks in a portfolio, one might be interested in estimating the
probability of all stocks crashing together,

$\displaystyle C_{(\min)}(\bm{x})$
$\displaystyle=\\{\bm{y}\in\mathbb{E}_{0}:y_{1}>x_{1}\text{ and }\ldots\text{
and }y_{d}>x_{d}\\}.$

Similarly to $C_{(\max)}$, we can calculate

$\displaystyle\nu(C_{(\min)}(\bm{x}))=\sum_{l=1}^{q}\min_{j=1,\ldots,d}\left(\frac{a_{jl}}{x_{j}}\right)^{\alpha}.$

Many simple regions (including the ones above when $x_{1}=\ldots=x_{d}$) can
be written as follows. Let
$C_{f}(x)=\\{\bm{y}\in\mathbb{E}_{0}:f(\bm{y})>x\\}$ denote a failure region
where $f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ is some function. For fixed $w$,
let $r_{*}=r_{*}(\bm{w})$ be such that $f(r_{*}\bm{w})=x$. Then

$\nu(C_{f}(x))=\sum_{l=1}^{q}\lVert{\bm{a}_{l}}\rVert^{\alpha}_{\alpha}\left\\{r_{*}(\bm{a}_{l}/\lVert{\bm{a}_{l}}\rVert_{\alpha})\right\\}^{-\alpha}.$

Hence, we can obtain an explicit expression for the exponent measure of any
failure region that can be written in the form of $C_{f}$. For example, if
$f(\bm{y})=\bm{v}^{T}\bm{y}$ for $\bm{v}\in\mathbb{E}_{0}$ such that
$v_{1}+\ldots+v_{d}=1$, we get

$\displaystyle C_{(\textnormal{sum})}(\bm{v},x)$
$\displaystyle=\\{\bm{y}\in\mathbb{E}_{0}:v_{1}y_{1}+\ldots+v_{d}y_{d}>x\\}.$

Such a region is both of interest for climate events (e.g. in attribution
studies such as Kiriliouk and Naveau (2020), among others) and in the
financial context (e.g. to calculate the Value-at-Risk for a portfolio), and
we calculate

$\displaystyle\nu\left(C_{(\textnormal{sum})}(\bm{v},x)\right)=x^{-\alpha}\sum_{l=1}^{q}\left(\bm{v}^{T}\bm{a}_{l}\right)^{\alpha}.$

Moreover, we find

$\nu\left(C_{(\textnormal{sum})}(\bm{v},x)\right)=\begin{dcases}x^{-1}\sum_{j=1}^{d}v_{j}\sigma_{jj}&\text{if
}\alpha=1,\\\ x^{-2}\left(\bm{v}^{T}\Sigma_{\bm{X}}\bm{v}\right)&\text{if
}\alpha=2.\end{dcases}$

In these two special cases, $\nu\left(C_{(\textnormal{sum})}(\bm{v},x)\right)$
can be calculated based on the TPDM only. In other words, the random variable
$\bm{v}^{T}\bm{Y}$ approximately follows a Pareto distribution whose scale is
a function of the TPDM. By contrast, when $\alpha\notin\\{1,2\\}$, the
parameter matrix $A$ needs to be recovered first.

##  3 Decomposing the pairwise tail dependence matrix

Proposition 2.2 shows that, for any random vector $\bm{X}$, a max-linear
vector $\bm{Y}$ exists such that $\Sigma_{\bm{X}}=\Sigma_{\bm{Y}}$. Obtaining
such a model requires a _completely positive decomposition_ of
$\Sigma_{\bm{X}}$, i.e., finding a $(d\times q)$ matrix $A_{*}$ with non-
negative entries such that $\Sigma_{\bm{X}}=A_{*}A_{*}^{T}$ (see also Section
2.2). Barioli and Berman (2003) show that the number of columns of such a
matrix $A_{*}$ is at most $d(d+1)/2-4$ for $d\geq 5$. Current algorithms
performing completely positive decompositions are computationally heavy
(Groetzner and Dür, 2020; Boţ and Nguyen, 2021), and hence they cannot be
applied repeatedly to obtain multiple factorizations when the dimension $d$ is
high. Notice that the factorization of a completely positive matrix is not
unique. Obtaining multiple different decompositions could be useful to
quantify the model uncertainty regarding estimations of the probability of
falling into a failure region.

In this section, we propose a simple algorithm to construct a $(d\times d)$
matrix $A_{*}$, which leads to an _approximate_ completely positive
decomposition in the following sense. Firstly, we define the notation
“$\preceq$”: for two TPDMs, we write $\Sigma_{\bm{X}}\preceq\Sigma_{\bm{Y}}$
if

$\sigma_{\bm{X}_{jk}}=\sigma_{\bm{Y}_{jk}}\,\,\,\,\,\text{for all
}j,k\in\\{1,\ldots,d\\},\,j\neq k,\quad\text{and
}\,\sigma_{\bm{X}_{jj}}\leq\sigma_{\bm{Y}_{jj}}\,\,\,\,\,\text{for all
}j\in\\{1,\ldots,d\\}.$

Then, for any random vector $\bm{X}$ with TPDM $\Sigma_{\bm{X}}$, we aim at
constructing a max-linear vector $\bm{Y}$ with a $(d\times d)$ coefficient
matrix $A_{*}$, such that $\Sigma_{\bm{X}}\preceq\Sigma_{\bm{Y}}$. The max-
linear model can be viewed as an approximate model for the original random
vector $\bm{X}$, “matching” the pairwise dependence structure of $\bm{X}$.

The details of the proposed algorithm are as follows. Define, for any TPDM
$\Sigma_{\bm{X}}=(\sigma_{{jk}})_{j,k=1,\ldots,d}$ and any
$i\in\\{1,\ldots,d\\}$,

$D_{i}(\Sigma_{\bm{X}}):=\max\left\\{j,k\in\\{1,\ldots,d\\}\setminus\\{i\\}:\,\frac{\sigma_{ji}\sigma_{ki}}{\sigma_{jk}\sigma_{ii}}\right\\}.$

We write $D_{i}=D_{i}(\Sigma_{\bm{X}})$ when no confusion can arise. Let
$\bm{\tau}_{i}=(\tau_{1,i},\ldots,\tau_{d,i})^{T}$ with

$\tau_{j,i}=\begin{dcases}\sigma_{ji}\left(\sigma_{ii}\max(D_{i},1)\right)^{-1/2}&\text{
if }j\neq i,\\\ \left(\sigma_{ii}\max(D_{i},1)\right)^{1/2}&\text{ if
}j=i.\end{dcases}$ (7)

Let
$\bm{\tau}_{-i}:=(\tau_{1,i},\ldots,\tau_{i-1,i},\tau_{i+1,i},\ldots,\tau_{d,i})^{T}\in\mathbb{R}^{d-1}$
and define

$\Sigma^{(i)}_{\bm{X}}:=\Sigma^{(-i,-i)}_{\bm{X}}-\bm{\tau}_{-i}\bm{\tau}_{-i}^{T}\quad\in\mathbb{R}^{(d-1)\times(d-1)},$
(8)

where $\Sigma^{(-i,-i)}_{\bm{X}}\in\mathbb{R}^{(d-1)\times(d-1)}$ is the
matrix obtained by removing the $i$-th row and column from $\Sigma_{\bm{X}}$.
The following lemma shows that $\Sigma^{(i)}_{\bm{X}}$ is still a valid TPDM.

######  Lemma 3.1.

For all $i\in\\{1,\ldots,d\\}$ with $D_{i}\neq 1$, the matrix
$\Sigma^{(i)}_{\bm{X}}$ in (8) is a TPDM since

  1. 1.

$\Sigma^{(i)}_{\bm{X}}>0$ componentwise;

  2. 2.

$\Sigma^{(i)}_{\bm{X}}$ is positive semi-definite.

The next proposition shows that if we can construct a $(d-1)$-dimensional
random vector with TPDM $\Sigma^{(i)}_{\bm{X}}$, we can then use this vector
to construct $\bm{Y}$ such that $\Sigma_{\bm{X}}\preceq\Sigma_{\bm{Y}}$. In
addition, the proposition illustrates the circumstances under which
$\Sigma_{\bm{X}}=\Sigma_{\bm{Y}}$ is achieved. The proofs of Lemma 3.1 and of
Proposition 3.2 are postponed to the appendix.

######  Proposition 3.2.

Let $\Sigma_{\bm{X}}=(\sigma_{jk})_{j,k=1,\ldots,d}\in\mathbb{R}^{d\times d}$
be a TPDM. Let $i\in\\{1,\ldots,d\\}$ and let $Z_{i}$ be a Fréchet random
variable with scale 1 and shape $\alpha$. Suppose there exists a random vector
$\bm{\tilde{Y}}_{-i}=(\tilde{Y}_{1},\ldots,\tilde{Y}_{i-1},\tilde{Y}_{i+1},\ldots,\tilde{Y}_{d})^{T}$
with TPDM $\Sigma^{(i)}_{\bm{X}}$. By defining

$Y_{j}:=\begin{dcases}\max\left(\tau_{j,i}^{2/\alpha}Z_{i},\tilde{Y}_{j}\right)&\text{
if }j\in\\{1,\ldots,d\\}\setminus\\{i\\},\\\ \tau_{i,i}^{2/\alpha}Z_{i}&\text{
if }j=i,\end{dcases}$

we have that $\Sigma_{\bm{X}}\preceq\Sigma_{\bm{Y}}$. In particular, if
$D_{i}(\Sigma_{\bm{X}})<1$, then $\Sigma_{\bm{Y}}=\Sigma_{\bm{X}}$.

Hence, we can repeatedly apply the procedure in expression (8) to
$\Sigma^{(i)}_{\bm{X}}$ to further decrease the dimension. This leads to an
iterative algorithm, throughout which we obtain the matrix $A_{*}$. First, let
$i_{1}\mapsto i_{2}\mapsto\ldots\mapsto i_{d}$ denote a _path_ , where
$(i_{1},\ldots,i_{d})$ is a permutation of $(1,\ldots,d)$. The order
determines which column of $\Sigma_{\bm{X}}$ will be treated first using
Proposition 3.2. For a given path, the iterative algorithm is as follows.

  * •

In the first step, we obtain the vector $\bm{\tau}_{i_{1}}$ from (7) by taking
$i=i_{1}$. We fill in the first column of the matrix $A_{*}$ with the obtained
vector $\bm{\tau}_{i_{1}}$ (as the order of the columns of $A_{*}$ does not
matter, we do not necessarily need to fill the $i_{1}$-th column). The
targeted TPDM is then reduced to a $(d-1)\times(d-1)$ matrix
$\Sigma_{\bm{X}}^{(-i_{1})}$ as defined in (8).

  * •

In the second step, we apply expressions (7) and (8) to
$\Sigma_{\bm{X}}^{(-i_{1})}$, with the focal dimension now indexed by $i_{2}$.
This step results in a $(d-1)$-dimensional vector. For the second column of
the matrix $A_{*}$, we set the $i_{1}$-th row to zero and fill in the rest
with the obtained $(d-1)$-dimensional vector. After the second step, the
targeted TPDM is reduced to a $(d-2)\times(d-2)$ matrix.

  * •

The previous step is iterated as follows. In each step $j$, we aim at filling
in the $j$-th column of the matrix $A_{*}$. First, we set the elements in the
$i_{1},i_{2},\cdots,i_{j-1}$-th row to zero. Then, we fill in the other
$(d-j+1)$ elements by the obtained $(d-j+1)$-dimensional vector of this
iteration. After each step, the dimension of the targeted TPDM is reduced by
one.

  * •

At step $d$, the targeted TPDM is reduced to $1\times 1$, i.e. a single
positive number. We simply fill in the square root of that value to the $d$-th
column and $i_{d}$-th row of $A_{*}$, while setting all other elements in the
$d$-th column to zero.

Proposition 3.2 guarantees that the matrix $A_{*}$ obtained with this
algorithm has only non-negative elements and satisfies $\Sigma_{\bm{X}}\preceq
A_{*}A_{*}^{T}$. In addition, define $A=A_{*}^{2/\alpha}$ (i.e. each element
of the matrix $A_{*}$ is raised to the power $2/\alpha$) and let
$\bm{Z}=(Z_{1},\cdots,Z_{d})^{T}$ where $Z_{1},\ldots,Z_{d}$ are i.i.d.
standard Fréchet distributed random variables with scale 1 and shape $\alpha$.
The max-linear vector $\bm{Y}=A\times_{\max}\bm{Z}$ has a TPDM that satisfies
$\Sigma_{\bm{X}}\preceq\Sigma_{\bm{Y}}$. In other words, we can construct a
max-linear vector $\bm{Y}$ based on only $d$ independent factors to
approximate the TPDM of $\bm{X}$.

Below is an example for this procedure when $d=3$ and the chosen path is
$1\mapsto 2\mapsto 3$.

######  Example 3.3 ($d=3,\alpha=2$, path $1\mapsto 2\mapsto 3$).

Let $\Sigma_{\bm{X}}=(\sigma_{jk})_{j,k=1,2,3}$ be a known (or estimated)
TPDM. We want to find a coefficient matrix $A$ (since $\alpha=2$, we have
$A=A_{*}$) such that the corresponding max-linear model
$\bm{Y}=A\times_{\max}(Z_{1},Z_{2},Z_{3})^{T}$ has TPDM $\Sigma_{\bm{Y}}$
satisfying $\Sigma_{\bm{X}}\preceq\Sigma_{\bm{Y}}$. Starting with $i=1$, we
obtain from (7) that

$\begin{pmatrix}\tau_{1,1}\\\ \tau_{2,1}\\\
\tau_{3,1}\end{pmatrix}=\frac{1}{\sqrt{\sigma_{11}\max(D_{1},1)}}\begin{pmatrix}\sigma_{11}\max(D_{1},1)\\\
\sigma_{12}\\\ \sigma_{13}\end{pmatrix}$

and

$\Sigma_{\bm{X}}^{(-1)}:=\begin{pmatrix}s_{22}&s_{23}\\\ s_{23}&s_{33}\\\
\end{pmatrix}:=\begin{pmatrix}\sigma_{22}-\frac{\sigma_{12}^{2}}{\sigma_{11}\max(D_{1},1)}&\sigma_{23}-\frac{\sigma_{12}\sigma_{13}}{\sigma_{11}\max(D_{1},1)}\\\
\sigma_{23}-\frac{\sigma_{12}\sigma_{13}}{\sigma_{11}\max(D_{1},1)}&\sigma_{33}-\frac{\sigma_{13}^{2}}{\sigma_{11}\max(D_{1},1)}\\\
\end{pmatrix}.$

Hence, the first column of $A$ is given by
$(\tau_{1,1},\tau_{2,1},\tau_{3,1})^{T}$.

To find the second column of $A$, we iterate the same algorithm to the newly
obtained TPDM $\Sigma_{\bm{X}}^{(-1)}$ to get that

$\begin{pmatrix}\tau_{2,2}\\\
\tau_{3,2}\end{pmatrix}=\frac{1}{\sqrt{s_{22}}}\begin{pmatrix}s_{22}\\\
s_{23}\end{pmatrix}.$

Note that since $\Sigma_{\bm{X}}^{(-1)}$ is positive semidefinite, in this
case we have
$D_{2}(\Sigma_{\bm{X}}^{(-1)})=\frac{s_{23}^{2}}{s_{33}s_{22}}\leq 1$. After
this iteration, the targeted TPDM is reduced to a single value
$s_{33}-s_{23}^{2}/s_{22}\geq 0$. Hence, we can define
$\tau_{3,3}:=\sqrt{s_{33}-s_{23}^{2}/s_{22}}$.

Eventually, we can construct the max-linear model as
$\bm{Y}=A\times_{\max}(Z_{1},Z_{2},Z_{3})$ where

$A=\begin{pmatrix}\tau_{1,1}&0&0\\\ \tau_{2,1}&\tau_{2,2}&0\\\
\tau_{3,1}&\tau_{3,2}&\tau_{3,3}\\\
\end{pmatrix},\qquad\text{and}\qquad\Sigma_{\bm{Y}}=AA^{T}=\begin{pmatrix}\sigma_{11}\max(D_{1},1)&\sigma_{12}&\sigma_{13}\\\
\sigma_{12}&\sigma_{22}&\sigma_{23}\\\
\sigma_{13}&\sigma_{23}&\sigma_{33}\end{pmatrix}.$

We have that $\Sigma_{\bm{Y}}=\Sigma_{\bm{X}}$ if $D_{1}<1$ and
$\Sigma_{\bm{X}}\preceq\Sigma_{\bm{Y}}$ otherwise.

If the matrix $\Sigma_{\bm{X}}$ has a zero eigenvalue, it means that the
components of $\bm{X}$ are linearly related when the radius $R$ is high. In
practice, an eigenvalue of $\Sigma_{\bm{X}}$ can indeed be close to zero when
the variables are strongly dependent given a high radius. This is comparable
with the concept of multicolinearity in multivariate statistics. In that case,
at least one column of $A_{*}$ will be zero: the column corresponding to the
first variable encountered that is close to a linear combination of previous
variables, given that the radius is high. We can remove that column in $A_{*}$
and correspondingly $A$. Eventually, the number of non-zero columns of $A_{*}$
corresponds to the rank of $\Sigma_{\bm{X}}$.

In Example 3.3, we applied the iterative algorithm implied by Proposition 3.2
with the path $1\mapsto 2\mapsto 3$. The resulting matrix $A_{*}$ (and
eventually $A$) is then lower triangular. By considering others paths we will
obtain different matrices $A_{*}$. Hence, in general, for different choices of
paths, a completely positive decomposition (approximate or not) is not unique.
It is also possible that one path leads to an exact decomposition while
another path leads to an approximation in the sense of $\preceq$. Proposition
3.4 shows that if the vector $\bm{X}$ indeed follows a max-linear model
constructed from a lower-triangular parameter matrix $A$, by choosing the path
_appropriately_ , it is possible to recover the exact max-linear model. Given
that the order of columns of $A$ does not matter, the proposition can of
course be generalized to parameter matrices $A$ in which each column $j$ has
at least one more zero than column $(j-1)$, for $j=2,\ldots,d$.

######  Proposition 3.4.

Suppose that $\bm{X}=A\times_{\max}\bm{Z}$ follows a max-linear model with
lower-triangular parameter matrix $A\in\mathbb{R}^{d\times d}$ and
$\bm{Z}=(Z_{1},\ldots,Z_{q})$, where $Z_{1},\ldots,Z_{q}$ are iid Fréchet with
scale 1 and shape $\alpha$. If we choose the canonical path, i.e.,
$1\mapsto\ldots\mapsto d$ in the decomposition of $\Sigma_{\bm{X}}$ to obtain
the matrix $A_{*}$, then $A_{*}A_{*}^{T}=\Sigma_{\bm{X}}$.

For a high-dimensional TPDM $\Sigma_{\bm{X}}$, it is not feasible to exhaust
all $d!$ possible paths in search for an exact decomposition. In practice, we
need an efficient searching algorithm to find a path that leads to a max-
linear model $\bm{Y}$ such that $\Sigma_{\bm{Y}}=\Sigma_{\bm{X}}$, if
possible. We consider three potential algorithms to achieve the goal.

  1. 1.

A simple approach: In step $j$ we pick the next dimension $i_{j}$
corresponding to the lowest value of $D_{i}$ among all remaining dimensions.

  2. 2.

An exhaustive approach: In step $j$, we denote $T_{j}$ as the _set_ of
remaining dimensions corresponding to a value of $D_{i}$ below 1, i.e.
$T_{j}:=\\{i:D_{i}<1,i\neq i_{1},i_{2},\cdots,i_{j-1}\\}$. For each $i\in
T_{j}$, we proceed with the iterative process for $\Sigma^{(-i)}_{\bm{X}}$.
This will result in a tree of possibilities after each step. For some
“branches” of the tree, the algorithm may end at a “leaf” where no $D_{i}$
lies below 1, i.e. $T_{j}=\varnothing$. If some “branch” can iterate $d$
steps, it leads to a path that corresponds to an exact decomposition. We
retain all paths (and hence all possible matrices $A_{*}$) that lead to an
exact decomposition.

  3. 3.

A pragmatic approach: In step $j$, we denote $T_{j}$ as the _set_ of remaining
dimensions corresponding to a value of $D_{i}$ below 1, i.e.
$T_{j}:=\\{i:D_{i}<1,i\neq i_{1},i_{2},\cdots,i_{j-1}\\}$. We randomly choose
a dimension from this set as $i_{j}$ and proceed with the iteration. If the
procedure stops after less than $d$ iterations and ends with a scenario where
no $D_{i}$ lies below 1, we restart the entire procedure from the first
iteration. We keep “restarting” the algorithm until finding an exact
decomposition.

Note that the simple approach may not lead to an exact decomposition because
there is no guarantee that the minimum of the $D_{i}$’s is below $1$ in all
steps. The exhaustive approach will find all exact decompositions if there
exists at least one. However, it is limited to moderate dimensions $d$ because
of its computational intensity. The pragmatic approach does not guarantee
finding an exact decomposition either. Nevertheless, for most examples in
Section 5, an exact decomposition is quickly found after a few “restarts”, and
hence it is the approach that we will most use in practice.

##  4 Estimation

###  4.1 Estimation of the TPDM

Let $\bm{X}_{i}=(X_{i1},\ldots,X_{id})$ for $i=1,\ldots,n$ be iid copies of
$\bm{X}$. First, we estimate the tail indices $\alpha_{j}$, or, equivalently,
the extreme value indices $\gamma_{j}=1/\alpha_{j}$, for each $j=1,\ldots,d$
using the sample $\left\\{X_{1j},\ldots,X_{nj}\right\\}$. Note that the
assumption of regular variation (1) implies that all margins have equal tail
indices, $\alpha_{1}=\ldots=\alpha_{d}:=\alpha$, which is a common assumption,
since different tail indices will lead to the largest one dominating the
others and masking tail dependence.

Let $X_{(1),j}\leq\ldots\leq X_{(n),j}$ denote the order statistics of
$X_{1j},\ldots,X_{nj}$ and let $k=k_{n}$ be an intermediate sequence such that
$k\rightarrow\infty$ and $k/n\rightarrow 0$ as $n\rightarrow\infty$. The _Hill
estimator_ ,

$\widehat{\gamma}_{j}=\frac{1}{k}\sum_{i=1}^{k}\log\left(\frac{X_{(n-i+1),j}}{X_{(n-k),j}}\right)$

is the best-known estimator of the extreme-value index (Hill, 1975), and will
be used to obtain the estimated tail indices
$\widehat{\alpha}_{j}=1/\widehat{\gamma}_{j}$. If the confidence intervals of
$\widehat{\alpha}_{1},\ldots,\widehat{\alpha}_{d}$ are overlapping, a single
tail index estimate $\widehat{\alpha}$ might be justified. In that case, the
entire sample
$X_{11},\ldots,X_{n1},X_{12},\ldots,X_{n2},\ldots,X_{1d},\ldots,X_{nd}$ of
length $nd$ can be used to obtain the estimate $\widehat{\alpha}$ (Einmahl et
al., 2022).

Next, for $i=1,\ldots,n$, set
$R_{i}=\lVert{\bm{X}_{i}}\rVert_{\widehat{\alpha}}$ and
$\bm{W}_{i}=\bm{X}_{i}/R_{i}$. Let $m=H_{\bm{X}}(\mathbb{S}_{d-1})$ denote the
total mass of the angular measure. If the data are standardized to ensure
equal scales for all $j$ (see subsection 4.2), then $m=d$. Otherwise, we
observe that

$m=\nu_{\bm{X}}(C_{1,\mathbb{S}_{d-1}})=\lim_{n\rightarrow\infty}n\,\mathbb{P}[b_{n}^{-1}\bm{X}_{i}\in
C_{1,\mathbb{S}_{d-1}}]=\lim_{n\rightarrow\infty}n\,\mathbb{P}[R_{i}>b_{n}]$

Following Cooley and Thibaud (2019), we estimate $b_{n}$ by a high quantile of
the empirical distribution of $R_{1},\ldots,R_{n}$, denoted as $r_{0}$. By
setting $n\approx r_{0}^{\widehat{\alpha}}$, $m$ can be estimated by
$\widehat{m}=r_{0}^{\widehat{\alpha}}(n_{\textnormal{exc}}/n)$ for
$n_{\textnormal{exc}}=\sum_{i=1}^{n}\mathbbm{1}\left\\{R_{i}>r_{0}\right\\}$.
A value for $r_{0}$ could be chosen by plotting different quantiles against
the estimates of $m$ and searching for a region where estimates are stable.
Finally, since
$\sigma_{jk}=m\int_{\mathbb{S}_{d-1}}w_{j}^{\alpha/2}w_{k}^{\alpha/2}\,\textnormal{d}\widetilde{H}_{\bm{X}}(\bm{w})$
where $\widetilde{H}_{\bm{X}}(\,\cdot\,):=m^{-1}H_{\bm{X}}(\,\cdot\,)$ is a
probability measure, the elements of the TPDM can be estimated by

$\widehat{\sigma}_{jk}=\frac{\widehat{m}}{n_{\textnormal{exc}}}\sum_{i=1}^{n}W_{ij}^{\widehat{\alpha}/2}W_{ik}^{\widehat{\alpha}/2}\mathbbm{1}\left\\{R_{i}>r_{0}\right\\}=\frac{r_{0}^{\widehat{\alpha}}}{n}\sum_{i:R_{i}>r_{0}}\left(W_{ij}W_{ik}\right)^{\widehat{\alpha}/2},$

giving $\widehat{\Sigma}_{\bm{X}}=(\widehat{\sigma}_{jk})_{j,k=1,\ldots,d}$.
Larsson and Resnick (2012) show the asymptotic normality of
$\widehat{\sigma}_{jk}$ when $\alpha=2$ is known.

###  4.2 Estimation of failure probabilities and standardization

Once the TPDM is estimated, we calculate the $d\times d$ matrix
$\widehat{A}_{*}$ obtained by applying the decomposition algorithm outlined in
Section 3 to $\widehat{\Sigma}_{\bm{X}}$, and the empirical $(d\times
n_{\textnormal{exc}})$ estimate

$\widetilde{A}_{*}=\left(\frac{\widehat{m}}{n_{\textnormal{exc}}}\right)^{1/2}\widehat{W}^{\widehat{\alpha}/2},$

where $\widehat{W}^{\widehat{\alpha}/2}$ is the matrix whose columns are the
vectors $\bm{W}_{i}^{\widehat{\alpha}/2}$ for which $R_{i}>r_{0}$. The matrix
$\widetilde{A}_{*}$ is “inefficient” in the sense that it has much more
columns than $\widehat{A}_{*}$, especially when $n_{\textnormal{exc}}$ is
large. However, it is easy to calculate and satisfies
$\widehat{\Sigma}_{\bm{X}}=\widetilde{A}_{*}\widetilde{A}_{*}^{T}$; we will
use it as a benchmark. We obtain two estimators of $A$ as

$\widetilde{A}=\left(\widetilde{A}_{*}\right)^{2/\widehat{\alpha}}\quad\textnormal{and}\quad\widehat{A}=\left(\widehat{A}_{*}\right)^{2/\widehat{\alpha}}.$

Now, $\widetilde{A}$ and $\widehat{A}$ can be used to estimate the probability
of falling into one of the failure regions presented in Section 2.3. We will
denote the estimators of these probabilities by $\widehat{p}$ for
$\widehat{A}$ and $\widetilde{p}$ for $\widetilde{A}$. In the applications, we
will usually take the pragmatic approach described in Section 3 to generate
$N$ different matrices $\widehat{A}_{*1},\ldots,\widehat{A}_{*N}$ that
correspond to _exact_ decompositions of $\widehat{\Sigma}_{\bm{X}}$. By an
exact decomposition, we mean that
$\lVert\widehat{\Sigma}_{\bm{X}}-\widehat{A}_{*l}\widehat{A}^{T}_{*l}\rVert_{F}\leq
10^{-12}$ for some $l=1,\ldots,N$, where $\lVert\,\cdot\,\rVert_{F}$ denotes
the Frobenius norm.

Finally, note that for certain failure regions we could standardize the
margins of $\bm{X}_{1},\ldots,\bm{X}_{n}$ without losing tractability of the
shape of the failure set. We set

$\displaystyle X_{ij}^{\\#}$
$\displaystyle=h_{j}(X_{ij})=\left(-\log\widehat{F}_{j}^{P}(X_{ij})\right)^{-1/2},\qquad
j=1,\ldots,d,i=1,\ldots,n,$

where $\widehat{F}_{j}^{P}$ denotes the semi-parametric estimate of $F_{j}$,

$\widehat{F}_{j}^{P}(y)=\begin{dcases}\widehat{F}_{j}(y)&\text{ if }y\leq
u_{j};\\\
1-\left(1-\widehat{F}_{j}(u_{j})\right)\overline{H}_{\widehat{\sigma}_{j},\widehat{\gamma}_{j}}(y-u_{j})&\text{
if }y>u_{j};\\\ \end{dcases}$

where $\widehat{F}_{j}$ denotes the empirical distribution function of
$X_{1j},\ldots,X_{nj}$ and
$\overline{H}_{\widehat{\sigma}_{j},{\widehat{\gamma}}_{j}}$ the survival
function of a generalized Pareto distribution (GPD) with estimated scale
$\widehat{\sigma}_{j}$ and shape $\widehat{\gamma}_{j}$, obtained by fitting
$(X_{ij}-u_{j}\mid X_{ij}>u_{j})_{i=1,\ldots,n}$ to the GPD for some high
threshold $u_{j}$. The distribution of $X^{\\#}_{1j},\ldots,X^{\\#}_{nj}$ is
Fréchet with scale $1$ and shape $\alpha=2$. Hence, we have $\sigma_{jj}=1$
for $j\in\\{1,\ldots,d\\}$ and $m=H_{\bm{X}}(\mathbb{S}_{d-1})=d$ does not
need to be estimated. Subsequently, we can decompose the TPDM
$\widehat{\Sigma}_{\bm{X}^{\\#}}$ and construct the associated max-linear
model.

##  5 Applications

###  5.1 Failure probabilities computed from repeated decompositions

The TPDM forms an incomplete characterization of tail dependence; hence, we
would like to assess by how much the exponent measure of a given failure
region varies if it is computed from repeated decompositions. We proceed as
follows, inspired by Section 2 of the supplementary material of Cooley and
Thibaud (2019). Let

$\displaystyle A_{1}$
$\displaystyle=\begin{pmatrix}1.00&0.50&1.75&0.00&0.50&0.75&1.00&0.25\\\
2.00&0.00&1.00&1.50&0.25&1.00&1.00&1.00\\\
1.75&1.25&0.50&0.75&2.00&1.75&0.25&0.25\\\
1.25&0.25&1.25&2.00&2.00&0.50&0.25&0.25\\\
1.75&0.50&0.25&0.75&0.50&1.75&0.00&1.25\\\ \end{pmatrix},$ $\displaystyle
A_{2}$ $\displaystyle=\begin{pmatrix}1.00&0.00&0.50&0.00&0.50\\\
2.00&1.50&0.00&1.00&0.00\\\ 1.75&0.75&1.25&0.25&0.75\\\
1.25&2.00&0.25&0.25&0.25\\\ 1.75&0.75&0.50&1.25&0.50\\\
\end{pmatrix},\qquad\textnormal{and }A_{3}=\begin{pmatrix}1.00&0.00&0.00\\\
2.00&1.50&0.00\\\ 1.75&0.75&0.25\\\ 1.25&2.00&1.00\\\ 1.75&0.75&0.50\\\
\end{pmatrix}$

Matrix $A_{1}$ is as in Cooley and Thibaud (2019) (randomly generated
elements). Matrices $A_{2}$ and $A_{3}$ are similar but contain less columns
and more zeroes; hence, we expect to obtain more exact decompositions for
$A_{2}$ and $A_{3}$ than for $A_{1}$. Let $\nu_{i}$ denote the exponent
measure of a max-linear random vector with coefficient matrix $A_{i}$,
$i=1,2,3$ and $\alpha=4$. Let $\Sigma_{i}=A_{i}A_{i}^{T}$; we obtain the
$d!=5!=120$ decompositions of $\Sigma_{i}$ and consider the 120 max-linear
vectors with corresponding coefficient matrices for each $i=1,2,3$. From the
120 possible paths, we use only 94 paths for $A_{1}$, 86 for $A_{2}$, and 88
for $A_{3}$, since the others lead to paths with $D_{i}=1$; see Section 3.

We assess the limiting exponent measure of
$C_{f}(x)=\\{\bm{y}\in\mathbb{E}_{0}:f(\bm{y})>x\\}$ where $f$ is some
function; see Section 2.3. We use the functions
$f_{1}(\bm{y})=\frac{1}{d}\sum_{j=1}^{d}y_{j}$,
$f_{2}(\bm{y})=\prod_{j=1}^{d}y_{j}$, $f_{3}(\bm{y})=\min(y_{1},\ldots,y_{d})$
and $f_{4}(\bm{y})=\max(y_{1},\ldots,y_{d})$. We take $\alpha=4$ and $x$ is
always chosen such that $\nu_{j}(C_{f})=0.1$. Figure 1 shows boxplots for the
measures $\nu_{i}$, $i=1,2,3$, of the sets $C_{f_{1}},C_{f_{2}}$, $C_{f_{3}}$
and $C_{f_{4}}$, obtained through the different paths. We show both the
results when taking only exact decompositions into account (12 for $A_{1}$, 16
for $A_{2}$, and 24 for $A_{3}$) and when taking all decompositions such that
$\lVert{\Sigma_{\bm{X}}-AA^{T}}\rVert_{F}\leq 5$ (58 for $A_{1}$, 72 for
$A_{2}$, and 76 for $A_{3}$). For $A_{3}$, 24 decompositions contain only 3
columns and are exactly equal to $A_{3}$, while 40 other decompositions
contain 4 columns. For the product-region ($f_{2}$) and the min-region
($f_{3}$), only the first column of the obtained coefficient matrices is used
in the calculation of $\nu_{i}$ (since all subsequent columns contain at least
one zero; see also Section 2.3). Overall, the TPDM seems a reasonable summary
of tail dependence. We see that the tail dependence properties of the failure
regions based on $f_{1}$ and $f_{4}$ are better summarized through the TPDM
than those of the failure regions based on $f_{2}$ and $f_{3}$. This can be
explained by the fact that the regions $C_{f_{1}}$ and $C_{f_{4}}$ “cut” the
axes and hence concern more scenarios that the regions $C_{f_{2}}$ and
$C_{f_{3}}$ that only concern scenarios involving simultaneous extremes. In
general, variability increases when using approximate decompositions rather
than exact decompositions.

Figure 1: Boxplots of limiting exponent measures of exact and approximate
decompositions of $\Sigma_{1}$, $\Sigma_{2}$, $\Sigma_{3}$ and $\Sigma_{4}$
for the four functions $f_{1}$, $f_{2}$, $f_{3}$, and $f_{4}$, and $\alpha=4$.
The numbers in parentheses represent the number of decompositions on which the
boxplot is based.

###  5.2 Wind gusts in the Netherlands

Extreme windstorms are one of the important natural hazards affecting Europe.
Impacts are not primarily caused by wind speed itself but by the gusts, i.e.,
the maxima of the wind speed during a few seconds. In February 2022, storm
Eunice caused massive damages in Western, Central, and Northern Europe. In the
Netherlands, the maximum wind gust measured during storm Eunice was 144 km/h
at station Cabauw, which set the record for harshest inland wind to have ever
been measured in the Netherlands.

We consider daily maximal speeds of wind gusts, measured in km/h, observed at
$d=35$ weather stations in the Netherlands during extended winter
(October–March), from October 2001 up to and including March 2022. The data
set is freely available from the Royal Netherlands Meteorological Institute
(KNMI), https://climexp.knmi.nl/. Please note that up to 50 stations are
available in the Netherlands, but 15 were removed because of a large amount of
missing data. The remaining 35 stations have between $0$ and $15$ missing data
points out of $n=3827$; these do not occur around clusters of extreme wind
gusts and hence they do not impact the following analysis. During the extended
winter season, the KNMI issues a code red alarm (for any Dutch province) when
there is at least 90 % certitude that wind gusts will exceed 120 km/h . The
$35$ weather stations are initially divided into coastal stations (17) and
inland stations (18) based on the same criteria used in Lenderink, Van
Meijgaard and Selten (2009), because most environmental variables are known to
behave differently over land or close to the sea; see Figure 3 (left) for a
map of the stations used in this study.

We first estimate the marginal tail indices $\alpha_{j}=1/\gamma_{j}$ using
the standard Hill estimator with automated threshold selection (the so-called
“automated eye-ball method” in Danielsson et al. (2016), where we consider
$k\geq 15$). However, the Hill estimator is very erratic because of the many
ties present in the data (due to rounding). To overcome this issue, we
calculate the smoothed Hill estimator (Resnick and Stărică, 1997) before
applying the automated eye-ball method to choose a threshold. Figure 2 (right)
shows the estimates, with upper and lower bootstrap confidence limits for a
$95\%$ confidence level, obtained through a block bootstrap using the R
package boot (Canty and Ripley, 2020). The first 17 bars correspond to coastal
stations, the following 18 bars to inland stations. We see no apparent
differences between the tail heaviness of maximum wind gusts at the coast and
inland. The horizontal line shows the smoothed Hill estimate
$\widehat{\alpha}=9.289$ (where $k$ is again chosen following the automated
eye-ball method), based on pooling all data from the 35 weather stations.

The next step concerns the choice of $r_{0}$. Figure 2 (right) show estimates
of $m$, the total mass of the angular measure, as a function of $r_{0}$. We
continue with $r_{0}$ corresponding to the $95\%$ quantile and hence
$\widehat{m}=0.175$, corresponding to using the 192 largest values.

Figure 2: Left: estimates $\widehat{\alpha}_{j}$ with upper and lower $95\%$
confidence limits; the horizontal line represents the estimate
$\widehat{\alpha}$. Right: estimated mass $\widehat{m}$ of the angular measure
as a function of $r_{0}$, a high quantile of $R_{1},\ldots,R_{n}$; the
horizontal line represents the estimate for $r_{0}=0.95$.

Next, we estimate pairwise tail dependence coefficients non-parametrically to
examine the difference in pairwise dependence structure between the coastal
and the inland stations. Figure 3 shows both the pairwise tail dependence
coefficients $\chi_{jk}$ (middle) and the tail pairwise dependence matrix
entries $\sigma_{jk}$ (right), for $j,k=1,\ldots,d$. The plot suggests that
there is no significant difference between pairs of inland or pairs of coastal
stations, and that the wind gusts are asymptotically dependent, even for large
distances.

Figure 3: Left: map of the Netherlands with 17 coastal stations in brown and
18 inland stations in green. Middle: empirical coefficients of tail dependence
$\chi_{jk}$ for all pairs (black), pairs of coastal stations (blue), and pairs
of inland stations (red), with loess curves. Right: same but for the TPDM
entries $\sigma_{jk}$.

Finally, we estimate $A$ as described in Section 4.2; in total, we obtain
$N=200$ exact decompositions of $\widehat{\Sigma}_{\bm{X}}$. Based on
$\widetilde{A}$ and on $\widehat{A}_{1},\ldots,\widehat{A}_{200}$, we estimate
the probability $p_{(\max)}=\mathbb{P}[\bm{X}\in C_{\max}(\bm{x})]$ for
$\bm{x}=x\bm{1}_{d}$ and $x\in\\{1,1.2,1.44\\}$; the first two values are the
thresholds for a code orange and a code red alarm respectively (issued by the
KNMI during winter season), while the last value corresponds to the record
wind gusts measured during storm Eunice.

Figure 4 shows boxplots of the estimates based on
$\widehat{A}_{1},\ldots,\widehat{A}_{200}$, together with estimates based on
$\widetilde{A}$ and empirical estimates, for the three thresholds (from left
to right), on the original scale and the transformed scale. Upper panels
represent the analysis using coastal stations only, middle panels are for
inland stations only, while lower panels use the full dataset ($d=35$). While
in the univariate and pairwise analysis, no apparent difference was present
between coastal and inland stations, when estimating failure probabilities
based on the $17$\- and $18$-dimensional datasets, we observe that values for
coastal stations are almost five times as high than for inland stations when
performing the analysis on the original scale. Smaller differences are found
when using the standardized data $\bm{X}_{1}^{\\#},\ldots,\bm{X}_{n}^{\\#}$.
The estimates based on the full $35-$dimensional dataset are closer to the
coastal estimates than to the inland estimates. Finally, we see larger
discrepancy with the empirical estimates for the 144 km/h threshold, which is
due to the estimation uncertainty associated to such a high threshold.

Figure 4: Boxplots of estimates of failure probabilities for the region
$C_{\textnormal{max}}$ for the 17 coastal stations (upper panels), the 18
inland stations (middle panels), and all 35 stations combined (lower panels).

###  5.3 Financial industry portfolios

We focus on the daily negative value-averaged returns of $d=30$ industry
portfolios, downloaded from the Kenneth French Data Library. Data is available
for the period January 1970 until December 2019, leading to a sample size
$n=12\,613$.

We start by estimating the marginal tail indices $\alpha_{j}=1/\gamma_{j}$
using the standard Hill estimator with automated threshold selection
(Danielsson et al., 2016). Figure 5 (left) show the estimates, with upper and
lower bootstrap confidence limits for a $95\%$ confidence level. Overlapping
confidence intervals suggest that the hypothesis of a single index of regular
variation $\alpha$ is not rejected. The marginal estimates are rather similar,
and the straight line corresponding to the estimate $\widehat{\alpha}$, based
on all data combined, is within all confidence intervals. Hence, in the
following, we use $\widehat{\alpha}=3.479$.

Figure 5: Left: estimates $\widehat{\alpha}_{j}$ with upper and lower $95\%$
confidence limits; the horizontal line represents the estimate
$\widehat{\alpha}$. Right: estimated mass $\widehat{m}$ of the angular measure
as a function of $r_{0}$, a high quantile of $R_{1},\ldots,R_{n}$; the
horizontal line represents the estimate for $r_{0}=0.975$.

The next step concerns the choice of $r_{0}$. Figure 5 (right) show estimates
of $m$, the total mass of the angular measure, as a function of $r_{0}$. We
continue with $r_{0}$ corresponding to the $97.5\%$ quantile and hence
$\widehat{m}=29.99$, corresponding to using the 316 largest values.

Finally, we estimate $A$ as described in Section 4.2; in total, we obtain
$N=200$ exact decompositions of $\widehat{\Sigma}_{\bm{X}}$. Based on
$\widetilde{A}$ and on $\widehat{A}_{1},\ldots,\widehat{A}_{200}$, we estimate
the following quantities:

  1. 1.

the probability $p_{\textnormal{sum}}=\mathbb{P}[\bm{X}\in
C_{(\textnormal{sum})}(\bm{v},x)]=\mathbb{P}[\bm{v}^{T}\bm{X}>x]$, where $x$
corresponds to the empirical $0.995$ quantile of
$(\bm{v}^{T}\bm{X}_{i})_{i=1,\ldots,n}$; in other words, we validate the
empirical $99.5\%$ value-at-risk (VaR) associated with a weighted sum of
industry portfolio returns.

  2. 2.

the probability $p_{\textnormal{minsum}}=\mathbb{P}[\bm{X}\in C_{f_{1}}(x)]$
associated with $C_{f_{1}}(x)=\\{\bm{y}\in\mathbb{E}_{0}:f_{1}(\bm{y})>x\\}$
and

$f_{1}(\bm{y})=\min\left(\sum_{j=1}^{10}v_{j}y_{j},\sum_{j=11}^{20}v_{j}y_{j},\sum_{j=21}^{30}v_{j}y_{j}\right),$

where $x$ is the empirical $0.995$ quantile of
$f_{1}(\bm{X}_{1}),\ldots,f_{1}(\bm{X}_{n})$; i.e., we validate the empirical
$99.5\%$ VaR associated with the minimum of three weighted sums of industry
portfolio returns.

  3. 3.

the probability $p_{\textnormal{maxsum}}=\mathbb{P}[\bm{X}\in C_{f_{2}}(x)]$
associated with $C_{f_{2}}(x)=\\{\bm{y}\in\mathbb{E}_{0}:f_{2}(\bm{y})>x\\}$
and

$f_{2}(\bm{y})=\max\left(\sum_{j=1}^{10}v_{j}y_{j},\sum_{j=11}^{20}v_{j}y_{j},\sum_{j=21}^{30}v_{j}y_{j}\right),$

where $x$ is the empirical $0.995$ quantile of
$f_{2}(\bm{X}_{1}),\ldots,f_{2}(\bm{X}_{n})$.

For all three failure regions, we compare equal weights,
$\bm{v}=(1/d,\ldots,1/d)$, and unequal weights,
$\bm{v}=(0.3,1.7,1,0.3,1.7,1,\ldots,0.3,1.7,1)$. All probabilities can be
calculated analytically based on the estimates $\widetilde{A}$,
$\widehat{A}_{1},\ldots,\widehat{A}_{200}$; see Section 2.3.

Figure 6 shows boxplots of the estimates based on
$\widehat{A}_{1},\ldots,\widehat{A}_{200}$ and $\widetilde{A}$. Some deviation
is visible with respect to the empirical VaR, which could also be due to the
empirical VaR suffering from estimation uncertainty. We note that variability
is lowest for the sum-region, and highest for the minsum-region.

Figure 6: Boxplots of estimates of failure probabilities for the regions
$C_{(\textnormal{sum})}$, $C_{f_{1}}$, and $C_{f_{2}}$ for the financial
industry portfolios.

##  6 Discussion

We proposed a general method for the estimation of failure probabilities based
on pairwise measures of tail dependence. Our approach is fast and convenient
for moderate to high dimensions. When $d\geq 40$, it becomes increasingly hard
to find exact decompositions; hence, the (asymptotic) scales of the variables
are overestimated, and the approach is appropriate for certain classes of
failure regions only.

When summarizing the tail dependence structure through bivariate summary
measures, it is natural that some information is lost. However, the ability to
obtain repeated decompositions helps in quantifying such uncertainty. Another
source of error stems from the estimation of the tail index $\alpha$ and the
tail pairwise dependence matrix $\Sigma$. Although a multitude of estimators
exist for $\alpha$, estimation of $\Sigma$ has yet to be investigated further
and a better estimator could potentially improve our failure probability
estimates.

##  Appendix A Proofs

###### Proof of Proposition 2.1.

For a vector $\bm{X}\in\mathbb{R}^{d}$, let $m=H_{\bm{X}}(\mathbb{S}_{d-1})$
denote the total mass of the angular measure. Then
$\widetilde{H}_{\bm{X}}(\,\cdot\,):=m^{-1}H_{\bm{X}}(\,\cdot\,)$ is a
probability measure. Let $\bm{W}$ be a random vector that satisfies
$\mathbb{P}[\bm{W}\in B]=\widetilde{H}_{\bm{X}}(B)$ for
$B\subset\mathbb{S}_{d-1}$. Then
$m\,\mathbb{E}\left[\bm{W}^{\alpha/2}\left(\bm{W}^{\alpha/2}\right)^{T}\right]=\Sigma_{\bm{X}}$.
Finally, for $\bm{y}\in\mathbb{R}^{d}\setminus\\{\bm{0}\\}$,

$\bm{y}^{T}\Sigma_{\bm{X}}\bm{y}=m\,\mathbb{E}\left[\bm{y}^{T}\bm{W}^{\alpha/2}\bm{W}^{\alpha/2}\bm{y}\right]=m\,\mathbb{E}\left[\left(\bm{y}^{T}\bm{W}^{\alpha/2}\right)^{2}\right]\geq
0.$

∎

###### Proof of Lemma 3.1.

We have $\Sigma^{(i)}>0$ if and only if

$\sigma_{jk}-\frac{\sigma_{ji}\sigma_{ki}}{\sigma_{ii}\max(D_{i},1)}>0\iff\max(D_{i},1)>\frac{\sigma_{ji}\sigma_{ki}}{\sigma_{ii}\sigma_{jk}},\qquad\text{for
all }j,k\in\\{1,\ldots,d\\}\setminus\\{i\\},$

which is true by definition of $D_{i}$. Next, for
$\bm{x}_{-i}=(x_{1},\ldots,x_{i-1},x_{i+1},\ldots,x_{d})^{T}\in\mathbb{R}^{d-1}$
we want to show that $\bm{x}^{T}_{-i}\Sigma^{(i)}\bm{x}_{-i}\geq 0$. Let
$\bm{y}=(y_{1},\ldots,y_{d})^{T}$ with

$y_{j}=\begin{dcases}x_{j}&\text{ if }j\in\\{1,\ldots,d\\}\setminus\\{i\\};\\\
\frac{-\bm{x}^{T}_{-i}\bm{\sigma}_{-i}}{\sigma_{ii}\max(D_{i},1)}&\text{ if
}j=i;\end{dcases}$

where
$\bm{\sigma}_{-i}=(\sigma_{1,i},\ldots,\sigma_{i-1,i},\sigma_{i+1,i},\ldots,\sigma_{d,i})^{T}$,
the $i$-th column of $\Sigma$ with the $i$-th element removed. Because
$\Sigma$ is positive semi-definite, we have $\bm{y}^{T}\Sigma\bm{y}\geq 0$.
Observing that

$\displaystyle 0\leq\bm{y}^{T}\Sigma\bm{y}$
$\displaystyle=\bm{x}^{T}_{-i}\Sigma_{(-i,-i)}\bm{x}_{-i}+2y_{i}\left(\bm{x}^{T}_{-i}\bm{\sigma}_{-i}\right)+y_{i}^{2}\sigma_{ii}$
$\displaystyle=\bm{x}^{T}_{-i}\Sigma_{(-i,-i)}\bm{x}_{-i}-\frac{2\left(\bm{x}^{T}_{-i}\bm{\sigma}_{-i}\right)^{2}}{\sigma_{ii}\max(D_{i},1)}+\frac{\left(\bm{x}^{T}_{-i}\bm{\sigma}_{-i}\right)^{2}}{\sigma_{ii}\max(D_{i},1)^{2}}$
$\displaystyle=\bm{x}^{T}_{-i}\Sigma_{(-i,-i)}\bm{x}_{-i}+\left(\bm{x}^{T}_{-i}\tau_{-i}\right)^{2}\left(\frac{1}{\max(D_{i},1)}-2\right)$
$\displaystyle\leq\bm{x}^{T}_{-i}\Sigma_{(-i,-i)}\bm{x}_{-i}-\left(\bm{x}^{T}_{-i}\tau_{-i}\right)^{2}=\bm{x}^{T}_{-i}\Sigma^{(i)}\bm{x}_{-i},$

where the third equality is obtained by observing that
$\tfrac{\left(\bm{x}_{-i}^{T}\bm{\sigma}_{-i}\right)^{2}}{\sigma_{ii}\max(D_{i},1)}=\left(\bm{x}_{-i}^{T}\bm{\tau}_{-i}\right)^{2}$.
Hence, we have shown that $\Sigma^{(i)}$ is positive semi-definite. ∎

###### Proof of Proposition 3.2.

Suppose without loss of generality that $i=1$. Since $\bm{Y}_{-1}$ is a max-
linear random vector with TPDM $\Sigma^{(1)}$, for some
$q\in\\{2,3,\ldots\\}$, there exists a coefficient matrix
$A_{(-1,-1)}=(a_{jk})_{j=2,\ldots,d,k=2,\ldots,q}\in\mathbb{R}^{(d-1)\times(q-1)}$
and Fréchet$(1,\alpha)$ random variables $Z_{2},\ldots,Z_{q}$ such that

$\bm{Y}_{-1}=A_{(-1,-1)}\times_{\max}(Z_{2},\ldots,Z_{q})^{T},\qquad\text{and
}\,\Sigma^{(1)}=A^{\alpha/2}_{(-1,-1)}\left(A^{\alpha/2}_{(-1,-1)}\right)^{T}$

Then

$\bm{X}=\begin{pmatrix}\tau_{1,1}^{2/\alpha}Z_{1}\\\
\max(\tau_{2,1}^{2/\alpha}Z_{1},Y_{2})\\\ \vdots\\\
\max(\tau_{d,1}^{2/\alpha}Z_{1},Y_{d})\end{pmatrix}=A\times_{\max}\bm{Z},\qquad\text{
for
}A:=\left(\begin{array}[]{c@{}|@{}c}\tau_{1,1}^{2/\alpha}&\begin{matrix}&0&\cdots&0\end{matrix}\\\
\hline\cr\begin{matrix}\tau_{2,1}^{2/\alpha}\\\ \vdots\\\
\tau_{d,1}^{2/\alpha}\end{matrix}&\,\,\,A_{(-1,-1)}\end{array}\right),$

and $\bm{Z}:=(Z_{1},\ldots,Z_{q})^{T}$. We have

$\displaystyle\Sigma_{\bm{X}}=A^{\alpha/2}\left(A^{\alpha/2}\right)^{T}=\left(\begin{array}[]{c@{}|@{}c}\tau_{1,1}^{2}&\begin{matrix}&\tau_{1,1}\tau_{2,1}&\cdots&\tau_{1,1}\tau_{d,1}\end{matrix}\\\
\hline\cr\begin{matrix}\tau_{1,1}\tau_{2,1}\\\ \vdots\\\
\tau_{1,1}\tau_{d,1}\end{matrix}&\,\,\bm{\tau}_{-1}\bm{\tau}_{-1}^{T}+\Sigma^{(1)}\end{array}\right)=\begin{pmatrix}\sigma_{11}{\textstyle\max(D_{1},1)}&\sigma_{12}&\ldots&\sigma_{1d}\\\
\sigma_{12}&\sigma_{22}&\ldots&\sigma_{2d}\\\ \vdots&\vdots&\ddots&\vdots\\\
\sigma_{1d}&\sigma_{2d}&\ldots&\sigma_{dd}\end{pmatrix},$

by the definitions of $\Sigma^{(1)}$, $\bm{\tau}_{1}$ and $\bm{\tau}_{-1}$.
Hence, $\Sigma=\Sigma_{\bm{X}}$ if $D_{1}<1$ and
$\sigma_{11}<\sigma_{\bm{X}_{11}}$ if $D_{1}>1$. ∎

###### Proof of Proposition 3.4.

We prove the Proposition by mathematical induction with respect to $d$. First,
we prove for $d=2$. Since $A$ is lower triangular, so is the corresponding
matrix $A^{\alpha/2}$. For $d=2$, assume that

$\displaystyle A^{\alpha/2}=\begin{pmatrix}a_{11}&0\\\
a_{21}&a_{22}\end{pmatrix}.$

Then

$\displaystyle\Sigma_{\bm{X}}=\begin{pmatrix}a^{2}_{11}&a_{11}a_{21}\\\
a_{11}a_{21}&a_{21}^{2}+a^{2}_{22}\end{pmatrix}.$

By taking the path $1\mapsto 2$, we calculate

$\displaystyle
D_{1}(\Sigma_{\bm{X}})=\frac{(a_{11}a_{21})^{2}}{a^{2}_{11}(a_{21}^{2}+a^{2}_{22})}=\frac{a_{21}^{2}}{a_{21}^{2}+a^{2}_{22}}.$

Since $D_{1}(\Sigma_{\bm{X}})<1$, we find $\bm{\tau}_{1}=(a_{11},a_{21})^{T}$.
Correspondingly, we calculate that

$\Sigma_{\bm{X}}^{(-1)}=a_{21}^{2}+a^{2}_{22}-a_{21}^{2}=a_{22}^{2}.$

Therefore, the iterative algorithm leads to the matrix

$\displaystyle A_{*}=\left(\bm{\tau}_{1},\begin{matrix}0\\\
a_{22}\end{matrix}\right)=A^{\alpha/2}.$

The proposition is thus proven for $d=2$.

Next, assume that the Proposition is valid for $d-1$ with $d\geq 3$, we show
that it is also valid for $d$. Denote the lower triangular matrix
$A^{\alpha/2}$ as

$\displaystyle A^{\alpha/2}=\begin{pmatrix}a_{11}&0&\cdots&0\\\
a_{21}&a_{22}&\cdots&0\\\ \vdots&\vdots&\ddots&\vdots\\\
a_{n1}&a_{n2}&\cdots&a_{nn}\end{pmatrix}=:\begin{pmatrix}a_{11}&\bm{0}^{T}\\\
\bm{a}_{-1}&\tilde{\bm{A}}^{(-1,-1)}\end{pmatrix},$

where $\bm{0}$ is a $(d-1)-$ dimensional vector with all elements zero,
$\bm{a}_{-1}=(a_{21},\cdots,a_{n1})^{T}$ and $\tilde{{A}}^{(-1,-1)}$ is a
$(d-1)\times(d-1)$ lower-triangular matrix. Then

$\displaystyle\Sigma_{\bm{X}}=\begin{pmatrix}a^{2}_{11}&a_{1,1}\bm{a}_{-1}^{T}\\\
a_{11}\bm{a}_{-1}&\bm{a}_{-1}\bm{a}_{-1}^{T}+\tilde{\bm{A}}^{(-1,-1)}(\tilde{A}^{(-1,-1)})^{T}\end{pmatrix}.$

Because $\tilde{A}^{(-1,-1)}$ contains only non-negative elements, we get
$\sigma_{jk}>a_{j1}a_{k1}$ for any $j,k\neq 1$. In addition,
$\sigma_{j1}=a_{11}a_{j1}$. Since we consider the path $1\mapsto
2\mapsto\cdots\mapsto d$, we first calculate $D_{1}$; for any $j,k\neq 1$,

$\frac{\sigma_{j1}\sigma_{k1}}{\sigma_{jk}\sigma_{11}}<\frac{a_{11}a_{j1}a_{11}a_{k1}}{a_{j1}a_{k1}a_{11}^{2}}=1,$

which implies that $D_{1}<1$ and thus
$\bm{\tau}_{1}=(a_{11},a_{21},\cdots,a_{d1})^{T}$ following the algorithm
specified in (7). Correspondingly, from (8) we get

$\displaystyle\Sigma_{\bm{X}}^{(1)}=\Sigma_{\bm{X}}^{(-1,-1)}-\bm{\tau}_{-1}\bm{\tau}_{-1}^{T}$
$\displaystyle=\bm{a}_{-1}\bm{a}_{-1}^{T}+\tilde{{A}}^{(-1,-1)}(\tilde{{A}}^{(-1,-1)})^{T}-\bm{a}_{-1}\bm{a}_{-1}^{T}$
$\displaystyle=\tilde{{A}}^{(-1,-1)}(\tilde{{A}}^{(-1,-1)})^{T}.$

In other words, the $(d-1)\times(d-1)$ matrix $\Sigma_{\bm{X}}^{(1)}$ is of
the form $\tilde{{A}}^{(-1,-1)}(\tilde{{A}}^{(-1,-1)})^{T}$ where
$\tilde{{A}}^{(-1,-1)}$ is a lower-triangular matrix with non-negative
elements only. By invoking the assumption that the proposition holds for
$d-1$, we complete the proof of the induction step and conclude that the
proposition hold for any positive integer $d$. ∎

##  Appendix B Extreme precipitation in Switzerland

We focus on daily precipitation amounts in mm for June, July and August for
the years 1962–2012 at $d=44$ weather stations near Zürich, Switzerland,
leading to a sample size $n=4691$. This data was analyzed in Thibaud and Opitz
(2015) and Cooley and Thibaud (2019) and found to be stationary, independent
in time, and asymptotically dependent in space. Estimates of marginal tail
indices tend to vary so that a common index of regular variation $\alpha$ is
not justified; we standardize the data to
$\bm{X}_{1}^{\\#},\ldots,\bm{X}_{n}^{\\#}$. We take $r_{0}$ corresponding to
the $95\%$ quantile, yielding $235$ “extreme” observations. Finally, we search
for estimates of $A$ corresponding to exact decompositions of
$\widehat{\Sigma}_{\bm{X}}$; a very large number of exact decompositions
exists (but with identical paths for the first 14 iterations); we keep
$N=2000$ of them. Based on $\widetilde{A}$ and
$\widehat{A}_{1},\ldots,\widehat{A}_{N}$, we estimate

  1. 1.

the probability $p_{(\min)}=\mathbb{P}[\bm{X}\in C_{\min}(\bm{x})]$, where
$\bm{x}$ is chosen as the 30 mm threshold for precipitation on the original
scale;

  2. 2.

the probability $p_{(\max)}=\mathbb{P}[\bm{X}\in C_{\max}(\bm{x})]$, where
$\bm{x}$ is chosen as the 60 and 80 mm threshold for precipitation on the
original scale.

The first probability has also been estimated in Cooley and Thibaud (2019),
who found a value of $4.8\times 10^{-4}$. For this region, the $N$
decompositions give the same estimate of $\widehat{p}_{(\min)}=3.3\times
10^{-3}$ since only the first column of
$\widehat{A}_{1},\ldots,\widehat{A}_{N}$ matters. We can compare these to the
empirical estimate of $4.3\times 10^{-4}$ and the estimate based on
$\tilde{A}$ which gives $1.7\times 10^{-3}$. Figure 7 shows boxplots of the
estimates based on $\widehat{A}_{1},\ldots,\widehat{A}_{N}$, together with
estimates based on $\widetilde{A}$ and empirical estimates for $p_{\max}$. We
see that the estimates based on $\widetilde{A}$ underestimate the probability
of a heavy rainfall event.

Figure 7: Boxplots of estimates of failure probabilities $p_{\max}$ for the
the Swiss precipitation data, for a 60 mm threshold (left) and a 80 mm
threshold (right) on the original scale.

[Acknowledgments] The authors would like to thank Emeric Thibaud for providing
the Swiss rainfall data used in Appendix B.

Anna Kiriliouk was supported in part by the Fonds de la Recherche Scientifique
- FNRS (grant number J.0124.20).

## REFERENCES

  * Barioli and Berman (2003) [author] Barioli, FrancescoF. and Berman, AbrahamA. (2003). The maximal cp-rank of rank k completely positive matrices. Linear Algebra and its Applications 363 17–33.
  * Beirlant et al. (2004) [author] Beirlant, J.J., Goegebeur, Y.Y., Segers, J.J. and Teugels, J.J. (2004). Statistics of Extremes: Theory and Applications. Wiley.
  * Boţ and Nguyen (2021) [author] Boţ, Radu IoanR. I. and Nguyen, Dang-KhoaD.-K. (2021). Factorization of completely positive matrices using iterative projected gradient steps. Numerical Linear Algebra with Applications 28 e2391.
  * Canty and Ripley (2020) [author] Canty, AngeloA. and Ripley, B. D.B. D. (2020). boot: Bootstrap R (S-Plus) Functions R package version 1.3-25.
  * Coles, Heffernan and Tawn (1999) [author] Coles, StuartS., Heffernan, JanetJ. and Tawn, JonathanJ. (1999). Dependence Measures for Extreme Value Analyses. Extremes 2 339–365.
  * Cooley and Thibaud (2019) [author] Cooley, DanielD. and Thibaud, EmericE. (2019). Decompositions of dependence for high-dimensional extremes. Biometrika 106 587-604.
  * Cui and Zhang (2018) [author] Cui, QiurongQ. and Zhang, ZhengjunZ. (2018). Max-linear competing factor models. Journal of Business & Economic Statistics 36 62–74.
  * Danielsson et al. (2016) [author] Danielsson, JonJ., Ergun, Lerby MuratL. M., de Haan, LaurensL. and de Vries, Casper GC. G. (2016). Tail index estimation: Quantile driven threshold selection. Available at SSRN 2717478.
  * de Haan and Ferreira (2006) [author] de Haan, LaurendL. and Ferreira, AnaA. (2006). Extreme Value Theory: an Introduction. Springer-Verlag Inc.
  * Drees and de Haan (2015) [author] Drees, HolgerH. and de Haan, LaurensL. (2015). Estimating failure probabilities. Bernoulli 21 957–1001.
  * Einmahl, Kiriliouk and Segers (2018) [author] Einmahl, John HJJ. H., Kiriliouk, AnnaA. and Segers, JohanJ. (2018). A continuous updating weighted least squares estimator of tail dependence in high dimensions. Extremes 21 205–233.
  * Einmahl et al. (2022) [author] Einmahl, John HJJ. H., Ferreira, AnaA., de Haan, LaurensL., Neves, ClaudiaC. and Zhou, ChenC. (2022). Spatial dependence and space–time trend in extreme events. The Annals of Statistics 50 30–52.
  * Fomichov and Ivanovs (2020) [author] Fomichov, VV. and Ivanovs, JJ. (2020). Detection of groups of concomitant extremes using clustering.
  * Fougeres, Mercadier and Nolan (2013) [author] Fougeres, Anne-LaureA.-L., Mercadier, CécileC. and Nolan, John PJ. P. (2013). Dense classes of multivariate extreme value distributions. Journal of Multivariate Analysis 116 109–129.
  * Gissibl and Klüppelberg (2018) [author] Gissibl, NadineN. and Klüppelberg, ClaudiaC. (2018). Max-linear models on directed acyclic graphs. Bernoulli 24 2693–2720.
  * Groetzner and Dür (2020) [author] Groetzner, PatrickP. and Dür, MirjamM. (2020). A factorization method for completely positive matrices. Linear Algebra and its Applications 591 1–24.
  * Hill (1975) [author] Hill, Bruce MB. M. (1975). A simple general approach to inference about the tail of a distribution. The Annals of Statistics 1163–1174.
  * Janßen and Wan (2020) [author] Janßen, AnjaA. and Wan, PhyllisP. (2020). $k$-means clustering of extremes. Electronic Journal of Statistics 14 1211–1233.
  * Kim and Kokoszka (2022) [author] Kim, MihyunM. and Kokoszka, PiotrP. (2022). Extremal dependence measure for functional data. Journal of Multivariate Analysis 189 104887\.
  * Kiriliouk (2020) [author] Kiriliouk, AnnaA. (2020). Hypothesis testing for tail dependence parameters on the boundary of the parameter space. Econometrics and Statistics 16 121–135.
  * Kiriliouk and Naveau (2020) [author] Kiriliouk, AnnaA. and Naveau, PhilippeP. (2020). Climate extreme event attribution using multivariate peaks-over-thresholds modeling and counterfactual theory. The Annals of Applied Statistics 14 1342–1358.
  * Klüppelberg and Krali (2021) [author] Klüppelberg, ClaudiaC. and Krali, MarioM. (2021). Estimating an extreme Bayesian network via scalings. Journal of Multivariate Analysis 181 104672\.
  * Larsson and Resnick (2012) [author] Larsson, MartinM. and Resnick, Sidney IS. I. (2012). Extremal dependence measure and extremogram: the regularly varying case. Extremes 15 231–256.
  * Lenderink, Van Meijgaard and Selten (2009) [author] Lenderink, GG., Van Meijgaard, EE. and Selten, FF. (2009). Intense coastal rainfall in the Netherlands in response to high sea surface temperatures: analysis of the event of August 2006 from the perspective of a changing climate. Climate dynamics 32 19–33.
  * Mainik and Embrechts (2013) [author] Mainik, GeorgG. and Embrechts, PaulP. (2013). Diversification in heavy-tailed portfolios: properties and pitfalls. Annals of Actuarial Science 7 26–45.
  * Resnick (2007) [author] Resnick, Sidney IS. I. (2007). Heavy-tail phenomena: probabilistic and statistical modeling. Springer Science & Business Media.
  * Resnick and Stărică (1997) [author] Resnick, SidneyS. and Stărică, CătălinC. (1997). Smoothing the Hill estimator. Advances in Applied Probability 29 271–293.
  * Samorodnitsky and Taqqu (1994) [author] Samorodnitsky, GennadyG. and Taqqu, Murad SM. S. (1994). Stable Non-Gaussian Random Processes: Stochastic Models with Infinite Variance.
  * Thibaud and Opitz (2015) [author] Thibaud, EmericE. and Opitz, ThomasT. (2015). Efficient inference and simulation for elliptical Pareto processes. Biometrika 102 855–870.
  * Valk (2016) [author] Valk, Cees deC. d. (2016). Approximation and estimation of very small probabilities of multivariate extreme events. Extremes 19 687–717.


# How to Use the graphicx Package

T. J. Mahoney

###### Abstract

This guide explains how to import and manipulate graphics into your article
for inclusion in a volume of proceedings in the ASP Conference Series.

Instituto de Astrofísica de Canarias, E-38205 La Laguna, Tenerife, Spain

##  1\. Introduction

Given the widely varying degrees of familiarity with
L$\scriptstyle\kern-2.10002ptA$TeX of our authors, we feel that the
`\plotone`, `\plottwo`, and `\plotfiddle` commands will serve most purposes.
However, those who are more au fait with L$\scriptstyle\kern-2.10002ptA$TeX
might wish to use the graphicx package. This package enables the user to wield
far greater control over the importation of graphics into a .tex file and,
apart from the usual sizing and rotational facilities, also enables the user
to crop or trim an image as desired (e.g., to get rid of surrounding blank
margins). The graphicx package is useful if you need to use only a part of a
complete image. This guide will take you through the various parameter
(“keys”) of the package and applying them to the sample figure demo.eps
(copied from Kopka & Daly 1999).

##  2\. Invoking graphicx

The first step to take is to ensure that you have the command

`\usepackage{graphicx}`

in the preamble of your article. Importing a graphic is then done by usingthe
command

`\includegraphics[`key1 = …, key2 = …, etc.`]{`filename`}`

Optional parameters—called “keys”—enable the figure to be resized, rotated,
cropped, trimmed, etc. These keys and their functions are listed below.

  * •

`scale =` number — a magnification factor

  * •

`width =` length — the width to which the figure should be scaled111If the
height is not given, it too is scaled by the same amount.

  * •

`height =` length — the height to which the figure should be scaled222If the
width is not given, it too is scaled by the same amount.

  * •

`totalheight =` length — height plus depth of figure (to be used if figure is
rotated)

  * •

`keepaspectratio = true/false` — maintains the height/width ratio

  * •

`angle =` number — angle (in degrees) by which the figure is to be rotated
counterclockwise

  * •

`origin =` location333bl = bottom left, c = center, t = top, r = right, B =
baseline, or any meaningful combination thereof. — the point about which
rotation is to occur

  * •

`draft = true/false` — prevents figure from being imported, but created a
named box with the dimensions of the figure (this option is used to speed up
processing)

  * •

`clip = true/false` — excludes whatever is outside the bounding box

  * •

`bb =` llx lly urx ury — enters the bounding box coordinates, which are given
by default in points (1/72 inch), manually (the bounding box might be missing
or you might want to alter it)

  * •

`viewport =` llx lly urx ury — specifies bounding box w.r.t. bottom left of
existing bounding box; used with clip to select a part of the image (or to
clear unwanted margins.

  * •

`trim =` dllx dlly durx dury — reduces the bounding box by the amount
specified

  * •

`hiresbb = true/false` — reads the bounding box information from the line
`%%HiResBoundingBox` in figure file

##  3\. Experimenting with graphicx

Figure 1 shows the default image produced by the graphics file demo.eps,
called using the command

`\includegraphics{demo}`

(note that none of the keys is used).

Figure 1.: Here is a demo.eps to play with.

Figure 2 shows the same figure scaled by a factor of three and rotated by 60
degrees anticlockwise and using the command

`\includegraphics[scale = 3, angle = 60]{demo}`

Figure 2.: Here is demo.eps again rescaled to three times its original size
and rotated counterclockwise through 60 degrees.

Figure 3 shows the default image with a different aspect ratio, created by
varying the height and width with the command

`\includegraphics[width = 400pt, height = 200pt]{demo}`

Figure 3.: Here is demo.eps with a different aspect ratio.

In Figure 4 we plot the file family demo2.eps, which has a larger bounding box
but using the default keys (`\includegraphics{demo2}`). We see that the
relevant part of the figure has been forced down to the lower left of the
figure in spite of our efforts to center it. There are three ways of
rectifying this situation. We can:

  1. 1.

Force a new bounding box manually with the bb key (keying in the original
coordinates)

  2. 2.

Use the viewport key

  3. 3.

Use the trim key

We can force a new bounding box by using:

`\includegraphics[bb = 169 158 233 242]{demo2}`

We have used the bounding box values for demo.eps for demo.eps2. The result is
shown in Figure 5. We chose these values for the new bounding box because we
already knew them from demo.eps; in practice we would need to experiment with
various values in order to find something that works.

Figure 4.: Here is demo2.eps, which is demo.eps with a different bounding box;
this creates wide margins at the top and to the right of the original figure
and causing the relevant part of the figure into the bottom left corner.
Figure 5.: We set a new bounding box manually with the key bb to eliminate the
unwanted margins and restore the figure to its central position.

We get identical results if we use the command

`\includegraphics[viewport = 0 0 64 84]{demo2.eps}`

(see Figure 6) or

`\includegraphics[trim = 0 0 300 300]{demo2.eps}`

(see Figure 7). (Note that the key takes the bottom left corner of the
existing bounding box as its origin.)

Figure 6.: Elimination of unwanted margins using the viewport key (note that
the origin, (0,0), is now the bottom left corner of the original bounding box.
Figure 7.: Elimination of unwanted margins using trim key.