#1 1

We obtain a functional analogue of the quantile of probability
measures admitting a continuous Lebesgue density on $\RR^d$, and use
it to characterize the class of non-trivial limit distributions of
radially recentered and rescaled multivariate exceedances in
geometric extremes. A new class of multivariate distributions is
identified, termed radially-stable generalized Pareto
distributions, and is shown to admit certain stability properties
that permit extrapolation to extremal sets along any direction in
$\RR^d$. Based on the limit Poisson point process likelihood of the
radially renormalized point process of exceedances, we develop
parsimonious statistical models that exploit theoretical links
between structural star-bodies and are amenable to Bayesian
inference. The star-bodies determine the mean measure of the limit
Poisson process through a hierarchical structure. Our framework
sharpens statistical inference by suitably including additional
information from the angular directions of the geometric exceedances
and facilitates efficient computations in dimensions $d=2$ and
$d=3$. Additionally, it naturally leads to the notion of the return
level-set, which is a canonical quantile set expressed in terms of
its average recurrence interval, and a geometric analogue of the
uni-dimensional return level. We illustrate our methods with a
simulation study showing superior predictive performance of
probabilities of rare events, and with two case studies, one
associated with river flow extremes, and the other with
oceanographic extremes.

Keywords: Bayesian inference, limit set, Poisson process,
quantile set, multivariate regular variation, radial function,
radially-stable distributions, return level-set, multivariate
threshold exceedances, starshaped set

§ INTRODUCTION

The realm of extreme events, by its multivariate nature, casts a
shadow of potentially devastating consequences upon ecosystems,
infrastructures, as well as financial, economic and insurance
sectors. Knowledge of the frequency and magnitude of extreme events
is crucial in enhancing planning strategies and adaptation efforts.The statistical properties of univariate extremes are
well-established <cit.>, but

statistical inference for multivariate random processes is much more intricate: one must analyse how random processes interact with and
influence each other.

A common way to describe the extremal dependence structure of a
real-valued random vector $\bm X=\left(X_1,\dots,X_d\right)$ is
through the coefficient of tail dependence defined by

\begin{equation}
\label{eq:chi-def}
\chi_q(A) = \PR\Big[\,\bigcap_{j\in A}\{F_j(X_{j})>q\}\,\Big]\big/(1-q),\quad q\in(0,1), \quad A\subseteq\left\{1,\dots,d\right\}, \left|A\right|>1,
\end{equation}

where $F_j$ is the cumulative distribution function of $X_j$. When
$\lim_{q\rightarrow1}\chi_q(A)=0$, the variables in $A$ are unlikely
to grow large together and we say that they exhibit asymptotic
independence. Conversely, when $\lim_{q\rightarrow1}\chi_q(A)>0$, the
variables in $A$ are likely to grow large together and the variables
in $A$ exhibit asymptotic dependence. The different dependence
structures that can be present within subgroups of the marginal
variables of $\bm X$

can make inference for rare events deeply challenging and render
extrapolation inacurrate.

Classical approaches to multivariate extreme events often rely on the
framework of multivariate regular variation (MRV) <cit.>,
which posits that the point processes of exceedances of a random
vector over a high threshold, when suitably renormalised, converge in
distribution to a non-degenerate non-homogenous Poisson point
process <cit.>. This provides a framework for understanding
the joint behavior of extreme observations and leads to meaningful
limit distributions that can be used for statistical inference of
multivariate extremes. A major development in this context is the
class of multivariate max-stable <cit.> and multivariate
generalized Pareto <cit.> distributions. These
distributions emerge respectively as the unique non-trivial limit
distributions of

renormalised componentwise maxima, and of renormalised exceedances.However, in practical implementations, MRV is applied in a way that
does not adequately describe relationships among asymptotically
independent random variables <cit.>, which is the
scenario that is of most practical interest. This is due to the type
of the renormalisation that is employed: given a random sample of $n$
$d$-dimentional observations, all components are normalised by the
same amount. This leads to considering the dependence structure only
in a single direction in $\RR^d$ (see Figure <ref> and
Section <ref>). If the extremal dependence
properties between the components of the random vector are strong,
that is, when $\chi_q({\{1,\dots,d\}})>0$ for all sufficiently large
$q$, then MRV provides with a rich description of the extremal
dependence structure of asymptotically dependent random
variables. Otherwise, it is not

Another drawback of statistical methods based on MRV is the limited
set of directions in the multivariate space in which one can
extrapolate the model. In particular, under MRV, the probability of
lying in an extremal set is estimated by shifting the extremal set and
performing empirical probability estimation on the translated set.When the translated set does not contain observations from the initial
dataset, then this leads to an inability to estimate the probability
of interest. To correct the joint rate of tail decay in the case of
weaker extremal dependence, the notion of HRV was introduced in
<cit.>. However, HRV
also suffers from the same extrapolation drawbacks as MRV since it
does not allow extrapolation along directions where not all variables
are simultaneously large. To extrapolate to a wider range of
directions, the frameworks of conditional extremes <cit.>
and angular dependence <cit.> have been introduced to
extrapolate to a wider range of extremal regions with a wider array of
dependence structures, but statistical methodology based on thes
frameworks comes with some drawbacks. Despite its wide applicability
and widespread adoption, the conditional extremal inference method
of <cit.> is based on composite likelihood methods and on
gluing post-fit separate models, making statistical inference and
computations unwiedly. While the angular dependence method
of <cit.> permits extrapolation in regions where variables
are not simultaneously extreme, it is only useful for joint survival
regions. The directions of extrapolation of the above described
frameworks in multivariate extreme value theory (MEVT) are shown in
Figure <ref>.

Directions along which MEVT frameworks allow extrapolation
to tail regions: $(a)$ MRV, $(b)$ and $(c)$ conditional extremes
given $X_2$ and $X_1$ are large, respectively, and $(d)$ geometric
extremes, with $\QS_q$ illustrating the posterior mean of the
quantile set ($q=0.95$). The support of the distribution of
multivariate exceedances is inscribed by the region containing the

Recently, a characterisation of extremal dependence through the
limiting geometry of suitably rescaled sample clouds—or observations
from $\bm X$—has become of interest. The limit set, with boundary
arising as the limiting hull of appropriately scaled sample clouds,
provides insight into the extremal dependence structure of $\bm X$.The so-called gauge function, whose unit level set is in
one-to-one correspondence with the boundary of the limit set, has been
shown to connect several known coefficients describing extremal
dependence of known
copulas <cit.>.<cit.> proposed a parametric approach to modelling the
gauge function associated with the limit set. Using a radial-angular
decomposition, this framework treats the gauge function evaluated at
angles as a rate parameter of a left-truncated gamma model for the
distribution of radii conditioned on angles on the unit simplex.Inference for the gauge function and its associated limit set is based
on parametric models derived from known copulas in $d$-dimensions and
in standard exponential margins. Also in exponential margins, but in
a bivariate setting, <cit.> model the conditional
distribution of the excess radii given angles on the simplex via the
generalized Pareto (GP) distribution <cit.>. In this
framework, the gauge function is seen as a rate parameter of the GP
distribution and is modelled via generalized additive
models <cit.>. <cit.>
also propose a statistical method to estimate the bivariate gauge
function and the shape of its associated limit set using Bézier

0.95 prediction intervals for the
return level-set on the Pang and Windrush river flow data set in original margins (left) for a return period of $T=10^3$ days and in Laplace margins (right) for return periods $T=10^3$, $10^5$, and $T=10^9$ days (from light to dark grey respectively).

In this paper, we develop a framework beyond the cone
$\RR_+^{d}$ by exploiting the structure of weak limits found in
suitably radially renormalised multivariate exceedances.Importantly, in our approach, the multivariate excedances are samples
that lie in the complement of a quantile set, a compact
subset of $\RR^d$ containing the origin. The boundary of the
quantile set can be interpreted as a canonical quantile in $\RR^d$,
allowing to bypass complications experienced in multivariate extreme
value theory due to the lack of natural ordering in $\RR^d$ <cit.>. By specifying an appropriate sequence of
quantile sets, we show under mild conditions that, not only can the
framework of MRV be suitably extended to the cone
$\RR^d\setminus \{0\}$, but that similar extensions are also
available for light-tailed and bounded multivariate distributions.A key property of our framework is that it is applicable to any
$d\geq 1$.

In the multivariate setting ($d>1$), our framework enables the
modelling of the entire joint tail, accommodating scenarios where
subsets of components are extreme. This approach can capture
behaviours across the entire spectrum of multivariate space and reveal
hidden dependencies, thereby bridging the gap between theory and
practice. Using radially recentred and rescaled exceedances over high
quantile sets, we characterise non-trivial limit distributions on
$\RR^d\setminus \{0\}$, termed radially-stable
distributions. The radial stability properties of these families
of distributions permit extrapolation beyond the range of observed
data along any direction (Figure <ref>, right) and lead
to the notion of a canonical return level-set in $\RR^d$. The
return level-set is a geometric extension of the uni-dimensional
(upper-tail) return level <cit.>. Specifically, it is a
quantile of the multivariate distribution expressed in terms of its
average recurrence interval. This gives a clear and interpretable way
to communicate the risk associated with extreme multivariate events.This allows decision-makers, policymakers, and the general public to
understand the likelihood of experiencing impact from joint
extreme events within, for example, a specified scale frame.Figure <ref> shows predictive intervals of
return level-set boundaries associated with different return periods
of daily river flow ($m^3/s$) at the rivers Pang and Windrush, in
England, previously analyzed by <cit.>.

Section <ref> provides background on univariate extreme
value theory and multivariate regular variation, and introduces the
main results, which include the weak convergence of radially
recentered and rescaled exceedances to a Poisson point process on
$\RR^d\setminus\{0\}$. Our findings lead to a novel class of limit
multivariate distributions that are presented in
Section <ref> and to novel return level-sets
presented in Section <ref>. Section <ref>
details methodology for statistical inference of $d=2$ and $d=3$
dimensional extremes using hierarchical Bayesian models with latent
Gaussian random effects on Euclidean spheres.Section <ref> discusses sampling from the joint
posterior distribution of the latent parameters and rare event
probability estimation. Last, in Section <ref> we
illustrate the merits of our approach on 2- and 3-dimensional data
sets of extreme flow in the rivers Pang and Windrush, both tributaries
to the Thames in England, and of extreme sea levels in Newlyn,

§ THEORY

§.§ Univariate extremes

When focusing on the upper right tail () of a random
variable $X \in \mathbb{R}$, prevailing approaches primarily rely on
the fundamental assumption of extended regular variation
<cit.>. This assumption postulates the existence of sequences
of norming scalars, $a_{n, \textsf{R}} \in \mathbb{R}_+$, and
$b_{n, \textsf{R}} \in \mathbb{R}$, $n = 1, 2, \dots$, satisfying

\begin{equation} n \mathbb{P}\left[\{X-b_{n, \textsf{R}}\}/a_{n,
\textsf{R}} > z\right] \vg \Lambda_{\textsf{R}}((-\infty, z]):=
\left(1+\xiR z\right)_+^{-1/\xiR} \quad \text{in }
M_+(\mathbb{R}_+\setminus\{0\}),\quad \text{as $n\to\infty$},
\label{eq:RV_univariate_R}
\end{equation}

where $\Lambda_{\textsf{R}}$ is a cumulative survival function, $\xiR \in \mathbb{R}$, $\left(x\right)_+=\max(x,0)$, and $\vg$ means
vague convergence in the space $M_+(\RR_+\setminus\{0\})$ of
nonnegative Radon measures on $\RR_+\setminus\{0\}$ <cit.> and $B'$ stands for the complement of the set
$B$ in $\RR$. The significance of this convergence becomes
particularly pronounced when examining its interplay with classical
limit theorems governing univariate extreme events. The convergence
outcomes extend to encompass two essential statistical methods that
have been developed to model the univariate extremes. The first
approach relates to the appropriately renormalised maxima, a concept
often harnessed within the statistical methodology based on block
maxima. Within this context, the convergence leads to the well-known
generalized extreme value (GEV) distribution, which stands as a
foundational model for the limiting behaviour of maxima over long time
intervals or blocks. The sequences of norming factors and translation
constants, which feature prominently in the analysis of block maxima,
are effectively absorbed into the limit GEV distribution as part of
the inference process. As a result, the GEV distribution acquires
distinct location, scale and shape parameters $\muR\in\RR$,
$\sigmaR>0$, $\xiR\in \RR$, respectively, and its distribution
function takes the form

\[
\GR(x) = \exp\left[-\left\{1 + \xiR \left(\frac{z -
\muR}{\sigmaR}\right)\right\}_+^{-1/\xi}\right] =
\exp\left\{-\Lambda_{\textsf{R}}(z; \thetaR)\right\},
\]

where $\thetaR=(\muR, \sigmaR, \xiR)$. The GEV distribution has
max-stability, a cornerstone trait that underpins its utility in
modelling extreme events across diverse contexts.

The second avenue emanating from convergence
(<ref>) relates to threshold exceedances. This
statistical strategy aligns with the peaks over threshold approach,
which focuses on modelling the distribution of extreme values that
exceed a predefined high threshold level. In this case, the
convergence leads to the generalized Pareto (GP) distribution function
given by

\[
\HR(y)= 1-\left\{1 + \xiR
\left(\frac{z}{\sigmaR}\right)\right\}_+^{-1/\xi},
\]

which offers a robust framework for analyzing extreme events that
exceed a prespecified threshold. By facilitating the modelling of
these exceedances, the GP provides a statistical method that is more
efficient than the block-maxima approach for estimating tail
probabilities and assessing risks associated with extreme events.

These limit distributions provide insights into the behaviour of
extreme values in the upper right tail. When interest is in the lower
left tail () of the distribution of $X$, the convergence
condition is rewritten with $X$ replaced by $-X$, and with a different
sequence of norming constants so as to allow flexibility on the shape
parameter. Thus, in this case, the convergence condition for the
lower left tail is typically given by

\begin{equation} n \mathbb{P}\left[\{-X-b_{n,
\textsf{L}}\}/a_{n,\textsf{L}} > z\right] \vg
\Lambda_{\textsf{L}}[(-\infty, z]'; (0, 1, \xiR)]:= \left(1+\xiL
z\right)_+^{-1/\xiL} \quad \text{in } M \left(\mathbb{R}_+
\setminus\{0\}\right),
\label{eq:RV_univariate_L}
\end{equation}

where $a_{n,\textsf{L}} >0$, $b_{n,\textsf{L}}\in \RR$ for all
$n\in\NN$ and $\xiL\in\RR$. When interest is in both left and right tails, the convergence
assumptions (<ref>) and (<ref>)
can be unified under a single theme. Perhaps not surprisingly, if
both convergences hold, then we can construct sequences of norming
functions $r_n^a\,:\,\mathbb{S}^0 \to \RR_+$ and
$r_n^b\,:\,\mathbb{S}^0\to \RR$, $n=1,2,\dots$, that are functions of
the direction $X/\lvert X\rvert\in\mathbb{S}^0$,
$\mathbb{S}^0=\{-1, 1\}$, along which extremes occur, such that for
$\zL < 0$ and $\zR>0$,

