# Regression-type analysis for block maxima on block maxima

Miguel de Carvalho  Gonçalo dos Reis  Alina Kumukova  School of Mathematics,
University of Edinburgh, The King’s Buildings, Edinburgh, EH9 3FD, UK  Centro
de Matemática e Aplicações (CMA), FCT, UNL, PT  Maxwell Institute for
Mathematical Sciences School of Mathematics, University of Edinburgh,
Edinburgh UK

###### Abstract

This paper devises a regression-type model for the situation where both the
response and covariates are extreme. The proposed approach is designed for the
setting where both the response and covariates are themselves block maxima,
and thus contrarily to standard regression methods it takes into account the
key fact that the limiting distribution of suitably standardized componentwise
maxima is an extreme value copula. An important target in the proposed
framework is the regression manifold, which consists of a family of regression
lines obeying the latter asymptotic result. To learn about the proposed model
from data, we employ a Bernstein polynomial prior on the space of angular
densities which leads to an induced prior on the space of regression
manifolds. Numerical studies suggest a good performance of the proposed
methods, and a finance real-data illustration reveals interesting aspects on
the comovements of extreme losses between two leading stock markets.

Key words and phrases: Bernstein polynomials, Block maxima, Extreme value
copula, Joint extremes, Multivariate extreme value distribution, Quantile
regression, Statistics of extremes.

##  1 Introduction

Block maxima data—such as annual maxima—are a mainstay of statistics of
extremes. Whereas classical statistical modeling is mostly concerned with
inferences surrounding the bulk of a distribution, the field of statistics of
extremes deals with the rather challenging situation of conducting inferences
about the tail of a distribution. The behavior of extreme values in large
samples is often mathematically tractable, and this tractability is often used
to build sound statistical methods for modeling risk and extreme values. As an
example of this asymptotic tractability, it is well known that if
$Y_{1},\dots,Y_{n}$ is a random sample with sample maximum
$M_{n}=\max(Y_{1},\dots,Y_{n})$ and if there exist sequences $\\{a_{n}>0\\}$
and $\\{b_{n}\\}$ such that
$(M_{n}-b_{n})/a_{n}\overset{\mathrm{d}}{\rightarrow}Z$, then $Z$ follows a
GEV (Generalized Extreme Value) distribution with location, scale, and shape
parameters $\mu\in\mathbb{R}$, $\sigma>0$, and $\xi\in\mathbb{R}$
respectively; see, for instance, (Coles, 2001, Theorem 3.1.1). Details on the
paradigm of statistics of extremes can be found, for instance, in the books of
Coles (2001), Beirlant et al. (2004), and de Haan and Ferreira (2006); for
applications in finance, insurance, and economics see, for instance, Embrechts
et al. (1997), (McNeil et al., 2015, Ch. 7), and Longin (2017).

In this paper, we devise a regression-type method for the situation where both
the response and the covariates are themselves block maxima. An important
result in the field of statistics of extremes—that will be fundamental for our
developments—is that the properly standardized vector of block maxima
converges in distribution to a so-called extreme value copula (Gudendorf and
Segers, 2010). Thus, a key target in the proposed framework is what we will
refer below as the regression manifold, that is, a family of regressions lines
that obeys the latter large sample result. Our methods thus take on board
information on the dependence structure between the extreme values so to
assess what effects block maxima covariates can have on a block maxima
response. To learn about the proposed model from data, we develop a prior in
the space of regression manifolds by resorting to a flexible Bernstein
polynomial prior on the space of angular densities as recently proposed by
Hanson et al. (2017).

Our approach contributes to the literature on conditional modeling given large
observed values (e.g. Wang and Stoev, 2011; Cooley et al., 2012), nonetheless,
our focus differs from the latter papers in a number of important ways as we
describe next. A main difference is that, as anticipated above, here the focus
is on devising a regression framework for a block maxima response on block
maxima covariate, whereas the latter papers focus mainly on using the
conditional density as a way to make probabilistic statements about the
likelihood of an extreme given the occurrence of another extreme. Since our
main target of analysis is regression, our method has some links with
statistical approaches for nonstationary extremes (e.g. Coles, 2001, Section
6); the most elementary version of approaches for nonstationary extremes aims
to learn about how the limiting law of a suitably standardized block maxima
response ($S_{\mathbf{x}}$) changes according to a covariate
$\mathbf{x}=(x_{1},\dots,x_{p})^{\mathrm{\scriptscriptstyle T}}$, via the
specification

$\displaystyle(S\mid\mathbf{X}=\mathbf{x})\sim\text{GEV}(\mu_{\mathbf{x}},\sigma_{\mathbf{x}},\xi_{\mathbf{x}}).$
(1)

Since the approach in (1) is built from the univariate theory of extremes it
is not tailored for conditioning on another variable being extreme as it fails
to take on board information from the dependence structure between the
extremes.

Additionally, the method proposed in this work has some relation to quantile
regression (Koenker, 2005), whose most elementary version consists of modeling
the conditional quantile of a response $Y$ given a covariate
$\mathbf{X}=(X_{1},\dots,X_{p})^{\mathrm{\scriptscriptstyle T}}$ in a linear
fashion, that is

$\displaystyle F^{-1}(q\mid\mathbf{x})=\mathbf{x}^{\mathrm{\scriptscriptstyle
T}}\bm{\beta}_{q},\qquad 0<q<1,$ (2)

where $F^{-1}(q\mid\mathbf{x})=\inf\\{y:F(y\mid\mathbf{x})\geq q\\}$ and
$F(y\mid\mathbf{x})$ is the distribution function of $Y\mid\mathbf{X}$.
Versions of quantile regression that aim to equip (2) with the ability to
extrapolate into the tail of $Y$ are often known as extremal quantile
regression methods (e.g. Chernozhukov, 2005). While flexible and sturdy, such
quantile regression-based approaches do not take into account information on
the fact that the limiting joint distribution of suitably standardized
componentwise maxima is an extreme value copula, and thus fail to be equipped
with the ability to extrapolate into the joint tail.

The remainder of the paper unfolds as follows. In Section 2 we introduce the
proposed model and devise an approach for learning about it from data. Section
3 reports the main findings of a Monte Carlo simulation study. We showcase the
proposed methodology in a real data application to stock market data in
Section 4. Finally, in Section 5 we present closings remarks. Proofs and
derivations can be found in the appendix, and further numerical experiments
and other technical details are presented in the supplementary material (E and
F).

##  2 Modelling limiting block maxima conditioned on block maxima

###  2.1 Conditional modeling based on multivariate extreme value
distributions

### Background on multivariate extremes

Prior to introducing a regression of block maxima on block maxima we need to
lay groundwork on multivariate extremes. Let
$\\{(S_{i},\mathbf{Z}_{i})\\}_{i=1}^{n}$ be a sequence of independent random
vectors with unit Fréchet marginal distributions, i.e. $\exp(-1/z)$, for
$z>0$. In our setup, $S_{i}$ should be understood as a response, whereas
$\mathbf{Z}_{i}=(Z_{1,i},\dots,Z_{p,i})$ should be understood as a
$p$-dimensional covariate. Let the componentwise block maxima be
$\mathbf{M}_{n}=(M_{n,s},M_{n,z_{1}},\dots,M_{n,z_{p}})$ with
$M_{n,s}=\max\\{S_{1},\dots,S_{n}\\}$ and
$M_{n,z_{j}}=\max(Z_{j,1},\dots,Z_{j,n})$, for $j=1,\dots,p$. Under this
setup, it is well-known that the vector of normalized componentwise maxima
$\mathbf{M}_{n}/n$ converges in distribution to a random vector
$(Y,\mathbf{X})$ which follows a multivariate extreme value distribution with
the joint distribution function

$\displaystyle G(y,\mathbf{x})=\exp\\{-V(y,\mathbf{x})\\},\quad
y,x_{1},\dots,x_{p}>0.$ (3)

Here,

$V(y,\mathbf{x})=d\int_{\Delta_{d}}\max\left(\dfrac{w_{1}}{y},\dfrac{w_{2}}{x_{1}},\dots,\dfrac{w_{d}}{x_{p}}\right)\,H(\mathrm{d}\mathbf{w}),$

is the so-called exponent measure; see, for instance, de Haan and Resnick
(1977), Pickands (1981), and (Coles, 2001, Theorem 8.1). In addition, $H$ is a
parameter of the multivariate extreme value distribution $G$ known as angular
measure, which controls the dependence between the extreme values;
specifically, $H$ is a probability measure on the unit simplex
$\Delta_{d}=\\{(w_{1},\dots,w_{d})\in[0,1]^{d},\sum_{i=1}^{d}w_{i}=1\\}\subset\mathbb{R}^{d}$,
with $d=p+1$, and obeying the mean constraint

$\displaystyle\int_{\Delta_{d}}\mathbf{w}\,H(\mathrm{d}\mathbf{w})=\dfrac{1}{d}1{1}_{d},$
(4)

where $1{1}_{d}$ is a vector of ones in $\mathbb{R}^{d}$. If $H$ is absolutely
continuous with respect to the Lebesgue measure then its density is given by
the Radon–Nikodym derivative $h=\mathrm{d}H/\mathrm{d}\mathbf{w}$, for
$\mathbf{w}\in\Delta_{d}$.

### Regression manifold for block maxima on block maxima

We are now ready to introduce our regression method for block maxima on block
maxima. We define the regression manifold as the family of regression lines,