\begin{equation}
n \PR\left[\frac{X -
r_{n}^a(X/|X|)}{ r_n^b(X/|X|)} \in[\zL, \zR]' \right] \vg
\Lambda([\zL,\zR]')\qquad \text{in $M(\RR \setminus \{0\})$,\qquad as
$n\to \infty$,}
\label{eq:RV_univariate_RL}
\end{equation}

\begin{equation} \Lambda([\zL, \zR]) = (1-\pi) (1-\xi_{\textsf{L}}
z)_+^{-1/\xi_{\textsf{L}}} + \pi (1+\xi_{\textsf{R}}
z)_+^{-1/\xi_{\textsf{R}}} ,\quad \pi = \PR(X/|X| = 1), \quad
\xi_{\textsf{L}} , \xi_{\textsf{R}} \in \RR.
\label{eq:Lambda_univariate}
\end{equation}

Convergence (<ref>) implies domain of attraction
properties for the distribution of univariate extremes. In
particular, as $n\to\infty$, the joint distribution of the
appropriately renormalised minima $m_n=\min_{i=1,\dots, n}(X_i)$ and
maxima $M_n=\max_{i=1,\dots, n}(X_i)$ of a random sample
$(X_1,\dots, X_n)$ from $X$, converges weakly to a non-degenerate
(m_n - b_n,/a_n, ≥, M_n - b_n,/a_n,
≤)    =    (X_i -
r_n^b(X/|X|)/r_n^a(X/|X|) ∈[,
])^n G([z_,]),

$G([z_\textsf{L},\zR]) = G([\zL,\zR]; \theta_{\textsf{L}},
\theta_{\textsf{R}})=\exp[-\Lambda([\zL, \zR]')]$ and $\Lambda$
defined in expression (<ref>). Similarly, the
distribution of the radially rescaled excesses leads to the radial
peaks-over-threshold model

\[
\PR\left[\frac{X-r_{t}^a(X/|X|)}{r_{t}^b(X/|X|)}\in[\zL,
\zR]' \mid |X|>r_{t}^a(X/|X|)\right] \wk \Lambda([\zL,
\zR]')\qquad \text{as $t\to\infty$}.
\]

Absorbing the sequences of norming functions as location and scale
parameters into the limit distribution leads to a class of
distributions satisfying a new notion of $\min-\max$ stability, that

\[
G([\zL, \zR]; \thetaL, \thetaR)^t =
G([\zL, \zR]; \thetaL^\star, \thetaR^\star),
\qquad \text{for all $t > 0$},
% G(z;
% , ([n^{\xi_{\textsf{R}}} - 1]/\xi_{\textsf{R}},
% n^{\xi_{\textsf{R}}}, \xi_{\textsf{R}}))^n
\]

$\thetaL^\star= ([n^{\xi_{\textsf{L}}} - 1]/\xi_{\textsf{L}},
n^{\xi_{\textsf{L}}}, \xi_{\textsf{L}})$ and $\thetaR^\star$ defined
similarly. This intriguing stability property, evident from the form
of $\Lambda$ in equation (<ref>), might be aptly
termed radial stability. While the term $\min-\max$
stable distribution could be tempting for $G$, it might not
comprehensively convey the unique characteristics of this
distribution, especially when extended to the multivariate setting,
which is the focus of this paper. In the multivariate context,
extremes cannot be solely described by component-wise maxima or
minima, but rather involve intricate dependencies among components.

This naming convention underscores the distinct behaviour of the limit
distribution in capturing extreme events across both left and right
tails while respecting the inherent relationships among components.Remarkably, the unifying convergence assumption
(<ref>), overlooked in the literature until now,
presents a substantial departure from traditional approaches to
extreme value modelling.

§.§ Multivariate regular variation

To understand radial stability properties in a $d$-dimensional setup,
we revisit classical multivariate extreme value theory through
consideration of multivariate regular variation (MRV), as it is still
a key extremal dependence assumption invoked in a wide range of
statistical methods including multivariate max-stable and multivariate
generalized Pareto distributions. Although often expressed in
standard margins, we here consider without loss the marginals to be
unit-exponentially distributed. With this choice for the margins, a
random vector $\bm X=\bm X_E$ is multivariate regularly varying on the
cone $\mathcal{E}=[-\infty,\infty]^d \setminus\{-\infty\}$ if there
exists a limit measure $\nu_{\textsf{MRV}}$ such that for any
relatively compact set $B\subset \mathcal{E}$,

\begin{equation}
r\PR[\bm X - (\log r)\bm 1 \in B] \to \nu_{\textsf{MRV}}(B), \quad \text{as $r\to \infty$},
\label{eq:MRV}
\end{equation}

with $\nu_{\textsf{MRV}} (\partial B) = 0$. A consequence of
convergence (<ref>) is that
$\nu_{\textsf{MRV}}(t \bm 1 + B)=\exp(-t)\nu_{\textsf{MRV}}(B)$. This
stability property implies domain of attraction properties for the
distribution of $\bm X$.  For example, the distribution
function of the componentwise maxima converges weakly after
appropriate recentering to

\[
\PR\left(\bm M_n - (\log n)\bm 1\leq \bm z\right)=\PR\left\{\bm X_i
- (\log n)\bm 1 \in [-\bm \infty, \bm z]\,:\,i=1,\dots,n\right\} \wk
\exp\left[-\Lambda_{\textsf{MRV}}([-\bm \infty, \bm z])\right],
\]

where $\Lambda_{\textsf{MRV}}(B)=\nu_{\textsf{MRV}}(B')$. The measure
$\Lambda_{\textsf{MRV}}$ is in one-to-one correspondence with a Radon
measure $H_{\textsf{MRV}}$ on the the unit-simplex
$\Delta_+^{d-1}=\{\bm \omega \in \RR^{d}_+ \,:\, \lVert \bm w \rVert_1=1\}$,
where $\lVert \bm x \rVert_p = (\sum_{i=1}^d \lvert x_i\rvert^p)^{1/p}$ for $\bm x = (x_1,\ldots,x_d)\in\RR^{d}$, termed the spectral
measure, satisfying the mass–moment constraint $H(\Delta_+^{d-1})=d$
and $\int_{\Delta_+^{d-1}} w_i H(\bm w) = 1$, for $i=1,\dots,d$. The
correspondence is given by

\[
\Lambda_{\textsf{MRV}}(B) = \int_{\Delta_+^{d-1}} \max\left[\bm w \exp(-\bm
x)\right] H_{\textsf{MRV}}(d\bm w),
\]

where vector algebra is interpreted component-wise. Similarly, the
vector of renormalised threshold exceedances converges in distribution
to a multivariate Pareto distribution, that is,

\[
\lim_{r\to \infty}\PR\left[\bm X -r\bm 1\in [-\bm \infty,\bm x] \mid \max \bm
X > r\right] = \frac{\Lambda_{\textsf{MRV}}(\bm x \wedge \bm
0)-\Lambda_{\textsf{MRV}}(\bm x)}{\Lambda_{\textsf{MRV}}(\bm 0)},
\]

at continuity points $\bm x \in\mathcal{E}$ of the limit distribution
function of the multivariate generalized Pareto distribution, where
$\Lambda_{\textsf{MRV}}(\bm x)=\nu_{\textsf{MRV}}([-\infty, \bm
x]')$. Multivariate generalized Pareto distributions arise as the
only non-trivial limit distributions when renormalised multivariate
exceedances are considered <cit.>. These distributions
possess several appealing theoretical and practical properties
<cit.>. Crucially, they are often thought to be the
multivariate analogue of the generalized Pareto distribution
<cit.> by being the only threshold-stable multivariate
distributions <cit.>. The practical importance of
their threshold-stability lies in the preservation of the model's
distributional type at higher levels, which is advantageous for
extrapolation into the joint tail region of the

Both component-wise maxima and multivariate threshold exceedance
approaches can be unified under a single theme. Let $r_q$ be the
$q$-th quantile of the distribution of $\lVert \bm X \rVert_1$. Then
the overarching assumption that unifies these approaches is that as
$q\to 1$, the point process of recentered exceedances
$\bm X_i- r_q\bm 1 \mid \lVert \bm X_i\rVert_1 > r_q$, converges in
distribution to Poisson point process on
$(0,\infty)\times \Delta_+^{d-1}$ with intensity measure

lim_q→1[‖X‖_1 - r_q>z, X/‖X‖_1 ∈B|‖X‖_1 >
r_q] →exp(-z) (B),  $z>0$,
$B\subset \Delta_+^{d-1}$.

In practice, the convergence (<ref>) is typically assumed
to hold with $\HMRV$ placing mass in the interior
$(\Delta_+^{d-1})^\circ:=\{\bm \omega\in\RR_+^d\,:\, \min_i\omega_i >
0, \lVert\bm \omega\rVert_1=1\}$ of the unit-simplex $\Delta_+^{d-1}$.Although useful for theoretical investigations, $\HMRV$ placing mass
in $(\Delta_+^{d-1})^\circ$ as a dependence assumption is
inadequate in practice since it imposes a strong form of dependence in
the joint tail, see for example Section <ref>. The strength
of this dependence, known as asymptotic dependence (see Section
<ref>), is necessary so that the mass of the distribution of
$\bm X$ does not rapidly diffuse along the sequence of joint tail
regions described by events of the form $\lVert\bm X\rVert_1> r_q$, as
$q$ increases. Although elegant as an assumption, this limitation has
been appreciated by a broad range of theoretical and practical

0 Proposition 2 of states that scaled random samples
converge in probability onto a limit set
$\G=\left\{\bm{x}\in\mathbb{R}_+^d\,:\,\gG(\bm{x})\leq1\right\}$ with
corresponding gauge function $\gG$ when the random vector's marginal
survivor functions have the form $1-F_j(x)\sim e^{-\psi_j(x)}$ for
$\psi_j(x)\sim\psi(x)$ as $x\rightarrow\infty$. Here, $\psi$ must
be $C^2$ and have a positive derivative such that
$\left({1}/{\psi^\prime (x)}\right)^\prime\rightarrow0$ for
$x\rightarrow\infty$. For our setting, recall the form of the
standard Laplace survival function

As we require $x\rightarrow\infty$, it is clear that we satisfy the
necessary conditions for the marginal distribution. The choice of
scaling factor $r_n=\log\left({n}/{2}\right)$ is also motivated by
the marginal survival function. Take $x>0$, then writing
$1-F(x)=e^{-x-\log(2)}$ gives $\psi(x)\sim x+\log(2)$. If we take
$\psi(r_n)\sim\log n$ then it follows that
$r_n\sim\log\left({n}/{2}\right)$ for $n\rightarrow\infty$. This
allows us to simply choose $r_n=\log\left({n}/{2}\right)$ in order
to satisfy all necessary conditions.

§.§ Geometric Poisson point process of exceedances

In this section, we consider a novel point process of multivariate
exceedances over suitably chosen starshaped sets. In what follows, we
say that $\A \in \bigstar \subset \mathcal{K}$ and call $\A$ a
star-body whenever $\A$ is starshaped at $0$ and
$0 \in (\text{ker}\,\A)^\circ$, see
Appendix <ref> for background information on
starshaped sets. Because $\A$ is compact and
$0 \in (\text{ker}\,G)^\circ$, then the radial function
$r _{\A}(\bm x)=\sup\{\lambda > 0 \,:\, \lambda \bm x \in \A\}$ of
$\A$ is a Lipschitz continuous $(-1)$-homogeneous function on
$\RR^{d}\setminus \{0\}\cong \SSS^{d-1}$, that is
$r_{\A}(\bm x)=\lVert \bm x\rVert^{-1}r_{\A}(\bm x/\lVert \bm
x\rVert)$, $\bm x\in\RR^{d}\setminus \{0\}$. The reciprocal of the
radial function $g_\A:=1/r_\A$ is known as the gauge function of
$\A$. The simplest star-body, perhaps, is the ball of radius $\ell$
centred about $0$, denoted $B_{\ell}(0)$: it radial function given by
$r_{B_{\ell}(0)}(\bm x) = \ell/\lVert \bm x\rVert$. The algebraic
structure on $\bigstar$ (see Appendix <ref>) enables
operations between star-bodies over the field of scalars in a coherent
and systematic manner, much like in a familiar algebraic structure
over a field. The space of compact star-bodies $\bigstar$ is endowed
with the operations $\radd$ and $\rmult$ of radial addition and
multiplication respectively. The radial addition operation combines
two compact star-bodies $G_1$ and $G_2$ by adding their radial
functions. It is commutative, associative, and has an identity
element (the empty set), which are properties akin to addition in a
ring. Similarly, the radial multiplication combines two compact star
bodies by multiplying their radial functions. The radial
multiplication $\rmult$ is distributive over $\radd$ and also
associative. This structure allows us to analyse and manipulate
star-bodies, revealing intuitive intrinsic properties of objects in
the geometric multivariate extremes framework.Figure <ref> illustrates operations commonly used
on star

Left: $\loc \radd \G$. Centre: $\loc \rdiv \G$. Right:
$\loc \radd B_{-\log(1-q)}(0)\rmult\G$, for $B_r(\bm x)$ the ball of radius $r$ centred at $\bm x$, and $q\in\{0.9,0.99,0.999\}$.  $\loc$ in lightest grey, $\G$ in mid grey, result of operations in
dark grey.

In the geometric setting, multivariate exceedances are no longer
defined with respect to $\lVert \bm X \rVert_p$ exceeding a high
quantile relative to the distribution of $\smash{F_{{\lVert \bm X\rVert}_p}}$.Instead, multivariate exceedances are understood as observations
$\bm X_i$ hitting the complement of a quantile set
$\QSq\in\bigstar$ in $\RR^d$. Henceforth, we define the polar transformation given by $R=\lVert \bm X \rVert$ and $\bm W = \bm X / R$, where $\lVert \bm X \rVert=\lVert \bm X \rVert_2$.

$\rquant(\bm w)=\inf\{r \in \RR\,:\, F_{R\mid \bm W}(r\mid \bm
w)=q\}$ be the conditional quantile function of the distribution of $R\mid \bm W=\bm w$.

The $q$-th quantile set of the distribution of
$\bm X=R\bm W$ is

\[
\QSq := \bigcup_{\bm w \in \SSS^{d-1}} [0\,:\, \rquant(\bm w)\,
\bm w],\qquad q \in (0,1),
\]

where $[0\,:\rquant(\bm w)\,\bm w)]$ denotes the closed line segment
from the origin to $\rquant(\bm w)\, \bm w$.

The conditional quantile function $\rquant$ is in one-to-one
correspondence with $\QS_q\in\bigstar$ in the sense that $\rquant$ is
the radial function of $\QSq$. From total probability, the quantile function satisfies the property

\begin{equation}
\PR[\bm X\in\QSq] = q.
\label{eq:quantile_set_property}
\end{equation}

Remarkably, the quantile function can be interpreted as
a quantile both conditionally and unconditionally, that is, for any
fixed $q\in(0,1)$, the quantile function $\QS_q$ is one of the class
of values of a variate that divides the total mass of a probability
measure on $\RR^d\setminus \{0\}$ into two parts, one that is visited
by any single observation $\bm X_i\in\RR^d$ with probability $q$ and
the other with probability $(1-q)$.  Thus, under random
sampling, the set $\QS_q'$ can be interpreted as a set that is
expected to be visited on average once in $1/(1-q)$ observations.This observation leads to the notion of return set $\QS_{q}'$, which
is a geometric analogue of the return level in univariate extreme
value analysis. Details on return level sets are given in
Section <ref>.
We now motivate our first main result via an assumption that

There exists $\xi \in \RR$ and a family of star-bodies
$\{\G_{q}\,:\, \G_q\in\bigstar, q\in(0, 1)\}$ such that

\[
\lim_{q\to 1} \frac{\partial}{\partial
z}\PR\left[\frac{R-\rquant(\bm W)}{r_{\G_q}(\bm W)}\leq z
~\Big|~\{R> \rquant(\bm W)\}, \bm W\right] = (1+\xi
z)_+^{-1-1/\xi} , \qquad \mathbb{P}-a.s.,
\]

where $(x)_+:=\max(x,0)$.

Let $\{\bm X_{i}\,:\,i=1,\dots, n\}$ be a random sample from $\bm X$
satisfying Assumption <ref>,
$q_{n,k}= 1-k/n$ where $1\leq k < n$, and consider the sequence of
point processes
P_n = { (‖X_i‖- (X_i/‖X_i‖;q_n,k)/r__q(X_i/‖X_i‖), X_i/‖X_i‖) :
i=1,…, n},     n=1,2….

Then, as $n$ approaches $\infty$,

\begin{equation}
P_n \wk P \quad \text{\normalfont in $M_p((0,\infty]\times \SSS^{d-1})$},
\label{eq:PntoP}
\end{equation}

where $P$ is a Poisson point process on
$(0,\infty]\times \SSS^{d-1}$ with intensity measure
$k\,(1+\xi r)^{-1-1/\xi}\,d r F_{\bm W}(d \bm w)$, and
$M_p((0,\infty]\times \SSS^{d-1})$ denotes the collection of all
Radon point measures on $(0,\infty]\times \SSS^{d-1}$.

Instead of relying on assumption (<ref>), we concentrate on the
convergence assumption

\begin{equation}
n \PR\left[\frac{\bm X-r_{\QS_{1-(1/n)}}(\bm X/\lVert\bm X\rVert)(\bm X/\lVert\bm X\rVert)}
{r_{\G_q}(\bm X/\lVert\bm X\rVert)}\in B\right] \vg \nu_\G(B),\qquad \text{in $M_p(\Rstar)$},
\label{eq:MERV}
\end{equation}

where $r_{\G_q},\rquant:\,\Rstar\to (0,\infty)$ are continuous radial
functions and $\nu_\G$ is a Radon measure on $\RR^d\setminus\{0\}$.Convergence (<ref>) implies domain of attraction
properties for the random vector $\bm X$. For example, the
analogue of convergence (<ref>) in the geometric setting
is that for any compact neighbourhood $B$ of the origin $o$, the
appropriately renormalised sample cloud converges weakly to a
non-degenerate limit distribution

\begin{equation}
\PR\left[\frac{\bm X_i-r_{\QS_{1-(1/n)}}(\bm X_i/\lVert\bm
X_i\rVert)(\bm X_i/\lVert\bm X_i\rVert)}{r_{\G_{1-1/n}}(\bm X_i/\lVert\bm
X_i\rVert)}\in B\,:\,i=1,\dots,n\right] \wk
\exp(-\Lambda_G(B)),
\label{eq:conv_max_stab}
\end{equation}

where $\Lambda_G$ is a survival measure and is defined by $\Lambda_G(B)=\nu_\G(B')$. The probability measure $\Lambda_G$
admits the spectral decomposition given by

\begin{equation}
\Lambda_G(B)=\int\limits_{\SSS^{d-1}}
\int\limits_{]0\,:\,\bm w) \cap B'} [1+\xi r]_+^{-1-1/\xi} \,\dd r \Hgauge(\dd
\bm w), \qquad B\in \Rstar,
\label{eq:spectral}
\end{equation}

where $\Hgauge$ is a probability measure on $\SSS^{d-1}$ termed the geometric spectral measure, and
$]0\,:\,\bm w)$ denotes the open half-line in the direction of
$\bm w$. The dependence of $\W$ on $\G$ is explained in Section <ref>. The limit probability in convergence (<ref>) is the void probability of a Poisson random variable with rate $\Lambda_G$. The void probability is interpreted as the probability that no point from a recentred and rescaled sample cloud fall outside a set $B$ as $n\to \infty$. This explains why the integral of the intensity measure is over the complement of $B$. The outer integral is with respect to
$\dd \bm w = \mu(\dd \bm w)$, where $\mu$ is the spherical Lebesgue
measure and the inner integral is with respect to the 1–dimensional
Lebesgue measure on $]0\,:\,\bm w)\cap A$.Representation (<ref>) suggests that $\Hgauge$ has similar
interpretation to the spectral measure $\HMRV$ <cit.>.The form of the spectral measure $\Hgauge$, however, enjoys remarkably
simpler structure than its $\HMRV$ counterpart as
$\Hgauge(\cdot)=F_{(\bm X/\lVert\bm X\rVert) \mid \G}(\cdot\mid \G)$.In particular, under the assumption that $\bm X$ has a Lebesgue

for arbitrary $S_B\subseteq \SSS^{d-1}$, directly implying that the
distribution of $\bm W\mid R>\rquant(\bm W)$ is equal to that of
$\bm W$. Further discussions about the implications on statistical
inference are found in Section <ref>. On
the other hand, when $\bm X$ has the copula of independence, $\HMRV$ places all of its mass on the corners of $\Delta_+^{d-1}$. In
this latter case, the measure $\HMRV$ bears no relation with
$F_{\bm X/\lVert\bm X\rVert }$.

Let $\mu_n$ be a sequence of measures on $\RR^d\setminus\{0\}$ and
$\mu$ a non-degenerate probability measure on
$\RR^d\setminus \{0\}$. Let $\A_n, \B_n\in\bigstar$ such that

\[
\mu_n(\A_n \rmult K \radd \B_n) \wk \mu(K) \qquad \text{for all
$K\in \bigstar$}.\
\]

Then for some non-degenerate probability measure $\mu^\star$ on
$\RR^d\setminus\{0\}$ and $A_n, B_n \in \bigstar$,

\[
\mu_n(A_n \rmult K \radd B_n) \wk \mu^\star(K) \qquad \text{for all
$K\in \bigstar$},
\]

if and only if $A_n\rdiv \A_n \to \mathcal{C}$ and
$(B_n \rsub \B_n)\rdiv\A_n \to \mathcal{D}$ for some
$\mathcal{C}, \mathcal{D}\in \bigstar$, and then
$\mu^\star(K) = \mu(\mathcal{C}\rmult K \radd \mathcal{D})$.

In Corollary <ref> we interpret multiplication
and addition of star-bodies radially. Conditioning on the event that
$\bm X$ exceeds the (continuous) boundary
$r_{\QSq}(\bm X/\lVert \bm X\rVert)$ of the quantile set $\QSq$, is
remarkably informative for the joint tail region $\QSq'$. An
equivalent formulation of convergence (<ref>) is through the
polar transform which gives

lim_q→1[‖X‖- (X/‖X‖)/r__q(X/‖X‖) >z,
X/‖X‖ ∈B |‖X‖>
(X/‖X‖)] = [1+ξz]_+^-1/ξ

where $z>0$, which imples that for any relatively compact set
$A \subset \Rstar$, convergence (<ref>) implies that the
appropriately renormalised multivariate threshold exceedances

[Z ∈A]    :=   lim_q→1 [X -
(X/‖X‖) (X/‖X‖)/r__q(X/‖X‖)  ∈A |
‖X‖> (X/‖X‖) ] =

The distribution of the limiting random vector
$\bm Z$ is defined through the geometric spectral measure $\Hgauge$
and due to the stability property of the generalized Pareto
distribution, it enjoys certain threshold stability properties. These
are described in Section <ref>.
Next, we give simple sufficient conditions subject to which
Assumption <ref> holds for the cases
of bounded ($\xi < 0$), light-tailed ($\xi = 0$) and heavy-tailed
($\xi > 0$) distributions. We remark that for the light-tailed case,
our assumption guarantees that a suitable sequence of scaling
constants $r_n$ can be found so that the random point-set
$\smash{N_n=\left\{\bm X^1/r_n, \dots, \bm

converges in probability onto a limit set $\G$
<cit.> whereas for the heavy tailed-case, our
assumption is in a similar spirit to <cit.> and
guarantees that the distribution of $\bm X$ is multivariate regularly
varying on $\RR^d\setminus\{0\}$. The corresponding assumption for
the case $\xi < 0$ is equivalent to the case $\xi > 0$ after a
suitable transformation, which is illustrated in
Figure <ref>.

Suppose that the random vector $\bm X$ has a Lebesgue density
$f_{\bm X}$ on $\RR^d$ and let $\G\in\bigstar$ be a star-body
described by continuous gauge function $\gG\,:\,\Rstar\to \RR_+$.  Suppose that one of the following conditions

$(i)$ The limit $\QS_{1}:=\lim_{q\to 1} \QS_q$ exists in
$\bigstar$, $f_{\bm X}(\bm x) > 0$ when $\bm x\in \G^\circ$,
$f_{\bm X}(\bm x) = 0$ when $\bm x\in \G'$, and there exists
$\psi_{\textsf{B}}\,:\,\RR_+\to\RR_+$, and a $\xi < 0$ such that
as $t\to \infty$,

\begin{equation}
\frac{f_{\bm X}[\{1-(t\lVert \bm y\rVert)^{-1}\}r_{\QS_1}(t \bm y)(t \bm y)]}{\psi_{\textsf{B}}(t)}\to
g_{\G}(\bm y)^{-1/\xi-1}, \qquad \bm y\in\Rstar.
\label{eq:RVneg}
\end{equation}

$(ii)$ There exists $\psi_{\textsf{L}}\,:\,\RR_+\to \RR_+$,
and a $\rho > 0$ such that as $t\to\infty$,

\begin{equation}
-\frac{\log f_{\bm X}(t \bm y)}{\psi_{\textsf{L}}(t)} \to \gG(\bm y)^\rho, \qquad \bm y\in\Rstar.
\label{eq:RV1}
\end{equation}

$(iii)$ There exists $\psi_{\textsf{H}}\,:\,\RR_+\to \RR_+$,
and $\xi > 0$, such that as $t\to\infty$,

\begin{equation}
\frac{f_{\bm X}(t \bm y)}{\psi_{\textsf{H}}(t)}\to
\gG(\bm y)^{-(d+\xi^{-1})}, \qquad \bm y\in\Rstar.
\label{eq:RVpos}
\end{equation}

If, in addition, the convergence is uniform on $\SSS^{d-1}$, then
Assumption <ref> holds.

$r_{\G_q}$ can be taken as

$r_{\G_q}(\bm w)=\xi \, [r_{\QS_1}(\bm w) - r_{\QS_q}(\bm w)]$
when $\xi < 0$;

$r_{\G_q}(\bm w)=[\psi_{\textsf{L}}'\{\rquant(\bm w)\}\rG(\bm
w)^\rho]^{-1}$ when $\xi = 0$;

$(iii)$ $r_{\G_q}(\bm w)=\rG(\bm w) (1-q)^{-\xi}/\xi$ when
$\xi > 0$.

The case $\xi < 0$ in Proposition <ref> we need to have a
condition for the rescaled vector $\bm Y = \bm X/t$ as $t\to \infty$.The density of the rescaled random variable is

Discuss interpretation for sequence of scaling sets $\G_q$. When
$\xi < 0$, case $\xi < 0$ in Proposition <ref> we need to
have a condition for the rescaled vector $\bm Y = \bm X/t$ as
$t\to \infty$. The density of the rescaled random variable is

$t^d f_{\bm X}(t \bm x)$, where $t^d$ is the Jacobian of the
transformation. This explains the condition for the term and the
exponent $\xi$. For the case $\xi < 0$, we are not interested in
$\bm X/t$ as $t \to \infty$, but instead we may care more for a
different transformation. Let
$\bm Y = [\rG (t \bm X)/\{\rG (t \bm X)-1\}] (t \bm X)/(\lVert t \bm
X\rVert)$. The transformation admits the inverse
$\bm X=(1-\lVert t \bm Y\rVert^{-1})\,\{\rG(\bm Y)\bm Y\}$. Here, the
Jacobian of the transformation

Illustration of the transformation
$\bm x=\rquant(t \bm y)(t \bm y - \bm y^\star)$ showing a sequence
$\bm x_1,\dots,\bm x_5$ in $\RR^d$ approaching the boundary point
$\rquant(\bm y) \bm y$ (arrows, $t=1,2,5,10, \infty$), from
(left) a region of positive density
($\bm y^\top \bm y^\star > 0$) and (right) a region of
zero density ($\bm y^\top \bm y^\star \leq 0$).

§.§ Probability density function of angles

The geometric framework has obvious advantage as it yields a limiting
distribution of angles that precisely mirrors the distribution of
angles found in the original variables, that is
$\Hgauge = F_{\bm X/ \lVert \bm X\rVert}$. For many real-world
applications, the assumption of an absolutely continuous distribution
for $F_{\bm X/ \lVert \bm X\rVert}$ is realistic, rendering our
framework highly suitable. Here, we describe further a link between
$\Hgauge$ and $\G$ which allows improving learning of the set $\G$
through pooling information from angles. This link explains the
substantial gains in uncertainty in our numerical simulations, even in
cases where $\Hgauge$ depends on $\G$ only parti    ally and not totally.Furthermore, we also see that in idealistic settings such as random
sampling, substantial gains can be made by also including in inference
the angles at which non-exceedances occur, that is, all observed
angles, see Section <ref>. This is a major benefit
when compared to other methods.

First, we describe a key connection between the distribution of angles
$\bm{W}=\bm X/ \lVert\bm X\rVert$ with star-bodies in
$\RR^d$. For simplicity, we will assume throughout that $F_{\bm
W}$ is absolutely continuous with respect to the spherical Lebesgue
measure. In this setting, we observe that any continuous probability
density function on the sphere is in one-to-one correspondence with a
strongly starshaped set
$\mathcal{W}$ defined by its radial function $r_{\Astar}(\bm
w):=f_{\bm W}(\bm w)$ as

\[
\Astar := \bigcup_{\bm w \in \SSS^{d-1}} [0\,:\, f_{\bm W}(\bm
w)\, \bm w] \in
\bigstar.% .\ \quad \color{red} \text{or} \quad \Astar^{1/d} :=
% \bigcup_{\bm w \in \SSS^{d-1}} [0\,:\, f_{\bm W}(\bm w)^{1/d}\,
% \bm w] \in \bigstar, \color{black}
% \qquad \text{and} \qquad \Astar_r :=
% \bigcup_{\bm w \in \SSS^{d-1}} [0\,:\, f_{\bm X}(r \bm w) \, \bm
% w] \in \bigstar, \quad \text{$\forall r>0$}.
\]

Because $\mathcal{W}\in\bigstar$ we see that
$r_{\Astar^{1/d}}(\bm w)^d = r_{\Astar}(\bm w) = f_{\bm W}(\bm w)$. Thus,
from the definition of the volume (<ref>) given in Appendix
<ref>, we can verify that for $f_{\bm W}$ to be a valid density, $\Astar$ must satisfy

\begin{equation}
1=\int_{\SSS^{d-1}}f_{\bm W}(\bm w)\,\dd \bm w = d\left[\frac{1}{d}\int_{\SSS^{d-1}}r_{\W^{1/d}}(\bm w)^d\,\dd \bm w\right] = d\vol{\mathcal{W}^{1/d}}.
\end{equation}

To construct any set $\mathcal{W}$ that satisfies the
geometric moment constraint, it suffices to consider
^1/d   =   (d )^1/d,

for some $\LL \in \bigstar$, a fact that follows again from the usual formula for the volume
(<ref>) in polar coordinates. The geometric moment constraint sets the stage for the subsequent
models we propose, each offering a unique perspective on this

One possibility for the form of the set $\LL$ in expression
(<ref>) is to have $\LL=\B$ where $\B\in\bigstar$ with
$\B$ independent of $\G$. In other words, this initial model holds
true for the class of limit distributions $\mathcal{W}$ that are
independent of the limit set $\G$. It is useful to consider a
concrete example which consists of the class homothetic densities on
$\RR^d$ with respect to the gauge function $g_{\G}$ of a starshaped
set $\G\in\bigstar$ <cit.>. Here, a density is called homothetic if
its has level-sets that are scaled copies of a set $\G$.
Proposition <ref> asserts that the density of
angles of homothetic densities is in one-to-one correspondence with
the gauge function $g_\G$ of $\G$.

Suppose that $f_{\bm X}(\bm x)=f_0\left(r_\G(\bm x)^{-1}\right)$ for a decreasing,
positive, a continuous function $f_0:[0,\infty)\to[0,\infty)$ and a
$({-}1)$–homogeneous function $r_{\G}$, describing a shape
$\G\in\bigstar$. Then $\mathcal{L}=\G$.

A proof of Proposition <ref>
based on the integration by parts formula of <cit.> is given in Appendix <ref>.Here we remark that Proposition <ref> provides
with a new probabilistic proof for the volume of a star-body in polar
coordinates <cit.>.  The usefulness
of Proposition <ref> can be seen from the ease
with which one can identify the density of angles in analytic form
when the density is homothetic and the joint density is known. For
concreteness, consider $\G\in\bigstar$ described by $r_{\G}(\bm x) =
(\bm x^\top \mathsf{Q} \bm
x^\top)^{-1/2}$ where $\mathsf{Q}=\mathsf{D}^{1/2}\,
\mathsf{\Sigma}^{-1} \,
\mathsf{D}^{1/2}$, with
$\mathsf{D}=\text{diag}(\mathsf{\Sigma})$ and
$\mathsf{\Sigma}$ a positive definite matrix. Homothetic densities
sharing the same homogeneous function
$g_{\G}$ but having different density generators
$f_0$ are given in Examples 1–3. As noted by <cit.> the
rate of decay of $f_0(r)$ as $r\to
\infty$ determines the properties of the margins.

$\bm X\sim {\normalfont\text{MVN}}(\bm 0,\mathsf{Q}^{-1})$. The
margins of $\bm X$ are standard normally distributed. The density
$f_{\bm X}$ is homothetic with respect to $r_{\G}$ and
$f_0(r) = C_1 \exp(-r^{-2}/2)$, $r >0$, where

Suppose $\bm X\sim \text{\normalfont MVL}(\bm 0,
\mathsf{Q}^{-1})$. The margins of $\bm
X$ are standard Laplace distributed. The density $f_{\bm X}(\bm
x)$ is homothetic with respect to $r_{\G}$ and $f_0(r) = C_2
$r>0$, where
$K_\nu$ is the modified Bessel function of the second kind
<cit.> and $\nu =

$\bm X\sim \text{t}_{\nu}(\bm 0, \mathsf{Q}^{-1})$. The margins of
$\bm X$ are standard Student-$t_{\nu}$ distributed. The density
$f_{\bm X}(\bm x)$ is homothetic with respect to $r_{\G}$ and
$f_0(r) = C_3 \left(1+r^{-2}/\nu\right)^{-(\nu + d)/2}$, $r>0$,

Proposition <ref> immediately gives the star
set $\mathcal{W}$ in analytic form for Examples 1–3. As expected
from Proposition <ref>, the set $\mathcal{W}$
is common to all examples and is defined by its radial function $r_{\G}^{d}/(d\vol{G})$.

Left to right: Circular histogram of
$10^6$ angles sampled from a standard bivariate normal
distribution (correlation
$\rho=0.5$) with standard normal margins, standard Laplace
margins, and from a bivariate max-stable logistic distribution
with Laplace margins (dependence parameter
$\theta=0.5$). A concentric circle corresponds to a leap of 0.1
in the density of angles. Solid lines correspond to the boundary of
$\G^{d}/(d\vol{\G})$, for
$\G$ of the respective distributions. Plot 4: Sets
\W^{-1}$ corresponding to Plots 1–3.

Standard practice in extreme value analysis consists of removing the
effect of marginal scale to studying the extremal dependence structure
of $\bm
X$. This is typically achieved by standardising the margins of $\bm
X$ to a common distribution, see Section <ref>.Consider the map $\RR^d \ni \bm x \mapsto T(\bm x) =
\big(\Psi^{\leftarrow}(F_{X_i}(x_i)):i=1,\ldots,d\big) \in
\RR^d$ performing a transformation of the marginal distributions of
$\bm X$, so that the $i$th element of $T(\bm X)$ follows
$\Psi$ where
$\Psi$ denotes a continuous cumulative distribution function on
$\RR$. The transformation preserves the cardinality of the set of
input vectors $\bm X/\lVert \bm
X\rVert$, ensuring a one-to-one correspondence with the output vector
$T(\bm X)/\lVert T(\bm
X)\rVert$. From a geometric perspective, this means that the
densities of $\bm X/\lVert \bm X\rVert$ and $T(\bm X)/\lVert T(\bm
X)\rVert$ can be defined over the same set of points. However, while
the cardinality is preserved, the geometric properties of the
distribution of angles may change, as for example when
$T$ introduces nonlinearities which distort the shape of the
distribution of $\bm X/\lVert \bm
X\rVert$. As a concrete example, suppose $\bm
X$ is a random vector distributed according to a multivariate normal
distribution with zero mean and covariance matrix
$\mathsf{Q}^{-1}$ where $\Q$ is such that each
$X_i$ follows the standard normal distribution. Let $\bm X_L =
F_L^{-1}[\Phi(\bm X)]$, were $F_L$ and
$\Phi$ denote the cumulative distribution functions of standard
Laplace and standard normal random variables, respectively. Although
the density of $\bm X$ is homothetic, the density of $\bm
X_L$ is not. However, convergence
$(ii)$ of Proposition <ref> holds true with
$\rho=0$ outside a set of Lebesgue measure zero with
$G=\G_L$ determined by its radial function

\begin{equation}
r_{\G_L}(\bm w)= \left[\left\{\text{sgn}(\bm w) \vert \bm
w\rvert^{1/2}\right\}^\top \mathsf{Q} \, \left\{\text{sgn}(\bm
w) \vert \bm w\rvert^{1/2}\right\}\right]^{-1}, \qquad \bm
\label{eq:MVN_Laplace}
\end{equation}

where $\text{sgn}(x)=x/\lvert
x\rvert$ denotes the signum function and $\text{sgn}(\bm w) |\bm
|w_i|^{1/2}\,:\,i=1,\dots,d)$. Then, the star set
$\mathcal{W}_L$ describing the density of $\bm X_L/\lVert \bm
X_L\rVert$ is no longer a constant scale multiple of
$\G_L$, but instead, it is a radial product of
$\G_L$ with another star-body. The second panel of
Figure <ref> shows the true density (histogram) and
w)^d/\{d\vol{\G_L}\}$ (solid curve). The star-body is the ratio
between $\W$ and
$\G_L^d/\{d\vol{\G_L}\}$, and for this example, it is shown in the
fourth panel of Figure <ref>.

The class of homothetic densities, although rich, serves at best as an
idealistic setting. The analysis of the marginal transformation,
however, reveals an intermediate case where $\mathcal{W}$ can depend
both on $\G$ and on some additional latent star-body $\B$ that is
independent of $\G$. More generally, the star-body $\B$ captures
residual angular variation that is not explained by a homothetic
density. For instance, assume that $\bm X$ follows a bivariate
max-stable logistic distribution in standard Laplace margins with
density given by expression (<ref>). It
is clear from Figure <ref> that
$\W\neq\G^d/(d\vol{\G})$. For the practical case where both
$r_{\Astar}$ and $r_{G}$ are positive on $\SSS^{d-1}$, then
$\LL=\B\rmult\G$, $ \B, \G, \in \bigstar$ with $\B$ independent of
$\G$. Here, the set $\B$ captures the deviation of the set associated
with a homothetic unnormalised density of angles, $\G^d$, from the set
$\W$ associated with a density $f_{\bm W}$. It follows that in the
specific context where $\fL$ is homothetic with respect to $\rG$,
$\B=B_{1}(0)$. Plot 4 in Figure <ref> displays the
sets $d\vol{\G}\W/\G$ associated with the first three plots.