$\mathscr{L}=\\{L_{q}:0<q<1\\}\quad\text{with}\quad
L_{q}=\\{y_{q\mid\mathbf{x}}:\mathbf{x}\in(0,\infty)^{p}\\},$ (5)

where

$\displaystyle
y_{q\mid\mathbf{x}}=\inf\left\\{y>0:G_{Y\mid\mathbf{X}}(y\mid\mathbf{x})\geq
q\right\\},$ (6)

is a conditional quantile of a multivariate extreme value distribution, with
$q\in(0,1)$ and $\mathbf{x}\in(0,\infty)^{p}$, and
$G_{Y|\mathbf{X}}(y\mid\mathbf{x})=\mathbb{P}\left\\{Y\leq
y\mid\mathbf{X}=\mathbf{x}\right\\}$ is a conditional multivariate extreme
value distribution function.

In higher dimensions $G_{Y\mid\mathbf{X}}$ can be expressed with the help of a
joint multivariate extreme value density $g_{Y,\mathbf{X}}$ and its expression
has been derived by Stephenson and Tawn (2005). By applying Bayes’ theorem, we
deduce
$G_{Y\mid\mathbf{X}}(y\mid\mathbf{x})=\int_{0}^{y}g_{Y\mid\mathbf{X}}(z\mid\mathbf{x})\,\mathrm{d}z$
from $g_{Y,\mathbf{X}}$ with $g_{Y|\mathbf{X}}$ given as follows:

$\displaystyle
g_{Y|\mathbf{X}}(y\mid\mathbf{x})=\dfrac{\exp\\{-V(y,\mathbf{x})\\}\sum\limits_{i=1}^{d}\sum\limits_{j=1}^{n_{i,d}}(-1)^{i}\prod\limits_{\Lambda\in
r_{ij,d}}V_{\Lambda}(y,\mathbf{x})}{\sum\limits_{i=1}^{d}\sum\limits_{j=1}^{n_{i,d}}(-1)^{i}\int\limits_{0}^{\infty}\exp\\{-V(y,\mathbf{x})\\}\prod\limits_{\Lambda\in
r_{ij,d}}V_{\Lambda}(y,\mathbf{x})\,\mathrm{d}y},\quad y,\mathbf{x}>0,$ (7)

where $V_{\Lambda}(y,\mathbf{x})$ corresponds to mixed partial derivative of
the exponent measure $V(y,\mathbf{x})$ with respect to the $i$th components of
$(y,\mathbf{x})$ such that $i\in\Lambda$, $n_{l,d}$ is the number of
partitions of $\\{1,\dots,d\\}$ of size $l=1,\dots,d$, and $r_{lm,d}$ is the
$m$th partition of $\\{1,\dots,d\\}$ of size $l$, with $1\leq m\leq n_{l,d}$.

In the particular case where we have a single covariate $(p=1)$, the
regression manifold $\mathscr{L}$ in (5) can be derived using properties of
bivariate copulas; see A. Accordingly, for an absolutely continuous angular
measure $H$ (with density $h)$, it follows that

$\displaystyle G_{Y\mid X}(y\mid
x)=2\exp\left\\{-2\int_{0}^{1}\max\left\\{\dfrac{w}{x},\dfrac{1-w}{y}\right\\}h(w)\,\mathrm{d}w+\frac{1}{x}\right\\}\int_{\omega(x,y)}^{1}wh(w)\,\mathrm{d}w,\quad
x,y>0,$ (8)

where $\omega(x,y)=x/(x+y)$, and $y_{q\mid x}$ is then calculated via (6). For
general $p\geq 1$, in the two special cases of independent and perfectly
dependent extremes the regression manifold $\mathscr{L}$ in (5) can be
obtained in a more direct way. For the former, we derive the regression lines
$L_{q}$ explicitly by first calculating $g_{Y\mid\mathbf{X}}$ and then using
that
$G_{Y\mid\mathbf{X}}(y\mid\mathbf{x})=\int_{0}^{y}g_{Y\mid\mathbf{X}}(z\mid\mathbf{x})\,\mathrm{d}z$
while for the latter, we provide approximations.

When extremes are independent, $H$ assigns equal mass to the boundaries of the
simplex, which also corresponds to asymptotic independence of $\mathbf{X}$ and
$Y$ (Hüsler and Li (2009)), resulting in

$\displaystyle L_{q}=\\{-1/\log q:\mathbf{x}\in(0,\infty)^{p}\\},$ (9)

with

$\displaystyle
G_{Y\mid\mathbf{X}}(q\mid\mathbf{x})=\int_{0}^{y}\dfrac{\exp(-z^{-1}-x^{-1}_{1}-\dots-
x_{p}^{-1})(-z^{-2})\prod\limits_{j=1}^{p}(-x_{j}^{-2})}{\exp(-x^{-1}_{1}-\dots-
x_{p}^{-1})\prod\limits_{j=1}^{p}(-x_{j}^{-2})}\,\mathrm{d}z=e^{-y^{-1}},\quad
y>0.$

When extremes are perfectly dependent, the angular measure $H$ assigns all its
mass to the barycenter of the simplex, $d^{-1}1{1}_{d}$, leading to
$G(y,\mathbf{x})=\exp\\{-\max(y^{-1},x_{1}^{-1},\dots,x_{p}^{-1})\\}$. Taking
derivatives of $G(y,\mathbf{x})$ in this case is non-trivial, and we replace
the maximum function with a soft maximum (Cook, 2011) so to obtain an
approximation for the shape of the regression lines for perfectly dependent
extremes. Thus, the soft maximum approximation for the regression lines for
perfectly dependent extremes is

$\displaystyle\tilde{L}_{q}=\\{\min(x_{1},\dots,x_{p}):\mathbf{x}\in(0,\infty)^{p}\\},$
(10)

with
$G_{Y\mid\mathbf{X}}(y\mid\mathbf{x})=I\\{y\geq\min(x_{1},\dots,x_{p})\\}$,
where $I$ is the indicator function; see B for the derivation.

Trivially, quantiles of a conditional extreme value distribution obey the
standard properties of quantile functions (see van der Vaart, 1998, Chap. 21).
Less trivial is however the fact that, monotone regression dependence of
bivariate extremes (Guillem, 2000, Theorem 1) implies that regression lines
$L_{q}$ in (5) are non-decreasing in $x$, for $p=1$, under some mild
assumptions.

######  Proposition 2.1 (Monotonicity of regression manifold).

Let $G_{Y|X}(y\mid x)=\mathbb{P}\left\\{Y\leq y\mid X=x\right\\}$ be a
conditional bivariate extreme value distribution function, which we assume to
be jointly continuously differentiable and strictly increasing in $y$ for any
fixed $x\in(0,\infty)$. Then, the regression lines for bivariate extremes
$(0,\infty)\ni x\mapsto y_{q\mid x}$ are non-decreasing for all $q\in(0,1)$.

###### Proof.

See C. ∎

An example of a bivariate extreme value distribution satisfying the
assumptions of Proposition 2.1 is the logistic model, whose regression
manifold is discussed in Example 2.2 below.

### Parametric examples of regression manifolds for block maxima on block
maxima

We now consider some parametric instances of regression manifolds as defined
in (5). Charts of regression manifolds for these parametric examples are
depicted in Figure 1. In D, we show that for sufficiently large $x$, the
following approximation holds for the regression manifold of the logistic
model,
$\tilde{L}_{q}=\\{\tilde{y}_{q\mid\mathbf{x}}:\mathbf{x}\in(0,\infty)^{p}\\}$,
with

$\tilde{y}_{q\mid
x}=\gamma_{q}+\beta_{q}x+\mathchoice{{\scriptstyle\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{{\scriptscriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptscriptstyle\mathcal{O}$}}{(x)}.$
(11)

Here $\gamma_{q}$ and $\beta_{q}$ are functions of $\alpha$ and $q$ (see (23)
and (22)), and
$\mathchoice{{\scriptstyle\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{{\scriptscriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptscriptstyle\mathcal{O}$}}{(x)}$
is
little-$\mathchoice{{\scriptstyle\mathcal{O}}}{{\scriptstyle\mathcal{O}}}{{\scriptscriptstyle\mathcal{O}}}{\scalebox{0.7}{$\scriptscriptstyle\mathcal{O}$}}$
of $x$ in Bachmann–Landau notation.

######  Example 2.2 (Logistic).

An instance of the logistic regression manifold can be found in Figure 1
(top). It stems from the logistic bivariate extreme value distribution
function given by

$\displaystyle G(x,y)=\exp\\{-(x^{-1/\alpha}+y^{-1/\alpha})^{\alpha}\\},\quad
x,y>0,$

where $\alpha\in(0,1]$ characterizes the dependence between extremes: The
closer $\alpha$ is to $0$, the stronger the dependence, with the limit
$\alpha\to 0$ corresponding to the case of perfect dependence. The conditional
distribution of $Y$ given $X$ is

$\displaystyle G_{Y\mid X}(y\mid
x)=G(x,y)(x^{-1/\alpha}+y^{-1/\alpha})^{\alpha-1}x^{1-1/\alpha}\exp(x^{-1}),\quad
x,y>0,$

thus leading to the following family of regression lines $L_{q}$ in (5) where

$\displaystyle y_{q\mid
x}=\left[\left\\{\dfrac{1-\alpha}{\alpha}xW\left(\dfrac{\alpha}{1-\alpha}x^{-1}e^{{\alpha}/(1-\alpha)x^{-1}}q^{{\alpha}/(\alpha-1)}\right)\right\\}^{1/\alpha}-1\right]^{-\alpha}x,$
(12)

and $x>0$. Here, $W$ is the so-called Lambert $W$ function, that is, the
multivalued analytic inverse of $f(z)=z\exp(z)$ with $z$ denoting a real or
complex number (Borwein and Lindstrom, 2016); see the supplementary material F
for further details. As it can be seen from Figure 1 (top), the regression
lines obey what is claimed in Proposition 2.1 in the sense that $(0,\infty)\ni
x\mapsto y_{q\mid x}$ are non-decreasing for all $q\in(0,1)$.

Logistic regression manifold

Husler–Reiss regression manifold

Coles–Tawn regression manifold

Figure 1: Regression manifold $\mathscr{L}$, as defined in (5), for bivariate
logistic, Husler–Reiss and Coles–Tawn models (top to bottom) with strong
dependence, intermediate and weak extremal dependence (left to right).

######  Example 2.3 (Husler–Reiss).

An instance of the Husler–Reiss regression manifold is depicted in Figure 1
(middle). It follows from the Husler–Reiss bivariate extreme value
distribution function which has the following form:

$\displaystyle
G(x,y)=\exp\left\\{-x^{-1}\Phi\left(\lambda+\dfrac{1}{2\lambda}\log\dfrac{y}{x}\right)-y^{-1}\Phi\left(\lambda+\dfrac{1}{2\lambda}\log\dfrac{x}{y}\right)\right\\},\quad
x,y>0,$

where $\Phi$ is the standard normal distribution function and
$\lambda\in(0,\infty]$ is the parameter regulating the dependence between
extremes: $\lambda\to 0$ corresponds to perfect dependence and the limit case
$\lambda\to\infty$ corresponds to complete independence. The family of
regression lines $L_{q}$ in (5) for this model does not have explicit
representations and is obtained using (6) with

$\displaystyle G_{Y\mid X}(y\mid x)$
$\displaystyle=\left[\Phi\left(\lambda+\dfrac{1}{2\lambda}\log\dfrac{y}{x}\right)+\dfrac{1}{2\lambda}\phi\left(\lambda+\log\dfrac{y}{x}\right)\dfrac{1}{2\lambda}-\dfrac{xy^{-1}}{2\lambda}\phi\left(\lambda+\dfrac{1}{2\lambda}\log\dfrac{x}{y}\right)\right]$
$\displaystyle\quad\times G(x,y)\exp\\{x^{-1}\\},\quad x,y>0,$

where $\phi$ is the standard normal density function.

######  Example 2.4 (Coles–Tawn).

An instance of the Coles–Tawn regression manifold is depicted in Figure 1
(bottom). It follows from the Coles–Tawn bivariate extreme value distribution
function which has the following form:

$\displaystyle
G(x,y)=\exp[-x^{-1}\\{1-\text{Be}(q;\alpha+1,\beta)\\}-y^{-1}\text{Be}(q;\alpha,\beta+1)],\quad
x,y>0,$

where $\text{Be}(q;a,b)$ is a cdf of a beta distribution function with
parameters $a,b>0$, $q=\alpha y^{-1}/(\alpha y^{-1}+\beta x^{-1})$ and
$\alpha,\beta>0$ are the parameters regulating dependence between extremes:
$\alpha=\beta=0$ corresponds to complete independence and
$\alpha=\beta\to\infty$ corresponds to perfect dependence; for fixed $\alpha$
($\beta$) the strength of dependence increases with $\beta$ ($\alpha$). The
family of regression lines $L_{q}$ in (5) for this model does not have an
explicit representation and is calculated using (6), for $x,y>0$, with

$\displaystyle G_{Y\mid X}(y\mid
x)=\Big{[}1-\text{Be}\left(q;\alpha+1,\beta\right)$
$\displaystyle+\dfrac{(\alpha+1)\beta}{\gamma}\text{be}\left(q;\alpha+2,\beta+1\right)$
$\displaystyle-\dfrac{x}{y}\dfrac{\alpha(\beta+1)}{\gamma}\text{be}\left(q;\alpha+1,\beta+2\right)\Big{]}x^{-2}G(x,y)\exp\\{x^{-1}\\},$

where $\text{be}(q;a,b)$ is the density function of the Beta distribution with
parameters $a,b>0$ and $\gamma=(\alpha+\beta+2)(\alpha+\beta+2)$.

This section introduced our key parameter of interest—regression manifolds for
block maxima on block maxima, i.e. $\mathscr{L}$ as in (5)—, it commented on
some of its properties, and gave examples of parametric instances. Next, we
discuss Bayesian inference for $\mathscr{L}$.

###  2.2 Learning about regression manifolds via Bernstein polynomials

### Induced prior on the space of regression manifolds

In this section we discuss how to learn about regression manifolds from data.
To achieve this, we will define an induced prior on the space of regression
manifolds by resorting to a flexible prior on the space of all angular
measures that was recently proposed by Hanson et al. (2017). To lay the
groundwork, we start by defining the setup of interest. Assume the
observations $\\{(S_{i},\mathbf{Z}_{i})\\}_{i=1}^{n}$, define the scalar
$R_{i}=S_{i}+\sum_{j=1}^{p}Z_{j,i}$ and the vector
$\mathbf{W}_{i}=(S_{i},\mathbf{Z}_{i})/R_{i}$, known as the pseudo-angular
decomposition of the observations. In de Haan and Resnick (1977) it is shown
the equivalence of the convergence of normalized componentwise maxima to $G$
to the following weak convergence of measures

$\displaystyle\mathbb{P}\left\\{\mathbf{W}\in\cdot\mid
R>u\right\\}\overset{\mathrm{d}}{\rightarrow}H(\cdot),\quad\text{as}\;u\rightarrow\infty.$

This means that when the radius $R$ is sufficiently large, the pseudo-angles
$\mathbf{W}$ are nearly independent of $R$ and follow approximately a
distribution associated with the angular measure $H$. Thus, to learn about
$L_{q}$ in (5), we first learn about $H$ based on
$\\{\mathbf{W}_{i}:R_{i}>u\\}_{i=1}^{n}$, for a large $u$, with a methodology
we describe next.

Following Hanson et al. (2017), we model the angular density $h$ via a
Bernstein polynomial defined on the unit simplex $\Delta_{d}$, and hence basis
polynomials are Dirichlet densities. Specifically, our specification for the
angular density is as follows

$\displaystyle
h_{\mathbf{v}}(\mathbf{w})=\sum\limits_{|\alpha|=J}v_{\alpha}\dfrac{\Gamma(|\alpha|)}{\prod\limits_{i=1}^{d}\Gamma(\alpha_{i})}\prod_{i=1}^{d}w_{i}^{\alpha_{i}-1},\quad\mathbf{w}\in\Delta_{d},$
(13)

where $v_{\alpha}>0$ are weights, $\alpha\in\mathbb{N}^{d}$,
$|\alpha|=\sum_{j=1}^{d}\alpha_{j}$, and where $J\in\mathbb{N}$ (with
$\mathbb{N}:=\\{1,2,3,\cdots\\}$). To ensure that the resulting
$h_{\mathbf{v}}(\mathbf{w})$ is a valid angular density (i.e. an actual
density satisfying the moment constraint (4)), the weights must obey

$\displaystyle\sum\limits_{|\alpha|=J}v_{\alpha}=1,\quad\sum\limits_{i=1}^{J-d+1}i\sum\limits_{|\alpha|=J,\alpha_{j}=i}v_{\alpha}=\dfrac{J}{d},$
(14)

for $j=1,\dots,d$.

To induce a prior in the space of regression manifolds we plug-in the angular
density in (13) into (7); subsequent integration with respect to $y$ and
inversion of $G_{Y|\mathbf{X}}(y|\mathbf{x})$ leads to an induced prior on the
space of regression lines $L_{q}$. In detail, to define a prior on the space
of regression manifolds we proceed as follows; to streamline the presentation
of ideas, we will first concentrate on the case $p=1$. The Bernstein
polynomial prior in (13) induces a prior on the space of regression lines
$L_{q}=\\{{y}_{q\mid x}:x\in(0,\infty)\\}$, where ${y}_{q\mid x}$ is a
solution to equation, ${G}_{Y\mid X}(y\mid x)=q$, for $q\in(0,1)$, and with

$\displaystyle{G}_{Y\mid X}(y\mid x)$
$\displaystyle=\dfrac{2}{J}\exp\Bigg{\\{}-\dfrac{2}{J}\sum\limits_{|\alpha|=J}v_{\alpha}\Big{(}\alpha_{1}x^{-1}\left[1-\text{Be}\left(\dfrac{x}{x+y};\alpha_{1}+1,\alpha_{2}\right)\right]+\alpha_{2}y^{-1}\text{Be}\left(\dfrac{x}{x+y};\alpha_{1},\alpha_{2}+1\right)\Big{)}\Bigg{\\}}$
$\displaystyle\quad\times\sum\limits_{|\alpha|=J}v_{\alpha}\alpha_{1}\left[1-\text{Be}\left(\dfrac{x}{x+y};\alpha_{1}+1,\alpha_{2}\right)\right]\exp\left\\{x^{-1}\right\\},\quad
x,y>0.$ (15)

When $p>2$ we proceed as above, that is our induced prior in the space of
regression lines is again induced by the Bernstein polynomial prior for the
angular density, and it follows by solving
${G}_{Y\mid\mathbf{X}}(y\mid\mathbf{x})=q$, with $h_{\mathbf{v}}(\mathbf{w})$
as in (13). The expression for the conditional multivariate extreme value
distribution ${G}_{Y\mid\mathbf{X}}(y\mid\mathbf{x})$ for $p>2$ is however not
as manageable as the one in (2).

### Induced prior on marginal effects

We now switch gears and focus on marginal effects. Generally speaking, one can
quantify the effect of a covariate on the quantiles of the response via so-
called marginal effects, expressed as the derivative of the regression lines
with respect to that covariate. For a given conditional quantile of the
response the marginal effect of the $i$th covariate is interpreted as the
change in this quantile of the response occurring due to a change in the $i$th
covariate, with other covariates being fixed.

The marginal effect of $x_{i}$ on $y_{q\mid\mathbf{x}}$, can be formally
defined by taking the derivative of the regression line $L_{q}$ with respect
to $x_{i}$, that is

$\beta_{q,x_{i}}=\frac{\partial y_{q\mid\mathbf{x}}}{\partial x_{i}},\quad
i=1,\dots,p.$

To ease notation we focus on the case $p=1$. Note first that $\beta_{q,x}\geq
0$ under mild assumptions. Indeed, recall that by Proposition 2.1, the
regression lines for bivariate extremes are non-decreasing in $x$ as
continuity of $H$ implies strict monotonicity of $G_{Y\mid X}$ in $y$, and
thus the marginal effects $\beta_{q,x}$ are non-negative for all
$x\in(0,\infty)$ and $q\in(0,1)$, provided that $h(w)>0$, for $w\in(0,1)$. To
learn about the marginal effect from data, we devise the following induced
prior

$\beta_{q,x}=\dfrac{A(x)\\{1-2A(x)\\}x^{-2}+B(x)y_{q\mid x}}{2A(x)C(x)y_{q\mid
x}^{-2}+B(x)x},$ (16)

with $q\in(0,1)$ and $x\in(0,\infty)$, and where

$\displaystyle A(x)$
$\displaystyle=J^{-1}\sum_{|\alpha|=J}v_{\alpha}\alpha_{1}\\{1-\text{Be}(x/(x+y_{q\mid
x});\alpha_{1}+1,\alpha_{2})\\},$ $\displaystyle B(x)$
$\displaystyle=x(x+y_{q\mid
x})^{-3}\sum_{|\alpha|=J}v_{\alpha}\text{be}(x/(x+y_{q\mid
x});\alpha_{1},\alpha_{2}),$ $\displaystyle C(x)$
$\displaystyle=J^{-1}\sum_{|\alpha|=J}v_{\alpha}\alpha_{2}\,\text{Be}(x/(x+y_{q\mid
x});\alpha_{1},\alpha_{2}+1).$

Some comments on the prior in (16) are in order. First, (16) is derived by
applying the implicit function theorem and calculating the derivative of
$y_{q\mid x}$ with respect to $x$, in a similar fashion as in (19). Second,
$\beta_{q,x}$ in (16) is a valid prior in the sense that $\beta_{q,x}\geq 0$;
the latter claim follows again from Proposition 2.1.

###  2.3 Induced prior for the $p$-covariate setting

The procedure discussed in Sections 2.1 and 2.2 becomes analytically
intractable and computationally expensive even for a relatively low-
dimensional covariate space. We thus propose an approach for learning about
the regression manifold $\mathscr{L}$, as defined in (5), via an approximation
to the conditional multivariate GEV density. Let
$\mathbf{u}=(y,\mathbf{x})\in(0,\infty)^{d}$ and
$\mathbf{u}(t)=(t,\mathbf{x})\in(0,\infty)^{d}$ with
$\lVert\mathbf{u}\rVert=y+\sum_{i=1}^{p}x_{i}>u$ for a large threshold $u$;
then, following Cooley et al. (2012, Proposition 1) the conditional density of
a multivariate extreme value distribution can be approximated, via a point
process representation for extremes, as follows

$g_{Y\mid\mathbf{X}}(y\mid\mathbf{x})=\dfrac{\lVert\mathbf{u}\rVert^{-d-1}h(\mathbf{u}/\lVert\mathbf{u}\rVert)}{\int_{0}^{\infty}\lVert\mathbf{u}(t)\rVert^{-d-1}h(\mathbf{u}(t)/\lVert\mathbf{u}(t)\rVert)\,\mathrm{d}t}.$
(17)

Scenario 1—strongly dependent extremes: Husler–Reiss model

Scenario 2—weakly dependent extremes: Logistic model

Scenario 3—asymmetric intermediate dependence: Coles–Tawn model

Figure 2: True regression manifold $\mathscr{L}$, as defined in (5), and its
estimate $\hat{\mathscr{L}}$ obtained using the methods from Section 2.2 for
Husler–Reiss, Logistic and Coles–Tawn bivariate extreme value models (top to
bottom) on a single-run experiment.

An induced prior for $g$ can be devised by plugging the approximation in (17)
with the specification from Sections 2.1, which leads to the following induced
prior for the conditional multivariate extreme value distribution function,

$\displaystyle
G_{Y\mid\mathbf{X}}(y\mid\mathbf{x})=\dfrac{\sum\limits_{|\alpha|=J}C_{\alpha}\prod\limits_{i=1}^{p}x_{i}^{\alpha_{i}-1}\int_{0}^{y}t^{\alpha_{d}-1}\lVert\mathbf{u}(t)\rVert^{-J-1}\,\mathrm{d}t}{\sum\limits_{|\alpha|=J}C_{\alpha}\prod\limits_{i=1}^{p}x_{i}^{\alpha_{i}-1}\int_{0}^{\infty}t^{\alpha_{d}-1}\lVert\mathbf{u}(t)\rVert^{-J-1}\,\mathrm{d}t},$
(18)

where
$C_{\alpha}=v_{\alpha}\Gamma(|\alpha|)/\prod_{i=1}^{d}\Gamma(\alpha_{i})$.
Hence, we can learn about the regression manifold $\mathscr{L}$ by estimating
$v_{\alpha}$ as described in Section 2.2, that is, by plugging in the
Bernstein polynomial estimates (13) into the approximation (17) and
numerically inverting it, (18), with respect to $y$. The merit of this
strategy is illustrated numerically in the supplementary materials (E).

##  3 Simulation study

###  3.1 Preliminary experiments

We study the finite sample performance of the proposed methods under three
data generating scenarios that were introduced in Section 2.1; see Examples
2.2—2.4. Specifically, we simulate data as follows:

  * 1.

Scenario 1—strongly dependent extremes: Husler–Reiss model with $\lambda=0.1$.

  * 2.

Scenario 2—weakly dependent extremes: Logistic model with $\alpha=0.9$.

  * 3.

Scenario 3—asymmetric intermediate dependence: Coles–Tawn model with
$\alpha=0.5$, $\beta=100$.

For now we focus on illustrating the methods in a single-run experiment; a
Monte Carlo simulation study will be reported in Section 3.2. To illustrate
how the resulted estimates compare with the true regression lines on a one-
shot experiment in each scenario we generate $n=5000$ samples
$\\{(Y_{i},X_{i})\\}_{i=1}^{n}$. For the analysis we use observations for
which $y_{i}+x_{i}>u$, where $u$ is a $95\%$ quantile of $Y+X$, providing
$k=250$ exceedances to fit the model. We take the number of basis functions
equal to the number of exceedances $k$. To learn about regression lines from
data, we exploit the single component adaptive Markov Chain Monte Carlo (MCMC)
with a wide Dirichlet prior, Dirichlet$(0.1\times 1{1}_{k})$, defined on a
generalized logit transformation of weights $v_{\alpha}$. The length of the
MCMC chain is $6000$ with a burn-in period of $2000$.

In Figure 2 we plot true and estimated regression manifolds under the three
scenarios above over the range $(x,y)\in(0,20]\times(0,20]$, where $20$
corresponds to the $95\%$ quantile of the unit Fréchet marginal distributions.
Figure 2 shows that the proposed estimator captures well the shape of
$\mathscr{L}$ for all three cases. However, for $q$ close to $0$ it
overestimates the dependence for weakly dependent extremes and underestimates
the dependence between strongly dependent extremes. This is even more evident
in Figure 3, where we depict cross sections of the regression manifolds for
the Husler–Reiss (top), logistic (central), and Coles–Tawn (bottom) model. The
plots on the left of Figure 3 represent regression lines for bivariate
extremes $L_{q}$ (for fixed $q$), and show that the stronger the dependence
between extremes, the closer the lines get to the identity $y=x$. Moreover,
one can notice that for strongly dependent extremes regression lines are
almost linear and for weakly dependent extremes they are curved for small $x$
and linear for large $x$. Indeed, for the logistic model for $\alpha<1$
regression lines become linear when $x$ increases; see D.

Scenario 1—strongly dependent extremes: Husler–Reiss model

Scenario 2—weakly dependent extremes: Logistic model

Scenario 3—asymmetric intermediate dependence: Coles–Tawn model

Figure 3: Posterior mean regression lines $L_{q}$ for $q=\\{0.1,0.5,0.9\\}$
and $x\in(0,20]$ (left) and conditional quantile curves $\\{y_{q\mid
x}:q\in(0,1)\\}$ along with credible bands, for $x=\\{1,10,20\\}$ (right) for
Husler–Reiss, Logistic and Coles–Tawn bivariate extreme value models (top to
bottom).

Figure 3 also illustrates that conditional quantiles for Scenarios $1$ and $3$
fall outside the credible intervals for extremal quantiles, i.e. $q\to 0\text{
or }1$, and for Scenario $2$ credible intervals do not cover the true values
for any $q$. The former drawback can be caused by less precise estimates of
the angular density close to the simplex boundaries, and the latter is
probably due to large bias that the estimators of angular densities often have
when applied to $U$-shaped angular densities. Nevertheless, since we work in
the extreme value setup, non-extremal quantiles are also useful when
considering extremal events.

Scenario 1—strongly dependent extremes: Husler–Reiss model

Scenario 2—weakly dependent extremes: Logistic model

Scenario 3—asymmetric intermediate dependence: Coles–Tawn model

Figure 4: Posterior mean regression lines $L_{q}$ for
$q=\\{0.1,0.45,0.55,0.9\\}$ and $x\in(0,20]$ estimated from $500$ Monte Carlo
simulations (gray lines) for sample size $k=500$ and true conditional
quantiles (black line) for Husler–Reiss, Logistic and Coles–Tawn bivariate
extreme value models (top to bottom).

###  3.2 Monte Carlo simulations

To conduct a simulation study we generate $500$ Monte Carlo samples of sizes
$n=5000$ and $n=10000$ resulting in $k=250$ and $k=500$ for the three
scenarios described in Section 3.1. We use the MCMC algorithm as described in
Section 3.1 with the same prior on $v_{\alpha}$’s. The performance of our
methods will be visualized via a comparison of posterior mean estimates of the
regression lines with the true regression lines $L_{q}$ for a few fixed
$q\in(0,1)$. We focus on the region $x\in(0,20]$ as the bivariate extreme
value concentrates most of its mass (at least $90$%) in the set
$(0,20]\times(0,20]$.

The regression lines corresponding to the described scenarios for $k=500$ are
shown in Figure 4 respectively; a similar chart for $k=250$ is available from
the supplementary materials (Figure 9). Figure 4 outlines that the model fits
the data from Scenario 1 reasonably well and, for weaker dependent extremes
(Scenario 2) and asymmetrically dependent extremes (Scenario 3), it provides
relatively precise estimates for middle values of $q$, but has expected bias
is visible for $q$ close to $0$ and $1$. Comparing different sample sizes
(i.e. comparing Figure 4 with Figure 9 in the supplementary material) we can
observe that increasing sample size reduces the variation of estimates for all
$q\in\\{0.1,0.45,0.55,0.9\\}$ and for all scenarios. Where the model performs
relatively well the estimates concentrate more around the true quantity
indicating improvement in the quality of fit, and where the fitted values
differ substantially from the true, raising sample size leads to larger bias.

##  4 Application to stock markets

###  4.1 Data, preprocessing, and applied rational for the analysis

We now apply the proposed method to the two of the world’s biggest stock
markets—the NASDAQ and NYSE (New York Stock Exchange). According to the
Statistics Portal of the World Federation of Exchanges
(https://statistics.world-exchanges.org), the total equity market
capitalization of the NASDAQ and NYSE are respectively 16 and 22 trillion US$,
as of 2020 / Oct, thus illustrating well the scale of these players in the
worldwide stock-exchange industry. The data were gathered from Yahoo Finance
(finance.yahoo.com), and consists of daily closing prices of the NASDAQ and
NYSE composite indices over the period from February 5, 1971 to November 23,
2020.

A key goal of the analysis will be to learn about the degree of association of
extreme losses in both of these markets through the lenses of our model, and
thus we focus on modeling negative log returns, which can be regarded as a
proxy for losses, and which consist of first differences of prices on a log-
scale; the resulting data for NASDAQ and NYSE is denoted below as
$\\{(X_{i},Y_{i})\\}_{i=1}^{n}$. The sample period under analysis is
sufficiently broad to cover a variety of major downturns and selloffs
including, for example, those related with the 2007–2010 subprime mortgage
crisis, the ongoing China–US trade war, and with the 2020 COVID-19 pandemic.
We take weekly maxima of negative log returns adjusted for heteroskedasticity
by fitting GARCH$(1,1)$ filter with normally distributed errors to both and
convert them to unit Fréchet margins via the transformation
$(\widehat{X}_{i},\widehat{Y}_{i})=(-1/\log\\{\widehat{F}_{X}(X_{i})\\},-1/\log\\{\widehat{F}_{Y}(Y_{i})\\}),$
where $\widehat{F}_{X}$ and $\widehat{F}_{Y}$ respectively denote the
empirical distribution functions of negative log return for NASDAQ ($X$) and
NYSE ($Y$). The raw data and resulting preprocessed data are depicted in
Figure 5. As can be seen from the latter figure the composite indices exhibit
a similar dynamics reacting to different economic shocks (9/11 attacks, 2001;
2008 financial crisis; China-US trade war started in 2018) alike. The shape of
the scatterplot of the negative log returns brought to unit Fréchet margins in
log-log scale above the boundary threshold evidences the extremal dependence
between log returns as they are quite concentrated around the identity line.

Figure 5:  Left: NASDAQ (red) and NYSE (blue) composite indices. Middle:
Scatterplot of negative log returns of NASDAQ and NYSE composite indices
converted to unit Fréchet margins; the solid line corresponds to the boundary
threshold in the log-log scale, with both axes being logarithmic. Right:
Angular density estimate with $95\%$ credible band along with a rug of pseudo-
angles.

###  4.2 Regression of extreme losses on extreme losses

We now apply our model so to learn about how the extreme losses on both
exchanges relate. To employ our model we start by fitting the angular density
via the Bernstein polynomial-based approach from Section 2 by using the
pseudo-angles based on thresholding the pseudo-radius at their 95% quantile.
Some comments on prior specification and on posterior inference are in order.
For our calculations we used a single component adaptive MCMC method with a
wide Dirichlet prior defined on a generalized logit transformation of weights
$v_{\alpha}$. We run a MCMC chain of length $10000$ with a burn-in period of
$2000$ and set the number of basis functions to be equal to the number of
exceedances.

Figure 6:  Left: Posterior mean regression manifold $\mathscr{L}$ for NYSE
given NASDAQ. Middle: Posterior mean regression lines $L_{q}$ for
$q=\\{0.1,0.5,0.9\\}$ for NYSE given NASDAQ along with $95\%$ credible bands
and plotted against filtered negative log returns on the indices converted to
unit Fréchet margins. Right: Posterior mean conditional quantile curves
$\\{y_{q\mid x}:q\in(0,1)\\}$ of filtered negative log returns on NYSE for
$x=\\{1,10,20\\}$, along with with $95\%$ credible bands, corresponding to
filtered negative log returns on NASDAQ converted to unit Fréchet margins.

The obtained fit for the angular density is reported in Figure 5 (right). As
is illustrated by this plot most of the observed pseudo-angles lie closer to
the middle of the interval $(0,1)$ and the estimate resembles a bell-shaped
density which suggests there is a dependence between extremal losses on NASDAQ
and NYSE composite indices. One can also observe that the estimated density
peaks at a point lying in the first half of the interval providing an
indication of an asymmetry in the dependence. Next we learn about the
regression manifold.

Figure 6 represents the resulting estimates of the regression manifold along
with cross-sections in $p$ and $x$ for negative log-returns on NASDAQ and NYSE
composite indices. The shape of the regression manifold suggests there is an
intermediate dependence between extremal losses on the specified indices and
the regression lines on the middle graph substantially differ from those
corresponding to independence and tend to be closer to the identity line.
Moreover, the cross-sections for different values of $x$ reveal considerable
variation in quantiles of $y$ supporting the conclusion about presence of the
dependence between negative log-returns.

Table 1: Marginal effects on $10\%$, $50\%$ and $90\%$ quantiles of filtered log-returns on NYSE, $Y$, evaluated at values $\\{1,10,20\\}$ of filtered log-returns on NASDAQ, $x$, in unit Fréchet margins; $95\%$ credible intervals in brackets quantiles | $x$
---|---
of $Y$ | $1$ | $10$ | $20$
$10\%$ | $0.165$ | $0.041$ | $0.040$
| $(0.141,0.193)$ | $(0.019,0.090)$ | $(0.016,0.093)$
$50\%$ | $1.133$ | $0.581$ | $0.567$
| $(1.080,1.191)$ | $(0.495,0.667)$ | $(0.473,0.657)$
$90\%$ | $3.677$ | $2.387$ | $2.317$
| $(3.193,4.539)$ | $(2.192,2.598)$ | $(2.110,2.516)$

Having assessed the degree of association between extreme values, we now
quantify the effect of extremal losses on NASDAQ index on extremal changes of
log-returns on NYSE index by calculating marginal effects. The marginal
effects on $10\%$, $50\%$ and $90\%$ quantiles of filtered log-returns on NYSE
calculated at values $\\{1,10,20\\}$ of filtered log-returns on NASDAQ (both
brought to unit Fréchet margins) are presented in Table 1. As can be seen from
the table the effect of NASDAQ on NYSE decreases when losses on NASDAQ
increase with more substantial fall in lower quantiles of losses on NYSE.
Moreover, for fixed losses on NASDAQ the effect on NYSE losses becomes larger
for higher quantiles.

##  5 Closing remarks

We propose a regression-type model for the setup where both the response and
the covariate are block maxima. The modeling starting point is the result that
the limiting behavior of the vector of properly standardized componentwise
maxima is given by a multivariate extreme value distribution. Conceptually,
the model is then constructed in a similar fashion as in quantile regression,
that is, by assessing how the conditional quantile of the response reacts to
changes in the covariate while it takes into account the latter asymptotic
result. An important target in the proposed framework is the regression
manifold, which consists of a family of regression lines obeying the proviso
of multivariate extreme value theory. A Bernstein polynomial prior on the
space of angular densities is used to learn about the model from data, with
numerical studies showcasing its flexibility.

One could wonder why not to resort to statistical models for nonstationary
extremes (e.g. Coles, 2001, Section 6) as an alternative to methods proposed
herein, as these can be used for assessing the effect of covariates on an
extreme-valued response, by indexing the parameters of the GEV distribution
with a covariate. Yet, since the latter models are built from the univariate
theory of extremes they are not tailored for conditioning on another variable
being extreme, as they fail to take on board information from the dependence
structure between the extremes. Other related approaches include extremal
quantile regression methods (Chernozhukov, 2005)—which similarly to the
statistical models for nonstationary extremes—have not been designed for
conditioning on another variable being extreme, as they do not take into
account the dependence structure between the extremes.

We close the manuscript with comments on future research. For regressions with
many predictors, it is likely that most covariates will have little effect on
the response and thus one could wonder how to devise a version of the proposed
method that shrinks towards zero the effect of such irrelevant covariates; the
development of a Lasso (Tibshirani, 1996) version of the proposed model would
thus seem natural for such situation, and is left as an open problem for
future research. Another natural avenue for future research would be to devise
regression-type methods for exceedances on exceedances by resorting to the so-
called multivariate generalized Pareto distribution (Kiriliouk et al., 2019),
rather than with the multivariate extreme value distribution as herein.
Finally, the development of a version of the model that could take into
account asymptotic independence by resorting to the hidden angular measure
(Ramos and Ledford, 2009), rather than the standard angular measure as herein,
would seem natural as well.

## Acknowledgement

We thank Johan Segers (Université catholique de Louvain) for insightful
discussions.

M. de Carvalho acknowledges support from the _Fundaç $\tilde{\text{a}}$o para
a Ci$\hat{e}$ncia e a Tecnologia_ (Portuguese Foundation for Science and
Technology) through the projects PTDC/MAT-STA/28649/2017 and
UID/MAT/00006/2020.

G. dos Reis acknowledges support from the _Fundaç $\tilde{\text{a}}$o para a
Ci$\hat{e}$ncia e a Tecnologia_ (Portuguese Foundation for Science and
Technology) through the project UIDB/00297/2020 (Centro de Matemática e
Aplicaç$\tilde{\text{o}}$es CMA/FCT/UNL).

A. Kumukova was supported by The Maxwell Institute Graduate School in Analysis
and its Applications, a Centre for Doctoral Training funded by the UK
Engineering and Physical Sciences Research Council (grant EP/L016508/01), the
Scottish Funding Council, Heriot-Watt University and the University of
Edinburgh.

##  Appendix A Conditional bivariate extreme value distribution

Here we derive the expression for the conditional bivariate extreme value
distribution function in (8). Recall that according to Sklar’s theorem
(Nelsen, 2007, Theorem 2.3.3), a joint bivariate distribution function
$G:\mathbb{R}^{2}\to[0,1]$ with continuous marginal distributions
$G_{X}:\mathbb{R}\to[0,1]$ and $G_{Y}:\mathbb{R}\to[0,1]$ can be uniquely
represented through a copula function $C$ for $(x,y)\in\mathbb{R}^{2}$
$G(x,y)=C(G_{X}(x),G_{Y}(y)),$ or, equivalently,
$C(u_{1},u_{2})=G(G^{-1}_{X}(u_{1}),G^{-1}_{Y}(u_{2}))$, for
$(u_{1},u_{2})\in[0,1]^{2}$, where $G^{-1}_{X}(q)=\inf\left\\{x:G_{X}(x)\geq
q\right\\}$. Using the following well-known property of copulas,

$\displaystyle C_{U_{2}\mid U_{1}}(u_{2}\mid
u_{1})\vcentcolon=\mathbb{P}\left(U_{2}\leq u_{2}\mid
U_{1}=u_{1}\right)=\frac{\partial C(u_{1},u_{2})}{\partial
u_{1}},\quad(u_{1},u_{2})\in[0,1]^{2},$

we calculate the conditional distribution $(Y\mid X)$ as

$\displaystyle G_{Y\mid X}(y\mid x)=C_{U_{2}\mid U_{1}}(e^{-y^{-1}}\mid
e^{-x^{-1}}).$

In our setting

$\displaystyle C(u_{1},u_{2})$
$\displaystyle=\exp\left\\{-2\int_{0}^{1}\max(-w\log u_{1},-(1-w)\log
u_{2})\,\mathrm{d}H(w)\right\\}.$

Assuming $H$ is absolutely continuous with density $h$, we have

$\displaystyle\frac{\partial}{\partial
x}2\int_{0}^{1}\max\left(\dfrac{w}{x},\dfrac{1-w}{y}\right)\,\mathrm{d}H(w)=2\frac{\partial}{\partial
x}\left(x^{-1}\int_{\omega(x,y)}^{1}w\,\mathrm{d}H(w)+y^{-1}\int_{0}^{\omega(x,y)}(1-w)\,\mathrm{d}H(w)\right)$
$\displaystyle\quad=2\left[-x^{-2}\int_{\omega(x,y)}^{1}w\,\mathrm{d}H(w)-(x^{-1}xy-y^{-1}y^{2})\dfrac{1}{(x+y)^{3}}h\\{\omega(x,y)\\}\right]=-2x^{-2}\int_{\omega(x,y)}^{1}wh(w)\,\mathrm{d}w,$

where $\omega(x,y)=x/(x+y)$. Then, the conditional copula has the following
form

$\displaystyle C_{U_{2}\mid U_{1}}(u_{2}\mid u_{1})$
$\displaystyle=2u^{-1}_{1}C(u_{1},u_{2})\int\limits_{\dfrac{\log u_{2}}{\log
u_{1}u_{2}}}^{1}wh(w)\,\mathrm{d}w,$

which in turn yields

$\displaystyle G_{Y|X}(y\mid
x)=2\exp\left\\{-2\int_{0}^{1}\max\left(\dfrac{w}{x},\dfrac{1-w}{y}\right)h(w)\,\mathrm{d}w+x^{-1}\right\\}\int_{w(x,y)}^{1}wh(w)\,\mathrm{d}w,\quad
x,y>0.$

##  Appendix B Soft-maximum approximation for regression manifold of perfectly
dependent extremes

Here we give details on the soft maximum approximation for the regression
lines for perfectly dependent extremes claimed in (10). We use a smooth
approximation of a maximum function called soft-maximum,

$\displaystyle
f(z_{1},\dots,z_{d};N)=\dfrac{1}{N}\log(e^{Nz_{1}}+\dots+e^{Nz_{d}}),$

which is infinitely differentiable everywhere and converges to a maximum
function as $N\to\infty$ (Cook, 2011). Then, the approximation of a
multivariate GEV distribution function for the case of perfect dependent
extremes, $G(y,\mathbf{x})=\max\\{y^{-1},x_{1}^{-1},\dots,x_{p}^{-1}\\}$ is

$\displaystyle\tilde{G}(y,\mathbf{x};N)=\exp\left\\{-\dfrac{1}{N}\log(e^{Ny^{-1}}+e^{Nx_{1}^{-1}}+\dots+e^{Nx_{p}^{-1}})\right\\}=(e^{Ny^{-1}}+e^{Nx_{1}^{-1}}+\dots+e^{Nx_{p}^{-1}})^{-1/N},$

and its partial derivative of order $d$ is

$\displaystyle\tilde{g}(y,\mathbf{x};N)$
$\displaystyle=\dfrac{\partial^{d}}{\partial y\partial x_{1}\cdots\partial
x_{p}}\tilde{G}(y,\mathbf{x};N)$
$\displaystyle=y^{-2}\prod\limits_{i=1}^{p}(1+iN)x_{i}^{-2}\exp\left(Ny^{-1}+\sum\limits_{i=1}^{p}x_{i}^{-1}\right)(e^{Ny^{-1}}+e^{Nx_{1}^{-1}}+\dots+e^{Nx_{p}^{-1}})^{-1/N-d}.$

This yields the following approximation of the conditional multivariate GEV
density

$\displaystyle\tilde{g}_{Y\mid\mathbf{X}}(y\mid\mathbf{x};N)$
$\displaystyle=\dfrac{y^{-2}\prod\limits_{i=1}^{p}(1+iN)x_{i}^{-2}\exp\left(Ny^{-1}+\sum\limits_{i=1}^{p}x_{i}^{-1}\right)(e^{Ny^{-1}}+e^{Nx_{1}^{-1}}+\dots+e^{Nx_{p}^{-1}})^{-1/N-d}}{\prod\limits_{i=1}^{p-1}(1+iN)\prod\limits_{i=1}^{p}x_{i}^{-2}\exp\left(\sum\limits_{i=1}^{p}x_{i}^{-1}\right)(e^{Nx_{1}^{-1}}+\dots+e^{Nx_{p}^{-1}})^{-1/N-p}}$
$\displaystyle=(1+pN)y^{-2}e^{Ny^{-1}}\frac{(e^{Ny^{-1}}+e^{Nx_{1}^{-1}}+\dots+e^{Nx_{p}^{-1}})^{-1/N-d}}{(e^{Nx_{1}^{-1}}\dots+e^{Nx_{p}^{-1}})^{-1/N-p}},$

and the conditional cumulative distribution function

$\displaystyle\tilde{G}_{Y\mid\mathbf{X}}(y\mid\mathbf{x};N)$
$\displaystyle=\int_{0}^{y}\tilde{g}_{Y\mid\mathbf{X}}(z\mid\mathbf{x})\,\mathrm{d}z$
$\displaystyle=(1+pN)(e^{Nx_{1}^{-1}}+\dots+e^{Nx_{p}^{-1}})^{1/N+p}\int_{0}^{y}z^{-2}e^{Nz^{-1}}(e^{Nz^{-1}}+e^{Nx_{1}^{-1}}+\dots+e^{Nx_{p}^{-1}})^{-1/N-d}\,\mathrm{d}z$
$\displaystyle=\dfrac{(1+pN)(-1/N)}{-1/N-d+1}(e^{Nx_{1}^{-1}}+\dots+e^{Nx_{p}^{-1}})^{1/N+p}(e^{Ny^{-1}}+e^{Nx_{1}^{-1}}+\dots+e^{Nx_{p}^{-1}})^{-1/N-d+1}$
$\displaystyle=\left(\dfrac{e^{Ny^{-1}}+e^{Nx_{1}^{-1}}+\dots+e^{Nx_{p}^{-1}}}{e^{Nx_{1}^{-1}}+\dots+e^{Nx_{p}^{-1}}}\right)^{-1/N-p}.$

Passing the last expression to the limit as $N\to\infty$ provides an ansatz
for the true conditional distribution function

$\displaystyle\tilde{G}_{Y\mid\mathbf{X}}(y\mid\mathbf{x})=\begin{cases}1,&y\geq\min(x_{1},\dots,x_{p})\\\
0,&y<\min(x_{1},\dots,x_{p})\end{cases}$

with $y,x_{1},\dots,x_{p}>0$, from where (10) follows.

##  Appendix C Proof of monotonicity of regression manifold

###### Proof of Proposition 2.1.

Since $y\mapsto G_{Y\mid X}(y\mid x)$ is continuous (strictly increasing) for
all $x\in(0,\infty)$, $y_{q\mid x}$ given by (6) is the solution to $G_{Y\mid
X}(y\mid x)=q$ for a fixed $q\in(0,1)$. Then $y$ satisfying $G_{Y\mid X}(y\mid
x)=q$ is an implicit function of $x$ parametrized by $q$. Under our
assumptions we apply the implicit function theorem and calculate the
derivative of $y_{q\mid x}$ with respect to $x$ via

$\displaystyle\frac{\partial}{\partial x}y_{q\mid
x}=-\dfrac{\frac{\partial}{\partial x}G_{Y\mid X}(y\mid
x)}{\frac{\partial}{\partial y}G_{Y\mid X}(y\mid x)}.$ (19)

Equation (19) combined with the monotone regression dependence property, i.e.
$x\mapsto G_{Y\mid X}(y\mid x)$ is non-increasing for all $y\in(0,\infty)$,
and the strict monotonicity of $y\mapsto G_{Y\mid X}(y\mid x)$ (increasing)
for all $x\in(0,\infty)$ gives

$\frac{\partial}{\partial x}y_{q\mid x}\geq 0.$

This completes the proof. ∎

##  Appendix D Exact and limiting regression manifolds for logistic model

Here we give details on how the exact (12) and approximated (11) regression
manifolds for the logistic model can be derived.

D.1. Exact regression manifold. Here we compute the conditional quantiles for
bivariate extreme value distribution and their linear approximation for large
$x$. Using (8), we calculate the conditional distribution function for the
logistic model: for $(x,y)\in(0,\infty)^{2}$

$\displaystyle G_{Y|X}(y\mid x)$
$\displaystyle=G(x,y)(x^{-1/\alpha}+y^{-1/\alpha})^{\alpha-1}x^{-1/\alpha-1}x^{2}\exp(x^{-1})$
$\displaystyle=\exp\\{-(x^{-1/\alpha}+y^{-1/\alpha})^{\alpha}+x^{-1}\\}(x^{-1/\alpha}+y^{-1/\alpha})^{\alpha-1}x^{1-1/\alpha}.$

The conditional quantiles behave differently depending on the strength of
dependence between extremes. The special case of the logistic model is for
$\alpha=1$, corresponding to independence between extremes for which the
family of regression lines are known to be given by (9). We now derive
conditional quantiles for bivariate dependent extremes, i.e. when
$\alpha\in[0,1)$. Since $y\mapsto G_{Y|X}$ is continuous, a conditional
quantile is a solution to

$\displaystyle\exp\\{-(x^{-1/\alpha}+y^{-1/\alpha})^{\alpha}\\}(x^{-1/\alpha}+y^{-1/\alpha})^{\alpha-1}x^{1-1/\alpha}\exp(x^{-1})=q,$
(20)

which can be written in terms of the Lambert $W$ function. Rewriting (20) as

$\displaystyle\exp\left\\{-\dfrac{\alpha}{\alpha-1}(x^{-1/\alpha}+y^{-1/\alpha})^{\alpha}\right\\}(x^{-1/\alpha}+y^{-1/\alpha})^{(\alpha-1){\alpha}/{(\alpha-1)}}x^{{(\alpha-1)}/{\alpha}{\alpha}/{(\alpha-1)}}\exp\left\\{\dfrac{\alpha}{\alpha-1}x^{-1}\right\\}=q^{{\alpha}/{(\alpha-1)}}$
$\displaystyle\Leftrightarrow$
$\displaystyle\exp\left\\{\dfrac{\alpha}{1-\alpha}(x^{-1/\alpha}+y^{-1/\alpha})^{\alpha}\right\\}\dfrac{\alpha}{1-\alpha}(x^{-1/\alpha}+y^{-1/\alpha})^{\alpha}=\dfrac{\alpha}{1-\alpha}x^{-1}\exp\left\\{\dfrac{\alpha}{1-\alpha}x^{-1}\right\\}q^{{\alpha}/(\alpha-1)}$
$\displaystyle\Leftrightarrow$
$\displaystyle\dfrac{\alpha}{1-\alpha}(x^{-1/\alpha}+y^{-1/\alpha})^{\alpha}=W\left(\dfrac{\alpha}{1-\alpha}x^{-1}e^{\frac{\alpha}{1-\alpha}x^{-1}}q^{{\alpha}/(\alpha-1)}\right),$

gives

$\displaystyle y_{q\mid
x}=\left[\left\\{\dfrac{1-\alpha}{\alpha}xW\left(\dfrac{\alpha}{1-\alpha}x^{-1}e^{\alpha/(1-\alpha)x^{-1}}q^{{\alpha}/(\alpha-1)}\right)\right\\}^{1/\alpha}-1\right]^{-\alpha}x,\quad
x>0.$ (21)

From properties of the Lambert $W$ function (see supplementary material for
details, F) it follows that

$\underset{x\rightarrow\infty}{\lim}x\,W\left(\dfrac{\alpha}{1-\alpha}\,x^{-1}\,e^{{\alpha}/(1-\alpha)x^{-1}}q^{{\alpha}/{(\alpha-1)}}\right)=\dfrac{\alpha}{1-\alpha}q^{{\alpha}/{(\alpha-1)}},$

and that the conditional quantiles tend to infinity as $x\rightarrow\infty$.

D.2. Limiting regression manifold. Below we show that (11) holds. To find
$\gamma_{q}$ and $\beta_{q}$ we use the expansion of the principal branch of
the Lambert $W$ function, a solution to $z=we^{w}$ when $z>0$, around $0$,
that is

$\displaystyle
W(z)=\sum\limits_{n=1}^{\infty}\dfrac{(-n)^{n-1}}{n!}z^{n}=z+\mathcal{O}(z^{2}),$

where $\mathcal{O}(z)$ is big-$\mathcal{O}$ of $z$. Substituting in (21) the
expansion of $W$ leads to

$\displaystyle y_{q\mid x}$
$\displaystyle=\left[\left\\{\dfrac{1-\alpha}{\alpha}x\left(\dfrac{\alpha}{1-\alpha}x^{-1}e^{\frac{\alpha}{1-\alpha}x^{-1}}q^{\dfrac{\alpha}{\alpha-1}}+\mathcal{O}(x^{-2}e^{{2\alpha}/{(1-\alpha)}x^{-1}})\right)\right\\}^{1/\alpha}-1\right]^{-\alpha}x$
$\displaystyle=\left[\left\\{e^{\frac{\alpha}{1-\alpha}x^{-1}}q^{\dfrac{\alpha}{\alpha-1}}+\mathcal{O}(x^{-1}e^{{2\alpha}/{(1-\alpha)}x^{-1}})\right\\}^{1/\alpha}-1\right]^{-\alpha}x,$

and taking $x\to\infty$ we find the linear asymptote of $x\mapsto y_{q\mid x}$
(we use L’Hospital in the computation of $\gamma_{q}$)

$\displaystyle\beta_{q}$
$\displaystyle=\lim\limits_{x\to\infty}\dfrac{y_{q\mid
x}}{x}=\\{q^{-\frac{1}{1-\alpha}}-1\\}^{-\alpha},$ (22)

and, the more involved calculation of $\gamma_{q}$. The derivative of the
asymptotic expansion of a function does not necessarily correspond to the
asymptotic expansion of the derivative of the function, hence, we use
L’Hospital with (21) and only then inject the asymptotic expansions. We have
then after some tedious derivations that

$\displaystyle\gamma_{q}$
$\displaystyle=\lim\limits_{x\to\infty}\left(y_{q\mid
x}-\beta_{q}x\right)=\dfrac{\alpha}{1-\alpha}\\{q^{{1}/{(\alpha-1)}}-1\\}^{-\alpha-1}\\{q^{{\alpha}/{(1-\alpha)}}-1\\}q^{\dfrac{1}{\alpha-1}}.$
(23)

The resulted linear approximation is

$\displaystyle\tilde{y}_{q\mid
x}=\dfrac{\alpha}{1-\alpha}\\{q^{{1}/{(\alpha-1)}}-1\\}^{-\alpha-1}\\{q^{{\alpha}/{(1-\alpha)}}-1\\}q^{{1}/{(\alpha-1)}}+\\{q^{-{1}/{(1-\alpha)}}-1\\}^{-\alpha}x,\quad
q\in(0,1),\;x\gg 1.$

Supplementary materials

##  Appendix E Additional numerical evidence

###  E.1 Induced prior for $p$-covariate setting

We report on two one-shot numerical experiments aimed at illustrating the
approach in Section 2.3 in the paper, that induces a prior on the space of all
regression manifolds by resorting to Bernstein polynomials and an
approximation of a multivariate GEV density due to Cooley et al. (2012). For
the numerical experiments in this supplementary material, we test our model by
taking a trivariate logistic extreme value distribution with dependence
parameter $\alpha=0.1$ (‘strongly’ dependent extremes) for the case $p=2$,
i.e. with the trivariate GEV distribution

$\displaystyle
G(y,x_{1},x_{2})=\exp\\{-(y^{-1/\alpha}+x_{1}^{-1/\alpha}+x_{2}^{-1/\alpha})^{\alpha}\\},\quad
y,x_{1},x_{2}>0.$

We generate two samples of sizes $n=10000$ and $n=20000$ which, after
thresholding at $95\%$ empirical quantiles of the observed pseudo angles,
yield $k=500$ and $k=1000$ data points to fit the model. Here, we use a
similar prior specification and MCMC setup as in Section 3.1 of the paper.

Figure 7 indicates that the proposed estimator of the angular density captures
reasonably well the dependence between extremes by concentrating around the
barycenter of the simplex, though in a less pronounced form than the true
density. As can be seen from Figure 8 the resulting regression lines estimates
resemble the true regression lines $L_{q}$, and increasing sample size
improves the fit as the lateral surfaces of the estimates become more
slanting, for $q=\\{0.3,0.5,0.7\\}$.

Figure 7:  Level plots of the true angular density (left) along with the
posterior mean estimate resulting from the methods from Section 2.2 (right) on
$n=10000$ observations for the trivariate logistic extreme value distribution,
on a single-run experiment, with the dependence parameter $\alpha=0.1$.

$q=0.3$

$q=0.5$

$q=0.7$

Figure 8: The true $L_{q}$ (left) along with the posterior mean estimate
resulting from the methods from Section 2.3 on $n=10000$ (middle) and
$n=20000$ (right) observations for $q=\\{0.3,0.5,0.7\\}$ (top to bottom) for
the trivariate logistic extreme value distribution, on a single-run
experiment, with the dependence parameter $\alpha=0.1$ over the domain
$\mathbf{x}=(x_{1},x_{2})\in(0,20]^{2}$.

###  E.2 Supporting Monte Carlo experiment

Figure 9 below complements Figure 4; the number of exceedances of the figure
reported here is $k=250$ while that on Figure 4 is $k=500$.

Scenario 1–strongly dependent extremes: Husler–Reiss model

Scenario 2–weakly dependent extremes: Logistic model

Scenario 3–asymmetric intermediate dependence: Coles–Tawn model

Figure 9: Regression lines $L_{q}$ for $q=\\{0.1,0.45,0.55,0.9\\}$ and
$x\in(0,20]$ estimated from $500$ Monte Carlo simulations (gray lines) for
sample size $k=250$ and true conditional quantiles (black line) for
Husler–Reiss, Logistic and Coles–Tawn bivariate extreme value models (top to
bottom).

##  Appendix F Details on the Lambert $W$ function

The Lambert $W$ function is used in the paper for deriving a linear
approximation to the logistic regression manifold (cf Appendix D), and thus we
offer here some details on it. Formally, the Lambert $W$ function is a set of
functions representing the inverse relation of the function $f(z)=ze^{z}$ for
any complex $z$. Since we deal only with positive real valued $z$, the
equation $f(z)=ze^{z}$ has only one solution $w=W(z)$, with $W$ being the
principal branch of the Lambert $W$ function. A useful property of this
function is that for any constant $a\in\mathbb{R}$ one has

$\underset{z\rightarrow\infty}{\lim}zW(a/z)=\underset{z\rightarrow\infty}{\lim}ae^{-W(a/z)}=a,$

which is derived from

$\underset{z\rightarrow\infty}{\lim}\frac{a}{z}=\underset{z\rightarrow\infty}{\lim}e^{W(a/z)}W(a/z)\quad\Rightarrow\quad\underset{z\rightarrow\infty}{\lim}W(a/z)=0.$

See Borwein and Lindstrom (2016) for further details.

## References

  * Beirlant et al. (2004) Beirlant, J., Goegebeur, Y., Segers, J., and Teugels, J. (2004), Statistics of Extremes: Theory and Applications, Wiley, Hoboken, NJ: Wiley.
  * Borwein and Lindstrom (2016) Borwein, J. M. and Lindstrom, S. B. (2016), “Meetings with Lambert $\mathcal{W}$ and other special functions in optimization and analysis,” Pure Appl. Funct. Anal., 1, 361–396.
  * Chernozhukov (2005) Chernozhukov, V. (2005), “Extremal Quantile Regression,” Annals of Statistics, 33, 806–839.
  * Coles (2001) Coles, S. (2001), An Introduction to Statistical Modeling of Extreme Values, Springer Series in Statistics, London: Springer London.
  * Cook (2011) Cook, J. D. (2011), “Basic properties of the soft maximum,” bepress,Working Paper Series 70.
  * Cooley et al. (2012) Cooley, D., Davis, R. A., and Naveau, P. (2012), “Approximating the conditional density given large observed values via a multivariate extremes framework, with application to environmental data,” Ann. Appl. Stat., 6, 1406–1429.
  * de Haan and Ferreira (2006) de Haan, L. and Ferreira, A. (2006), Extreme Value Theory: An Introduction, New York: Springer.
  * de Haan and Resnick (1977) de Haan, L. and Resnick, S. I. (1977), “Limit Theory for Multivariate Sample Extremes,” Zeitschrift für Wahrscheinlichkeitstheorie und Verwandte Gebiete, 40, 317–337.
  * Embrechts et al. (1997) Embrechts, P., Klüppelberg, C., and Mikosch, T. (1997), Modelling Extremal Events for Insurance and Finance, New York: Springer.
  * Gudendorf and Segers (2010) Gudendorf, G. and Segers, J. (2010), “Extreme-Value Copulas,” in Copula Theory and Its Applications, Springer, Berlin, Heidelberg, pp. 127–145.
  * Guillem (2000) Guillem, A. I. G. (2000), “Structure de dépendance des lois de valeurs extrêmes bivariées,” Comptes Rendus de l’Académie des Sciences-Series I-Mathematics, 330, 593–596.
  * Hanson et al. (2017) Hanson, T. E., de Carvalho, M., and Chen, Y. (2017), “Bernstein Polynomial Angular Densities of Multivariate Extreme Value Distributions,” Statistics and Probability Letters, 128, 60–66.
  * Hüsler and Li (2009) Hüsler, J. and Li, D. (2009), “Testing asymptotic independence in bivariate extremes,” Journal of Statistical Planning and Inference, 139, 990–998.
  * Kiriliouk et al. (2019) Kiriliouk, A., Rootzén, H., Segers, J., and Wadsworth, J. L. (2019), “Peaks over thresholds modeling with multivariate generalized Pareto distributions,” Technometrics, 61, 123–135.
  * Koenker (2005) Koenker, R. (2005), Quantile Regression, Cambridge, MA: Cambridge University Press.
  * Longin (2017) Longin, F. (ed.) (2017), Extreme Events in Finance: A Handbook of Extreme Value Theory and Its Applications, New York: Wiley.
  * McNeil et al. (2015) McNeil, A. J., Frey, R., and Embrechts, P. (2015), Quantitative risk management: concepts, techniques and tools-revised edition, Princeton: Princeton University Press.
  * Nelsen (2007) Nelsen, R. B. (2007), An introduction to copulas, Springer Science & Business Media.
  * Pickands (1981) Pickands, J. (1981), “Multivariate Extreme Value Distributions,” in Proc. 43th Sess. Int. Statist. Inst, pp. 859–878.
  * Ramos and Ledford (2009) Ramos, A. and Ledford, A. (2009), “A New Class of Models for Bivariate Joint Tails,” Journal of the Royal Statistical Society. Series B (Statistical Methodology), 71, 219–241.
  * Stephenson and Tawn (2005) Stephenson, A. and Tawn, J. (2005), “Exploiting Occurrence Times in Likelihood Inference for Componentwise Maxima,” Biometrika, 92, 213–227.
  * Tibshirani (1996) Tibshirani, R. (1996), “Regression shrinkage and selection via the lasso,” Journal of the Royal Statistical Society: Series B, 58, 267–288.
  * van der Vaart (1998) van der Vaart, A. W. (1998), Asymptotic Statistics, Cambridge series in statistical and probabilistic mathematics, Cambridge, UK ; New York, NY, USA: Cambridge University Press.
  * Wang and Stoev (2011) Wang, Y. and Stoev, S. A. (2011), “Conditional sampling for spectrally discrete max-stable random fields,” Advances in Applied Probability, 43, 461–483.