§.§ Multivariate radially-stable distributions

In this Section, we present the class of multivariate radially-stable (MRS)
distributions. The elements of this class are the only non-trivial
limits of radially renormalised exceedances above
$\QS_q$. Due to their stability properties, MRS distributions make them suitable for extrapolation
in the joint tail region. For the remainder of this section, we note
that $\RR^{d-1} \cong \SSS^{d-1}$ and $\RR_+^{d-1} \cong \SSS^{d-1}_+$
so that integration in (<ref>) is interpreted over the
spherical domain that is congruent to the domain of the distribution.Also, $K,K_1,K_2\in\bigstar$ all contain 0 and are such that
$\loc \subseteq K$ and $\loc \subseteq K_1 \subseteq K_2$.

A random vector $\bm Z$ is said to follow the multivariate
radially-stable exponential distribution with location
$\loc \in\bigstar$, scale $\Sigma\in\bigstar$, and shape
$\lvert \Astar^{1/d}\rvert=1/d$, if its distribution function is
given by

\begin{equation}              %IP 2 OCT
\PR(\bm Z\leq \bm z) = 1-\int\limits_{\SSS^{d-1}_+} \exp\left\{-\frac{\min(\bm z/\bm w)-r_{\loc }(\bm w)}{r_{\Sigma}(\bm w)}\right\}\,\rW(\bm w)~\dd \bm w, \quad \bm z\in\RR_+^{d}\setminus\loc , %\RR_+^{d-1}\cap
\label{eq:CDF_exp_pos_orthant}
\end{equation}

and if

\begin{equation}
\PR\left[\bm Z\in K^\prime\right] = \int_{\SSS^{d-1}}\exp\left\{-\frac{r_{K}(\bm w)-r_{\loc }(\bm w)}{r_{\Sigma}(\bm w)}\right\}\,\rW(\bm w)~\dd \bm w, \quad K'\subset\RR^{d}\setminus\loc .
\end{equation}

A notable property of random vectors $\bm Z$ following a MRS exponential distribution is that their quantile set $\QSq$ is given compactly by $\QSq = \loc \radd B_{-\log(1-q)}(0)\rmult\Sigma$. The distribution of $\bm Z$ also satisfies the stability property

\begin{equation}
\label{eq:stab_sets_exp}
\PR\left[\bm Z\in K_2^\prime\mid \bm Z \in K_1^\prime\right]= \PR\left[\bm Z\in (\loc \radd K_2\rsub K_1)^\prime\right].
\end{equation}

A random vector $\bm Z$ is said to follow the multivariate
radially-stable generalized Pareto distribution with location
$\loc \in\bigstar$, scale $\Sigma\in\bigstar$, and tail-index
$\xi \in \RR$, and shape $\Astar\in\bigstar$,
$\lvert \Astar^{1/d}\rvert=1/d$, if its distribution function is
given by

\begin{equation}
\PR(\bm Z\leq \bm z) = 1-\int\limits_{\SSS^{d-1}_+}
\left\{1+\xi\frac{\min(\bm z/\bm w)-r_{\loc }(\bm w)}{r_{\Sigma}(\bm w)}\right\}_+^{-1/\xi}\,\rW(\bm w) \,\dd \bm w, \quad \bm z\in\RR_+^{d}\setminus\loc ,
\label{eq:CDF_GP_pos_orthant}
\end{equation}

and if

\begin{equation}
\PR\left[\bm Z\in K^\prime\right]= \int\limits_{\SSS^{d-1}}
\left\{1+\xi \frac{r_{K}(\bm w)-r_{\loc }(\bm w)}{r_{\Sigma}(\bm w)}\right\}_+^{-1/\xi}\,\rW(\bm w) \,\dd \bm w, \quad K^\prime\in\RR^{d}\setminus\loc .
\label{eq:RingStab_exp}
\end{equation}

Analogously to the MRS exponential
distribution, we can obtain the quantile set $\QSq$ of a random vector
$\bm Z$ following a MRS generalized Pareto
distribution through
$\QSq = \loc \radd B_{\xi^{-1}\{(1-q)^{-\xi}-1\}}(0)\rmult\Sigma$. The
distribution of $\bm Z$ also satisfies the stability property

\begin{equation}
\label{eq:stab_sets_gP}
\PR\left[\bm Z\in K_2^\prime\mid \bm Z \in K_1^\prime\right] = \PR\left[\bm Z \in \left\{\frac{\loc \radd K_2\rsub K_1}{B_1(0)\radd B_{\xi}(0)\rmult(K_1\rsub\loc )\rdiv\Sigma}\right\}^\prime\right].
\end{equation}

The stability properties given in equations (<ref>) and (<ref>) and the expressions for the quantile set $\QS_q$ of both MRS distributions in terms of their starshaped parameters provide a powerful mean of extrapolating in tail regions. These properties are pivotal for the following section, introducing the notion of return level-sets, as well as Section <ref>, on probability estimation.

§.§ Return level-sets

The MRS developed in Section <ref> are central to the concept of return level-sets introduced in this section. Analogously to a univariate return level corresponding to a value expected to be exceeded by a random variable $X$ once in a certain return period, return level-sets correspond to a specific set expected to be attained by a random vector $\bm X$ once in a certain return period. Below, we detail the important properties of return level-sets, which give insight into the extremes of a random vector of interest.

For any Bernoulli experiment with probability $1-q$ of success, a
return period $T=1/(1-q)$ defines the expected number of experiments
needed to obtain one success. Given a random vector $\bm X\in\RR^d$
and an arbitrary set $B\in\RR^d$, the event $\bm X \in B$ can be
interpreted as a Bernoulli experiment with some return period $T$.Perhaps unsurprisingly, depending on the properties of $\bm X$, there
can be infinitely many distinct sets in $\RR^d$ satisfying this
property. If $\bm X$ has a Lebesgue density on $\RR^d$, then we can
always uniquely determine a canonical return set as in
Definition <ref>. In particular, two open subsets
$U_{1/(1-q)}$ and $V_{1/(1-q)}$ of $\RR^d\setminus\{0\}$ are said to
be equivalent with respect to a random vector $\bm Z$ and a measure
$\mathbb{P}$, denoted $U_{1/(1-q)}{\sim}^{\PR} V_{1/(1-q)}$, if there
exists $q\in (0,1)$ such that $\PR[\bm Z \in U_{1/(1-q)}] = 1-q$ and
$\PR[\bm Z \in V_{1/(1-q)}] = 1-q$.

The relation $\sim^\PR$ is an equivalence relation and the collection
of subsets of $\RR^d\setminus\{0\}$,
$\{U_{1/(1-q)}\}_{\sim^{\PR}} = \{\text{$V_{1/(1-q)}$ open in
$\RR^d\setminus\{0\}$,} \,:\, V_{1/(1-q)} \sim^\PR U_{1/(1-q)}\}$
is an equivalence class of $U_{1/(1-q)}$. An equivalence class
$\{U_{1/(1-q)}\}_{\sim^{\PR}}$ always contains a unique set
$\RS_{1/(1-q)}$ defined by

\begin{equation}
\RS_{1/(1-q)} = \RR^d \setminus \QSq, \quad q\in(0,1),
\label{eq:Rq}
\end{equation}

where $\QSq$ is the quantile set introduced in
Definition <ref>.

The set $\RS_T$ defined by (<ref>) is called a
return level-set with return-period $T$.

By Definition <ref>, a sequence of return
level-sets $\RS_{T_1},\ldots,\RS_{T_k}$ with $T_1\leq \cdots\leq T_k$ forms a nested sequence of open
sets. It thus follows that for return periods $T_i\leq T_j$ with $1\leq i\leq j \leq k$,
[Z∈_T_j|Z∈_T_i]   =    [Z∈_T_j]/[Z∈_T_i]=T_i/T_j.

From Section <ref>, we know that the MRS distributions described in Section <ref> arise as limiting distributions of suitably scaled excesses of a quantile set. We here describe the geometry of the return level-sets corresponding to these distributions.

If the random vector $\bm Z$ follows a MRS exponential distribution with location parameter $\loc$ and scale parameter $\Sigma$, the return level-set $\RS_{T_j}$ can be written in terms of another $\RS_{T_i}$ with $T_j\geq T_i$ via

_T_j    =    [_T_i^'B_{log(T_j/T_i)}(0)Σ]^'.

Similarly, if the random vector $\bm Z$ follows a MRS generalized Pareto distribution with location parameter $\loc$ and scale parameter $\Sigma$ and tail-index $\xi$, the return level-set $\RS_{T_j}$ can be written in terms of another return level-set $\RS_{T_i}$ with $T_j\geq T_i$ as

_T_j    =    [_T_i^' B_{(T_j^ξ-T_i^ξ)/ξ}(0)Σ]^'.

A classical
assumption in extreme value theory, phrased in our setting, is that
there exists $q\in(0,1)$ such that suitably scaled exceedances of the
quantile set $\QS_q$ follow the MRS distributions exactly. Under this
assumption, expressions (<ref>)
and (<ref>) for the return sets provide a
well-posed way of extrapolating far beyond
$\RS_{1/(1-q)}=\QS_q^\prime$. This assumption is central to the
upcoming Section <ref>, which details the statistical
inference procedure based on the limiting results provided in
Section <ref>.

§ STATISTICAL INFERENCE

§.§ Standardisation of margins

Statistical inference for real data, encoded into realisations from a
random sample $\{\bm X_i\,:\, i=1,\ldots,n\}$, is messy as the support
and tail properties of the marginal distributions of $\bm X_i$ may
behave in wildly differing manners. It is therefore convenient, and
frequently done in practice, to standardise observed data to common
marginal distributions prior to statistical modelling. In practice,
this is often achieved through the probability integral transform. In
this work, we model exceedances of a high quantile set using the MRS
exponential distribution after transforming the marginal distributions
the standard Laplace distributions. This leads to a relatively simple
and parsimonious model which covers a wide range of extremal
dependence structures <cit.>.

§.§ Quantile regression

Quantile regression methods are typically implemented using a pinball
loss function <cit.>, and not on a distributional
model for the density. Although <cit.> propose the use of the
asymmetric Laplace distribution for the model density due to the
equivalence of the negative log density of the asymmetric Laplace with
the pinball loss function, naively treating the asymmetric Laplace as
an adequate model for the data is precarious in a Bayesian setting.For example <cit.> show that the resulting
posterior prediction intervals have poor frequentist calibration
properties, and this is especially pronounced for tail quantiles,
which are essential in our setup. Second, the scale parameter of the
asymmetric Laplace distribution is arbitrary in a Bayesian framework
and even maximum likelihood based estimators are known to lead to
inacurrate quantile estimates, see <cit.>.
We adopt a generalized linear model based approach for
quantile regression, requiring an adequate distributional model for the density of $R\mid \bm W$. In particular, we
assume that $R\mid \bm W=\bm w$ follows a Gamma distribution and model
the logarithm of conditional $q$-quantile $\log \rquant(\bm w)$ using
a finite-dimensional continuously specified Gaussian process prior on
$\SSS^{d-1}$ (see Section <ref>). The latent function
$\rquant(\bm w)$ satisfies

\[
\PR[R \leq \rquant(\bm w)\mid W = \bm w] = q, \quad\text{for all
$\bm w\in \SSS^{d-1}$}.
\]

From total probability, we have

\begin{equation}
\PR[R\bm W \in \QSq]=\PR[R \leq \rquant(\bm W)] = q \int_{\SSS^{d-1}} f_{\bm W}(\bm w)\,\dd \bm w = q
\label{eq:prob_exc}
\end{equation}

There are at least two major strengths behind this choice. First, due
to Proposition <ref>, the Gamma distribution serves as a good
approximation for the density in the tail region of the distribution
of $R\mid \bm W$. Since estimators of high quantiles are not
influenced by the bulk of the distribution, a likely misspecification
between our choice and the true density of $R\mid \bm W$ in the body
of the distribution is not of concern. Second, our choice
exploits the form of the decay of the conditional density and allows
for Bayesian inference for the conditional quantile under a model that can adequately describe the behaviour in the tail.

In this paper, we adopt a Gamma quantile regression models due to the
simplicity that it affords, but remark that further improvements may
be achieved by using, for example, a truncated Gamma distribution, so
as to eliminate any effect from the body of the distribution of
$R\mid \bm W$.

§.§ Likelihood function

We assume that given a sample
$\dat =\{(r_i,\bm w_i)\,:\,i=1,\dots,n\}$ of $n$ observations, there
exists an integer $k_0 < n$ such that the limit expression
(<ref>) holds exactly $\forall k>k_0$. This type of
assumption is the standard practical assumption adopted in extreme
value statistics where the distribution of observed extremes is
approximated by the limiting theoretical distribution. In what
follows, we use this assumption to motivate a likelihood based
approach based on the Poisson process
convergence (<ref>). Given that $\bm X$ has standard Laplace margins, this approach leads to a MRS exponential limiting distribution with location
parameter $\QSq$ and scale parameter $\G$ as described in Definition <ref>.We prefer this approach over the analogous one based on
Definition <ref> for its parsimony and convenient memoryless stability
properties. In this section, the quantities of interest
are the limit set $\G$ and the set $\W$ with respective radial
functions $\rG$ and $f_{\bm W}$.

We recall from equation (<ref>) that the set $\W$ can be
expressed in terms of a starshaped set $\LL$ such that
$f_{\bm W}(\bm w)= r_{\LL}(\bm w)^d/(d\vol{\LL})$. Based on
Section <ref>, we motivate three models for $\W$
given by $\Mh:\LL=\B$, $\Mg:\LL=\G$, and $\Mgh:\LL=\B\rmult\G$, for
$\B,\G\in\bigstar$. Without loss, we write the model likelihood in terms of $\LL$.
Conditionally on the quantile set $\QSq$ and the data $\dat$, the
likelihood of $\bm \theta = (\G,\LL)$, is given by
$\smash{L}(\bm \theta \mid \QSq, \dat):= \prod_{i=1}^n f_{R, \bm
W}(r_i, \bm w_i\mid \bm \theta,\QSq)$. Letting
$\mathcal{S}_q:=\{i\in \{1,\dots,n\}\,:\, r_i > \rquant(\bm w_i)\}$
denote the set of random indices corresponding to exceedances of
$\QSq$, and
$\mathcal{S}_q^\prime:=\{1,\ldots,n\}\backslash \mathcal{S}_q$,
$\smash{L}$ can be expressed in terms of contributions of the radii
$\{r_1,\ldots,r_n\}$ above and below $\rquant$ via
L(θ|, ) = [∏_i∈𝒮_q^'f_R|R≤(W),W(r_i)f_W(w_i)][∏_i ∈𝒮_qf_R|R> (W),W(r_i)f_W(w_i)].

Based on Theorem <ref>, we assume that all
non-exceedance radii $\{r_i:i\in\mathcal{S}_q^{C}\}$ do not carry
information about $\G$ and $\LL$, and hence that $L$ is constant with
respect to them. If $\bm X$ has a Lebesgue density on $\RR^d$, then
equation (<ref>) suggests that under random
sampling, all angles $\{\bm w_1,\ldots,\bm w_n\}$ may be used in the likelihood. Hence, let $\mathcal{S}_{\bm w}$ denote the set of indices containing at least those exceedances and at most those of all observations— or
$\mathcal{S}_{q}\subseteq \mathcal{S}_{\bm w}\subseteq\{1,\ldots,n\}$. Then, we write
L(θ|,)    ∝
∏_i∈𝒮_q f_R|R> (W),W(r_i|w_i) ∏_i∈𝒮_w f_W(w_i)

=    exp{-|𝒮_w |log
}∏_i∈𝒮_qexp[-{r_i - (w_i)}/(w_i)]/(w_i)∏_i∈𝒮_wr_(w_j)^d.

In practice, a key complication arises
when using the likelihood in expression (<ref>).Evaluating the likelihood function (<ref>) requires
computing the constant $d\vol{\LL}$ which ensures the density of
$\bm W$ integrates to one. This normalising constant is in most cases impossible to compute exactly, which makes
inference difficult. However, using the Poisson transform
<cit.>, we can map the likelihood into an equivalent
likelihood function $L( \bm \theta,\beta \mid \QSq,\dat)$ of a Poisson
point process defined in an expanded space given by

\begin{align}\label{eq:likelihood-v2}
L(\bm \theta,\beta)
= \exp\left[-\lvert \mathcal{S}_{\bm w} \rvert
e^\beta (d\vol{\LL})\right] \prod_{i\in
\mathcal{S}_q}\frac{1}{\rG(\bm w_i)}\exp\left(-\frac{r_i - \rquant(\bm w_i)}{\rG(\bm
w_i)}\right) \prod_{i\in\mathcal{S}_{\bm w}}e^\beta r_{\LL}(\bm w_i)^d,
\end{align}

where the latent variable $\beta$ estimates the normalising
constant $d\vol{\LL}$. The normalising constant is inferred as just
another parameter at no loss of information
<cit.>, whilst inference can be performed either
using frequentist methods or in a fully Bayesian manner by assigning
suitable prior distributions on $\G$, $\LL$, and $\beta$
<cit.>. We adopt the latter approach due to the
optimality of Bayes as an information processing rule
<cit.>. The volume $d\lvert \LL\rvert$ is
computed using the computationally efficient numerical integration
scheme of <cit.>.

Hierarchical structure of
models $\Mh$ (left), $\Mg$ (centre), and $\Mgh$ (right). The
vectors $\alpha_{\QSq}$, $\alpha_{\G}$, and $\alpha_{\LL}$ are hyperparameters for the prior distributions of the latent
variables $\QSq$, $\G$ and $\LL$. The variable $R_E$ denotes the excesses $R-\rquant(\bm W) \mid R>\rquant(\bm W)$.

Based on the results from Section <ref>, we
motivate the following models for the density of angles
$f_{\bm W}(\bm w)$: $\Mh: \LL = \B$, $\Mg: \LL = \G$, and
$\Mgh: \LL = \B\rmult\G$, for $\B \in \bigstar$ independent of
$\G\in\bigstar$. We note that the latent star-bodies $\B$ in models
$\Mh$ and $\Mgh$ do not have the same interpretation, but the slight
abuse of notation allows to simplify the notation for the parameter
space in later sections without loss of interpretation.Figure <ref> shows the hierarchical structure of
models $\Mh$, $\Mg$ and $\Mgh$ via their directed acyclic conditional
independence graphs.
It is apparent from Figure <ref> that the graph of
$\Mg$ forms a sub-graph of $\Mgh$, and that $\Mgh$ is the graph sum of
$\Mh$ and $\Mg$. The nested structure of model $\Mg$ within the
parameter space of model $\Mgh$ translates into a bias-variance
trade-off as the latter offers additional flexibility at the cost of a
possibly increased variance for the latent set $\G$. A similar
trade-off occurs for $\G$ between models $\Mh$ and $\Mgh$ since the
former ignores possible information contained in the observed angles.Inference for the latent set $\B$ of model $\Mgh$ as well as
calibration checks for the posterior predictive distribution obtained
from each model yield possible methods for model selection. This is
discussed in Section <ref>. Next, we delve into the
prior distribution for a specific subclass of star-shaped objects in
both $\RR^2$ and $\RR^3$, along with the corresponding prior
distribution for its hyperparameters.

§.§ Inference for latent variables

Proposition <ref> and

Given $n$ independent
observations $\dat:=\{r_1\bm w_1,\ldots,r_n\bm w_n\}$ from
$\bm X_L=R\bm W\in \RR^d$, the latent quantities of interest are the
limit set $\G$, the quantile set $\QS_q$ and the set $\LL$. In what
follows, we detail the procedure to obtain realisations from the joint
posterior distribution of these quantities.
We model the logarithms of the random radial functions $\rquant$,
$\rG$ and $r_{\LL}$ as Matérn (Gaussian) fields on $\SSS^{d-1}$
using the stochastic partial differential equation (SPDE) approach by
<cit.>, using $\alpha=2$ (see their Equation (2)) which is
also the default option in the  package
(). We denote the respective random intercepts
$\beta_{\QS},\beta_{\G},\beta_{\LL}\in\RR$ as well as their stochastic
weights in the finite element representation
<cit.> as
$\bm z_{\QS},\bm z_{\G},\bm z_{\LL}\in\RR^p$.

We denote by $\bm \theta=(\beta_{\QS},\bm z_{\QS},\beta_{\G},\bm z_{\G},\beta_{\LL},\bm z_{\LL})\in\RR^{3p+3}$ the set of all latent variables for our models, and it follows that the joint posterior distribution of $\bm \theta$ fully determines that of $\QSq$, $\G$ and $\LL$. Due to the hierarchical structure of all proposed models, as detailed in Sections <ref> and <ref>, the joint posterior
density of $\bm \theta$ factorises according to $\pi[\bm \theta \mid \dat]= \pi[\beta_{\G},\bm z_{\G}, \beta_{\LL},\bm z_{\LL} \mid \beta_{\QS},\bm z_{\QS},\dat]\pi[\beta_{\QS},\bm z_{\QS}\mid \dat]$.We first fit the Bayesian Gamma log-linear quantile regression model
described in Section <ref> for a fixed probability $q$
to all observations $\dat$. Samples
$\{(\beta_{\QS,i},\bm z_{\QS,i}):i=1,\ldots,n_{\QS}\}$ from the posterior density
$\pi[\beta_{\QS},\bm z_{\QS}\mid \dat]$ map to a set of radial functions $\{\rquants{i}:i=1,\ldots,n_{\QS}\}$ which are then interpreted as
candidate radial functions of a $q$-th return set.
Given a radial
function $\rquants{i}$, we define a set of exceedances through

\begin{equation}
\label{eq:set_Exc}
\setExc_i = \left\{\left(r_j,\bm w_j\right)\in (0,\infty)\times\SSS^{d-1}\;:\; r_j>\rquants{i}(\bm w_j),\, r_j\bm w_j \in \RR^{d},\,j=1,\ldots,n\right\}.
\end{equation}

We then fit model $\Mh$, $\Mg$, or $\Mgh$ to each collection of
exceedances $\setExc_i$, as detailed in Section <ref>,
and possibly include angles from non-exceedance data in the
likelihood. This procedure yields a conditional posterior density
$\pi[\beta_{\G},\bm z_{\G}, \beta_{\LL},\bm z_{\LL} \mid
\beta_{\QS,i},\bm z_{\QS,i},\dat]$ for each ${i=1,\ldots,n_{\QS}}$.Sampling $n_{\G \LL}$ realisations jointly from each
$\pi[\beta_{\G},\bm z_{\G}, \beta_{\LL},\bm z_{\LL} \mid
\beta_{\QS,i},\bm z_{\QS,i},\dat]$ provides an assembled sample of
$n_{\QS}\cdot n_{\G \LL}$ realisations from the joint posterior
distribution of $\bm \theta$,

\begin{equation}
\label{eq:samp_theta_ij}
\left\{\bm \theta_{i,j}=\left(\beta_{\QS,i},\bm z_{\QS,i},\beta_{\G,(i,j)},\bm z_{\G,(i,j)},\beta_{\LL,(i,j)},\bm z_{\LL,(i,j)}\right)\in\bm \Theta\;:\;i=1,\ldots,n_{\QS},
j=1,\ldots,n_{\G \LL}\right\}.
\end{equation}

For simplicity and without loss, we re-index the
sample (<ref>) from the posterior distribution of
$\bm \theta$ to $\{\bm \theta_{i}\,:\,i=1,\ldots,n_{\bm \theta}\}$
with $n_{\bm \theta}:=n_{\QS}\cdot n_{\G \LL}$ and use this notation
in the next sections. We shall also refer to the sampled latent
functions $\rquants{i}$, $\rGs{i}$, and $r_{\LL,i}$ constructed from
$\bm \theta_i$ using a piecewise linear basis function
representation of the logarithms of the radial functions.

Interval estimation for the latent fields $\rquant$, $\rG$, and $r_{\LL}$ is
accomplished via prediction
intervals <cit.>. A $(1-\alpha)$
prediction-interval for the value of a random field
$Y\,:\,\Omega \to \smash{\RR^{\SSS^{d-1}}}$ at an angle $\bm w\in\SSS^{d-1}$
is the closed line segment
$[q_\alpha(\bm w) \bm w : q_{1-\alpha}(\bm w) \bm w]$
where $q_{\alpha}(\bm w)$ is the $\alpha$-quantile of the distribution
of $Y(\bm w)$. A $(1-\alpha)$ prediction interval for the process $Y$
defined on $\SSS^{d-1}$ consists of the strip
$R_{1-\alpha}:=\cup_{\bm w\in \SSS^{d-1}}R_{1-\alpha}(\bm w)$ defined through $R_{1-\alpha}(\bm w):=[q_\rho(\bm w) \bm w :
q_{1-\rho}(\bm w) \bm w]$
for some $\rho$ such that $q_{\rho}(\bm w)$ and $q_{1-\rho}(\bm w)$ satisfy

\[
\PR\left[ q_\rho(\bm w) \leq Y(\bm w)\leq q_{1-\rho}(\bm w),\bm w\in\SSS^{d-1}\right] =
\]

In the context of our latent fields $\rquant$, $\rG$, and $r_{\LL}$,
$(1-\alpha)$ prediction interval consist of sets withing within which
the true functions lie entirely with probability $1-\alpha$.Prediction intervals can be obtained from a sample from the posterior
distribution of a parameter of interest using the
package in .<cit.>.

§.§ Rare event probability estimation

Given $n$ independent observations
$\dat:=\{r_1\bm w_1,\ldots,r_n\bm w_n\}$ from
$\bm X_L=R\bm W\in \RR^d$, interest lies in the estimation of the
probability that a new draw from $R\bm W$ falls within some Borel set
$B\in\RR^d\setminus\{0\}$ far enough from the origin ($0\in\RR^d$).We provide a framework for inference for
${\PR_{B\mid\dat}:=\PR[R\bm W\in B\mid \dat]}$. For practical
reasons, we restrict ourselves to sets $B$ that are starshaped at $0$
(see Appendix <ref>).

This is useful in our setting as it allows for an exact probability
calculation with respect to our model specification for exceedances
(see expression (<ref>) below). In practice, common
types of sets of interest are boxes
$\{\bm x \in \RR^d: \bm a \leq \bm x \leq \bm b,\, \bm a,\bm
b\in\RR^d\}$ and
sets of the form
$\{r\bm w :r>h(\bm w)>0,r\in(0,\infty),\bm w \in \SSS^{d-1}\}$ for some positive function $h$ defined on $\SSS^{d-1}$, both
starshaped at $0$, and our choice is hence not too restrictive.Inference for more complex sets is amenable through simple
modifications of what follows.

For a set $B\in\RR^d\setminus\{0\}$ starshaped at $0$, define $S_B:=\{\bm w \in \SSS^{d-1}: r\bm w\in B, r\in(0,\infty]\} \subseteq
\SSS^{d-1}$ and consider, for any $\bm w \in S_B$, the
partition $I_{\inf}(\bm w)\cup I_B(\bm w)\cup I_{\sup} (\bm w)$ of
$(0,\infty)$, where

$I_{\inf}(\bm w):=(0,r_{B_{\inf}}(\bm w))$,
$I_B(\bm w):=[r_{B_{\inf}}(\bm w),r_{B}(\bm w)]$, and
$I_{\sup}(\bm w):=(r_{B}(\bm w),\infty)$ for the radial function
$r_{B}(\bm w)$ of $B$ and the function $r_{B_{\inf}}(\bm w)=\inf\{r>0:r\bm w \in B\}$. Then,

_B|   =    [R∈I_B(W),W ∈S_B|],

The posterior predictive distribution in expression (<ref>) is given from the posterior density of $\bm \theta$ via

_B|   =    ∫_^3p+3[R∈I_B(W),W ∈S_B|θ]π[θ|] dθ,

where $\pi[\bm \theta\mid\dat]$ is obtained following the procedure described in Section <ref>.We now express $\PR_{B\mid \bm \theta}:=\PR[R\in I_{B}(\bm W),\bm W \in S_B\mid \bm \theta]$ in terms of our model formulation and assumptions.

Unlike standard
approaches that base inference on a fixed threshold and do not take into
account the uncertainty in its estimation, all asymptotically motivated models proposed in Section <ref> are defined in the random tail-region $R>\rquant(\bm W)$ where $\rquant(\bm w)=\exp\{\beta_{\QS}+\bm \psi(\bm w)^\top\bm z_{\QS}\}$ is a latent radial function capturing the uncertainty for the threshold set $\QSq$. It follows that $\rquant(\bm w)$ intersects with any of the intervals $I_{\inf}(\bm w)$, $I_B(\bm w)$, or $I_{\sup}(\bm w)$ with non-zero probability for $\bm w\in S_B$. Using total probability, we thus split the calculation of $\PR_{B\mid \bm \theta}$
into the sum

_B|θ= [RW ∈B, R>(W)|θ] + [RW ∈B, R≤(W)|θ],

and our model is best suited for sets $B$ such that the second term in the
equation (<ref>) is small. The first term in equation (<ref>) decomposes into

[RW ∈B, R>(W)|θ]= [R∈I_B(W)|W ∈S_B, R>(W),θ]×

×  [W ∈S_B |R>(W),θ][R>(W)|θ].
By the assumptions detailed in Section <ref>, the last term $\PR[R>\rquant(\bm W)\mid \bm \theta]$ in
equation (<ref>) equals $1-q$ for all $\bm \theta\sim\pi[\bm \theta\mid\dat]$.
Following our model formulation for $\bm W$, the second term in
equation (<ref>) corresponds to the predictive
distribution of the angles given $\bm \theta$ and is given by

[W ∈S_B |R>(W),θ]   =
where $\widetilde{f}_{\bm W}$ is chosen from $\Mh$, $\Mg$, or $\Mgh$
(see Section <ref>), and
$\rG(\bm w) = \exp\{\beta_{\G}+\bm \psi(\bm w)^\top\bm z_{\G}\}$,
$r_{\LL}(\bm w) = \beta_{\LL}+\bm \psi(\bm w)^\top\bm z_{\LL}$. In practice,
the integrals in equation (<ref>) are obtained
efficiently via numerical integration for $d=2$ or $3$. The first
term in equation (<ref>) is obtained through

[R∈I_B(W)|W ∈S_B, R>(W),θ]=_W|W∈S_B,([R∈I_B(W)|R>(W),θ]),

where the term $\PR[R\in I_B(\bm w)\mid \bm W = \bm w, R>\rquant(\bm w),\bm \theta]$, needed for the expectation calculation in equation (<ref>), is equal to

Equation (<ref>) follows from our model formulation for
exponential excesses with scale $\G$ and the memoryless property of
the exponential distribution. In practice, we approximate
expectation (<ref>) through Monte-Carlo integration by
sampling angles
$\{\widetilde{\bm w}_1,\ldots,\widetilde{\bm w}_{n_w}\}$ from
$f_{\bm W}$ restricted to the set $S_B$ via an
acceptance-rejection scheme and calculating

n_w^-1∑_i=1^n_w[R∈I_B(w_i)|R>(w_i),θ]    →    [R∈I_B(W)|W ∈S_B, R>(W),θ],

as $n_w\to\infty$. The second term in equation (<ref>),
$p:=\PR[R\bm W \in B, R\leq \rquant(\bm W)\mid \bm \theta]$,
corresponds to a probability over a region in which our model is not
specified. The probability in question is thus assumed to be the
success probability of a binomial experiment with size $n$ for the
number of observations falling in the set $\QS_q \cap B$. In practice, for each
$i\in\{1,\ldots,n_{\bm \theta}\}$, we impose a binomial likelihood and
beta prior for the parameter $p$, leading to a closed form beta
posterior distribution. When $\QS_q \cap B\neq \varnothing$, we sample a probability from the posterior distribution of
$p$ and add it accordingly in equation (<ref>).
In practice, inference for $\PR_{B\mid \dat}$ is performed by first collecting a sample $\{\bm \theta_1,\ldots,\bm \theta_{n_{\bm\theta}}\}$ from $\pi[\bm \theta\mid\dat]$ through the procedure described in Section <ref>, and then calculating $\PR_{B\mid \bm \theta_i}$ for each $i\in\{1,\ldots,n_{\bm \theta}\}$. A posteriori mean estimates corresponding to the Monte-Carlo integration

n_θ^-1∑_i=1^n_θ_B|θ_i   →    _B|,   as  n_θ →∞,

approximating equation (<ref>), and equal-tailed prediction intervals from the sample $\{\PR_{B\mid \bm\theta_i}:i=1,\ldots,n_{\theta}\}$ hence complete the inference procedure for $\PR_{B\mid \dat}$.

§.§ Model selection and validation

Our model formulation gives rise to various modelling choices needing to be assessed and validated. For instance, as discussed in Section <ref>, selecting the most appropriate model within the set of candidate models $\Mh$, $\Mg$, and $\Mgh$ given observed data amounts to analysing the properties of the distribution of $R,\bm W\mid R>\rquant(\bm W)$ with respect to the association between $\gG$ and $f_{\bm W}$. Another modelling choice requiring assessment is that of the set of angles contributing to likelihood (<ref>) since it translates into a bias-variance trade-off for $\G$ and $\LL$. Further, hyper-parameters values imposed on the prior distributions of $\QSq$, $\G$, and $\LL$ imply a priori information on the differentiability properties of their boundaries (see Section <ref>), and their impact on the posterior may need assessment. Finally, the usual concern of the sensitivity of the posterior distribution of the latent variables to the return period for the latent threshold function $\rquant$ (Section <ref>) remains. We detail possible methods for model selection and validation below.

The nested structure of model $\Mg$
within the parameter space of model $\Mgh$ raises the question of
whether model $\Mg$ can serve as a sensible simplification of model
$\Mgh$ given observed data
$\dat = \{r_1\bm w_1,\ldots, r_n\bm w_n\}$. As discussed in
Section <ref>, evidence for a constant $r_{\B}$ in model $\Mgh$ points to evidence that $\Mg$ can serve as an
sensible simplification of $\Mgh$, or equivalently that the observed
data $\dat$ may come from a process following a homothetic density with respect to its gauge function $\gG$.Given a sample of random functions $\{\psi(\bm w)^\top\bm z_{\B,i}:i=1,\ldots,n_{\theta}\}$
from the posterior distribution of the latent random field
$\psi(\bm w)^\top\bm z_{\B}$, we obtain a $0.95$
prediction interval, $R_{0.95}$, for $\psi(\bm W)^\top\bm z_{\B}$ through the procedure described in Section <ref>.
We say there is evidence against a homothetic underlying desnity whenever there is some non-empty subset $S_0\subset\SSS^{d-1}$ such that $0\notin R_{0.95}(\bm w)$ for $\bm w \in S_0$.
The sensitivity of the posterior distributions of the latent
parameters to the remaining modelling choices discussed in the
beginning of this section can be assessed through the quality of the
calibration of the posterior predictive distribution to the observed
data. Given a set of observed exceedances $\setExc_i$ of a sampled
latent function $\rquants{i}$ (see expression (<ref>)), we
wish to assess the calibration of the predictive distribution for the
excess variable $R-\rquants{i}(\bm W)\mid R>\rquants{i}(\bm W)$. We thus
compare each observed excess $r_j-\rquants{i}(\bm w_j)$ with the
quantile of the Exponential$(\rG(\bm w_j))$, with
$\log\rG(\bm w)=\beta_{\G} + \bm \psi(\bm w)^\top \bm z_{\G}$ and
$(\beta_{\G},\bm z_{\G})\sim\pi[\beta_{\G}, \bm
z_{\G}\mid\beta_{\QS,i},\bm z_{\QS,i},\dat]$. We draw exceedance
radii $\{\widetilde r_1,\ldots,\widetilde r_{n_s}\}$ from their
predictive distribution Exp$(\rG(\bm w_j))$, and define the empirical
distribution function
$F_{\bm w_j,n_s}(r):={n_s}^{-1}\sum_{k=1}^{n_s}\mathbbm{1}[\widetilde
r_k\leq r]$, for $r\in(0,\infty)$. A probability-probability (PP)
plot for the exceedances $\setExc_i$ is then given by

{(j/|_i |+1 ,p_(j)):j=1,…,|_i|},

where $p_{(j)}$ denotes the $j$-th order statistic of the sample $\{F_{\bm w_j,n_s}(r_j-\rquants{i}(\bm w_j)):j=1,\ldots,\lvert\setExc_i\vert\}$. A quantile-quantile (QQ) plot in unit exponential margins for the exceedances $\setExc_i$ is then easily obtained via
probability-probability (PP) plot for the exceedances $\setExc_i$ via

{(-log(1-j/|_i |+1) ,-log(1-p_(j))):j=1,…,|_i|}.

In a similar fashion, we can obtain PP plots for the predictive distribution of $\bm W\mid R>\rquants{i}(\bm W)$. This is achieved by sampling angles $\{\widetilde{\bm w}_1,\ldots,\widetilde{\bm w}_{n_s}\}$ from the predictive density $f_{\bm W\mid \beta_{\G},\bm z_{\G}, \beta_{\B},\bm z_{\B},\dat}$ where $(\beta_{\G},\bm z_{\G}, \beta_{\B},\bm z_{\B})\sim\pi[\beta_{\G},\bm z_{\G}, \beta_{\B},\bm z_{\B} \mid \beta_{\QS,i},\bm z_{\QS,i},\dat]$. We transform the sampled angles and the observed angles from $\setExc$ to spherical coordinates, respectively denoted $\{\widetilde{\bm \varphi}_j=(\widetilde{\varphi}_{1,j},\ldots,\widetilde{\varphi}_{d-1,j})\in\Phi:j=1,\ldots,n_s\}$ and $\{\bm \varphi_j=(\varphi_{1,j},\ldots,\varphi_{d-1,j})\in\Phi:j=1,\ldots,\lvert\setExc_i\rvert\}$, and consider the empirical distribution function $F_{n_s}(\bm \varphi):={n_s}^{-1}\sum_{k=1}^{n_s}\mathbbm{1}[\widetilde{\bm \varphi}_k\leq \bm \varphi]$. A PP plot for the calibration of the predictive distribution of angles to the observed angles in $\setExc_i$ is then given by expression (<ref>) with $p_{(j)}$ denoting the $j$-th order statistic of the sample $\{F_{n_s}(\bm \varphi_j):j=1,\ldots,\lvert\setExc_i\vert\}$.

§ DATA ANALYSES

§.§ Thames tributary river flow

We test our model-fitting on a bivariate dataset consisting of river flow
measurements for the rivers Pang and Windrush, both tributaries of the
River Thames in England, measured in cubic meters per second, $m^{3}s^{-1}$.
For each site, have daily measurements from 1968 to 2008 are recorded, and we note that the data is
already provided in Laplace margins, bypassing the need for marginal transformations. This dataset has been previously analyzed by <cit.> who advocated the use of the
Laplace distribution in the context of modelling extreme values using
the conditional approach to multivariate extremes first introduced in <cit.>. Given that there are a total of $n=14545$ bivariate
measurements, we have sufficient data to fit a hierarchical Bayesian Gamma quantile regression model for the radial function $\rquant$ at the high quantile of $q=0.95$. Once fitted, we
sample 20 radial functions $\rquants{i}$, $i=1,\dots,20$, from the posterior distribution of the threshold $\pi[\beta_{\QS},\bm z_{\QS}\mid \dat]$. Figure <ref>
shows the fitted posterior radial quantile set $\QSq$ and it is observed that our posterior mean estimate of the anglewise threshold capture the quantile of radii at given angle values with little uncertainty.

Left to right: Posterior mean of $\partial\QS_q$ with $0.95$
predictive intervals. Posterior mean of the boundary of the limit set, $\partial\G$, with $0.95$
predictive intervals. Histogram of observed angles with
posterior mean density of angles and $0.95$ prediction
intervals. Empirical (solid line) and posterior mean (dashed line) of $\chi_q$ values with 0.95 confidence intervals (light grey) and model predictive intervals (dark grey).

At each of these sampled threshold curves $\rquants{i}$ we fit the joint model for the gauge function
and the angular density via the
fitting of the likelihood function (<ref>), and subsequently
sample 50 posterior values of the gauge function
$\gG$ evaluated at a mesh for the unit sphere $\SSS^1$ for ${n_{\bm{\theta}}=1000}$ posterior samples in total. In all, we fit 6 different models for the density of angles, $\Mh$, $\Mg$, and $\Mgh$, where each model is fitted using either exceedance angles and the entire dataset's angles.

From Figure <ref> in Appendix <ref>, model $\Mg$ fitted on exceedance angles only produces a posterior estimate of the gauge function whose unit level set best captures the shape of the $\log({n}/{2})$-scaled sample cloud of data in Laplace margins with the least amount of uncertainty. Furthermore, Figure <ref> shows that this model provides a good approximation of the density of angles, but not the best approximation. Despite this, we say that model $\Mg$ fitted on exceedance angles is best overall. We proceed by only considering this model.

It is of interest to assess how well our joint modeling procedure approximates the distribution of angles in the dataset. As shown in Figure <ref>, the posterior mean curve overlays nicely with the empirical histogram of angles from the dataset, suggesting our joint modeling of exceedances and angles provide well-fitting posteriors.
The PP plot in Figure <ref> of Section <ref> shows good agreement between the posterior predictive distribution and the observed data, further justifying the performance of our model.

Figure <ref> shows posterior estimates of $\chi_q = \chi_q(1,2)$ for values of $q$ approaching 1. The agreement with empirical and posterior mean model-based estimates  as $q$ varies, as well as the overlapping confidence intervals, show the ability of our posterior model to capture the extremal dependence of this dataset.

Furthermore, the return level boundary curves in Figure <ref> in Section <ref> accurately portray the intricacies of the dataset as angles vary along the unit simplex
$\SSS^1$. Above each curve corresponding to a value $T$, we are given an estimated region in which we expect to observe joint river flow measurements with probability ${1}/{T}$. Based on the size of the dataset, we expect to see a small number of points from the initial dataset beyond the return level curve corresponding to $T=10^{3}$, which is confirmed by Figure <ref>. The increased uncertainty in the further return level sets reflect the further extrapolation needed to describe events occurring in these regions.
Lastly, the QQ plots for exponential exceedances in Section <ref> of the supplementary material shows there is good agreement with the empirical and model-based estimates in their abilities to extrapolate extreme values.

The intricate extremal dependence structure displayed in this data set illustrates the
flexibility of our approach which
provides with a reasonable fit. A crucial remark, however, is that the
assumption of independence between observations does not hold here.Estimates of the probability that extremes of water flow occur
simultaneously at the rivers Pang and Windrush based on this
assumption are hence likely to suffer from severe bias. In what
follows, we discuss a different modelling question to which our method
is perhaps better suited; we investigate the use of our models to
characterise the temporal extremal dependence properties of river
flows. We analyse the pair $(x_t,x_{t+\ell})$ where $x_t$ denotes the
measured river flow at river Pang during day $t$, with $\ell$ denoting
a time lag in days.

Posterior mean (black line) and $0.95$-simultaneous predictive intervals (grey bands) for the limit set $\G$ associated with the pair $(X_{t},X_{t+\ell})$. From left to right: $\ell=1,2,7,30$.

From the $n=14545$ river flow measurements at river Pang, we form the $n-\ell$ pairs $(x_t,x_{t+\ell})$. We treat the pairs as observations from $(X_t,X_{t+\ell})$ fit the bivariate model $\Mgh$ to obtain samples from the posterior distribution of the latent parameters as described in Sections <ref> to <ref>.
Figure <ref> displays the posterior mean and the $0.95$-simultaneous predictive interval for the limit set $\G_{\ell}$ associated with the pair $(X_{t},X_{t+\ell})$, $\ell=1,2,7,30$.
The increasing distance between the posterior distribution of $\partial \G_\ell$ along the angle $(1/\sqrt{2},1/\sqrt{2})$ and the point $(1,1)$ as $\ell$ increases is expected, and corresponds to weakening positive extremal dependence with time. This behaviour can also be observed in terms of the decreasing sub-asymptotic
coefficient of tail dependence $\chi_{t,t+\ell}(q)$ as a function of time lags between measurements, as seen in Figure <ref>. It displays the mean and $0.95$ prediction
intervals based on the posterior distribution of the coefficient
$\chi_{t,t+\ell}(u)$ for
the different threshold values $u=0.95$, $0.9875$, and $0.9975$. It is also observed that extremal dependence pertains for
as long as a month, implying that extremes of river flows may last or
cluster for long periods of time.

Coefficient $\chi_{t,t+\ell}(q)$ for $u=0.95$, $0.9875$,
$0.9975$ from dark to light grey. Left: Empirical estimates and
$0.95$ confidence intervals. Right: Model based posterior mean
and $0.95$-simultaneous predictive intervals.

§.§ Newlyn wave heights

We now proceed to fit a model on a dataset with dimension $d=3$. The
dataset consists of hourly measurements of wave characteristics: wave
height ($H$, in meters), wave period ($P$, in seconds), and surge
($S$, in meters). measured over the period 1971–1977 at the Newlyn
port in south-west England. The Newlyn wave height dataset was first
analysed in <cit.> in the case where asymptotic dependence
was assumed between all 3 variables. Like in the previous literature,
we take componentwise maxima over 15-hour periods, resulting in a
dataset of size $n=2,894$. The margins $X_H, X_P, X_S$ of the data
are unknown, and are therefore transformed to standard Laplace
$ \smash{X_{L,j} = F^{-1}(\widehat{F}(X_j)) }$ for $j=H,S,P$ where $F^{-1}$
is the quantile function of the univariate standard Laplace
distribution, and

$\smash{\widehat{F}}$ is given by an empirical cdf below a high marginal threshold and by a GP distribution above a high marginal threshold. This marginal operation was previously used in <cit.> for transforming to unit Fréchet margins, and was implemented in the context of the Newlyn wave dataset in <cit.> for transforming the margins to standard exponential.
We note that in practice, the GP parameters were estimated using maximum likelihood in a frequentist manner. In future work, a fully Bayesian technique can be used to account for the uncertainty from the estimation of these marginal parameters. For the time being, this additional layer of complexity is neglected.

Left: Estimated posterior mean unit level set of the gauge function with $\log({n}/{2})$-scaled data in Laplace margins. Centre: Posterior mean return set with $T=10^{3}$ and data in Laplace margins. Right: Posterior mean return set with $T=10^{4}$ and data in Laplace margins.

Subsequently, we fit the hierarchical Bayesian gamma quantile regression model described in
Section <ref> at the $q=0.9$ quantile to obtain the posterior distribution $\pi[\beta_{\QS},\bm z_{\QS}\mid\dat]$.

For the subsequent joint model of exceedances and angles, it is deemed through our model checks that $\Mg$ fitted on exceedance angles only provides the best overall joint fit (see Section <ref> of the supplementary material), and we are left with the joint posterior model $\pi[\beta_{\G},\bm z_{\G}\mid\beta_{\QS},\bm z_{\QS}, \dat]$.
As in the fiver flow analysis of Section <ref>, we sample 20 observations from the posterior threshold model $\pi[\beta_{\QS},\bm z_{\QS}\mid\dat]$, and at each threshold, we sample 50 observations from the joint model of exceedances and angles $\pi[\beta_{\G},\bm z_{\G}\mid\beta_{\QS},\bm z_{\QS}, \dat]$, resulting in $n_{\bm{\theta}}=1000$.
Figure <ref>
shows the unit level set for the model $\Mg$ with the posterior mean boundaries of the $T=10^3 , 10^4$ return level sets, as described in Section <ref>. We see
that we are able to capture the
extremal dependence structure of our dataset quite well by the clear
correspondence between the posterior mean unit level set of $\gG$ and
the scaled sample cloud.

$\chi_q$ plots for all pairwise combinations of marginal variables and for all three marginal variables. Solid line is an empirical estimate with 95% bootstrap confidence intervals in light grey, and the dashed line is a model-based estimate with 95% posterior confidence intervals in dark grey.

Figure <ref> shows posterior mean estimates of $\chi_q(A)$ for $A\in\left\{HP, HS,PS, HPS\right\}$ and for $q\in(0.9,1)$ with posterior 95% confidence intervals.

<cit.> establish that the variable groups $HP$, $PS$, and $HPS$ are asymptotically independent, and that $HS$ is asymptotically dependent. With this in mind, we conclude that our posterior joint model is able to accurately describe the extremal dependence structure presented in this dataset. Figure <ref> correctly shows posterior estimates of $\chi_q(A)$ tending to 0 as $q$ tends to $1$ for the asymptotically independent groups, and tending to a nonzero value close to the empirical estimate for the asymptotically dependent pair $HS$.

Furthermore, the QQ plots in Figures <ref> in Section <ref> of the supplementary material show that there is good agreement with the empirical and model-based estimates in their abilities to extrapolate exceedances to extreme values.
Figure <ref> in Section <ref> of the supplementary material shows PP plots for angles generated from the posterior density of angles from $\pi[\beta_{\LL},\bm z_{\LL}\mid\beta_{\QS},\bm z_{\QS}, \dat]$ have good agreement with the dataset's empirical angular distribution.

Lastly, for random vectors in the dataset's original margins $\bm{X}=\left(X_H,X_S,X_P\right)^\top$, it is of interest to estimate the return level for the overtopping discharge rate $Q(v;\bm{X})$ over a wall of height $v$. First introduced in <cit.>, $Q$ is interpreted as the volume of water (in cubic meters, $m^3)$ overtopping the sea-wall per unit length (in meters, $m$) over a fixed duration (in seconds, $s$) and is thus measured in $m^3 s^{-1} m^{-1}$. For the application at hand, $Q(v;\bm{X}) := a_1 X_S X_P \exp\{{a_2\left(v-X_S-l\right)}/{(X_P {X_H^\ast}^{{1}/{2}})}\}$. The wave height component $X_H^\ast$ is a calibration of the wave height marginal variable $X_H$ to approximate the actual off-shore wave height, since measurements are taken on-shore.

Estimated return levels for sea-wall height. Presented are results from empirical fits (solid black line), GP (dotted black line), and our semi-parametric method (dashed black line). Grey regions are 95% confidence intervals, with the boundary line types agreeing with their corresponding black lines.

We estimate the sea-wall height structure variable $v_p$ (in meters, $m$) for which the discharge rate value $Q(v_p;\bm{X})$ is expected to exceed the design standard of $0.002 m^3 s^{-1} m ^{-1}$ with probability $p$. That is, setting $V=Q^{-1}(0.002; \bm{X})$, we solve $\Pr(V>v_p) = p$ for $v_p$ using realizations of $V$ generated through simulation and reverse marginal transformation. As in <cit.>, we fix the sea-wall design feature constants to $a_1=0.25$, $a_2=26$, and tidal level relative to the seabed to $l=4.3$. We obtain three separate estimates for $v_p$ at a range of $p$ values. First, an empirical estimate of $v_p$ is obtained by using the empirical quantile function on values of the structure variable $V$ computed from the dataset. Second, a GP model is fitted using the dataset's computed values of $V$ obtained from the dataset, and quantiles are then obtained for $v_p$.
Lastly, we use our estimated joint posterior distribution for the exceedances and angles $\pi[\beta_{\G},\bm z_{\G}\mid\beta_{\QS},\bm z_{\QS}, \dat]$, to generate exceedance an observed dataset of $\bm{X}^\prime$ length $n=2,894$. We sample a posterior threshold $\rquants{i}(\bm{w})$ from the fitted distribution $\pi[\beta_{\QS},\bm z_{\QS}\mid \dat]$. Based on this threshold, sample angles $\bm{w}$ gauge function values $\gG(\bm{w})$ from the joint distribution $\pi[\beta_{\G},\bm z_{\G}\mid\beta_{\QS},\bm z_{\QS}, \dat]$ and radial exceedance observations $r-\rquants{i}(\bm{w})$ from the $\text{Exponential}(\gG(\bm{w}))$ distribution. The values $r\bm{w}$ is added to our new dataset, joined with sub-asymptotic values from the initial dataset. With this new dataset $\bm{X}^\prime$, we first convert to original margins via
for $j=H,S,P$, and we compute values of the structure variable, $v_p$, for $p$ values that correspond to extremal observations of sea-wall height, $p\in(0,0.10)$. This procedure is repeated 200 times, and we compute $p$-wise sample means and 95% confidence intervals for $v_p$. Figure <ref> shows the resulting estimates of $v_p$ for a range of $p$ values. In this figure, we follow convention and plot $-\log(-\log(1-p))$ against $v_p$, and it can be concluded that, compared to the empirical and GP distribution fit approach, our method accurately estimates the sea-wall height variable $v_p$ across all values of $p\in(0,0.10)$. The larger confidence intervals corresponding to our semi-parametric method are due to accounting for sampling variation. For each of the 200 repetitions, we generate a new dataset $\bm{X}^\prime$. Each repetition will lead to different structure variable estimates. Note that this variability is not accounted for in the empirical and GPD fits, where bootstrap confidence intervals are computed from the original dataset.

§.§ Acknowledgments

We thank Jennifer L. Wadsworth for insightful discussions on the framework of geometric extremes, and Finn Lindgren for guidance with the  and  packages in the  computing language. Part of this paper was written during a visit of IP and
LDM at King Abdullah University of Science and Technolodgy and during
a visit of RC at the University of Edinburgh.
Ryan Campbell acknowledges funding by the EPSRC DTP EP/W523811/1 fund at Lancaster University, with additional funding provided by NSERC PGS D (Canada), and the FRQNT doctoral fund (Québec).

§.§ Disclaimer

The main results of this paper, including
Theorem <ref> and the concepts presented in Section <ref> were first presented
at the 13th International Conference on Extreme Value Analysis (EVA), June 2023, Milan. After EVA, we became aware through personal communication of parallel independent work related to the notion of return level-sets, as described in Section <ref>, done by Emma S. Simpson at University College London and Jonathan A. Tawn at Lancaster University. As of now, there is no publicly available document to reference their work.

Baker, S. G.  1994, `The
multinomial-Poisson transformation', Journal of the Royal Statistical
Society: Series D (The Statistician) 43(4), 495–504.

Balkema  de Haan1974balkema1974residual
Balkema, A. A.  de Haan, L.  1974, `Residual life time at great age', The Annals of
probability 2(5), 792–804.

Balkema  Nolde2010balknold10
Balkema, G.  Nolde, N.  2010,
`Asymptotic independence for unimodal densities', Advances in Applied
Probability 42(2), 411–432.

Barnett, V.  1976, `The ordering of
multivariate data (with discussion)', J. Roy. Statist. Soc., A 139, 318–355.

Barthelmé  Chopin2015barthelme2015poisson
Barthelmé, S.  Chopin, N.  2015, `The Poisson transform for unnormalised statistical
models', Statistics and Computing 25(4), 767–780.

Bolin  Lindgren2015Bolin_excursions_2015
Bolin, D.  Lindgren, F.  2015,
`Excursion and contour uncertainty regions for latent Gaussian models',
Journal of the Royal Statistical Society, Series B (Statistical
Methodology) 77(1), 85–106.

Bolin  Lindgren2017Bolin_excursions_2017
Bolin, D.  Lindgren, F.  2017,
`Quantifying the uncertainty of contour maps', Journal of Computational
and Graphical Statistics 26(3), 513–524.

Bolin  Lindgren2018Bolin_excursions_2018
Bolin, D.  Lindgren, F.  2018,
`Calculating probabilistic excursion sets and related quantities using
excursions', Journal of Statistical Software 86(5), 1–20.

[Bortot et al.]Bortot, Coles   Tawn2000Bortotetal01
Bortot, P., Coles, S.  Tawn, J.  2000, `The multivariate Gaussian tail model: An
application to oceanographic data', JRSSC 49(1), 31–049.

Coles, S. G.  2001, An Introduction to
Statistical Modeling of Extreme Values, Springer–Verlag, London.

Coles  Tawn1991coletawn91
Coles, S. G.  Tawn, J. A.  1991,
`Modelling extreme multivariate events', J. Roy. Statist. Soc., B 53, 377–392.

Coles  Tawn1994coletawn94
Coles, S. G.  Tawn, J. A.  1994,
`Statistical methods for multivariate extremes: an application to structural
design (with discussion)', Appl. Statist. 43, 1–48.

Davison  Smith1990davismit90
Davison, A. C.  Smith, R. L.  1990, `Models for exceedances over high thresholds', J. Roy. Statist. Soc., B 52, 393–442.

de Haan1970haan70
de Haan, L.  1970, On Regular
Variation and its Application to the Weak Convergence of Sample
Extremes, Vol. 32 of Mathematical Centre Tracts, Mathematisch
Centrum, Amsterdam.

de Haan1984haan84
de Haan, L.  1984, `A spectral
representation for max-stable processes', Ann. Probab. 12, 1194–1204.

de Haan  Resnick1987de1987regular
de Haan, L.  Resnick, S.  1987, `On regular variation of probability densities', Stochastic processes
and their applications 25, 83–93.

de Haan  Resnick1977haanresn77
de Haan, L.  Resnick, S. I.  1977, `Limit theory for multivariate sample extremes', Z. Wahrsch. Theor.
40, 317–337.

de Haan  Stadtmüller1996de1996generalized
de Haan, L.  Stadtmüller, U.  1996, `Generalized regular variation of second order', Journal of the Australian Mathematical Society 61(3), 381–395.

Engelke  Hitz2020engehitz19
Engelke, S.  Hitz, A. S.  2020,
`Graphical models for extremes (with discussion)', J. Roy. Statist. Soc.
Ser. B .

[Fasiolo et al.]Fasiolo, Wood, Zaffran, Nedellec   Goude2021fasiolo2021fast
Fasiolo, M., Wood, S. N., Zaffran, M., Nedellec, R.  Goude, Y.
2021, `Fast calibrated additive quantile
regression', Journal of the American Statistical Association 116(535), 1402–1412.

Gradshteyn  Ryzhik2014gradshteyn2014table
Gradshteyn, I. S.  Ryzhik, I. M.  2014, Table of Integrals, Series, and Products,
Academic press.

[Guillen et al.]Guillen, Lindgren, Muff, Glass, Breed   Schlägel2023Lindgren_ExpLik_2023
Guillen, R. A., Lindgren, F., Muff, S., Glass, T. W., Breed, G. A.   Schlägel, U. E.  2023, `Accounting for
unobserved spatial variation in step selection analyses of animal movement
via spatial random effects', bioRxiv .

[Hansen et al.]Hansen, Herburt, Martini   Moszyńska2020hansen2020starshaped
Hansen, G., Herburt, I., Martini, H.  Moszyńska, M.
2020, `Starshaped sets', Aequationes
mathematicae 94, 1001–1092.

Heffernan  Tawn2004hefftawn04
Heffernan, J. E.  Tawn, J. A.  2004, `A conditional approach for multivariate extreme
values (with discussion)', J. Roy. Statist. Soc., B 66(3), 1–34.

[Keef et al.]Keef, Papastathopoulos   Tawn2013keefpaptawn13
Keef, C., Papastathopoulos, I.  Tawn, J. A.  2013, `Estimation of the conditional distribution of a
multivariate variable given that one of its components is large: additional
constraints for the Heffernan and Tawn model', J. Mult. Anal 115, 396–404.

[Kiriliouk et al.]Kiriliouk, Rootzén, Segers   Wadsworth2019kiriliouk2019peaks
Kiriliouk, A., Rootzén, H., Segers, J.  Wadsworth, J. L.
2019, `Peaks over thresholds modeling with
multivariate generalized Pareto distributions', Technometrics 61(1), 123–135.

Klain, D. A.  1997, `Invariant valuations on
star-shaped sets', advances in mathematics 125(1), 95–113.

Koenker, R.  2005, Quantile
regression, Vol. 38, Cambridge university press.

[Kotz et al.]Kotz, Kozubowski   Podgorski2001kotz2001laplace
Kotz, S., Kozubowski, T. J.  Podgorski, K.  2001, `The laplace distribution and generalizations', (No Title) .

Ledford  Tawn1996ledtawn96
Ledford, A. W.  Tawn, J. A.  1996, `Statistics for near independence in multivariate extreme values', Biometrika 83, 169–187.

Ledford  Tawn1997ledtawn97
Ledford, A. W.  Tawn, J. A.  1997, `Modelling dependence within joint tail regions', J. Roy. Statist.
Soc., B 59, 475–499.

[Lindgren et al.]Lindgren, Rue   Lindström2011lindetal11
Lindgren, F., Rue, H.  Lindström, J.  2011, `An explicit link between Gaussian fields and
Gaussian Markov random fields: the stochastic partial differential
equation approach (with discussion)', J. Roy. Statist. Soc., B 73, 423–498.

Lutwak, E.  1988, `Intersection bodies and
dual mixed volumes', Advances in Mathematics 71(2), 232–261.

[Majumder et al.]Majumder, Shaby, Reich   Cooley2023majumder2023semiparametric
Majumder, R., Shaby, B. A., Reich, B. J.  Cooley, D.
2023, `Semiparametric estimation of the
shape of the limiting multivariate point cloud'.

Maulik  Resnick2004maulik2004characterizations
Maulik, K.  Resnick, S.  2004,
`Characterizations and examples of hidden regular variation', Extremes
7(1), 31–67.

Nolde, N.  2014, `Geometric interpretation
of the residual dependence coefficient', Journal of Multivariate
Analysis 123, 85–95.

Nolde  Wadsworth2021nolde2021linking
Nolde, N.  Wadsworth, J.  2021,
`Linking representations for multivariate extremes via a limit set', Advances in Applied Probability .

Papastathopoulos  Paulin2023pappau23
Papastathopoulos, I.  Paulin, D.  2023, `Markov random tail fields and extended multivariate
regularly varying functions'.
Unpublished research.

Pickands, J.  1975, `Statistical inference
using extreme order statistics', Ann. Statist. 3, 119–131.

Resnick, S. I.  2007, Heavy-tail
phenomena: probabilistic and statistical modeling, Springer Science &
Business Media.

Rootzén  Tajvidi2006roo+ta:06
Rootzén, H.  Tajvidi  2006,
`The multivariate generalized Pareto distribution', Bernoulli 12, 917–930.

Seneta, E.  2006, Regularly Varying
Functions, Vol. 508, Springer, New York.

[Simpson et al.]Simpson, Illian, Lindgren, Sørbye   Rue2016simpsetal16
Simpson, D., Illian, J. B., Lindgren, F., Sørbye, S. H.  Rue, H.
2016, `Going off grid: computationally
efficient inference for log-Gaussian Cox processes', Biometrika
103(1), 49–70.

Simpson  Tawn2022simptawn22
Simpson, E. S.  Tawn, J. A.  2022, `Estimating the limiting shape of bivariate scaled sample clouds for
self-consistent inference of extremal dependence properties', arXiv
preprint arXiv:2207.02626 .

Soms, A. P.  1976, `An asymptotic expansion
for the tail area of the t-distribution', Journal of the American
Statistical Association 71(355), 728–730.

Wadsworth  Campbell2022wadscamp22
Wadsworth, J.  Campbell, R.  2022, `Statistical inference for multivariate extremes via a geometric approach',
arXiv preprint arXiv:2208.14951 .

Wadsworth  Tawn2013wadtawn13
Wadsworth, J.  Tawn, J.  2013,
`A new representation for multivariate tail probabilities', Bernoulli
19(5B), 2689–2714.

[Waldmann et al.]Waldmann, Kneib, Yue, Lang   Flexeder2013waldmann2013bayesian
Waldmann, E., Kneib, T., Yue, Y. R., Lang, S.  Flexeder, C.
2013, `Bayesian semiparametric additive
quantile regression', Statistical Modelling 13(3), 223–252.

Wood, S. N.  2017, Generalized
Additive Models: an Introduction with R, CRC press.

Yu  Moyeed2001yumoy01
Yu, K.  Moyeed, R. A.  2001,
`Bayesian quantile regression', Stat. Probabil. Lett. 54(4), 437–447.

Zellner, A.  1988, `Optimal information
processing and Bayes's theorem', The American Statistician 42(4), 278–280.

§ STARSHAPED SETS

§.§ Background

An expository review and survey of material about the fundamental
mathematical notion of starshaped sets, emphasising their
geometric, analytical, combinatorial, and topological properties, is
given by <cit.> and references therein.

If $a$ and $b$ are different elements of $\RR^d$, then
$[a:b]=\{(1-t)a + t b\,:\, 0 \leq t\leq 1\}$ denotes the segment with
endpoints $a$ and $b$, $[a:b)=\{(1-t)a + t b\,:\, t \geq 0\}$ denotes
the half-line with origin $a$ through $b$,
$(a:b)=\{(1-t)a + t b\,:\, 0 \leq t\in \RR\}$ denotes the line through
$a$ and $b$, $]a:b]=\{(1-t)a + t b\,:\, 0 < t\leq 1\}$ denotes the
segment $[a:b]$ without its endpoint $a$,
$[a:b[=\{(1-t)a + t b\,:\, 0 \leq t < 1\}$ denotes the segment
$[a: b]$ without its endpoint $b$ and
$]a:b)=\{(1-t)a + t b\,:\, t > 0\}$ denotes the ray $[a, b)$ without
its origin $a$. A set $G\subset \RR^d$ is called starshaped if there
exists $x\in G$ such that every point $y\in G$ is visible from $x$ in
the sense that $[x:y]\subset G$. If $G$ is a set and $x\in G$, the
star of $x$ in $G$ is the set of all points in $G$ that are visible
from $x$, that is, $\text{st}(x:G)=\{y\in G\,:\, [x:y]\subset G\}$. A
set $G$ is starshaped if and only if there exists $x\in G$ such that
$\text{st}(x:G)=G$. The kernel of a set $G$, denoted by
$\text{ker}\,G$, is the convex set of all of its star centres, that
is, $\text{ker}\,G= \{x\in G\,:\, \text{st}(x:G)=G\}$. A star set is
called a star-body if it is a compact regular domain, that is if $G$
is a compact starshaped set such that $G^\circ$ is connected and
$G=\overline{G^{\circ}}$. Let $G$ be a closed proper subset of
$\RR^d$ with non-empty interior. Then $G$ is said to be strongly
starshaped at $a$ whenever there exists $a\in\RR^d$ such that for
every $\bm w\in\SSS^{d-1}$, the halfline $a+[0: \bm w)$ does not
intersect the boundary $\partial G$ of $G$ more than once. A strongly
starshaped set at $a$ is a starshaped set with
$a \in \text{ker}\, G$. A set $G$ is said to be strongly starshaped
if it is strongly starshaped at $0$. We write $\Delta_{\bm w}$ for
the half-line through the origin $[0:\bm w)$ passing through
$\bm w\in\SSS^{d-1}$. If $G$ is starshaped at $m\in \text{ker}\,G$,
we say that $z\in G$ is the last point of a ray $[m:x)$ in $G$ if
$z\in[m:x)$ and there is no point $y\in[m:x)\cap G$ such that
$y\notin [m:z]$. Let $G$ be a starshaped set such that
$(\text{ker}\, G)^\circ \not= \varnothing$. If
$m\in (\text{ker}\, G)^\circ$ and $x\in \overline{G}$, then
$[m:x[\subset G$. If $G$ is a closed starshaped set such that
$(\text{ker}\, G)^\circ\not= \varnothing$, then $G^\circ$ is connected
and $G=\overline{(G^\circ)}$, whence a regular domain.

§.§ Radial addition and gauges of starshaped sets

If $G$ is a closed starshaped set, $m\in \text{ker}\,G$ and
$\bm w\in\SSS^{d-1}$, then there are two cases. There exists a last
point $p$ of $m+\Delta_{\bm w}$ in $G$ or
$m + \Delta_{\bm w} \subset G$. In the first case, let

\[
\rho_{m,G}(\bm w) = \sup\{\lambda \in \RR\,:\, m + \lambda \bm w \in G\}
\]

The radial function of $G$ at $m$ is the function
$r_{m,G}\,:\, \SSS^{d-1}\to \RR_+$ defined by

\[
r_{m,G} (\bm w)=
\begin{cases}
\rho_{m,G}(\bm w)& m+ \Delta_{\bm w} \not\subseteq G\\
+\infty& \text{$m+\Delta_{\bm w}\subseteq G$}
\end{cases}
\]

We have the equivalence
$m \in (\text{ker}\,G)^\circ \Leftrightarrow r_{m, G} > 0$. Also, if
$G$ is compact and $m \in (\text{ker}\,G)^\circ$, then $r_{m,G}$ is a
Lipschitz function.

Let $\bigstar$ denote the space of star-bodies in $\RR^d$ and
$\radd,\rmult\,:\,\bigstar \times \bigstar \to \bigstar$ be the
operators of radial addition and radial
multiplication respectively, defined by means of radial functions
as as
ρ_G_1 G_2 (w)    =    ρ_G_1(w) + ρ_G_2(w),     G_1, G_2 ∈

ρ_G_1 G_2 (w)    =    ρ_G_1(w) ρ_G_2(w),
G_1, G_2 ∈.

Multiplication of a set $G \in \bigstar$ by a nonnegative scalar $c$
is equivalent to radial multiplication of $G$ by $B_c(0)$, that is,
$c G:=B_c(0)\rmult G$, $c \geq 0$. The set $\bigstar$ equiped with the binary
operation $\radd$ is a commutative monoid. In other words,
$(\bigstar, \radd)$ satisfies the following.

Identity element: There exists an element $0\in\bigstar$ such
that $0 \radd G = G \radd 0 = G$ $\forall G \in\bigstar$.

Commutative law: For $G_1$ and $G_2$ in $\bigstar$,
$G_1 \radd G_2 = G_2 \radd G_1$.

Additionally, we have

Scalar Multiplication Closure: The result of $\rmult$
between a compact star set and a non-negative scalar remains within
the space of compact star sets.

Distributive Laws: The operations $\radd$ and $\rmult$ obey
distributive laws similar to those in a ring. For example,
$(G_1 + G_2) \rmult B_c(0) = G_1 \rmult B_c(0) + G_2 \rmult B_c(0)$.

Associativity: Both $\radd$ and
$\rmult$ are associative operations. For example, $(G_1 + G_2)
+ G_3 = G_1 + (G_2 + G_3)$.

The monoid $(\bigstar, \radd)$ is endowed with its algebraic
preordering $\leq$ defined by $G_1 \leq G_2$ whenever there exists
$H\in\bigstar$ such that $G_1 \radd H = G_2$. In this case, we define
radial subtraction $H = G_2 - G_1$ through ordinary difference of
radial functions

The algebraic structure on the space of compact star sets with the
operations $\radd$ and $\rmult$ provides a mathematical
framework to explore relationships, transformations, and compositions
of these sets. This algebraic approach facilitates the development of
consistent theories and methodologies when dealing with operations
involving compact star sets.

Last, let $G_1$ and $G_2$ be two star-bodies such that $\rho_{S_i}$ is
continuous on $\SSS^{d-1}$. Then the dual Brunn-Minkowski inequality
<cit.> yields

\[
\lvert G_1 \radd G_2\rvert^{1/d} \leq \lvert G_1\rvert^{1/d} + \lvert G_2\rvert^{1/d},
\]

and equality holds if $d=1$ or $d\geq 2$ and $G_2 = r G_1$, for some
$r>0$. If $G$ is compact and strongly starshaped, then the volume of
$G$ can be expressed by the formula

\begin{equation}
\vol{G}=\frac{1}{d}\int\limits_{\SSS^{d-1}} \rG(\bm w)^d d\bm w,
\label{eq:vol}
\end{equation}

see <cit.>. The reciprocal of the radial function
of a set $G$ starshaped at $0$ is called the gauge function
$g_G\,:\,\Rstar \to \RR_+$, defined by

\[
\gG(\bm x) = \inf\{t\geq 0\,:\,\bm x\in t G\}, \quad \text{for
$\bm x\in\Rstar$}.
\]

§ EXAMPLES

A random vector $\bm X\in\RR^d$ is said to follow a multivariate normal
distribution with precision matrix $\mathsf{Q}=\Sigma^{-1}$
for a positive definite covariance $\mathsf{\Sigma}$ if its density $f_{\bm X}$ satisfies

\[
f_{\bm X}(\bm x) = \frac{\lvert\mathsf{Q}\rvert^{1/2}}{(2\pi)^{d/2}}\exp\left(-\bm x^\top \mathsf{Q} \bm x^\top/2\right).
\]

A random vector $\bm X\in \RR^d$ is said to follow a multivariate standard normal distribution with Laplace margins if its joint density $\fL$ is given by

\[
\fL(\bm x) = \frac{\lvert \QQ \rvert^{1/2}\exp\{-\lVert \lvert\bm x\rvert^{1/2} \rVert_1\}}{2^{3}\pi}\exp\left(-\frac{1}{2} [\Phi^{-1}\{F( \bm x)\}]^\top \QQ \Phi^{-1}\{F(\bm x)\}\right)\prod_{i=1}^2\frac{f(x_i)}{\varphi\left[ \Phi^{-1}\{F(x_i)\}\right]},
\]

where $f$ and $F$ respectively denote the density and
distribution functions of a univariate Laplace distribution with function evaluation $F(\bm{x})$ and $f(\bm{x})$ performed componentwise. Furthermore,
$\QQ=\Sigma^{-1}$ with $\Sigma$ as before. In what follows, we use the
bivariate normal distribution as an instance of asymptotically
independent distribution and draw samples using the value $\rho=0.8$.A random vector $\bm X\in \RR^2$ has a bivariate max-stable logistic
distribution with Laplace margins if its density $\fL$ is given
(x)   =    [∏_j=1^df(x_j)/F(x_j)[logF(x_j)]^2 ](∑_π∈Π∏_s∈π

where the exponent function is given by
x_j^{-{1}/{\theta}}\right)^\theta$ and the dependence parameter
$\theta\in(0,1)$ fully determines a member of the family.

§ REMAINDER TERM AND QUALITY OF CONVERGENCE

§.§ Rates of convergence

\begin{equation}
\fL(t \bm x)=\exp\{-g_{\G_L}(t\bm x) + \text{remainder}\}(1+\oo(1)), \quad \text{as $t\to \infty$},
\label{eq:rof}
\end{equation}

where the remainder is interpreted as the evaluation $\roc(r \bm w)$
of a function $\roc\,:\, \mathbb{K}\to \RR$ at
$r\bm w \in \mathbb{K}$. The domain of the function $\roc$ is assumed
to be a solid cone $\mathbb{K}$ in $\Rstar$. A characterisation of
the map $\roc$ under second-order conditions <cit.>
is given in Section <ref>. A functional
$\roc\,:\,\Rstar\to\RR$ is said to determine the rate of convergence
to $\gG$ in expression (<ref>) if it satisfies

\begin{equation}
\fL(t\bm x) = \exp\left(-\left[t \gG(\bm x)
\{1+\oo(1)\}\right]\right)=\exp\left[-\left\{\gG(t\bm x) +
\roc(t\bm x)\right\}\right]\{1+\oo(1)\},
\quad\bm x \in \Rstar,
% \{-\log f(t\bm x)/t\} - \gG(\bm x) = \roc(t\bm x)\{1+\oo(1)\}, \quad
% \text{as $r\to \infty$}
\label{eq:rof_defn}
\end{equation}

as $t\to \infty$. Two functionals $\roc_1$ and $\roc_2$ that
determine the rate of convergence to $\gG$ are asymptotically
equivalent, that is,
$\lim_{t\to \infty} \roc_1(t \bm x)/\roc_2(t \bm x) = 1$, $\forall \bm
x\in\Rstar$. As such, the set of all functionals
$\roc$ satisfying expression (<ref>) forms an equivalence
class. To better understand the properties of the members of this
class, it is helpful to consider the special case where the leading
order term of
$\roc$ is equal to the product of a positive scaling function
$\lfirst(t)$ and a function $\qfirst(\bm x)$, that is, when $\roc(t
\bm x) = \OO(\lfirst(t)\qfirst(\bm x))$ as $t\to \infty$, for $\bm
x\in\Rstar$. Assumption <ref> delineates a sufficient
condition subject to which the rate of convergence has such a

Suppose that condition $(ii)$ of Proposition <ref> holds.  There exists a non-decreasing function
$\lfirst\,:\,(0,\infty) \to (0,\infty)$ and a function
$\qfirst\,:\,\Rstar\to \RR$ which is not everywhere zero, such that

\begin{equation}
\lim_{r\to\infty}\frac{\{-\log \fL(t \bm x)/t\} - \gG(\bm x)}{\afirst(t)} =
\qfirst(\bm x),
\label{eq:RV2}
\end{equation}

where $\afirst(t)=\lfirst(t)/t$ for $t>0$.

Due to condition $(ii)$ of Proposition <ref> we have that the
scale sequence $\afirst(t) = \oo(1)$ as $t\to\infty$. This implies
that $\lfirst(t) = \oo(t)$ as $t\to \infty$. Also
Assumption (<ref>) requires that $-\log \fL$ is a
multivariate extended regularly varying function in the sense
of Definition 1 of <cit.>. To be precise,
convergence (<ref>) is equivalent to assuming that there exists
$\pfirst\,:\,\Rstar\to \RR$ which is not everywhere zero, such that

\begin{equation}
\lim_{t\to\infty}\frac{[-\log \fL(t\bm x)/\gG(t\bm x)] - [-\log \fL(t\bm
1)/\gG(t\bm 1)]}{\afirst(t)} = \pfirst(\bm x), \qquad \bm x\in\Rstar.
\label{eq:ERVd}
\end{equation}

The correspondence between the two limit functions in (<ref>)
and (<ref>) is
$\pfirst(\bm x)=[\qfirst(\bm x)/\gG(\bm x)]-[\qfirst(\bm 1)/\gG(\bm
1)]$. From this, it immediately follows that $p_1(\bm 1)=0$.
Also, it is known that the only possible non-trivial limit functions
that can arise from convergence (<ref>) are those that satisfy
the functional equation

\[
[\pfirst(t \bm x)-\pfirst(t \bm 1)]/t^{\gamma - 1} =\pfirst(\bm x) ,
\qquad t > 0,\qquad \bm x\in\Rstar,
\]

for some $\gamma \in \RR$, see Theorem 2 of <cit.> for a
proof. These observations lead to Proposition <ref> which
identifies the key properties of the scaling function
$\afirst(t)=\lfirst(t)/t$ and of the limit function $\qfirst(\bm x)$
under Assumption (<ref>).

Suppose that Assumption <ref> holds. There exists
$\gamma < 1$ such that $\afirst(t)\in\textsf{RV}_{\gamma-1}^\infty$
and $\qfirst$ is $\gamma$–homogeneous, that is,
$\qfirst(t \bm x)=t^\gamma \qfirst(\bm x)$ $\forall t>0$ and
$\bm x\in\Rstar$.

From Proposition <ref> we obtain
$\lfirst(t)\in \text{RV}_{\gamma}^\infty$. Because Assumption
<ref> requires $\lfirst(t)$ to be a non-decreasing function,
it suffices to consider the following two cases. First, the case when
$\lfirst(t)$ converges to a constant $K\in (0,\infty)$ and second, the
case when $\lfirst(t)$ grows without bound. For the former case we
consider without loss of generality the case $\lfirst(t)=1$
$\forall t>0$. This is because positive real constants can be
absorbed in the limit function $\qfirst$. Below, we treat the two
cases in turn.

First, when $\lfirst(t)=1$, the leading order term of $\roc(t\bm x)$
is constant in $t$ which implies that the rate of convergence is
determined by the function $\qfirst$, that is,
$\roc(t \bm x) = \qfirst(t \bm x) + \oo(1)$, as $t\to \infty$. Also,
because $\lfirst\in \textsf{RV}_{0}^\infty$, Proposition <ref>
yields that $\qfirst$ is a $0$–homogeneous function. Consequently,

\begin{equation}
\fL(t \bm x) = \exp\left[-\left\{t \gG(\bm x) + \qfirst(\bm x)\right\}\right] \{1+\oo(1)\}, \qquad
\text{as $t\to \infty$}.
\label{eq:approx_g0}
\end{equation}

This special case gives a characterisation of the rate of convergence
which is independent of the distance from the origin $t$. This
behaviour is specific to the case $l(t)=1$ and is waived when
$\lfirst(t)$ grows without bound because in the latter case, the rate
of convergence may have non-negligible remainder terms.
Using the same working principle as in Assumption <ref>, below
in Assumption <ref> we outline a sufficient condition subject
to which the rate of convergence exhibits a so-called
second-order behaviour <cit.>.

Suppose that Assumption <ref> holds with
$\lim_{t\to\infty} \lfirst(t) = \infty$. There exists a
non-decreasing function $\lsecond$ and a function
$\qsecond\,:\,\Rstar\to \RR$ which is not equal to the product of
$\qfirst$ and a $0$–homogeneous function on
$\RR^{d}\setminus\{0\}$, such that

\begin{equation}
\lim_{t\to\infty} \frac{([\{-\log f_{L}(t \bm x)/t\} - \gG(\bm x)]/\afirst(t)) - \qfirst(\bm x)}{\asecond(t)} =
\qsecond(\bm x),
\label{eq:RV3}
\end{equation}

where $\asecond(t)=\lsecond(t)/\lfirst(t)$.

It is clear that when Assumption <ref> holds with
$\lsecond(t)=1$, then

\begin{equation*}
\fL(t \bm x) = \exp[-\{t \gG(\bm x) + \lfirst(t)\,\qfirst(\bm x) + \qsecond(\bm x)\} \{1+\oo(1)], \qquad \text{as $t\to \infty$}.
\end{equation*}

More generally, we can construct an asymptotic expansion of
$-\log \fL(t\bm x)$ of the form

\begin{equation}
-\log \fL(t\bm x) = \gG(t \bm x) + \sum_{j=1}^{\infty}l_j(t)
u_j(\bm x), \quad\text{as $t\to\infty$}.
\label{eq:expansion}
\end{equation}

for a sequence of positive scale functions $\{l_j\}_{j=1}^M$ with
$M\in\mathbb{N}\cup\{\infty\}$, satisfying $l_j(t)=\oo(l_{j+1}(t))$ as
$t\to \infty$, and a sequence of functions $\{u_j\}_{j=1}^M$
mapping $\Rstar$ into $\RR$. The function $u_j$ captures the
contribution of $-\log\fL(t\bm x)$ at a given order of the series and
$l_{j}$ is the scaling factor associated with this term. The $j$-th
term $l_j(t)u_j(\bm x)$ in the sum represents the $j$-th
approximation to $-\log \fL(t\bm x)$ as $t\to \infty$, and the scaling
factors $l_j$ determine the rate of convergence to $u_{j-1}$,
with $u_0(\bm x):=\gG(\bm x)$.

In Sections
<ref>–<ref> of the
Supplementary Material, we derive the gauge $\gG$ of the limit set
$\G$ and, the limit functions $\qfirst$ and $\qsecond$ of the rate of
convergence $\roc$, for three examples of $d$-dimensional copulas in
standard Laplace margins, the multivariate normal copula (MVN), the
multivariate Laplace (MVL) and the max-stable logistic. The forms of
the functions are summarised in Table <ref>. We remark
that in Assumptions <ref> and <ref>, the limit
functions $\qfirst$ and $\qsecond$ may be discontinuous. For example,
Figure <ref> shows the gauge function $\gG$ and the limit
functions $\qfirst$ and $\qsecond$, plotted on a subset of
$\RR^2\setminus \{0\}$, for the bivariate normal copula in Laplace
margins (Supplementary Material <ref>). Here,
the limit functions $\qfirst$ and $\qsecond$ exhibit discontinuities
of the second kind at the axes $\{(x_1, 0)\,:\, x_1 > 0\}$ and
$\{(0, x_2)\,:\, x_2 > 0\}$.

Limit functions in asymptotic expansion of the logarithm of
the density of the multivariate normal copula with Laplace
marginals. Left: gauge function $\gG$ of limit set
$\G$. Centre: limit function $\qfirst$. Right:
Limit function $\qsecond$. Bottom: alternative viewing
angle obtained by rotating the images shown in top row by
$60^\circ$ about $x$ axis

A natural question arises: can we determine the properties of the
limit functions $\lambda_1,\dots,\lambda_M$ in some detail? (can we
use induction on $M$?)

Suppose that Assumption <ref> holds for the
random vector $\bm X_L$. There exists $\gamma < 1$ and $\rho\leq 0$
such that

\[
\qsecond (t \bm x) = t^\gamma \qsecond(\bm x) + c
t^\gamma\frac{t^\rho - 1}{ \rho} \qfirst(\bm x).
\]

and the solution of this functional equation is

\begin{equation}
\qsecond (\bm x) = \lVert \bm x \rVert^{\gamma+\rho}
h\left(\frac{\bm x}{\lVert \bm x\rVert}\right) + c \lVert \bm x
\rVert^{\gamma} \log \lVert \bm x
\rVert \qfirst\left(\frac{\bm x} {\lVert \bm x
\rVert}\right), \qquad \bm x\in\RR^{d}\setminus\{0\}.
\label{eq:g0star_sol}
\end{equation}

for a measurable function $h$ that maps $\SSS^{d-1}$ into $\RR$.

Copula     $\gG(\bm{x})$

5*Gaussian, $\QQ=\Sigma^{-1}$
5* $\left\{\text{sgn}(\bm x) \vert \bm
x\rvert^{1/2}\right\}^\top \QQ \, \left\{\text{sgn}(\bm x)
\vert
\bm x\rvert^{1/2}\right\}$

4*Laplace, $\QQ=\Sigma^{-1}$
4* $\left(\bm{x}^\top \QQ\bm{x}\right)^{{1}/{2}}$

11*$\begin{cases}\frac{1}{\theta}\sum\limits_{j=1}^d x_j + \left(1-\frac{d}{\theta}\right)\underset{k=1,\dots,d}{\min}x_k & ; x_j>0\forall j \\ \frac{1}{\theta}\sum\limits_{j:x_j>0}x_j + \left(\sum\limits_{k:x_k<0} \left(-x_k\right)^{{1}/{\theta}}\right)^\theta&;\text{o.w.}\end{cases}$

logistic dependence,

3*Student-$t$, $\nu>0$     3* $-\frac{1}{\nu}\sum\limits_{j=1}^d \left|x_j\right| + \left(1+\frac{d}{\nu}\right)\max_{j=1,\dots,d}\left|x_j\right|$

Gauge functions for a suite of $d$-dimensional copulas.

Copula      $\qfirst(\bm{x})$     $\qsecond(\bm{x})$

5*Gaussian, $\QQ=\Sigma^{-1}$
5*$\frac{d}{2} - \frac{1}{2} \left\{\text{sgn}(\bm x)
\lvert \bm x\rvert^{1/2}\right\}^\top \, \QQ \,
\left\{\text{sgn}(\bm x) \lvert \bm x\rvert^{-1/2}\right\}$

$\left\{\text{sgn}(\bm x) \vert \bm
x\rvert^{1/2}\right\}^\top \QQ \, \left\{\text{sgn}(\bm x)
\log2\lvert \bm x \rvert^{-1/2}\right\}$

$\hspace{0.3cm}- \frac{1}{2}
\left\{\text{sgn}(\bm x) \lvert \bm x\rvert^{1/2}\right\}^\top \,
\QQ \, \left\{\text{sgn}(\bm x) \lvert \bm x\rvert^{-1/2}
\log (4\pi \lvert \bm x \rvert)\right\}$

$\hspace{0.3cm}+ \frac{1}{2}\log \lvert \QQ \rvert + \sum_{i=1}^d \frac{1}{2} \log 4\pi \lvert x_i\rvert$

4*Laplace, $\QQ=\Sigma^{-1}$

$-d\log 2 + \frac{d}{2}\log 2\pi -\frac{1}{2}\log\left|\QQ\right| $

$\hspace{0.3cm}- \frac{1}{2}\log\left(\frac{\pi}{2}\right) + \frac{d}{4}\log\left(\bm{x}^\top \QQ\bm{x}\right)$

-\log\left(2^{-1} (-1)^{d+1} \left(\prod\limits_{\ell=1}^{d-1} \left(1-\frac{\ell}{\theta}\right)\right)\right)&; x_j>0\forall j\\
d-1&; x_j<0\forall j \\
-\left(\frac{1}{\theta}-1\right)|B| -1 +\frac{d}{\theta}&; \text{o.w.}
\end{cases}$
\begin{cases}
0&; x_j>0\forall j\\
-\log\left[(-1)^{d+1} \left(\prod\limits_{\ell=0}^{d-1} \left(1-\frac{\ell}{\theta}\right)\right)\right.\ \\
\left.\hspace{0.8cm}\times\left(\prod\limits_{j=1}^d \left(-x_j\right)\right)^{\frac{1}{\theta}-1}\left(\sum\limits_{j=1}^d \left(-x_j\right)^{{1}/{\theta}}\right)^{\theta-d}\right]&; x_j<0\forall j
\\
-\log\left[(-1)^{d+1} \left(\prod\limits_{\ell=0}^{d-1} \left(1-\frac{\ell}{\theta}\right)\right) 2^{-\frac{|A|}{\theta}}(\log 2)^{-2|C|}\right] \\ \hspace{0.5cm}-\left(\frac{1}{\theta}-1\right)\sum\limits_{k\in B}\log(-x_k) \\ \hspace{0.5cm}- (\theta-d)\log\left(\sum\limits_{k\in B} \left(-x_k\right)^{{1}/{\theta}}\right)&;\text{o.w.}
\end{cases}$

logistic dependence,

3*Student-$t$, $\nu>0$     3*$-\log\left( \frac{c_\nu^{d+2} k_{\nu,\bm{Q}}q_{j^\star j^\star}}{\nu^{d+1}}\right)$ (see Section <ref> of the Supplementary Material)     3*0

Higher order terms for the negative logarithm of the joint density in Laplace margins for a suite of $d$-dimensional copulas.